FEI Version of the Day
Loading...
Searching...
No Matches
fei_Factory_Trilinos.cpp
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43
44#include "fei_trilinos_macros.hpp"
45
46#include <fei_Factory_Trilinos.hpp>
47
48#include <fei_VectorReducer.hpp>
49#include <fei_MatrixReducer.hpp>
50#include <fei_Matrix_Local.hpp>
51#include <fei_Vector_Local.hpp>
52
53#ifdef HAVE_FEI_AZTECOO
54#include <fei_Solver_AztecOO.hpp>
55#endif
56#ifdef HAVE_FEI_BELOS
57#include <fei_Solver_Belos.hpp>
58#endif
59#ifdef HAVE_FEI_AMESOS
60#include <fei_Solver_Amesos.hpp>
61#endif
62
63#ifdef HAVE_FEI_EPETRA
64
65Factory_Trilinos::Factory_Trilinos(MPI_Comm comm)
66 : fei::Factory(comm),
67 comm_(comm),
68 reducer_(),
69 lpm_epetrabasic_(),
70 use_lpm_epetrabasic_(false),
71 useAmesos_(false),
72 useBelos_(false),
73 use_feiMatrixLocal_(false),
74 blockEntryMatrix_(false),
75 orderRowsWithLocalColsFirst_(false),
76 outputLevel_(0)
77{
78}
79
80Factory_Trilinos::~Factory_Trilinos()
81{
82}
83
84int Factory_Trilinos::parameters(int numParams,
85 const char* const* paramStrings)
86{
87 std::vector<std::string> stdstrings;
88 fei::utils::char_ptrs_to_strings(numParams, paramStrings, stdstrings);
89
90 fei::ParameterSet paramset;
91 fei::utils::parse_strings(stdstrings, " ", paramset);
92
93 parameters(paramset);
94 return(0);
95}
96
97void Factory_Trilinos::parameters(const fei::ParameterSet& parameterset)
98{
99 fei::Factory::parameters(parameterset);
100
101 parameterset.getIntParamValue("outputLevel", outputLevel_);
102
103 bool blkGraph = false;
104 bool blkMatrix = false;
105
106 parameterset.getBoolParamValue("BLOCK_GRAPH", blkGraph);
107 parameterset.getBoolParamValue("BLOCK_MATRIX", blkMatrix);
108
109 blockEntryMatrix_ = (blkGraph || blkMatrix);
110
111 const fei::Param* param = parameterset.get("Trilinos_Solver");
112 if (param != 0) {
113 if (param->getType() == fei::Param::STRING) {
114 std::string strval = param->getStringValue();
115 std::string::size_type ii = strval.find("Amesos");
116 if (ii != std::string::npos) {
117 useAmesos_ = true;
118 }
119 ii = strval.find("Belos");
120 if (ii != std::string::npos) {
121 useBelos_ = true;
122 }
123 }
124 }
125
126 parameterset.getBoolParamValue("LPM_EpetraBasic", use_lpm_epetrabasic_);
127 if (use_lpm_epetrabasic_) {
128 create_LinProbMgr();
129 }
130
131 parameterset.getBoolParamValue("USE_FEI_MATRIX_LOCAL", use_feiMatrixLocal_);
132
133 parameterset.getBoolParamValue("ORDER_ROWS_WITH_LOCAL_COLS_FIRST",
134 orderRowsWithLocalColsFirst_);
135}
136
138Factory_Trilinos::createMatrixGraph(fei::SharedPtr<fei::VectorSpace> rowSpace,
140 const char* name)
141{
142 static fei::MatrixGraph_Impl2::Factory factory2;
143 if (!use_lpm_epetrabasic_) {
144 return(factory2.createMatrixGraph(rowSpace, colSpace, name));
145 }
146
147 return(factory2.createMatrixGraph(rowSpace, colSpace, name));
148}
149
151Factory_Trilinos::wrapVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
153{
155 if (vecSpace->getNumIndices_Owned() != multiVec->Map().NumMyPoints()) {
156 //vecSpace isn't compatible with multiVec's map, so return empty vecptr
157 return(vecptr);
158 }
159
161 multiVec.get(),
162 multiVec->Map().NumMyPoints(), false));
163 return(vecptr);
164}
165
167Factory_Trilinos::wrapVector(fei::SharedPtr<fei::MatrixGraph> matGraph,
169{
171 if (matGraph->getRowSpace()->getNumIndices_Owned() !=
172 multiVec->Map().NumMyPoints()) {
173 //vector-space isn't compatible with multiVec's map, so return empty vecptr
174 return(vecptr);
175 }
176
177 vecptr.reset(new fei::Vector_Impl<Epetra_MultiVector>(matGraph->getRowSpace(),
178 multiVec.get(),
179 multiVec->Map().NumMyPoints(), false));
180 return(vecptr);
181}
182
184Factory_Trilinos::createVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
185 bool isSolutionVector,
186 int numVectors)
187{
188 if (use_feiMatrixLocal_) {
189 fei::SharedPtr<fei::Vector> feivec(new fei::Vector_Local(vecSpace));
190 return(feivec);
191 }
192
193// if (blockEntryMatrix_) {
194// throw std::runtime_error("Factory_Trilinos: fei ERROR, block-entry matrices/vectors not currently supported.");
195// }
196
197 std::vector<int> indices;
198 int err = 0, localSize;
199 if (reducer_.get() != NULL) {
200 indices = reducer_->getLocalReducedEqns();
201 localSize = indices.size();
202 }
203 else {
204 if (blockEntryMatrix_) {
205 localSize = vecSpace->getNumBlkIndices_Owned();
206 indices.resize(localSize*2);
207 err = vecSpace->getBlkIndices_Owned(localSize, &indices[0], &indices[localSize], localSize);
208 }
209 else {
210 localSize = vecSpace->getNumIndices_Owned();
211 indices.resize(localSize);
212 err = vecSpace->getIndices_Owned(indices);
213 }
214 }
215 if (err != 0) {
216 throw std::runtime_error("error in vecSpace->getIndices_Owned");
217 }
218
219 fei::SharedPtr<fei::Vector> feivec, tmpvec;
220 if (!use_lpm_epetrabasic_) {
221 try {
222 Epetra_BlockMap emap = blockEntryMatrix_ ?
223 Trilinos_Helpers::create_Epetra_BlockMap(vecSpace) :
224 Trilinos_Helpers::create_Epetra_Map(comm_, indices);
225
226 Epetra_MultiVector* emvec = new Epetra_MultiVector(emap, numVectors);
227
229 emvec, localSize,
230 isSolutionVector, true));
231 }
232 catch(std::runtime_error& exc) {
233 fei::console_out() << "Factory_Trilinos::createVector: caught exception '"
234 << exc.what() << "', re-throwing..." << FEI_ENDL;
235 throw exc;
236 }
237 }
238 else {
240 lpm_epetrabasic_.get(), localSize,
241 isSolutionVector));
242 }
243
244 if (reducer_.get() != NULL) {
245 feivec.reset(new fei::VectorReducer(reducer_,
246 tmpvec, isSolutionVector));
247 }
248 else {
249 feivec = tmpvec;
250 }
251
252 return(feivec);
253}
254
256Factory_Trilinos::createVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
257 int numVectors)
258{
259 bool isSolnVector = false;
260 return(createVector(vecSpace, isSolnVector, numVectors));
261}
262
264Factory_Trilinos::createVector(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
265 int numVectors)
266{
267 bool isSolnVector = false;
268 return(createVector(matrixGraph, isSolnVector, numVectors));
269}
270
272Factory_Trilinos::createVector(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
273 bool isSolutionVector,
274 int numVectors)
275{
276 if (use_feiMatrixLocal_) {
277 fei::SharedPtr<fei::Vector> feivec(new fei::Vector_Local(matrixGraph->getRowSpace()));
278 return(feivec);
279 }
280
281 int globalNumSlaves = matrixGraph->getGlobalNumSlaveConstraints();
282
283 if (globalNumSlaves > 0 && reducer_.get()==NULL) {
284 reducer_ = matrixGraph->getReducer();
285 }
286
287 fei::SharedPtr<fei::Vector> feivec, tmpvec;
288
289 std::vector<int> indices;
290 int err = 0, localSize;
291 fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
292 if (reducer_.get() != NULL) {
293 indices = reducer_->getLocalReducedEqns();
294 localSize = indices.size();
295 }
296 else {
297 localSize = vecSpace->getNumIndices_Owned();
298 indices.resize(localSize);
299 err = vecSpace->getIndices_Owned(indices);
300 }
301 if (err != 0) {
302 throw std::runtime_error("error in vecSpace->getIndices_Owned");
303 }
304
305 if (!use_lpm_epetrabasic_) {
306#ifdef HAVE_FEI_EPETRA
307 try {
308 Epetra_BlockMap emap = blockEntryMatrix_ ?
309 Trilinos_Helpers::create_Epetra_BlockMap(vecSpace) :
310 Trilinos_Helpers::create_Epetra_Map(comm_, indices);
311
312 Epetra_MultiVector* emvec = new Epetra_MultiVector(emap, numVectors);
313
314 tmpvec.reset(new fei::Vector_Impl<Epetra_MultiVector>(matrixGraph->getRowSpace(), emvec,
315 localSize, isSolutionVector, true));
316 }
317 catch(std::runtime_error& exc) {
318 fei::console_out() << "Factory_Trilinos::createVector: caught exception '"
319 << exc.what() << "', re-throwing..." << FEI_ENDL;
320 throw exc;
321 }
322#else
323 fei::console_out() << "fei_Factory_Trilinos::createVector ERROR, HAVE_FEI_EPETRA not defined."
324 << FEI_ENDL;
325#endif
326 }
327 else {
328#ifdef HAVE_FEI_EPETRA
329 vecSpace = matrixGraph->getRowSpace();
330
331 lpm_epetrabasic_->setRowDistribution(indices);
333 lpm_epetrabasic_.get(), localSize, isSolutionVector));
334#endif
335 }
336
337 if (reducer_.get() != NULL) {
338 feivec.reset(new fei::VectorReducer(reducer_, tmpvec, isSolutionVector));
339 }
340 else {
341 feivec = tmpvec;
342 }
343
344 return(feivec);
345}
346
348Factory_Trilinos::createMatrix(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
349{
351#ifdef HAVE_FEI_EPETRA
352 int globalNumSlaves = matrixGraph->getGlobalNumSlaveConstraints();
353
354 if (globalNumSlaves > 0 && reducer_.get()==NULL) {
355 reducer_ = matrixGraph->getReducer();
356 }
357
358 if (use_lpm_epetrabasic_) {
359 return(
360 Trilinos_Helpers::create_from_LPM_EpetraBasic(matrixGraph,
361 blockEntryMatrix_,
362 reducer_,
363 lpm_epetrabasic_)
364 );
365 }
366
367 if (use_feiMatrixLocal_) {
368 return(fei::Matrix_Local::create_Matrix_Local(matrixGraph, blockEntryMatrix_));
369 }
370
371 return(
372 Trilinos_Helpers::create_from_Epetra_Matrix(matrixGraph, blockEntryMatrix_,
373 reducer_, orderRowsWithLocalColsFirst_)
374 );
375#else
376 fei::console_out() << "fei_Factory_Trilinos::createMatrix ERROR, HAVE_FEI_EPETRA "
377 << "not defined."<<FEI_ENDL;
378 return feimat;
379#endif
380}
381
383Factory_Trilinos::createSolver(const char* name)
384{
386
387 std::string::size_type ii_amesos = std::string::npos;
388 std::string::size_type ii_belos = std::string::npos;
389
390 if (name != 0) {
391 std::string sname(name);
392 ii_amesos = sname.find("Amesos");
393 ii_belos = sname.find("Belos");
394 }
395
396 if (useAmesos_ || ii_amesos != std::string::npos) {
397#ifdef HAVE_FEI_AMESOS
398 solver.reset(new Solver_Amesos);
399 return solver;
400#else
401 fei::console_out() << "fei_Factory_Trilinos::createSolver: ERROR, Amesos not available (not enabled at compile-time?)"<<FEI_ENDL;
402#endif
403 }
404
405 if (useBelos_ || ii_belos != std::string::npos) {
406#ifdef HAVE_FEI_BELOS
407 solver.reset(new Solver_Belos);
408 return solver;
409#else
410 fei::console_out() << "fei_Factory_Trilinos::createSolver: ERROR, Belos not available (not enabled at compile-time?)"<<FEI_ENDL;
411#endif
412 }
413
414#ifdef HAVE_FEI_AZTECOO
415 solver.reset(new Solver_AztecOO);
416 return solver;
417#else
418 fei::console_out() << "fei_Factory_Trilinos::createSolver: ERROR, AztecOO not "
419 << "available." << FEI_ENDL;
420 return solver;
421#endif
422}
423
424void Factory_Trilinos::create_LinProbMgr(bool replace_if_already_created)
425{
426 if (!use_lpm_epetrabasic_) {
427 return;
428 }
429
430 bool need_to_create_lpm = false;
431
432 if (lpm_epetrabasic_.get() == NULL) {
433 need_to_create_lpm = true;
434 }
435
436 if (replace_if_already_created) {
437 need_to_create_lpm = true;
438 }
439
440 if (need_to_create_lpm) {
441#ifdef HAVE_FEI_EPETRA
443 newlpm(new LinProbMgr_EpetraBasic(comm_));
444
445 lpm_epetrabasic_ = newlpm;
446#else
447 fei::console_out() << "fei_Factory_Trilinos::create_LinProbMgr ERROR, HAVE_FEI_EPETRA"
448 <<" not defined."<<FEI_ENDL;
449#endif
450 }
451}
452
453//HAVE_FEI_EPETRA
454#endif
455
virtual void parameters(const fei::ParameterSet &paramset)
virtual void setRowDistribution(const std::vector< int > &ownedGlobalRows)=0
virtual fei::SharedPtr< fei::MatrixGraph > createMatrixGraph(fei::SharedPtr< fei::VectorSpace > rowSpace, fei::SharedPtr< fei::VectorSpace > columnSpace, const char *name)
const std::string & getStringValue() const
ParamType getType() const
Definition fei_Param.hpp:98
int getIntParamValue(const char *name, int &paramValue) const
const Param * get(const char *name) const
int getBoolParamValue(const char *name, bool &paramValue) const
void reset(T *p=0)
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
void char_ptrs_to_strings(int numStrings, const char *const *charstrings, std::vector< std::string > &stdstrings)
std::ostream & console_out()