FEI Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
fei_Solver_Amesos.cpp
Go to the documentation of this file.
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_sstream.hpp>
45
47
48#ifdef HAVE_FEI_AMESOS
49#ifdef HAVE_FEI_EPETRA
50
51#include <fei_Solver_Amesos.hpp>
52#include <fei_ParameterSet.hpp>
53#include <snl_fei_Utils.hpp>
54#include <fei_utils.hpp>
55
56//fei_Include_Trilinos.hpp includes the actual Trilinos headers (epetra, aztecoo, ml ...)
61
62#include <fei_Vector.hpp>
63#include <fei_Matrix.hpp>
64#include <fei_LinearSystem.hpp>
65
66//---------------------------------------------------------------------------
67Solver_Amesos::Solver_Amesos()
68 : tolerance_(1.e-6),
69 maxIters_(500),
70 amesos_factory_(new Amesos()),
71 amesos_solver_(NULL),
72 epetra_linearproblem_(NULL)
73{
74 paramlist_ = new Teuchos::ParameterList;
75}
76
77//---------------------------------------------------------------------------
78Solver_Amesos::~Solver_Amesos()
79{
80 delete amesos_solver_;
81 delete amesos_factory_;
82 delete epetra_linearproblem_;
83 delete paramlist_;
84}
85
86//---------------------------------------------------------------------------
87Teuchos::ParameterList& Solver_Amesos::get_ParameterList()
88{
89 if (paramlist_ == NULL) {
90 paramlist_ = new Teuchos::ParameterList;
91 }
92
93 return( *paramlist_ );
94}
95
96//---------------------------------------------------------------------------
97int Solver_Amesos::solve(fei::LinearSystem* linearSystem,
98 fei::Matrix* preconditioningMatrix,
99 const fei::ParameterSet& parameterSet,
100 int& iterationsTaken,
101 int& status)
102{
103 Trilinos_Helpers::copy_parameterset(parameterSet, *paramlist_);
104
105 int numParams = 0;
106 const char** paramStrings = NULL;
107 std::vector<std::string> stdstrings;
108 fei::utils::convert_ParameterSet_to_strings(&parameterSet, stdstrings);
109 fei::utils::strings_to_char_ptrs(stdstrings, numParams, paramStrings);
110
111 int err = solve(linearSystem, preconditioningMatrix, numParams, paramStrings,
112 iterationsTaken, status);
113
114 int olevel = 0;
115 parameterSet.getIntParamValue("outputLevel", olevel);
116
117 std::string param2;
118 parameterSet.getStringParamValue("FEI_OUTPUT_LEVEL", param2);
119
120 if (olevel >= 3 || param2 == "MATRIX_FILES" || param2 == "ALL") {
121 std::string param1;
122 parameterSet.getStringParamValue("debugOutput", param1);
123
124 FEI_OSTRINGSTREAM osstr;
125 if (!param1.empty()) {
126 osstr << param1 << "/";
127 }
128 else osstr << "./";
129
130 static int counter = 1;
131 osstr << "x_Amesos.vec.slv"<<counter++;
132 fei::SharedPtr<fei::Vector> feix = linearSystem->getSolutionVector();
133 feix->writeToFile(osstr.str().c_str());
134 }
135
136 delete [] paramStrings;
137
138 return(err);
139}
140
141//---------------------------------------------------------------------------
142int Solver_Amesos::solve(fei::LinearSystem* linearSystem,
143 fei::Matrix* preconditioningMatrix,
144 int numParams,
145 const char* const* solverParams,
146 int& iterationsTaken,
147 int& status)
148{
149 Epetra_MultiVector* x = NULL;
150 Epetra_MultiVector* b = NULL;
151 Epetra_CrsMatrix* A = NULL;
152 Epetra_Operator* opA = NULL;
153
154 fei::SharedPtr<fei::Matrix> feiA = linearSystem->getMatrix();
155 fei::SharedPtr<fei::Vector> feix = linearSystem->getSolutionVector();
156 fei::SharedPtr<fei::Vector> feib = linearSystem->getRHS();
157
158 Trilinos_Helpers::get_Epetra_pointers(feiA, feix, feib,
159 A, opA, x, b);
160
161 if (opA == 0 || x == 0 || b == 0) {
162 fei::console_out() << "Error, couldn't obtain Epetra objects from "
163 << "fei container-objects."<<FEI_ENDL;
164 return(-1);
165 }
166
167 if (epetra_linearproblem_ == NULL) {
168 epetra_linearproblem_ = new Epetra_LinearProblem;
169 }
170
171 epetra_linearproblem_->SetOperator(A);
172 epetra_linearproblem_->SetLHS(x);
173 epetra_linearproblem_->SetRHS(b);
174
175 const char* param = snl_fei::getParamValue("Trilinos_Solver",
176 numParams, solverParams);
177 if (param != NULL) {
178 if (amesos_solver_ == 0) {
179 amesos_solver_ = amesos_factory_->Create(param, *epetra_linearproblem_);
180 }
181 if (amesos_solver_ == 0) {
182 std::cerr << "Solver_Amesos::solve ERROR, couldn't create Amesos solver named "
183 << param << ", amesos_factory::Create returned NULL." << std::endl;
184 status = -999;
185 return(-1);
186 }
187 }
188 else {
189 static char amesosklu[] = "Amesos_Klu";
190 if (amesos_solver_ == 0) {
191 amesos_solver_ = amesos_factory_->Create( amesosklu,
192 *epetra_linearproblem_);
193 }
194 if (amesos_solver_ == 0) {
195 std::cerr << "Solver_Amesos::solve ERROR, couldn't create Amesos solver named "
196 << amesosklu << ", it's apparently not supported." << std::endl;
197 status = -999;
198 return(-1);
199 }
200 }
201
202 amesos_solver_->SetParameters(*paramlist_);
203
204 if (feiA->changedSinceMark()) {
205 amesos_solver_->SymbolicFactorization();
206 amesos_solver_->NumericFactorization();
207 feiA->markState();
208 }
209
210 amesos_solver_->Solve();
211 status = 0;
212
213 return(0);
214}
215
216//---------------------------------------------------------------------------
217int Solver_Amesos::parseParameters(int numParams,
218 const char*const* params)
219{
220
221 return(0);
222}
223
224#endif //HAVE_FEI_EPETRA
225#endif //HAVE_FEI_AMESOS
void SetOperator(Epetra_RowMatrix *A)
virtual fei::SharedPtr< fei::Vector > getRHS()
virtual fei::SharedPtr< fei::Matrix > getMatrix()
virtual fei::SharedPtr< fei::Vector > getSolutionVector()
int getIntParamValue(const char *name, int &paramValue) const
int getStringParamValue(const char *name, std::string &paramValue) const
#define FEI_ENDL
#define FEI_OSTRINGSTREAM
void copy_parameterset(const fei::ParameterSet &paramset, Teuchos::ParameterList &paramlist)
void convert_ParameterSet_to_strings(const fei::ParameterSet *paramset, std::vector< std::string > &paramStrings)
void strings_to_char_ptrs(std::vector< std::string > &stdstrings, int &numStrings, const char **&charPtrs)
std::ostream & console_out()
const char * getParamValue(const char *key, int numParams, const char *const *paramStrings, char separator=' ')