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