FEI Version of the Day
Loading...
Searching...
No Matches
fei_Solver_AztecOO.cpp
1/*
2// @HEADER
3// ************************************************************************
4// FEI: Finite Element Interface to Linear Solvers
5// Copyright (2005) Sandia Corporation.
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
8// U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Alan Williams (william@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41*/
42
43
44#include "fei_trilinos_macros.hpp"
45
46#ifdef HAVE_FEI_AZTECOO
47
48//fei_Include_Trilinos.h includes the actual Trilinos headers (epetra, aztecoo, ...)
49#include <fei_Include_Trilinos.hpp>
50#include <fei_Trilinos_Helpers.hpp>
51
52#include <fei_Solver_AztecOO.hpp>
53#include <fei_LinProbMgr_EpetraBasic.hpp>
54#include <fei_ParameterSet.hpp>
55#include <fei_utils.hpp>
56
57#include <fei_VectorTraits_Epetra.hpp>
58#include <fei_Vector_Impl.hpp>
59#include <fei_MatrixTraits_Epetra.hpp>
60#include <fei_Matrix_Impl.hpp>
61#include <fei_LinearSystem.hpp>
62
63#ifdef HAVE_FEI_ML
64#include <ml_LevelWrap.h>
65#endif
66
67//---------------------------------------------------------------------------
68Solver_AztecOO::Solver_AztecOO()
69 : tolerance_(1.e-6), maxIters_(500), useTranspose_(false), paramlist_(NULL),
70 linProb(NULL),
71 useML_(false),
72#ifdef HAVE_FEI_ML
73 ml_prec_(NULL), ml_defaults_set_(false),
74 ml_aztec_options_(NULL), ml_aztec_params_(NULL),
75#endif
76 name_(),
77 dbgprefix_("SlvAzoo: ")
78{
79 azoo_ = new AztecOO;
80 paramlist_ = new Teuchos::ParameterList;
81
82#ifdef HAVE_FEI_ML
83 ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
84 ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
85#endif
86}
87
88//---------------------------------------------------------------------------
89Solver_AztecOO::~Solver_AztecOO()
90{
91 delete azoo_;
92 delete paramlist_;
93 delete linProb;
94#ifdef HAVE_FEI_ML
95 delete [] ml_aztec_options_;
96 delete [] ml_aztec_params_;
97 delete ml_prec_;
98#endif
99
100 AZ_manage_memory(0, AZ_CLEAR_ALL, 0, NULL, NULL);
101}
102
103//---------------------------------------------------------------------------
104AztecOO& Solver_AztecOO::getAztecOO()
105{
106 return( *azoo_ );
107}
108
109//---------------------------------------------------------------------------
110void Solver_AztecOO::setUseML(bool useml)
111{
112 useML_ = useml;
113}
114
115//---------------------------------------------------------------------------
116Teuchos::ParameterList& Solver_AztecOO::get_ParameterList()
117{
118 if (paramlist_ == NULL) {
119 paramlist_ = new Teuchos::ParameterList;
120 }
121
122 return( *paramlist_ );
123}
124
125//---------------------------------------------------------------------------
126int Solver_AztecOO::solve(fei::LinearSystem* linearSystem,
127 fei::Matrix* preconditioningMatrix,
128 const fei::ParameterSet& parameterSet,
129 int& iterationsTaken,
130 int& status)
131{
132 std::string pcstring;
133 parameterSet.getStringParamValue("AZ_precond", pcstring);
134 if (pcstring == "ML_Op") {
135 useML_ = true;
136 }
137
138 Teuchos::ParameterList& paramlist = get_ParameterList();
139
140#ifdef HAVE_FEI_ML
141 if (ml_aztec_options_ == NULL)
142 ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
143 if (ml_aztec_params_ == NULL)
144 ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
145
146 if (!ml_defaults_set_ && useML_) {
147 Teuchos::ParameterList mlparams;
148 ML_Epetra::SetDefaults("SA", mlparams, ml_aztec_options_,ml_aztec_params_);
149 mlparams.setParameters(paramlist);
150 paramlist = mlparams;
151 ml_defaults_set_ = true;
152 }
153#endif
154
155 Trilinos_Helpers::copy_parameterset(parameterSet, paramlist);
156
157 fei::SharedPtr<fei::Matrix> feiA = linearSystem->getMatrix();
158 fei::SharedPtr<fei::Vector> feix = linearSystem->getSolutionVector();
159 fei::SharedPtr<fei::Vector> feib = linearSystem->getRHS();
160
161 Epetra_MultiVector* x = NULL;
162 Epetra_MultiVector* b = NULL;
163 Epetra_Operator* epetra_op = 0;
164 Epetra_CrsMatrix* crsA = NULL;
165
166 Trilinos_Helpers::get_Epetra_pointers(feiA, feix, feib,
167 crsA, epetra_op, x, b);
168
169 if (epetra_op == 0 || x == 0 || b == 0) {
170 fei::console_out() << "Solver_AztecOO::solve Error, couldn't obtain Epetra objects"
171 << " from fei container-objects."<<FEI_ENDL;
172 return(-1);
173 }
174
175 //when we call azoo_->SetProblem, it will set some options. So we will
176 //first take a copy of all options and params, then reset them after the
177 //call to SetProblem. That way we preserve any options that have already
178 //been set.
179
180 std::vector<int> azoptions(AZ_OPTIONS_SIZE);
181 std::vector<double> azparams(AZ_PARAMS_SIZE);
182
183 const int* azoptionsptr = azoo_->GetAllAztecOptions();
184 const double* azparamsptr = azoo_->GetAllAztecParams();
185
186 int i;
187 for(i=0; i<AZ_OPTIONS_SIZE; ++i) {
188 azoptions[i] = azoptionsptr[i];
189 }
190 for(i=0; i<AZ_PARAMS_SIZE; ++i) {
191 azparams[i] = azparamsptr[i];
192 }
193
194 Epetra_RowMatrix* precond = NULL;
195 if (preconditioningMatrix != NULL) {
196 fei::Matrix_Impl<Epetra_CrsMatrix>* snl_epetra_crs =
197 dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(preconditioningMatrix);
198 fei::Matrix_Impl<Epetra_VbrMatrix>* snl_epetra_vbr =
199 dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(preconditioningMatrix);
200 if (snl_epetra_crs != NULL) {
201 precond = snl_epetra_crs->getMatrix().get();
202 }
203 else if (snl_epetra_vbr != NULL) {
204 precond = snl_epetra_vbr->getMatrix().get();
205 }
206 else {
207 fei::console_out() << "Solver_AztecOO::solve: ERROR getting epetra row matrix"
208 << " from preconditioningMatrix."<<FEI_ENDL;
209 return(-1);
210 }
211 }
212
213 if (precond != NULL) {
214 Epetra_LinearProblem * newlinProb = new Epetra_LinearProblem(epetra_op,x,b);
215 azoo_->SetProblem(*newlinProb);
216 delete linProb;
217 linProb = newlinProb;
218
219 azoo_->SetAllAztecOptions(&(azoptions[0]));
220 azoo_->SetAllAztecParams(&(azparams[0]));
221
222 azoo_->SetUseAdaptiveDefaultsTrue();
223
224 azoo_->SetPrecMatrix(precond);
225 }
226
227 bool needNewPreconditioner = false;
228
229 if (feiA->changedSinceMark()) {
230 feiA->markState();
231 needNewPreconditioner = true;
232 }
233
234 if (needNewPreconditioner) {
235 Epetra_LinearProblem * newlinProb = new Epetra_LinearProblem(epetra_op,x,b);
236 azoo_->SetProblem(*newlinProb);
237 delete linProb;
238 linProb = newlinProb;
239
240 azoo_->SetAllAztecOptions(&(azoptions[0]));
241 azoo_->SetAllAztecParams(&(azparams[0]));
242
243 azoo_->SetUseAdaptiveDefaultsTrue();
244
245 if (useML_) {
246#ifdef HAVE_FEI_ML
247 setup_ml_operator(*azoo_, crsA);
248#else
249 fei::console_out() <<"Solver_AztecOO::solve ERROR, ML requested but HAVE_FEI_ML not defined."
250 << FEI_ENDL;
251 return(-1);
252#endif
253 }
254 else {
255 azoo_->SetAztecOption(AZ_pre_calc, AZ_calc);
256 azoo_->SetAztecOption(AZ_keep_info, 1);
257 }
258 }
259 else {
260 if (!useML_) {
261 azoo_->SetAztecOption(AZ_pre_calc, AZ_reuse);
262 }
263 }
264
265 epetra_op->SetUseTranspose(useTranspose_);
266
267 azoo_->SetParameters(paramlist);
268
269 maxIters_ = azoptionsptr[AZ_max_iter];
270 tolerance_= azparamsptr[AZ_tol];
271
272 status = azoo_->Iterate(maxIters_,
273 //2,//maxSolveAttempts
274 tolerance_);
275
276 iterationsTaken = azoo_->NumIters();
277
278 int olevel = 0;
279 parameterSet.getIntParamValue("outputLevel", olevel);
280
281 std::string param2;
282 parameterSet.getStringParamValue("FEI_OUTPUT_LEVEL", param2);
283
284 if (olevel >= 3 || param2 == "MATRIX_FILES" || param2 == "ALL") {
285 std::string param1;
286 parameterSet.getStringParamValue("debugOutput", param1);
287
288 FEI_OSTRINGSTREAM osstr;
289 if (!param1.empty()) {
290 osstr << param1 << "/";
291 }
292 else osstr << "./";
293
294 osstr << "x_AztecOO.vec";
295 feix->writeToFile(osstr.str().c_str());
296 }
297
298 return(0);
299}
300
301//---------------------------------------------------------------------------
302int Solver_AztecOO::solve(fei::LinearSystem* linearSystem,
303 fei::Matrix* preconditioningMatrix,
304 int numParams,
305 const char* const* solverParams,
306 int& iterationsTaken,
307 int& status)
308{
309 std::vector<std::string> stdstrings;
310 fei::utils::char_ptrs_to_strings(numParams, solverParams, stdstrings);
311
312 fei::ParameterSet paramset;
313 fei::utils::parse_strings(stdstrings, " ", paramset);
314
315 return( solve(linearSystem, preconditioningMatrix, paramset,
316 iterationsTaken, status) );
317}
318
319//---------------------------------------------------------------------------
320int Solver_AztecOO::setup_ml_operator(AztecOO& azoo, Epetra_CrsMatrix* A)
321{
322#ifdef HAVE_FEI_ML
323 if (ml_aztec_options_ == NULL) {
324 ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
325 }
326 if (ml_aztec_params_ == NULL) {
327 ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
328 }
329
330 if (!ml_defaults_set_) {
331 Teuchos::ParameterList mlparams;
332 ML_Epetra::SetDefaults("SA", mlparams,ml_aztec_options_,ml_aztec_params_);
333 mlparams.setParameters(*paramlist_);
334 *paramlist_ = mlparams;
335 ml_defaults_set_ = true;
336 }
337
338 if (ml_prec_ != NULL) {
339 delete ml_prec_; ml_prec_ = NULL;
340 }
341
342 const bool variableDof = paramlist_->get("ML variable DOF",false);
343 const bool levelWrap = paramlist_->get("ML LevelWrap",false);
344 assert(!variableDof || !levelWrap);
345 if (variableDof)
346 {
347 bool * dofPresent = paramlist_->get("DOF present",(bool*)0);
348 int nOwnedEntities = paramlist_->get("num owned entities", 0);
349 assert( (nullptr != dofPresent) || (nOwnedEntities == 0));
350 int maxDofPerEntity = paramlist_->get("max DOF per entity", 0);
351
352 Epetra_Vector *dummy = nullptr;
353 ml_prec_ = new ML_Epetra::MultiLevelPreconditioner(*A, *paramlist_,
354 nOwnedEntities,maxDofPerEntity, dofPresent,
355 *dummy, *dummy, false,true);
356 }
357 else if (levelWrap)
358 {
359 // Set the LevelWrap default parameters
360 ML_Epetra::SetDefaultsLevelWrap(*paramlist_, false);
361 Teuchos::ParameterList &crsparams = paramlist_->sublist("levelwrap: coarse list");
362 // Set the smoother parameters to the provided arrays
363 Teuchos::ParameterList &ifparams = crsparams.sublist("smoother: ifpack list");
364 int N_indices = paramlist_->get("LevelWrap number of local smoothing indices",(int)0);
365 int * smth_indices = paramlist_->get("LevelWrap local smoothing indices",(int*)0);
366 ifparams.set("relaxation: number of local smoothing indices",N_indices);
367 ifparams.set("relaxation: local smoothing indices",smth_indices);
368
369 Epetra_CrsMatrix * volAvgMtx = paramlist_->get("LevelWrap volume average matrix",(Epetra_CrsMatrix*)0);
370 ml_prec_ = new ML_Epetra::LevelWrap(Teuchos::rcp(A, false), Teuchos::rcp(volAvgMtx, false), *paramlist_);
371 }
372 else
373 {
374 ml_prec_ = new ML_Epetra::MultiLevelPreconditioner(*A, *paramlist_, true);
375 }
376
377 azoo_->SetPrecOperator(ml_prec_);
378#endif
379
380 return(0);
381}
382
383#endif
384//HAVE_FEI_AZTECOO
385
virtual int SetUseTranspose(bool UseTranspose)=0
virtual fei::SharedPtr< fei::Vector > getRHS()
virtual fei::SharedPtr< fei::Matrix > getMatrix()
virtual fei::SharedPtr< fei::Vector > getSolutionVector()
fei::SharedPtr< T > getMatrix()
int getIntParamValue(const char *name, int &paramValue) const
int getStringParamValue(const char *name, std::string &paramValue) const
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
void char_ptrs_to_strings(int numStrings, const char *const *charstrings, std::vector< std::string > &stdstrings)
std::ostream & console_out()