FEI  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 //---------------------------------------------------------------------------
64 Solver_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 //---------------------------------------------------------------------------
83 Solver_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 //---------------------------------------------------------------------------
93 void Solver_Belos::setUseML(bool useml)
94 {
95  useML_ = useml;
96 }
97 
98 //---------------------------------------------------------------------------
99 Teuchos::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 //---------------------------------------------------------------------------
109 int 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);
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 //---------------------------------------------------------------------------
255 int 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 //---------------------------------------------------------------------------
273 int 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 
void char_ptrs_to_strings(int numStrings, const char *const *charstrings, std::vector< std::string > &stdstrings)
Definition: fei_utils.cpp:164
virtual int writeToFile(const char *filename, bool matrixMarketFormat=true)=0
virtual fei::SharedPtr< fei::Matrix > getMatrix()
virtual int SetUseTranspose(bool UseTranspose)=0
virtual void markState()=0
virtual fei::SharedPtr< fei::Vector > getRHS()
int getStringParamValue(const char *name, std::string &paramValue) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T * get() const
fei::SharedPtr< T > getMatrix()
virtual bool changedSinceMark()=0
Ptr< T > release()
ParameterList & setParameters(const ParameterList &source)
std::ostream & console_out()
void parse_strings(std::vector< std::string > &stdstrings, const char *separator_string, fei::ParameterSet &paramset)
Definition: fei_utils.cpp:191
virtual fei::SharedPtr< fei::Vector > getSolutionVector()
int getIntParamValue(const char *name, int &paramValue) const