fei_Solver_AztecOO.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_AZTECOO
00012 
00013 //fei_Include_Trilinos.h includes the actual Trilinos headers (epetra, aztecoo, ...)
00014 #include <fei_Include_Trilinos.hpp>
00015 #include <fei_Trilinos_Helpers.hpp>
00016 
00017 #include <fei_Solver_AztecOO.hpp>
00018 #include <fei_LinProbMgr_EpetraBasic.hpp>
00019 #include <fei_ParameterSet.hpp>
00020 #include <fei_utils.hpp>
00021 
00022 #include <fei_VectorTraits_Epetra.hpp>
00023 #include <fei_Vector_Impl.hpp>
00024 #include <fei_MatrixTraits_Epetra.hpp>
00025 #include <fei_Matrix_Impl.hpp>
00026 #include <fei_LinearSystem.hpp>
00027 //---------------------------------------------------------------------------
00028 Solver_AztecOO::Solver_AztecOO()
00029   : tolerance_(1.e-6), maxIters_(500), useTranspose_(false), paramlist_(NULL),
00030     useML_(false),
00031 #ifdef HAVE_FEI_ML
00032    ml_prec_(NULL), ml_defaults_set_(false),
00033    ml_aztec_options_(NULL), ml_aztec_params_(NULL),
00034 #endif
00035   name_(),
00036   dbgprefix_("SlvAzoo: ")
00037 {
00038   azoo_ = new AztecOO;
00039   paramlist_ = new Teuchos::ParameterList;
00040 
00041 #ifdef HAVE_FEI_ML
00042   ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
00043   ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
00044 #endif
00045 }
00046 
00047 //---------------------------------------------------------------------------
00048 Solver_AztecOO::~Solver_AztecOO()
00049 {
00050   delete azoo_;
00051   delete paramlist_;
00052 #ifdef HAVE_FEI_ML
00053   delete [] ml_aztec_options_;
00054   delete [] ml_aztec_params_;
00055   delete ml_prec_;
00056 #endif
00057 
00058   AZ_manage_memory(0, AZ_CLEAR_ALL, 0, NULL, NULL);
00059 }
00060 
00061 //---------------------------------------------------------------------------
00062 AztecOO& Solver_AztecOO::getAztecOO()
00063 {
00064   return( *azoo_ );
00065 }
00066 
00067 //---------------------------------------------------------------------------
00068 void Solver_AztecOO::setUseML(bool useml)
00069 {
00070   useML_ = useml;
00071 }
00072 
00073 //---------------------------------------------------------------------------
00074 Teuchos::ParameterList& Solver_AztecOO::get_ParameterList()
00075 {
00076   if (paramlist_ == NULL) {
00077     paramlist_ = new Teuchos::ParameterList;
00078   }
00079 
00080   return( *paramlist_ );
00081 }
00082 
00083 //---------------------------------------------------------------------------
00084 int Solver_AztecOO::solve(fei::LinearSystem* linearSystem,
00085         fei::Matrix* preconditioningMatrix,
00086         const fei::ParameterSet& parameterSet,
00087         int& iterationsTaken,
00088         int& status)
00089 {
00090   std::string pcstring;
00091   parameterSet.getStringParamValue("AZ_precond", pcstring);
00092   if (pcstring == "ML_Op") {
00093     useML_ = true;
00094   }
00095 
00096   Teuchos::ParameterList& paramlist = get_ParameterList();
00097 
00098 #ifdef HAVE_FEI_ML
00099   if (ml_aztec_options_ == NULL)
00100     ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
00101   if (ml_aztec_params_ == NULL)
00102     ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
00103 
00104   if (!ml_defaults_set_ && useML_) {
00105     Teuchos::ParameterList mlparams;
00106     ML_Epetra::SetDefaults("SA", mlparams, ml_aztec_options_,ml_aztec_params_);
00107     mlparams.setParameters(paramlist);
00108     paramlist = mlparams;
00109     ml_defaults_set_ = true;
00110   }
00111 #endif
00112 
00113   Trilinos_Helpers::copy_parameterset(parameterSet, paramlist);
00114 
00115   fei::SharedPtr<fei::Matrix> feiA = linearSystem->getMatrix();
00116   fei::SharedPtr<fei::Vector> feix = linearSystem->getSolutionVector();
00117   fei::SharedPtr<fei::Vector> feib = linearSystem->getRHS();
00118 
00119   Epetra_MultiVector*    x = NULL;
00120   Epetra_MultiVector*    b = NULL;
00121   Epetra_Operator* epetra_op = 0;
00122   Epetra_CrsMatrix* crsA = NULL;
00123 
00124   Trilinos_Helpers::get_Epetra_pointers(feiA, feix, feib,
00125                                         crsA, epetra_op, x, b);
00126 
00127   if (epetra_op == 0 || x == 0 || b == 0) {
00128     FEI_CERR << "Solver_AztecOO::solve Error, couldn't obtain Epetra objects"
00129      << " from fei container-objects."<<FEI_ENDL;
00130     return(-1);
00131   }
00132 
00133   Epetra_LinearProblem linProb(epetra_op,x,b);
00134 
00135   //when we call azoo_->SetProblem, it will set some options. So we will
00136   //first take a copy of all options and params, then reset them after the
00137   //call to SetProblem. That way we preserve any options that have already
00138   //been set.
00139 
00140   std::vector<int> azoptions(AZ_OPTIONS_SIZE);
00141   std::vector<double> azparams(AZ_PARAMS_SIZE);
00142 
00143   const int* azoptionsptr = azoo_->GetAllAztecOptions();
00144   const double* azparamsptr = azoo_->GetAllAztecParams();
00145 
00146   int i;
00147   for(i=0; i<AZ_OPTIONS_SIZE; ++i) {
00148     azoptions[i] = azoptionsptr[i];
00149   }
00150   for(i=0; i<AZ_PARAMS_SIZE; ++i) {
00151     azparams[i] = azparamsptr[i];
00152   }
00153 
00154   Epetra_RowMatrix* precond = NULL;
00155   if (preconditioningMatrix != NULL) {
00156     fei::Matrix_Impl<Epetra_CrsMatrix>* snl_epetra_crs =
00157       dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(preconditioningMatrix);
00158     fei::Matrix_Impl<Epetra_VbrMatrix>* snl_epetra_vbr =
00159       dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(preconditioningMatrix);
00160     if (snl_epetra_crs != NULL) {
00161       precond = snl_epetra_crs->getMatrix().get();
00162     }
00163     else if (snl_epetra_vbr != NULL) {
00164       precond = snl_epetra_vbr->getMatrix().get();
00165     }
00166     else {
00167       FEI_CERR << "Solver_AztecOO::solve: ERROR getting epetra row matrix"
00168          << " from preconditioningMatrix."<<FEI_ENDL;
00169       return(-1);
00170     }
00171   }
00172 
00173   if (precond != NULL) {
00174     azoo_->SetProblem(linProb);
00175 
00176     azoo_->SetAllAztecOptions(&(azoptions[0]));
00177     azoo_->SetAllAztecParams(&(azparams[0]));
00178 
00179     azoo_->SetUseAdaptiveDefaultsTrue();
00180 
00181     azoo_->SetPrecMatrix(precond);
00182   }
00183 
00184   bool needNewPreconditioner = false;
00185 
00186   if (feiA->changedSinceMark()) {
00187     feiA->markState();
00188     needNewPreconditioner = true;
00189   }
00190 
00191   if (needNewPreconditioner) {
00192     azoo_->SetProblem(linProb);
00193 
00194     azoo_->SetAllAztecOptions(&(azoptions[0]));
00195     azoo_->SetAllAztecParams(&(azparams[0]));
00196 
00197     azoo_->SetUseAdaptiveDefaultsTrue();
00198 
00199     if (useML_) {
00200 #ifdef HAVE_FEI_ML
00201       setup_ml_operator(*azoo_, crsA);
00202 #else
00203       FEI_CERR <<"Solver_AztecOO::solve ERROR, ML requested but HAVE_FEI_ML not defined."
00204          << FEI_ENDL;
00205       return(-1);
00206 #endif
00207     }
00208     else {
00209       azoo_->SetAztecOption(AZ_pre_calc, AZ_calc);
00210       azoo_->SetAztecOption(AZ_keep_info, 1);
00211     }
00212   }
00213   else {
00214     if (!useML_) {
00215       azoo_->SetAztecOption(AZ_pre_calc, AZ_reuse);
00216     }
00217   }
00218 
00219   epetra_op->SetUseTranspose(useTranspose_);
00220 
00221   azoo_->SetParameters(paramlist);
00222 
00223   maxIters_ = azoptionsptr[AZ_max_iter];
00224   tolerance_= azparamsptr[AZ_tol];
00225 
00226   status = azoo_->Iterate(maxIters_,
00227        //2,//maxSolveAttempts
00228        tolerance_);
00229 
00230   iterationsTaken = azoo_->NumIters();
00231 
00232   int olevel = 0;
00233   parameterSet.getIntParamValue("outputLevel", olevel);
00234 
00235   std::string param2;
00236   parameterSet.getStringParamValue("FEI_OUTPUT_LEVEL", param2);
00237 
00238   if (olevel >= 3 || param2 == "MATRIX_FILES") {
00239     std::string param1;
00240     parameterSet.getStringParamValue("debugOutput", param1);
00241 
00242     FEI_OSTRINGSTREAM osstr;
00243     if (!param1.empty()) {
00244       osstr << param1 << "/";
00245     }
00246     else osstr << "./";
00247 
00248     osstr << "x_AztecOO.vec";
00249     feix->writeToFile(osstr.str().c_str());
00250   }
00251 
00252   return(0);
00253 }
00254 
00255 //---------------------------------------------------------------------------
00256 int Solver_AztecOO::solve(fei::LinearSystem* linearSystem,
00257         fei::Matrix* preconditioningMatrix,
00258         int numParams,
00259         const char* const* solverParams,
00260         int& iterationsTaken,
00261         int& status)
00262 {
00263   std::vector<std::string> stdstrings;
00264   fei::utils::char_ptrs_to_strings(numParams, solverParams, stdstrings);
00265 
00266   fei::ParameterSet paramset;
00267   fei::utils::parse_strings(stdstrings, " ", paramset);
00268 
00269   return( solve(linearSystem, preconditioningMatrix, paramset,
00270     iterationsTaken, status) );
00271 }
00272 
00273 //---------------------------------------------------------------------------
00274 int Solver_AztecOO::setup_ml_operator(AztecOO& azoo, Epetra_CrsMatrix* A)
00275 {
00276 #ifdef HAVE_FEI_ML
00277   if (ml_aztec_options_ == NULL) {
00278     ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
00279   }
00280   if (ml_aztec_params_ == NULL) {
00281     ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
00282   }
00283 
00284   if (!ml_defaults_set_) {
00285     Teuchos::ParameterList mlparams;
00286     ML_Epetra::SetDefaults("SA", mlparams,ml_aztec_options_,ml_aztec_params_);
00287     mlparams.setParameters(*paramlist_);
00288     *paramlist_ = mlparams;
00289     ml_defaults_set_ = true;
00290   }
00291 
00292   if (ml_prec_ != NULL) {
00293     delete ml_prec_; ml_prec_ = NULL;
00294   }
00295 
00296   ml_prec_ = new ML_Epetra::MultiLevelPreconditioner(*A, *paramlist_, true);
00297   azoo_->SetPrecOperator(ml_prec_);
00298 
00299 #endif
00300 
00301   return(0);
00302 }
00303 
00304 #endif
00305 //HAVE_FEI_AZTECOO
00306 

Generated on Tue Jul 13 09:27:46 2010 for FEI by  doxygen 1.4.7