FEI Version of the Day
fei_Factory_Trilinos.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #include "fei_trilinos_macros.hpp"
00045 
00046 #include <fei_Factory_Trilinos.hpp>
00047 
00048 #include <fei_VectorReducer.hpp>
00049 #include <fei_MatrixReducer.hpp>
00050 #include <fei_Matrix_Local.hpp>
00051 #include <fei_Vector_Local.hpp>
00052 
00053 #ifdef HAVE_FEI_AZTECOO
00054 #include <fei_Solver_AztecOO.hpp>
00055 #endif
00056 #ifdef HAVE_FEI_BELOS
00057 #include <fei_Solver_Belos.hpp>
00058 #endif
00059 #ifdef HAVE_FEI_AMESOS
00060 #include <fei_Solver_Amesos.hpp>
00061 #endif
00062 
00063 #ifdef HAVE_FEI_EPETRA
00064 
00065 Factory_Trilinos::Factory_Trilinos(MPI_Comm comm)
00066   : fei::Factory(comm),
00067     comm_(comm),
00068     reducer_(),
00069     lpm_epetrabasic_(),
00070     use_lpm_epetrabasic_(false),
00071     useAmesos_(false),
00072     useBelos_(false),
00073     use_feiMatrixLocal_(false),
00074     blockEntryMatrix_(false),
00075     orderRowsWithLocalColsFirst_(false),
00076     outputLevel_(0)
00077 {
00078 }
00079 
00080 Factory_Trilinos::~Factory_Trilinos()
00081 {
00082 }
00083 
00084 int Factory_Trilinos::parameters(int numParams,
00085                                   const char* const* paramStrings)
00086 {
00087   std::vector<std::string> stdstrings;
00088   fei::utils::char_ptrs_to_strings(numParams, paramStrings, stdstrings);
00089 
00090   fei::ParameterSet paramset;
00091   fei::utils::parse_strings(stdstrings, " ", paramset);
00092 
00093   parameters(paramset);
00094   return(0);
00095 }
00096 
00097 void Factory_Trilinos::parameters(const fei::ParameterSet& parameterset)
00098 {
00099   fei::Factory::parameters(parameterset);
00100 
00101   parameterset.getIntParamValue("outputLevel", outputLevel_);
00102 
00103   bool blkGraph = false;
00104   bool blkMatrix = false;
00105 
00106   parameterset.getBoolParamValue("BLOCK_GRAPH", blkGraph);
00107   parameterset.getBoolParamValue("BLOCK_MATRIX", blkMatrix);
00108 
00109   blockEntryMatrix_ = (blkGraph || blkMatrix);
00110 
00111   const fei::Param* param = parameterset.get("Trilinos_Solver");
00112   if (param != 0) {
00113     if (param->getType() == fei::Param::STRING) {
00114       std::string strval = param->getStringValue();
00115       std::string::size_type ii = strval.find("Amesos");
00116       if (ii != std::string::npos) {
00117         useAmesos_ = true;
00118       }
00119       ii = strval.find("Belos");
00120       if (ii != std::string::npos) {
00121         useBelos_ = true;
00122       }
00123     }
00124   }
00125 
00126   parameterset.getBoolParamValue("LPM_EpetraBasic", use_lpm_epetrabasic_);
00127   if (use_lpm_epetrabasic_) {
00128     create_LinProbMgr();
00129   }
00130 
00131   parameterset.getBoolParamValue("USE_FEI_MATRIX_LOCAL", use_feiMatrixLocal_);
00132 
00133   parameterset.getBoolParamValue("ORDER_ROWS_WITH_LOCAL_COLS_FIRST",
00134                                  orderRowsWithLocalColsFirst_);
00135 }
00136 
00137 fei::SharedPtr<fei::MatrixGraph>
00138 Factory_Trilinos::createMatrixGraph(fei::SharedPtr<fei::VectorSpace> rowSpace,
00139                                    fei::SharedPtr<fei::VectorSpace> colSpace,
00140                                    const char* name)
00141 {
00142   static fei::MatrixGraph_Impl2::Factory factory2;
00143   if (!use_lpm_epetrabasic_) {
00144     return(factory2.createMatrixGraph(rowSpace, colSpace, name));
00145   }
00146 
00147   return(factory2.createMatrixGraph(rowSpace, colSpace, name));
00148 }
00149 
00150 fei::SharedPtr<fei::Vector>
00151 Factory_Trilinos::wrapVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
00152                              fei::SharedPtr<Epetra_MultiVector> multiVec)
00153 {
00154   fei::SharedPtr<fei::Vector> vecptr;
00155   if (vecSpace->getNumIndices_Owned() != multiVec->Map().NumMyPoints()) {
00156     //vecSpace isn't compatible with multiVec's map, so return empty vecptr
00157     return(vecptr);
00158   }
00159 
00160   vecptr.reset(new fei::Vector_Impl<Epetra_MultiVector>(vecSpace,
00161                                                        multiVec.get(),
00162                                      multiVec->Map().NumMyPoints(), false));
00163   return(vecptr);
00164 }
00165 
00166 fei::SharedPtr<fei::Vector>
00167 Factory_Trilinos::wrapVector(fei::SharedPtr<fei::MatrixGraph> matGraph,
00168                              fei::SharedPtr<Epetra_MultiVector> multiVec)
00169 {
00170   fei::SharedPtr<fei::Vector> vecptr;
00171   if (matGraph->getRowSpace()->getNumIndices_Owned() !=
00172       multiVec->Map().NumMyPoints()) {
00173     //vector-space isn't compatible with multiVec's map, so return empty vecptr
00174     return(vecptr);
00175   }
00176 
00177   vecptr.reset(new fei::Vector_Impl<Epetra_MultiVector>(matGraph->getRowSpace(),
00178                                                        multiVec.get(),
00179                                         multiVec->Map().NumMyPoints(), false));
00180   return(vecptr);
00181 }
00182 
00183 fei::SharedPtr<fei::Vector>
00184 Factory_Trilinos::createVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
00185                                bool isSolutionVector,
00186                                int numVectors)
00187 {
00188   if (use_feiMatrixLocal_) {
00189     fei::SharedPtr<fei::Vector> feivec(new fei::Vector_Local(vecSpace));
00190     return(feivec);
00191   }
00192 
00193 //  if (blockEntryMatrix_) {
00194 //    throw std::runtime_error("Factory_Trilinos: fei ERROR, block-entry matrices/vectors not currently supported.");
00195 //  }
00196 
00197   std::vector<int> indices;
00198   int err = 0, localSize;
00199   if (reducer_.get() != NULL) {
00200     indices = reducer_->getLocalReducedEqns();
00201     localSize = indices.size();
00202   }
00203   else {
00204     if (blockEntryMatrix_) {
00205       localSize = vecSpace->getNumBlkIndices_Owned();
00206       indices.resize(localSize*2);
00207       err = vecSpace->getBlkIndices_Owned(localSize, &indices[0], &indices[localSize], localSize);
00208     }
00209     else {
00210       localSize = vecSpace->getNumIndices_Owned();
00211       indices.resize(localSize);
00212       err = vecSpace->getIndices_Owned(indices);
00213     }
00214   }
00215   if (err != 0) {
00216     throw std::runtime_error("error in vecSpace->getIndices_Owned");
00217   }
00218 
00219   fei::SharedPtr<fei::Vector> feivec, tmpvec;
00220   if (!use_lpm_epetrabasic_) {
00221     try {
00222       Epetra_BlockMap emap = blockEntryMatrix_ ? 
00223         Trilinos_Helpers::create_Epetra_BlockMap(vecSpace) :
00224         Trilinos_Helpers::create_Epetra_Map(comm_, indices);
00225 
00226       Epetra_MultiVector* emvec = new Epetra_MultiVector(emap, numVectors);
00227 
00228       tmpvec.reset(new fei::Vector_Impl<Epetra_MultiVector>(vecSpace,
00229                                               emvec, localSize,
00230                                         isSolutionVector, true));
00231     }
00232     catch(std::runtime_error& exc) {
00233       fei::console_out() << "Factory_Trilinos::createVector: caught exception '"
00234                << exc.what() << "', re-throwing..." << FEI_ENDL;
00235       throw exc;
00236     }
00237   }
00238   else {
00239     tmpvec.reset(new fei::Vector_Impl<fei::LinearProblemManager>(vecSpace,
00240                                     lpm_epetrabasic_.get(), localSize,
00241                                         isSolutionVector));
00242   }
00243 
00244   if (reducer_.get() != NULL) {
00245     feivec.reset(new fei::VectorReducer(reducer_,
00246                                         tmpvec, isSolutionVector));
00247   }
00248   else {
00249     feivec = tmpvec;
00250   }
00251 
00252   return(feivec);
00253 }
00254 
00255 fei::SharedPtr<fei::Vector>
00256 Factory_Trilinos::createVector(fei::SharedPtr<fei::VectorSpace> vecSpace,
00257                                int numVectors)
00258 {
00259   bool isSolnVector = false;
00260   return(createVector(vecSpace, isSolnVector, numVectors));
00261 }
00262 
00263 fei::SharedPtr<fei::Vector>
00264 Factory_Trilinos::createVector(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00265                               int numVectors)
00266 {
00267   bool isSolnVector = false;
00268   return(createVector(matrixGraph, isSolnVector, numVectors));
00269 }
00270 
00271 fei::SharedPtr<fei::Vector>
00272 Factory_Trilinos::createVector(fei::SharedPtr<fei::MatrixGraph> matrixGraph,
00273                                bool isSolutionVector,
00274                                int numVectors)
00275 {
00276   if (use_feiMatrixLocal_) {
00277     fei::SharedPtr<fei::Vector> feivec(new fei::Vector_Local(matrixGraph->getRowSpace()));
00278     return(feivec);
00279   }
00280 
00281   int globalNumSlaves = matrixGraph->getGlobalNumSlaveConstraints();
00282 
00283   if (globalNumSlaves > 0 && reducer_.get()==NULL) {
00284     reducer_ = matrixGraph->getReducer();
00285   }
00286 
00287   fei::SharedPtr<fei::Vector> feivec, tmpvec;
00288 
00289   std::vector<int> indices;
00290   int err = 0, localSize;
00291   fei::SharedPtr<fei::VectorSpace> vecSpace = matrixGraph->getRowSpace();
00292   if (reducer_.get() != NULL) {
00293     indices = reducer_->getLocalReducedEqns();
00294     localSize = indices.size();
00295   }
00296   else {
00297     localSize = vecSpace->getNumIndices_Owned();
00298     indices.resize(localSize);
00299     err = vecSpace->getIndices_Owned(indices);
00300   }
00301   if (err != 0) {
00302     throw std::runtime_error("error in vecSpace->getIndices_Owned");
00303   }
00304 
00305   if (!use_lpm_epetrabasic_) {
00306 #ifdef HAVE_FEI_EPETRA
00307     try {
00308       Epetra_BlockMap emap = blockEntryMatrix_ ?
00309         Trilinos_Helpers::create_Epetra_BlockMap(vecSpace) :
00310         Trilinos_Helpers::create_Epetra_Map(comm_, indices);
00311 
00312       Epetra_MultiVector* emvec = new Epetra_MultiVector(emap, numVectors);
00313 
00314       tmpvec.reset(new fei::Vector_Impl<Epetra_MultiVector>(matrixGraph->getRowSpace(), emvec,
00315                                                   localSize, isSolutionVector, true));
00316     }
00317     catch(std::runtime_error& exc) {
00318       fei::console_out() << "Factory_Trilinos::createVector: caught exception '"
00319                << exc.what() << "', re-throwing..." << FEI_ENDL;
00320       throw exc;
00321     }
00322 #else
00323     fei::console_out() << "fei_Factory_Trilinos::createVector ERROR, HAVE_FEI_EPETRA not defined."
00324       << FEI_ENDL;
00325 #endif
00326   }
00327   else {
00328 #ifdef HAVE_FEI_EPETRA
00329     vecSpace = matrixGraph->getRowSpace();
00330 
00331     lpm_epetrabasic_->setRowDistribution(indices);
00332     tmpvec.reset(new fei::Vector_Impl<fei::LinearProblemManager>(vecSpace,
00333                                    lpm_epetrabasic_.get(), localSize, isSolutionVector));
00334 #endif
00335   }
00336 
00337   if (reducer_.get() != NULL) {
00338     feivec.reset(new fei::VectorReducer(reducer_, tmpvec, isSolutionVector));
00339   }
00340   else {
00341     feivec = tmpvec;
00342   }
00343 
00344   return(feivec);
00345 }
00346 
00347 fei::SharedPtr<fei::Matrix>
00348 Factory_Trilinos::createMatrix(fei::SharedPtr<fei::MatrixGraph> matrixGraph)
00349 {
00350   fei::SharedPtr<fei::Matrix> feimat;
00351 #ifdef HAVE_FEI_EPETRA
00352   int globalNumSlaves = matrixGraph->getGlobalNumSlaveConstraints();
00353 
00354   if (globalNumSlaves > 0 && reducer_.get()==NULL) {
00355     reducer_ = matrixGraph->getReducer();
00356   }
00357 
00358   if (use_lpm_epetrabasic_) {
00359     return(
00360       Trilinos_Helpers::create_from_LPM_EpetraBasic(matrixGraph,
00361                                                     blockEntryMatrix_,
00362                                                     reducer_,
00363                                                     lpm_epetrabasic_)
00364     );
00365   }
00366 
00367   if (use_feiMatrixLocal_) {
00368     return(fei::Matrix_Local::create_Matrix_Local(matrixGraph, blockEntryMatrix_));
00369   }
00370 
00371   return(
00372       Trilinos_Helpers::create_from_Epetra_Matrix(matrixGraph, blockEntryMatrix_,
00373                                                   reducer_, orderRowsWithLocalColsFirst_)
00374   );
00375 #else
00376   fei::console_out() << "fei_Factory_Trilinos::createMatrix ERROR, HAVE_FEI_EPETRA "
00377      << "not defined."<<FEI_ENDL;
00378   return feimat;
00379 #endif
00380 }
00381 
00382 fei::SharedPtr<fei::Solver>
00383 Factory_Trilinos::createSolver(const char* name)
00384 {
00385   fei::SharedPtr<fei::Solver> solver;
00386 
00387   std::string::size_type ii_amesos = std::string::npos;
00388   std::string::size_type ii_belos = std::string::npos;
00389 
00390   if (name != 0) {
00391     std::string sname(name);
00392     ii_amesos = sname.find("Amesos");
00393     ii_belos = sname.find("Belos");
00394   }
00395 
00396   if (useAmesos_ || ii_amesos != std::string::npos) {
00397 #ifdef HAVE_FEI_AMESOS
00398     solver.reset(new Solver_Amesos);
00399     return solver;
00400 #else
00401     fei::console_out() << "fei_Factory_Trilinos::createSolver: ERROR, Amesos not available (not enabled at compile-time?)"<<FEI_ENDL;
00402 #endif
00403   }
00404 
00405   if (useBelos_ || ii_belos != std::string::npos) {
00406 #ifdef HAVE_FEI_BELOS
00407     solver.reset(new Solver_Belos);
00408     return solver;
00409 #else
00410     fei::console_out() << "fei_Factory_Trilinos::createSolver: ERROR, Belos not available (not enabled at compile-time?)"<<FEI_ENDL;
00411 #endif
00412   }
00413 
00414 #ifdef HAVE_FEI_AZTECOO
00415   solver.reset(new Solver_AztecOO);
00416   return solver;
00417 #else
00418   fei::console_out() << "fei_Factory_Trilinos::createSolver: ERROR, AztecOO not "
00419          << "available." << FEI_ENDL; 
00420   return solver;
00421 #endif
00422 }
00423 
00424 void Factory_Trilinos::create_LinProbMgr(bool replace_if_already_created)
00425 {
00426   if (!use_lpm_epetrabasic_) {
00427     return;
00428   }
00429 
00430   bool need_to_create_lpm = false;
00431 
00432   if (lpm_epetrabasic_.get() == NULL) {
00433     need_to_create_lpm = true;
00434   }
00435 
00436   if (replace_if_already_created) {
00437     need_to_create_lpm = true;
00438   }
00439 
00440   if (need_to_create_lpm) {
00441 #ifdef HAVE_FEI_EPETRA
00442     fei::SharedPtr<fei::LinearProblemManager>
00443       newlpm(new LinProbMgr_EpetraBasic(comm_));
00444 
00445     lpm_epetrabasic_ = newlpm;
00446 #else
00447     fei::console_out() << "fei_Factory_Trilinos::create_LinProbMgr ERROR, HAVE_FEI_EPETRA"
00448        <<" not defined."<<FEI_ENDL;
00449 #endif
00450   }
00451 }
00452 
00453 //HAVE_FEI_EPETRA
00454 #endif
00455 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends