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

Generated on Wed May 12 21:30:41 2010 for FEI by  doxygen 1.4.7