00001
00002
00003
00004
00005
00006 #include "PlayaAztecSolver.hpp"
00007 #include "PlayaEpetraVector.hpp"
00008 #include "PlayaEpetraMatrix.hpp"
00009 #include "Ifpack_Preconditioner.h"
00010 #include "Ifpack.h"
00011 #include "EpetraPlayaOperator.hpp"
00012 #include "Teuchos_basic_oblackholestream.hpp"
00013
00014
00015
00016 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00017 #include "PlayaVectorImpl.hpp"
00018 #include "PlayaLinearOperatorImpl.hpp"
00019 #include "PlayaLinearSolverImpl.hpp"
00020 #endif
00021
00022 #ifdef HAVE_ML
00023 #include "ml_include.h"
00024 #include "ml_epetra_utils.h"
00025 #include "ml_epetra_operator.h"
00026 #include "ml_aztec_utils.h"
00027 #include "ml_MultiLevelPreconditioner.h"
00028 using namespace ML_Epetra;
00029 #else
00030 #error blarf
00031 #endif
00032
00033 using namespace Playa;
00034 using namespace Teuchos;
00035
00036
00037 AztecSolver::AztecSolver(const ParameterList& params)
00038 : LinearSolverBase<double>(params),
00039 options_(AZ_OPTIONS_SIZE),
00040 parameters_(AZ_PARAMS_SIZE),
00041 useML_(false),
00042 useIfpack_(false),
00043 useUserPrec_(false),
00044 aztec_recursive_iterate_(false),
00045 precParams_(),
00046 userPrec_(),
00047 aztec_status(AZ_STATUS_SIZE),
00048 aztec_proc_config(AZ_PROC_SIZE)
00049 {
00050 initParamMap();
00051
00052
00053 AZ_defaults((int*) &(options_[0]), (double*) &(parameters_[0]));
00054
00055
00056
00057 ParameterList::ConstIterator iter;
00058 for (iter=params.begin(); iter != params.end(); ++iter)
00059 {
00060 const std::string& name = params.name(iter);
00061 const ParameterEntry& entry = params.entry(iter);
00062
00063 if (entry.isList())
00064 {
00065 if (name=="Preconditioner")
00066 {
00067 precParams_ = params.sublist("Preconditioner");
00068 TEUCHOS_TEST_FOR_EXCEPTION(!precParams_.isParameter("Type"), std::runtime_error,
00069 "preconditioner type not specified in parameter list "
00070 << precParams_);
00071 if (precParams_.get<string>("Type") == "ML")
00072 {
00073 useML_ = true;
00074 }
00075 else if (precParams_.get<string>("Type") == "Ifpack")
00076 {
00077 useIfpack_ = true;
00078 }
00079 else if (precParams_.get<string>("Type") == "User")
00080 {
00081 useUserPrec_ = true;
00082 }
00083 continue;
00084 }
00085 }
00086
00087
00088 if (paramMap().find(name) == paramMap().end()) continue;
00089
00090
00091 int aztecCode = paramMap()[name];
00092
00093 if (name=="Recursive Iterate")
00094 {
00095 if (getValue<int>(entry))
00096 aztec_recursive_iterate_ = true;
00097 else
00098 aztec_recursive_iterate_ = false;
00099 continue;
00100 }
00101
00102
00103
00104
00105
00106
00107
00108
00109 if (entry.isType<string>())
00110 {
00111 std::string val = getValue<string>(entry);
00112 TEUCHOS_TEST_FOR_EXCEPTION(paramMap().find(val) == paramMap().end(),
00113 std::runtime_error,
00114 "Aztec solver ctor: [" << val << "] is not a "
00115 "valid Aztec option value");
00116 int optionVal = paramMap()[val];
00117 options_[aztecCode] = optionVal;
00118 }
00119 else if (entry.isType<int>())
00120 {
00121 int val = getValue<int>(entry);
00122 options_[aztecCode] = val;
00123 if (name=="Verbosity") setVerb(val);
00124 }
00125 else if (entry.isType<double>())
00126 {
00127 double val = getValue<double>(entry);
00128 parameters_[aztecCode] = val;
00129 }
00130 }
00131 }
00132
00133
00134 AztecSolver::AztecSolver(const Teuchos::map<int, int>& aztecOptions,
00135 const Teuchos::map<int, double>& aztecParameters)
00136 : LinearSolverBase<double>(ParameterList()),
00137 options_(AZ_OPTIONS_SIZE),
00138 parameters_(AZ_PARAMS_SIZE),
00139 useML_(false),
00140 useIfpack_(false),
00141 useUserPrec_(false),
00142 aztec_recursive_iterate_(false),
00143 precParams_(),
00144 userPrec_(),
00145 aztec_status(AZ_STATUS_SIZE),
00146 aztec_proc_config(AZ_PROC_SIZE)
00147 {
00148 if (aztecOptions.find(AZ_recursive_iterate) != aztecOptions.end())
00149 {
00150 aztec_recursive_iterate_ = true;
00151 }
00152
00153
00154 AZ_defaults((int*) &(options_[0]), (double*) &(parameters_[0]));
00155
00156
00157 map<int, int>::const_iterator opIter;
00158 for (opIter=aztecOptions.begin(); opIter!=aztecOptions.end(); opIter++)
00159 {
00160 int opKey = opIter->first;
00161 if (opKey==AZ_recursive_iterate) continue;
00162 int opValue = opIter->second;
00163 options_[opKey] = opValue;
00164 }
00165
00166
00167 map<int, double>::const_iterator parIter;
00168 for (parIter=aztecParameters.begin(); parIter!=aztecParameters.end();
00169 parIter++)
00170 {
00171 int parKey = parIter->first;
00172 double parValue = parIter->second;
00173 parameters_[parKey] = parValue;
00174 }
00175 }
00176
00177
00178 void AztecSolver::updateTolerance(const double& tol)
00179 {
00180 parameters_[AZ_tol] = tol;
00181 }
00182
00183 SolverState<double> AztecSolver::solve(const LinearOperator<double>& op,
00184 const Vector<double>& rhs,
00185 Vector<double>& soln) const
00186 {
00187 RCP<ostream> out;
00188 if (verb()==0)
00189 {
00190 out = rcp(new oblackholestream());
00191 }
00192 else
00193 {
00194 out = rcp(&Out::os(), false);
00195 }
00196 RCP<MultiLevelPreconditioner> mlPrec;
00197 RCP<Ifpack_Preconditioner> ifpackPrec;
00198 RCP<Epetra_Operator> prec;
00199
00200 Playa::Vector<double> bCopy = rhs.copy();
00201 Playa::Vector<double> xCopy = rhs.copy();
00202
00203 if (verb() > 4)
00204 {
00205 *out << "rhs=" << bCopy << std::endl;
00206 }
00207
00208 Epetra_Vector* b = EpetraVector::getConcretePtr(bCopy);
00209 Epetra_Vector* x = EpetraVector::getConcretePtr(xCopy);
00210
00211 Epetra_CrsMatrix& A = EpetraMatrix::getConcrete(op);
00212
00213 AztecOO aztec(&A, x, b);
00214
00215 aztec.SetAllAztecOptions((int*) &(options_[0]));
00216 aztec.SetAllAztecParams((double*) &(parameters_[0]));
00217 aztec.SetOutputStream(*out);
00218 aztec.SetErrorStream(*out);
00219
00220 int maxIters = options_[AZ_max_iter];
00221 double tol = parameters_[AZ_tol];
00222
00223 if (useML_)
00224 {
00225 std::string precType = precParams_.get<string>("Problem Type");
00226 ParameterList mlParams;
00227 ML_Epetra::SetDefaults(precType, mlParams);
00228
00229
00230
00231 ParameterList::ConstIterator iter;
00232 ParameterList mlSettings = precParams_.sublist("ML Settings");
00233 for (iter=mlSettings.begin(); iter!=mlSettings.end(); ++iter)
00234 {
00235 const std::string& name = mlSettings.name(iter);
00236 const ParameterEntry& entry = mlSettings.entry(iter);
00237 mlParams.setEntry(name, entry);
00238 }
00239
00240 mlPrec = rcp(new ML_Epetra::MultiLevelPreconditioner(A, mlParams));
00241 prec = rcp_dynamic_cast<Epetra_Operator>(mlPrec);
00242 }
00243 else if (useIfpack_)
00244 {
00245 Ifpack precFactory;
00246 int overlap = precParams_.get<int>("Overlap");
00247 std::string precType = precParams_.get<string>("Prec Type");
00248
00249 ParameterList ifpackParams = precParams_.sublist("Ifpack Settings");
00250
00251 ifpackPrec = rcp(precFactory.Create(precType, &A, overlap));
00252 prec = rcp_dynamic_cast<Epetra_Operator>(ifpackPrec);
00253 ifpackPrec->SetParameters(ifpackParams);
00254 ifpackPrec->Initialize();
00255 ifpackPrec->Compute();
00256 }
00257 else if (useUserPrec_)
00258 {
00259 TEUCHOS_TEST_FOR_EXCEPT(userPrec_.get() == 0);
00260 prec = userPrec_;
00261 }
00262
00263
00264 if (prec.get() != 0) aztec.SetPrecOperator(prec.get());
00265
00266 aztec.CheckInput();
00267
00268
00269
00270 if (aztec_recursive_iterate_)
00271 {
00272 aztec.recursiveIterate(maxIters, tol);
00273 }
00274 else
00275 {
00276 aztec.Iterate(maxIters, tol);
00277 }
00278
00279 soln = xCopy;
00280
00281 const double* status = aztec.GetAztecStatus();
00282 SolverStatusCode state = SolveCrashed;
00283
00284 std::string msg;
00285 switch((int) status[AZ_why])
00286 {
00287 case AZ_normal:
00288 state = SolveConverged;
00289 msg = "converged";
00290 break;
00291 case AZ_param:
00292 state = SolveCrashed;
00293 msg = "failed: parameter not available";
00294 break;
00295 case AZ_breakdown:
00296 state = SolveCrashed;
00297 msg = "failed: numerical breakdown";
00298 break;
00299 case AZ_loss:
00300 state = SolveCrashed;
00301 msg = "failed: numerical loss of precision";
00302 break;
00303 case AZ_ill_cond:
00304 state = SolveCrashed;
00305 msg = "failed: ill-conditioned Hessenberg matrix in GMRES";
00306 break;
00307 case AZ_maxits:
00308 state = SolveFailedToConverge;
00309 msg = "failed: maxiters reached without converged";
00310 break;
00311 }
00312 SolverState<double> rtn(state, "Aztec solver " + msg, (int) status[AZ_its],
00313 status[AZ_r]);
00314 return rtn;
00315 }
00316
00317
00318 void AztecSolver::setUserPrec(const LinearOperator<double>& P,
00319 const LinearSolver<double>& pSolver)
00320 {
00321 if (useUserPrec_)
00322 {
00323 userPrec_ = rcp(new Epetra::Epetra_PlayaOperator(P, pSolver));
00324 }
00325 else
00326 {
00327 TEUCHOS_TEST_FOR_EXCEPTION(!useUserPrec_, std::runtime_error,
00328 "Attempt to set user-defined preconditioner "
00329 "after another preconditioner has been specified");
00330 }
00331 }
00332
00333 void AztecSolver::initParamMap()
00334 {
00335 static bool first = true;
00336 if (first)
00337 {
00338 paramMap()["Method"]=AZ_solver;
00339 paramMap()["CG"]=AZ_cg;
00340 paramMap()["GMRES"]=AZ_gmres;
00341 paramMap()["CGS"]=AZ_cgs;
00342 paramMap()["TFQMR"]=AZ_tfqmr;
00343 paramMap()["BICGSTAB"]=AZ_bicgstab;
00344 paramMap()["Direct"]=AZ_lu;
00345 paramMap()["Precond"]=AZ_precond;
00346 paramMap()["None"]=AZ_none;
00347 paramMap()["Jacobi"]=AZ_Jacobi;
00348 paramMap()["Neumann Series"]=AZ_Neumann;
00349 paramMap()["Symmetric Gauss-Seidel"]=AZ_sym_GS;
00350 paramMap()["Least-Squares Polynomial"]=AZ_ls;
00351 paramMap()["Recursive Iterate"]=AZ_recursive_iterate;
00352 paramMap()["Domain Decomposition"]=AZ_dom_decomp;
00353 paramMap()["Subdomain Solver"]=AZ_subdomain_solve;
00354 paramMap()["Approximate Sparse LU"]=AZ_lu;
00355 paramMap()["Saad ILUT"]=AZ_ilut;
00356 paramMap()["ILU"]=AZ_ilu;
00357 paramMap()["RILU"]=AZ_rilu;
00358 paramMap()["Block ILU"]=AZ_bilu;
00359 paramMap()["Incomplete Cholesky"]=AZ_icc;
00360 paramMap()["Residual Scaling"]=AZ_conv;
00361 paramMap()["Initial"]=AZ_r0;
00362 paramMap()["RHS"]=AZ_rhs;
00363 paramMap()["Matrix"]=AZ_Anorm;
00364 paramMap()["Solution"]=AZ_sol;
00365 paramMap()["No Scaling"]=AZ_noscaled;
00366 paramMap()["Verbosity"]=AZ_output;
00367 paramMap()["All"]=AZ_all;
00368 paramMap()["Silent"]=AZ_none;
00369 paramMap()["Warnings"]=AZ_warnings;
00370 paramMap()["Final Residual"]=AZ_last;
00371 paramMap()["Graph Fill"]=AZ_graph_fill;
00372 paramMap()["Max Iterations"]=AZ_max_iter;
00373 paramMap()["Polynomial Order"]=AZ_poly_ord;
00374 paramMap()["Overlap"]=AZ_overlap;
00375 paramMap()["Overlap Type"]=AZ_type_overlap;
00376 paramMap()["Standard"]=AZ_standard;
00377 paramMap()["Symmetric"]=AZ_symmetric;
00378 paramMap()["Restart Size"]=AZ_kspace;
00379 paramMap()["Reorder ILU"]=AZ_reorder;
00380 paramMap()["Keep Factorization"]=AZ_keep_info;
00381 paramMap()["GMRES Orthogonalization"]=AZ_orthog;
00382 paramMap()["Classical Gram-Schmidt"]=AZ_classic;
00383 paramMap()["Modified Gram-Schmidt"]=AZ_modified;
00384 paramMap()["Auxiliary Vector"]=AZ_aux_vec;
00385 paramMap()["Residual"]=AZ_resid;
00386 paramMap()["Random"]=AZ_rand;
00387 paramMap()["Tolerance"]=AZ_tol;
00388 paramMap()["Drop Tolerance"]=AZ_drop;
00389 paramMap()["Fill Ratio"]=AZ_ilut_fill;
00390 paramMap()["Damping"]=AZ_omega;
00391
00392 first = false;
00393 }
00394 }
00395
00396