FEI Version of the Day
fei_Factory_Trilinos.cpp
00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include "fei_trilinos_macros.hpp"
00010 
00011 #include <fei_Factory_Trilinos.hpp>
00012 
00013 #include <fei_VectorReducer.hpp>
00014 #include <fei_MatrixReducer.hpp>
00015 #include <fei_Matrix_Local.hpp>
00016 #include <fei_Vector_Local.hpp>
00017 
00018 #ifdef HAVE_FEI_AZTECOO
00019 #include <fei_Solver_AztecOO.hpp>
00020 #endif
00021 #ifdef HAVE_FEI_AMESOS
00022 #include <fei_Solver_Amesos.hpp>
00023 #endif
00024 
00025 #ifdef HAVE_FEI_EPETRA
00026 
00027 Factory_Trilinos::Factory_Trilinos(MPI_Comm comm)
00028   : fei::Factory(comm),
00029     comm_(comm),
00030     reducer_(),
00031     lpm_epetrabasic_(),
00032     use_lpm_epetrabasic_(false),
00033     useAmesos_(false),
00034     use_feiMatrixLocal_(false),
00035     blockEntryMatrix_(false),
00036     orderRowsWithLocalColsFirst_(false),
00037     outputLevel_(0)
00038 {
00039 }
00040 
00041 Factory_Trilinos::~Factory_Trilinos()
00042 {
00043 }
00044 
00045 int Factory_Trilinos::parameters(int numParams,
00046                                   const char* const* paramStrings)
00047 {
00048   std::vector<std::string> stdstrings;
00049   fei::utils::char_ptrs_to_strings(numParams, paramStrings, stdstrings);
00050 
00051   fei::ParameterSet paramset;
00052   fei::utils::parse_strings(stdstrings, " ", paramset);
00053 
00054   parameters(paramset);
00055   return(0);
00056 }
00057 
00058 void Factory_Trilinos::parameters(const fei::ParameterSet& parameterset)
00059 {
00060   fei::Factory::parameters(parameterset);
00061 
00062   parameterset.getIntParamValue("outputLevel", outputLevel_);
00063 
00064   bool blkGraph = false;
00065   bool blkMatrix = false;
00066 
00067   parameterset.getBoolParamValue("BLOCK_GRAPH", blkGraph);
00068   parameterset.getBoolParamValue("BLOCK_MATRIX", blkMatrix);
00069 
00070   blockEntryMatrix_ = (blkGraph || blkMatrix);
00071 
00072   const fei::Param* param = parameterset.get("Trilinos_Solver");
00073   if (param != 0) {
00074     if (param->getType() == fei::Param::STRING) {
00075       std::string strval = param->getStringValue();
00076       std::string::size_type ii = strval.find("Amesos");
00077       if (ii != std::string::npos) {
00078         useAmesos_ = true;
00079       }
00080     }
00081   }
00082 
00083   parameterset.getBoolParamValue("LPM_EpetraBasic", use_lpm_epetrabasic_);
00084   if (use_lpm_epetrabasic_) {
00085     create_LinProbMgr();
00086   }
00087 
00088   parameterset.getBoolParamValue("USE_FEI_MATRIX_LOCAL", use_feiMatrixLocal_);
00089 
00090   parameterset.getBoolParamValue("ORDER_ROWS_WITH_LOCAL_COLS_FIRST",
00091                                  orderRowsWithLocalColsFirst_);
00092 }
00093 
00094 fei::SharedPtr<fei::MatrixGraph>
00095 Factory_Trilinos::createMatrixGraph(fei::SharedPtr<fei::VectorSpace> rowSpace,
00096                                    fei::SharedPtr<fei::VectorSpace> colSpace,
00097                                    const char* name)
00098 {
00099   static fei::MatrixGraph_Impl2::Factory factory2;
00100   if (!use_lpm_epetrabasic_) {
00101     return(factory2.createMatrixGraph(rowSpace, colSpace, name));
00102   }
00103 
00104   return(factory2.createMatrixGraph(rowSpace, colSpace, name));
00105 }
00106 
00107 fei::SharedPtr<fei::Vector>
00108 Factory_Trilinos::wrapVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
00109                              fei::SharedPtr<Epetra_MultiVector> multiVec)
00110 {
00111   fei::SharedPtr<fei::Vector> vecptr;
00112   if (vecSpace->getNumIndices_Owned() != multiVec->Map().NumMyPoints()) {
00113     //vecSpace isn't compatible with multiVec's map, so return empty vecptr
00114     return(vecptr);
00115   }
00116 
00117   vecptr.reset(new fei::Vector_Impl<Epetra_MultiVector>(vecSpace,
00118                                                        multiVec.get(),
00119                                      multiVec->Map().NumMyPoints(), false));
00120   return(vecptr);
00121 }
00122 
00123 fei::SharedPtr<fei::Vector>
00124 Factory_Trilinos::wrapVector(fei::SharedPtr<fei::MatrixGraph> matGraph,
00125                              fei::SharedPtr<Epetra_MultiVector> multiVec)
00126 {
00127   fei::SharedPtr<fei::Vector> vecptr;
00128   if (matGraph->getRowSpace()->getNumIndices_Owned() !=
00129       multiVec->Map().NumMyPoints()) {
00130     //vector-space isn't compatible with multiVec's map, so return empty vecptr
00131     return(vecptr);
00132   }
00133 
00134   vecptr.reset(new fei::Vector_Impl<Epetra_MultiVector>(matGraph->getRowSpace(),
00135                                                        multiVec.get(),
00136                                         multiVec->Map().NumMyPoints(), false));
00137   return(vecptr);
00138 }
00139 
00140 fei::SharedPtr<fei::Vector>
00141 Factory_Trilinos::createVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
00142                                bool isSolutionVector,
00143                                int numVectors)
00144 {
00145   if (use_feiMatrixLocal_) {
00146     fei::SharedPtr<fei::Vector> feivec(new fei::Vector_Local(vecSpace));
00147     return(feivec);
00148   }
00149 
00150 //  if (blockEntryMatrix_) {
00151 //    throw std::runtime_error("Factory_Trilinos: fei ERROR, block-entry matrices/vectors not currently supported.");
00152 //  }
00153 
00154   std::vector<int> indices;
00155   int err = 0, localSize;
00156   if (reducer_.get() != NULL) {
00157     indices = reducer_->getLocalReducedEqns();
00158     localSize = indices.size();
00159   }
00160   else {
00161     if (blockEntryMatrix_) {
00162       localSize = vecSpace->getNumBlkIndices_Owned();
00163       indices.resize(localSize*2);
00164       err = vecSpace->getBlkIndices_Owned(localSize, &indices[0], &indices[localSize], localSize);
00165     }
00166     else {
00167       localSize = vecSpace->getNumIndices_Owned();
00168       indices.resize(localSize);
00169       err = vecSpace->getIndices_Owned(indices);
00170     }
00171   }
00172   if (err != 0) {
00173     throw std::runtime_error("error in vecSpace->getIndices_Owned");
00174   }
00175 
00176   fei::SharedPtr<fei::Vector> feivec, tmpvec;
00177   if (!use_lpm_epetrabasic_) {
00178     try {
00179       Epetra_BlockMap emap = blockEntryMatrix_ ? 
00180         Trilinos_Helpers::create_Epetra_BlockMap(vecSpace) :
00181         Trilinos_Helpers::create_Epetra_Map(comm_, indices);
00182 
00183       Epetra_MultiVector* emvec = new Epetra_MultiVector(emap, numVectors);
00184 
00185       tmpvec.reset(new fei::Vector_Impl<Epetra_MultiVector>(vecSpace,
00186                                               emvec, localSize,
00187                                         isSolutionVector, true));
00188     }
00189     catch(std::runtime_error& exc) {
00190       fei::console_out() << "Factory_Trilinos::createVector: caught exception '"
00191                << exc.what() << "', re-throwing..." << FEI_ENDL;
00192       throw exc;
00193     }
00194   }
00195   else {
00196     tmpvec.reset(new fei::Vector_Impl<fei::LinearProblemManager>(vecSpace,
00197                                     lpm_epetrabasic_.get(), localSize,
00198                                         isSolutionVector));
00199   }
00200 
00201   if (reducer_.get() != NULL) {
00202     feivec.reset(new fei::VectorReducer(reducer_,
00203                                         tmpvec, isSolutionVector));
00204   }
00205   else {
00206     feivec = tmpvec;
00207   }
00208 
00209   return(feivec);
00210 }
00211 
00212 fei::SharedPtr<fei::Vector>
00213 Factory_Trilinos::createVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
00214                                int numVectors)
00215 {
00216   bool isSolnVector = false;
00217   return(createVector(vecSpace, isSolnVector, numVectors));
00218 }
00219 
00220 fei::SharedPtr<fei::Vector>
00221 Factory_Trilinos::createVector(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00222                               int numVectors)
00223 {
00224   bool isSolnVector = false;
00225   return(createVector(matrixGraph, isSolnVector, numVectors));
00226 }
00227 
00228 fei::SharedPtr<fei::Vector>
00229 Factory_Trilinos::createVector(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00230                                bool isSolutionVector,
00231                                int numVectors)
00232 {
00233   if (use_feiMatrixLocal_) {
00234     fei::SharedPtr<fei::Vector> feivec(new fei::Vector_Local(matrixGraph->getRowSpace()));
00235     return(feivec);
00236   }
00237 
00238   int globalNumSlaves = matrixGraph->getGlobalNumSlaveConstraints();
00239 
00240   if (globalNumSlaves > 0 && reducer_.get()==NULL) {
00241     reducer_ = matrixGraph->getReducer();
00242   }
00243 
00244   fei::SharedPtr<fei::Vector> feivec, tmpvec;
00245 
00246   std::vector<int> indices;
00247   int err = 0, localSize;
00248   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
00249   if (reducer_.get() != NULL) {
00250     indices = reducer_->getLocalReducedEqns();
00251     localSize = indices.size();
00252   }
00253   else {
00254     localSize = vecSpace->getNumIndices_Owned();
00255     indices.resize(localSize);
00256     err = vecSpace->getIndices_Owned(indices);
00257   }
00258   if (err != 0) {
00259     throw std::runtime_error("error in vecSpace->getIndices_Owned");
00260   }
00261 
00262   if (!use_lpm_epetrabasic_) {
00263 #ifdef HAVE_FEI_EPETRA
00264     try {
00265       Epetra_BlockMap emap = blockEntryMatrix_ ?
00266         Trilinos_Helpers::create_Epetra_BlockMap(vecSpace) :
00267         Trilinos_Helpers::create_Epetra_Map(comm_, indices);
00268 
00269       Epetra_MultiVector* emvec = new Epetra_MultiVector(emap, numVectors);
00270 
00271       tmpvec.reset(new fei::Vector_Impl<Epetra_MultiVector>(matrixGraph->getRowSpace(), emvec,
00272                                                   localSize, isSolutionVector, true));
00273     }
00274     catch(std::runtime_error& exc) {
00275       fei::console_out() << "Factory_Trilinos::createVector: caught exception '"
00276                << exc.what() << "', re-throwing..." << FEI_ENDL;
00277       throw exc;
00278     }
00279 #else
00280     fei::console_out() << "fei_Factory_Trilinos::createVector ERROR, HAVE_FEI_EPETRA not defined."
00281       << FEI_ENDL;
00282 #endif
00283   }
00284   else {
00285 #ifdef HAVE_FEI_EPETRA
00286     vecSpace = matrixGraph->getRowSpace();
00287 
00288     lpm_epetrabasic_->setRowDistribution(indices);
00289     tmpvec.reset(new fei::Vector_Impl<fei::LinearProblemManager>(vecSpace,
00290                                    lpm_epetrabasic_.get(), localSize, isSolutionVector));
00291 #endif
00292   }
00293 
00294   if (reducer_.get() != NULL) {
00295     feivec.reset(new fei::VectorReducer(reducer_, tmpvec, isSolutionVector));
00296   }
00297   else {
00298     feivec = tmpvec;
00299   }
00300 
00301   return(feivec);
00302 }
00303 
00304 fei::SharedPtr<fei::Matrix>
00305 Factory_Trilinos::createMatrix(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00306 {
00307   fei::SharedPtr<fei::Matrix> feimat;
00308 #ifdef HAVE_FEI_EPETRA
00309   int globalNumSlaves = matrixGraph->getGlobalNumSlaveConstraints();
00310 
00311   if (globalNumSlaves > 0 && reducer_.get()==NULL) {
00312     reducer_ = matrixGraph->getReducer();
00313   }
00314 
00315   if (use_lpm_epetrabasic_) {
00316     return(
00317       Trilinos_Helpers::create_from_LPM_EpetraBasic(matrixGraph,
00318                                                     blockEntryMatrix_,
00319                                                     reducer_,
00320                                                     lpm_epetrabasic_)
00321     );
00322   }
00323 
00324   if (use_feiMatrixLocal_) {
00325     return(fei::Matrix_Local::create_Matrix_Local(matrixGraph, blockEntryMatrix_));
00326   }
00327 
00328   return(
00329       Trilinos_Helpers::create_from_Epetra_Matrix(matrixGraph, blockEntryMatrix_,
00330                                                   reducer_, orderRowsWithLocalColsFirst_)
00331   );
00332 #else
00333   fei::console_out() << "fei_Factory_Trilinos::createMatrix ERROR, HAVE_FEI_EPETRA "
00334      << "not defined."<<FEI_ENDL;
00335   return feimat;
00336 #endif
00337 }
00338 
00339 fei::SharedPtr<fei::Solver>
00340 Factory_Trilinos::createSolver(const char* name)
00341 {
00342   fei::SharedPtr<fei::Solver> solver;
00343 
00344   if (useAmesos_ || name != 0) {
00345     if (name != 0) {
00346       std::string sname(name);
00347       std::string::size_type ii = sname.find("Amesos");
00348       if (ii == std::string::npos) {
00349 #ifdef HAVE_FEI_AZTECOO
00350         solver.reset(new Solver_AztecOO);
00351         return(solver);
00352 #else
00353         fei::console_out() << "fei_Factory_Trilinos::createSolver: ERROR, AztecOO not "
00354            << "available." << FEI_ENDL; 
00355         return(solver);
00356 #endif
00357       }
00358     }
00359 #ifdef HAVE_FEI_AMESOS
00360     solver.reset(new Solver_Amesos);
00361 #else
00362     fei::console_out() << "fei_Factory_Trilinos::createSolver: ERROR, Amesos requested,"
00363       << " but HAVE_FEI_AMESOS is not defined so Amesos is not available."
00364       <<FEI_ENDL;
00365     return(solver);
00366 #endif
00367   }
00368   else {
00369 #ifdef HAVE_FEI_AZTECOO
00370     solver.reset(new Solver_AztecOO);
00371 #endif
00372   }
00373 
00374   return(solver);
00375 }
00376 
00377 void Factory_Trilinos::create_LinProbMgr(bool replace_if_already_created)
00378 {
00379   if (!use_lpm_epetrabasic_) {
00380     return;
00381   }
00382 
00383   bool need_to_create_lpm = false;
00384 
00385   if (lpm_epetrabasic_.get() == NULL) {
00386     need_to_create_lpm = true;
00387   }
00388 
00389   if (replace_if_already_created) {
00390     need_to_create_lpm = true;
00391   }
00392 
00393   if (need_to_create_lpm) {
00394 #ifdef HAVE_FEI_EPETRA
00395     fei::SharedPtr<fei::LinearProblemManager>
00396       newlpm(new LinProbMgr_EpetraBasic(comm_));
00397 
00398     lpm_epetrabasic_ = newlpm;
00399 #else
00400     fei::console_out() << "fei_Factory_Trilinos::create_LinProbMgr ERROR, HAVE_FEI_EPETRA"
00401        <<" not defined."<<FEI_ENDL;
00402 #endif
00403   }
00404 }
00405 
00406 //HAVE_FEI_EPETRA
00407 #endif
00408 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends