FEI Version of the Day
Loading...
Searching...
No Matches
fei_Solver_Belos.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_BELOS
47
48//fei_Include_Trilinos.h includes the Trilinos headers (epetra, belos, ...)
49#include <fei_Include_Trilinos.hpp>
50#include <fei_Trilinos_Helpers.hpp>
51
52#include <fei_Solver_Belos.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//---------------------------------------------------------------------------
64Solver_Belos::Solver_Belos()
65 : tolerance_(1.e-6), maxIters_(500), useTranspose_(false), paramlist_(),
66 useML_(false),
67#ifdef HAVE_FEI_ML
68 ml_prec_(NULL), ml_defaults_set_(false),
69 ml_aztec_options_(NULL), ml_aztec_params_(NULL),
70#endif
71 name_(),
72 dbgprefix_("SlvBelos: ")
73{
74 paramlist_ = Teuchos::rcp(new Teuchos::ParameterList);
75
76#ifdef HAVE_FEI_ML
77 ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
78 ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
79#endif
80}
81
82//---------------------------------------------------------------------------
83Solver_Belos::~Solver_Belos()
84{
85#ifdef HAVE_FEI_ML
86 delete [] ml_aztec_options_;
87 delete [] ml_aztec_params_;
88 delete ml_prec_;
89#endif
90}
91
92//---------------------------------------------------------------------------
93void Solver_Belos::setUseML(bool useml)
94{
95 useML_ = useml;
96}
97
98//---------------------------------------------------------------------------
99Teuchos::ParameterList& Solver_Belos::get_ParameterList()
100{
101 if (paramlist_.get() == NULL) {
102 paramlist_ = Teuchos::rcp(new Teuchos::ParameterList);
103 }
104
105 return( *paramlist_ );
106}
107
108//---------------------------------------------------------------------------
109int Solver_Belos::solve(fei::LinearSystem* linearSystem,
110 fei::Matrix* preconditioningMatrix,
111 const fei::ParameterSet& parameterSet,
112 int& iterationsTaken,
113 int& status)
114{
115 std::string krylov_solver_name;
116 parameterSet.getStringParamValue("krylov_solver", krylov_solver_name);
117
118 Teuchos::RCP<Teuchos::ParameterList>& paramlist = paramlist_;
119
120#ifdef HAVE_FEI_ML
121 if (ml_aztec_options_ == NULL)
122 ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
123 if (ml_aztec_params_ == NULL)
124 ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
125
126 if (!ml_defaults_set_ && useML_) {
127 Teuchos::ParameterList mlparams;
128 ML_Epetra::SetDefaults("SA", mlparams, ml_aztec_options_,ml_aztec_params_);
129 mlparams.setParameters(*paramlist);
130 *paramlist = mlparams;
131 ml_defaults_set_ = true;
132 }
133#endif
134
135 Trilinos_Helpers::copy_parameterset(parameterSet, *paramlist);
136
137 fei::SharedPtr<fei::Matrix> feiA = linearSystem->getMatrix();
138 fei::SharedPtr<fei::Vector> feix = linearSystem->getSolutionVector();
139 fei::SharedPtr<fei::Vector> feib = linearSystem->getRHS();
140
141 Epetra_MultiVector* x = NULL;
142 Epetra_MultiVector* b = NULL;
143 Epetra_Operator* epetra_op = 0;
144 Epetra_CrsMatrix* crsA = NULL;
145
146 Trilinos_Helpers::get_Epetra_pointers(feiA, feix, feib,
147 crsA, epetra_op, x, b);
148
149 Teuchos::RCP<Epetra_CrsMatrix> rcp_A(crsA);
150 Teuchos::RCP<Epetra_MultiVector> rcp_x(x);
151 Teuchos::RCP<Epetra_MultiVector> rcp_b(b);
152
153 if (epetra_op == 0 || x == 0 || b == 0) {
154 fei::console_out() << "Solver_Belos::solve Error, couldn't obtain Epetra objects"
155 << " from fei container-objects."<<FEI_ENDL;
156 return(-1);
157 }
158
159 Epetra_RowMatrix* precond = NULL;
160 if (preconditioningMatrix != NULL) {
161 fei::Matrix_Impl<Epetra_CrsMatrix>* snl_epetra_crs =
162 dynamic_cast<fei::Matrix_Impl<Epetra_CrsMatrix>*>(preconditioningMatrix);
163 fei::Matrix_Impl<Epetra_VbrMatrix>* snl_epetra_vbr =
164 dynamic_cast<fei::Matrix_Impl<Epetra_VbrMatrix>*>(preconditioningMatrix);
165 if (snl_epetra_crs != NULL) {
166 precond = snl_epetra_crs->getMatrix().get();
167 }
168 else if (snl_epetra_vbr != NULL) {
169 precond = snl_epetra_vbr->getMatrix().get();
170 }
171 else {
172 fei::console_out() << "Solver_Belos::solve: ERROR getting epetra row matrix"
173 << " from preconditioningMatrix."<<FEI_ENDL;
174 return(-1);
175 }
176 }
177
178 if (precond != NULL) {
179//TODO: set up preconditioner for Belos here
180 }
181
182 bool needNewPreconditioner = false;
183
184 if (feiA->changedSinceMark()) {
185 feiA->markState();
186 needNewPreconditioner = true;
187 }
188
189 if (needNewPreconditioner) {
190//
191// if (useML_) {
192#ifdef HAVE_FEI_ML
193// setup_ml_operator(*azoo_, crsA);
194#else
195// fei::console_out() <<"Solver_Belos::solve ERROR, ML requested but HAVE_FEI_ML not defined."
196// << FEI_ENDL;
197// return(-1);
198#endif
199// }
200// else {
201// azoo_->SetAztecOption(AZ_pre_calc, AZ_calc);
202// azoo_->SetAztecOption(AZ_keep_info, 1);
203// }
204 }
205 else {
206// if (!useML_) {
207// azoo_->SetAztecOption(AZ_pre_calc, AZ_reuse);
208// }
209 }
210
211 epetra_op->SetUseTranspose(useTranspose_);
212
213 Belos::SolverFactory<double,Epetra_MultiVector,Epetra_Operator> belos_factory;
214 belos_solver_manager_ = belos_factory.create(krylov_solver_name, paramlist);
215
216 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));
217
218 belos_lin_prob->setProblem();
219
220 belos_solver_manager_->setProblem(belos_lin_prob);
221
222 belos_solver_manager_->solve();
223 status = 0;
224
225 iterationsTaken = belos_solver_manager_->getNumIters();
226
227 rcp_A.release();
228 rcp_x.release();
229 rcp_b.release();
230
231 int olevel = 0;
232 parameterSet.getIntParamValue("outputLevel", olevel);
233
234 std::string param2;
235 parameterSet.getStringParamValue("FEI_OUTPUT_LEVEL", param2);
236
237 if (olevel >= 3 || param2 == "MATRIX_FILES" || param2 == "ALL") {
238 std::string param1;
239 parameterSet.getStringParamValue("debugOutput", param1);
240
241 FEI_OSTRINGSTREAM osstr;
242 if (!param1.empty()) {
243 osstr << param1 << "/";
244 }
245 else osstr << "./";
246
247 osstr << "x_Belos.vec";
248 feix->writeToFile(osstr.str().c_str());
249 }
250
251 return(0);
252}
253
254//---------------------------------------------------------------------------
255int Solver_Belos::solve(fei::LinearSystem* linearSystem,
256 fei::Matrix* preconditioningMatrix,
257 int numParams,
258 const char* const* solverParams,
259 int& iterationsTaken,
260 int& status)
261{
262 std::vector<std::string> stdstrings;
263 fei::utils::char_ptrs_to_strings(numParams, solverParams, stdstrings);
264
265 fei::ParameterSet paramset;
266 fei::utils::parse_strings(stdstrings, " ", paramset);
267
268 return( solve(linearSystem, preconditioningMatrix, paramset,
269 iterationsTaken, status) );
270}
271
272//---------------------------------------------------------------------------
273int Solver_Belos::setup_ml_operator(AztecOO& azoo, Epetra_CrsMatrix* A)
274{
275#ifdef HAVE_FEI_ML
276 if (ml_aztec_options_ == NULL) {
277 ml_aztec_options_ = new int[AZ_OPTIONS_SIZE];
278 }
279 if (ml_aztec_params_ == NULL) {
280 ml_aztec_params_ = new double[AZ_PARAMS_SIZE];
281 }
282
283 if (!ml_defaults_set_) {
284 Teuchos::ParameterList mlparams;
285 ML_Epetra::SetDefaults("SA", mlparams,ml_aztec_options_,ml_aztec_params_);
286 mlparams.setParameters(*paramlist_);
287 *paramlist_ = mlparams;
288 ml_defaults_set_ = true;
289 }
290
291 if (ml_prec_ != NULL) {
292 delete ml_prec_; ml_prec_ = NULL;
293 }
294
295 ml_prec_ = new ML_Epetra::MultiLevelPreconditioner(*A, *paramlist_, true);
296 //azoo_->SetPrecOperator(ml_prec_);
297
298#endif
299
300 return(0);
301}
302
303#endif
304//HAVE_FEI_BELOS
305
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()