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