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