fei_Solver_Amesos.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 #ifdef HAVE_FEI_AMESOS
00012 #ifdef HAVE_FEI_EPETRA
00013 
00014 #include <fei_Solver_Amesos.hpp>
00015 #include <fei_ParameterSet.hpp>
00016 #include <snl_fei_Utils.hpp>
00017 #include <fei_utils.hpp>
00018 
00019 //fei_Include_Trilinos.hpp includes the actual Trilinos headers (epetra, aztecoo, ml ...)
00020 #include <fei_Include_Trilinos.hpp>
00021 #include <fei_Trilinos_Helpers.hpp>
00022 #include <fei_VectorTraits_Epetra.hpp>
00023 #include <fei_MatrixTraits_Epetra.hpp>
00024 
00025 #include <fei_Vector.hpp>
00026 #include <fei_Matrix.hpp>
00027 #include <fei_LinearSystem.hpp>
00028 
00029 //---------------------------------------------------------------------------
00030 Solver_Amesos::Solver_Amesos()
00031   : tolerance_(1.e-6),
00032     maxIters_(500),
00033     amesos_factory_(new Amesos()),
00034     amesos_solver_(NULL),
00035     epetra_linearproblem_(NULL)
00036 {
00037   paramlist_ = new Teuchos::ParameterList;
00038 }
00039 
00040 //---------------------------------------------------------------------------
00041 Solver_Amesos::~Solver_Amesos()
00042 {
00043   delete amesos_solver_;
00044   delete amesos_factory_;
00045   delete epetra_linearproblem_;
00046   delete paramlist_;
00047 }
00048 
00049 //---------------------------------------------------------------------------
00050 Teuchos::ParameterList& Solver_Amesos::get_ParameterList()
00051 {
00052   if (paramlist_ == NULL) {
00053     paramlist_ = new Teuchos::ParameterList;
00054   }
00055 
00056   return( *paramlist_ );
00057 }
00058 
00059 //---------------------------------------------------------------------------
00060 int Solver_Amesos::solve(fei::LinearSystem* linearSystem,
00061         fei::Matrix* preconditioningMatrix,
00062         const fei::ParameterSet& parameterSet,
00063         int& iterationsTaken,
00064         int& status)
00065 {
00066   Trilinos_Helpers::copy_parameterset(parameterSet, *paramlist_);
00067 
00068   int numParams = 0;
00069   const char** paramStrings = NULL;
00070   std::vector<std::string> stdstrings;
00071   fei::utils::convert_ParameterSet_to_strings(&parameterSet, stdstrings);
00072   fei::utils::strings_to_char_ptrs(stdstrings, numParams, paramStrings);
00073 
00074   int err = solve(linearSystem, preconditioningMatrix, numParams, paramStrings,
00075       iterationsTaken, status);
00076 
00077   delete [] paramStrings;
00078   
00079   return(err);
00080 }
00081 
00082 //---------------------------------------------------------------------------
00083 int Solver_Amesos::solve(fei::LinearSystem* linearSystem,
00084         fei::Matrix* preconditioningMatrix,
00085         int numParams,
00086         const char* const* solverParams,
00087         int& iterationsTaken,
00088         int& status)
00089 {
00090   Epetra_MultiVector*    x = NULL;
00091   Epetra_MultiVector*    b = NULL;
00092   Epetra_CrsMatrix* A = NULL;
00093   Epetra_Operator* opA = NULL;
00094 
00095   fei::SharedPtr<fei::Matrix> feiA = linearSystem->getMatrix();
00096   fei::SharedPtr<fei::Vector> feix = linearSystem->getSolutionVector();
00097   fei::SharedPtr<fei::Vector> feib = linearSystem->getRHS();
00098 
00099   Trilinos_Helpers::get_Epetra_pointers(feiA, feix, feib,
00100                                         A, opA, x, b);
00101 
00102   if (opA == 0 || x == 0 || b == 0) {
00103     FEI_CERR << "Error, couldn't obtain Epetra objects from "
00104       << "fei container-objects."<<FEI_ENDL;
00105     return(-1);
00106   }
00107 
00108   if (epetra_linearproblem_ == NULL) {
00109     epetra_linearproblem_ = new Epetra_LinearProblem;
00110   }
00111 
00112   epetra_linearproblem_->SetOperator(A);
00113   epetra_linearproblem_->SetLHS(x);
00114   epetra_linearproblem_->SetRHS(b);
00115 
00116   const char* param = snl_fei::getParamValue("Trilinos_Solver",
00117                numParams, solverParams);
00118   if (param != NULL) {
00119     if (amesos_solver_ == 0) {
00120       amesos_solver_ = amesos_factory_->Create(param, *epetra_linearproblem_);
00121     }
00122     if (amesos_solver_ == 0) {
00123       cerr << "Solver_Amesos::solve ERROR, couldn't create Amesos solver named "
00124      << param << ", amesos_factory::Create returned NULL." << endl;
00125       status = -999;
00126       return(-1);
00127     }
00128   }
00129   else {
00130     static char amesosklu[] = "Amesos_Klu";
00131     if (amesos_solver_ == 0) {
00132       amesos_solver_ = amesos_factory_->Create( amesosklu,
00133             *epetra_linearproblem_);
00134     }
00135     if (amesos_solver_ == 0) {
00136       cerr << "Solver_Amesos::solve ERROR, couldn't create Amesos solver named "
00137      << amesosklu << ", it's apparently not supported." << endl;
00138       status = -999;
00139       return(-1);
00140     }
00141   }
00142 
00143   amesos_solver_->SetParameters(*paramlist_);
00144 
00145   if (feiA->changedSinceMark()) {
00146     amesos_solver_->SymbolicFactorization();
00147     amesos_solver_->NumericFactorization();
00148     feiA->markState();
00149   }
00150 
00151   amesos_solver_->Solve();
00152   status = 0;
00153 
00154   return(0);
00155 }
00156 
00157 //---------------------------------------------------------------------------
00158 int Solver_Amesos::parseParameters(int numParams,
00159            const char*const* params)
00160 {
00161  
00162   return(0);
00163 }
00164 
00165 #endif //HAVE_FEI_EPETRA
00166 #endif //HAVE_FEI_AMESOS

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