FEI Version of the Day
fei_Solver_Belos.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 #ifdef HAVE_FEI_BELOS
00047 
00048 //fei_Include_Trilinos.h includes the Trilinos headers (epetra, belos, ...)
00049 #include <fei_Include_Trilinos.hpp>
00050 #include <fei_Trilinos_Helpers.hpp>
00051 
00052 #include <fei_Solver_Belos.hpp>
00053 #include <fei_LinProbMgr_EpetraBasic.hpp>
00054 #include <fei_ParameterSet.hpp>
00055 #include <fei_utils.hpp>
00056 
00057 #include <fei_VectorTraits_Epetra.hpp>
00058 #include <fei_Vector_Impl.hpp>
00059 #include <fei_MatrixTraits_Epetra.hpp>
00060 #include <fei_Matrix_Impl.hpp>
00061 #include <fei_LinearSystem.hpp>
00062 
00063 //---------------------------------------------------------------------------
00064 Solver_Belos::Solver_Belos()
00065   : tolerance_(1.e-6), maxIters_(500), useTranspose_(false), paramlist_(),
00066     useML_(false),
00067 #ifdef HAVE_FEI_ML
00068    ml_prec_(NULL), ml_defaults_set_(false),
00069    ml_aztec_options_(NULL), ml_aztec_params_(NULL),
00070 #endif
00071   name_(),
00072   dbgprefix_("SlvBelos: ")
00073 {
00074   paramlist_ = Teuchos::rcp(new Teuchos::ParameterList);
00075 
00076 #ifdef HAVE_FEI_ML
00077   ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
00078   ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
00079 #endif
00080 }
00081 
00082 //---------------------------------------------------------------------------
00083 Solver_Belos::~Solver_Belos()
00084 {
00085 #ifdef HAVE_FEI_ML
00086   delete [] ml_aztec_options_;
00087   delete [] ml_aztec_params_;
00088   delete ml_prec_;
00089 #endif
00090 }
00091 
00092 //---------------------------------------------------------------------------
00093 void Solver_Belos::setUseML(bool useml)
00094 {
00095   useML_ = useml;
00096 }
00097 
00098 //---------------------------------------------------------------------------
00099 Teuchos::ParameterList& Solver_Belos::get_ParameterList()
00100 {
00101   if (paramlist_.get() == NULL) {
00102     paramlist_ = Teuchos::rcp(new Teuchos::ParameterList);
00103   }
00104 
00105   return( *paramlist_ );
00106 }
00107 
00108 //---------------------------------------------------------------------------
00109 int Solver_Belos::solve(fei::LinearSystem* linearSystem,
00110         fei::Matrix* preconditioningMatrix,
00111         const fei::ParameterSet& parameterSet,
00112         int& iterationsTaken,
00113         int& status)
00114 {
00115   std::string krylov_solver_name;
00116   parameterSet.getStringParamValue("krylov_solver", krylov_solver_name);
00117 
00118   Teuchos::RCP<Teuchos::ParameterList>& paramlist = paramlist_;
00119 
00120 #ifdef HAVE_FEI_ML
00121   if (ml_aztec_options_ == NULL)
00122     ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
00123   if (ml_aztec_params_ == NULL)
00124     ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
00125 
00126   if (!ml_defaults_set_ && useML_) {
00127     Teuchos::ParameterList mlparams;
00128     ML_Epetra::SetDefaults("SA", mlparams, ml_aztec_options_,ml_aztec_params_);
00129     mlparams.setParameters(*paramlist);
00130     *paramlist = mlparams;
00131     ml_defaults_set_ = true;
00132   }
00133 #endif
00134 
00135   Trilinos_Helpers::copy_parameterset(parameterSet, *paramlist);
00136 
00137   fei::SharedPtr<fei::Matrix> feiA = linearSystem->getMatrix();
00138   fei::SharedPtr<fei::Vector> feix = linearSystem->getSolutionVector();
00139   fei::SharedPtr<fei::Vector> feib = linearSystem->getRHS();
00140 
00141   Epetra_MultiVector*    x = NULL;
00142   Epetra_MultiVector*    b = NULL;
00143   Epetra_Operator* epetra_op = 0;
00144   Epetra_CrsMatrix* crsA = NULL;
00145 
00146   Trilinos_Helpers::get_Epetra_pointers(feiA, feix, feib,
00147                                         crsA, epetra_op, x, b);
00148 
00149   Teuchos::RCP<Epetra_CrsMatrix> rcp_A(crsA);
00150   Teuchos::RCP<Epetra_MultiVector> rcp_x(x);
00151   Teuchos::RCP<Epetra_MultiVector> rcp_b(b);
00152 
00153   if (epetra_op == 0 || x == 0 || b == 0) {
00154     fei::console_out() << "Solver_Belos::solve Error, couldn't obtain Epetra objects"
00155      << " from fei container-objects."<<FEI_ENDL;
00156     return(-1);
00157   }
00158 
00159   Epetra_RowMatrix* precond = NULL;
00160   if (preconditioningMatrix != NULL) {
00161     fei::Matrix_Impl<Epetra_CrsMatrix>* snl_epetra_crs =
00162       dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(preconditioningMatrix);
00163     fei::Matrix_Impl<Epetra_VbrMatrix>* snl_epetra_vbr =
00164       dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(preconditioningMatrix);
00165     if (snl_epetra_crs != NULL) {
00166       precond = snl_epetra_crs->getMatrix().get();
00167     }
00168     else if (snl_epetra_vbr != NULL) {
00169       precond = snl_epetra_vbr->getMatrix().get();
00170     }
00171     else {
00172       fei::console_out() << "Solver_Belos::solve: ERROR getting epetra row matrix"
00173          << " from preconditioningMatrix."<<FEI_ENDL;
00174       return(-1);
00175     }
00176   }
00177 
00178   if (precond != NULL) {
00179 //TODO: set up preconditioner for Belos here
00180   }
00181 
00182   bool needNewPreconditioner = false;
00183 
00184   if (feiA->changedSinceMark()) {
00185     feiA->markState();
00186     needNewPreconditioner = true;
00187   }
00188 
00189   if (needNewPreconditioner) {
00190 //
00191 //    if (useML_) {
00192 #ifdef HAVE_FEI_ML
00193 //      setup_ml_operator(*azoo_, crsA);
00194 #else
00195 //      fei::console_out() <<"Solver_Belos::solve ERROR, ML requested but HAVE_FEI_ML not defined."
00196 //         << FEI_ENDL;
00197 //      return(-1);
00198 #endif
00199 //    }
00200 //    else {
00201 //      azoo_->SetAztecOption(AZ_pre_calc, AZ_calc);
00202 //      azoo_->SetAztecOption(AZ_keep_info, 1);
00203 //    }
00204   }
00205   else {
00206 //    if (!useML_) {
00207 //      azoo_->SetAztecOption(AZ_pre_calc, AZ_reuse);
00208 //    }
00209   }
00210 
00211   epetra_op->SetUseTranspose(useTranspose_);
00212 
00213   Belos::SolverFactory<double,Epetra_MultiVector,Epetra_Operator> belos_factory;
00214   belos_solver_manager_ = belos_factory.create(krylov_solver_name, paramlist);
00215 
00216   Teuchos::RCP<Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator> > belos_lin_prob = Teuchos::rcp(new Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator>(rcp_A, rcp_x, rcp_b));
00217 
00218   belos_lin_prob->setProblem();
00219 
00220   belos_solver_manager_->setProblem(belos_lin_prob);
00221 
00222   belos_solver_manager_->solve();
00223   status = 0;
00224 
00225   iterationsTaken = belos_solver_manager_->getNumIters();
00226 
00227   rcp_A.release();
00228   rcp_x.release();
00229   rcp_b.release();
00230 
00231   int olevel = 0;
00232   parameterSet.getIntParamValue("outputLevel", olevel);
00233 
00234   std::string param2;
00235   parameterSet.getStringParamValue("FEI_OUTPUT_LEVEL", param2);
00236 
00237   if (olevel >= 3 || param2 == "MATRIX_FILES") {
00238     std::string param1;
00239     parameterSet.getStringParamValue("debugOutput", param1);
00240 
00241     FEI_OSTRINGSTREAM osstr;
00242     if (!param1.empty()) {
00243       osstr << param1 << "/";
00244     }
00245     else osstr << "./";
00246 
00247     osstr << "x_AztecOO.vec";
00248     feix->writeToFile(osstr.str().c_str());
00249   }
00250 
00251   return(0);
00252 }
00253 
00254 //---------------------------------------------------------------------------
00255 int Solver_Belos::solve(fei::LinearSystem* linearSystem,
00256         fei::Matrix* preconditioningMatrix,
00257         int numParams,
00258         const char* const* solverParams,
00259         int& iterationsTaken,
00260         int& status)
00261 {
00262   std::vector<std::string> stdstrings;
00263   fei::utils::char_ptrs_to_strings(numParams, solverParams, stdstrings);
00264 
00265   fei::ParameterSet paramset;
00266   fei::utils::parse_strings(stdstrings, " ", paramset);
00267 
00268   return( solve(linearSystem, preconditioningMatrix, paramset,
00269     iterationsTaken, status) );
00270 }
00271 
00272 //---------------------------------------------------------------------------
00273 int Solver_Belos::setup_ml_operator(AztecOO& azoo, Epetra_CrsMatrix* A)
00274 {
00275 #ifdef HAVE_FEI_ML
00276   if (ml_aztec_options_ == NULL) {
00277     ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
00278   }
00279   if (ml_aztec_params_ == NULL) {
00280     ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
00281   }
00282 
00283   if (!ml_defaults_set_) {
00284     Teuchos::ParameterList mlparams;
00285     ML_Epetra::SetDefaults("SA", mlparams,ml_aztec_options_,ml_aztec_params_);
00286     mlparams.setParameters(*paramlist_);
00287     *paramlist_ = mlparams;
00288     ml_defaults_set_ = true;
00289   }
00290 
00291   if (ml_prec_ != NULL) {
00292     delete ml_prec_; ml_prec_ = NULL;
00293   }
00294 
00295   ml_prec_ = new ML_Epetra::MultiLevelPreconditioner(*A, *paramlist_, true);
00296   //azoo_->SetPrecOperator(ml_prec_);
00297 
00298 #endif
00299 
00300   return(0);
00301 }
00302 
00303 #endif
00304 //HAVE_FEI_BELOS
00305 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends