Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_Amesos.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
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 Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include "Ifpack_ConfigDefs.h"
45#include "Ifpack_Amesos.h"
46#include "Ifpack_Condest.h"
47#include "Epetra_MultiVector.h"
48#include "Epetra_Map.h"
49#include "Epetra_Comm.h"
50#include "Amesos.h"
52#include "Epetra_RowMatrix.h"
53#include "Epetra_Time.h"
54#include "Teuchos_ParameterList.hpp"
55
56static bool FirstTime = true;
57
58//==============================================================================
60 Matrix_(Teuchos::rcp( Matrix_in, false )),
61 Label_("Amesos_Klu"),
62 IsEmpty_(false),
63 IsInitialized_(false),
64 IsComputed_(false),
65 UseTranspose_(false),
66 NumInitialize_(0),
67 NumCompute_(0),
68 NumApplyInverse_(0),
69 InitializeTime_(0.0),
70 ComputeTime_(0.0),
71 ApplyInverseTime_(0.0),
72 ComputeFlops_(0),
73 ApplyInverseFlops_(0),
74 Condest_(-1.0)
75{
76 Problem_ = Teuchos::rcp( new Epetra_LinearProblem );
77}
78
79//==============================================================================
81 Matrix_(Teuchos::rcp( &rhs.Matrix(), false )),
82 Label_(rhs.Label()),
83 IsEmpty_(false),
84 IsInitialized_(false),
85 IsComputed_(false),
86 NumInitialize_(rhs.NumInitialize()),
87 NumCompute_(rhs.NumCompute()),
88 NumApplyInverse_(rhs.NumApplyInverse()),
89 InitializeTime_(rhs.InitializeTime()),
90 ComputeTime_(rhs.ComputeTime()),
91 ApplyInverseTime_(rhs.ApplyInverseTime()),
92 ComputeFlops_(rhs.ComputeFlops()),
93 ApplyInverseFlops_(rhs.ApplyInverseFlops()),
94 Condest_(rhs.Condest())
95{
96
97 Problem_ = Teuchos::rcp( new Epetra_LinearProblem );
98
99 // copy the RHS list in *this.List
100 Teuchos::ParameterList RHSList(rhs.List());
101 List_ = RHSList;
102
103 // I do not have a copy constructor for Amesos,
104 // so Initialize() and Compute() of this object
105 // are called if the rhs did so
106 if (rhs.IsInitialized()) {
107 IsInitialized_ = true;
108 Initialize();
109 }
110 if (rhs.IsComputed()) {
111 IsComputed_ = true;
112 Compute();
113 }
114
115}
116//==============================================================================
117int Ifpack_Amesos::SetParameters(Teuchos::ParameterList& List_in)
118{
119
120 List_ = List_in;
121 Label_ = List_in.get("amesos: solver type", Label_);
122 return(0);
123}
124
125//==============================================================================
127{
128 using std::cerr;
129 using std::endl;
130
131 IsEmpty_ = false;
132 IsInitialized_ = false;
133 IsComputed_ = false;
134
135 if (Matrix_ == Teuchos::null)
136 IFPACK_CHK_ERR(-1);
137
138#if 0
139 using std::cout;
140
141 // better to avoid strange games with maps, this class should be
142 // used for Ifpack_LocalFilter'd matrices only
143 if (Comm().NumProc() != 1) {
144 cout << "Class Ifpack_Amesos must be used for serial runs;" << endl;
145 cout << "for parallel runs you should declare objects as:" << endl;
146 cout << "Ifpack_AdditiveSchwarz<Ifpack_Amesos> APrec(Matrix)" << endl;
147 exit(EXIT_FAILURE);
148 }
149#endif
150
151 // only square matrices
152 if (Matrix_->NumGlobalRows64() != Matrix_->NumGlobalCols64())
153 IFPACK_CHK_ERR(-1);
154
155 // if the matrix has a dimension of 0, this is an empty preconditioning object.
156 if (Matrix_->NumGlobalRows64() == 0) {
157 IsEmpty_ = true;
158 IsInitialized_ = true;
160 return(0);
161 }
162
163 Problem_->SetOperator(const_cast<Epetra_RowMatrix*>(Matrix_.get()));
164
165 // create timer, which also starts it.
166 if (Time_ == Teuchos::null)
167 Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
168
169 Amesos Factory;
170 Solver_ = Teuchos::rcp( Factory.Create((char*)Label_.c_str(),*Problem_) );
171
172 if (Solver_ == Teuchos::null)
173 {
174 // try to create KLU, it is generally enabled
175 Label_ = "Amesos_Klu";
176 Solver_ = Teuchos::rcp( Factory.Create("Amesos_Klu",*Problem_) );
177 }
178 if (Solver_ == Teuchos::null)
179 {
180 // finally try to create LAPACK, it is generally enabled
181 // NOTE: the matrix is dense, so this should only be for
182 // small problems...
183 if (FirstTime)
184 {
185 cerr << "IFPACK WARNING: In class Ifpack_Amesos:" << endl;
186 cerr << "IFPACK WARNING: Using LAPACK because other Amesos" << endl;
187 cerr << "IFPACK WARNING: solvers are not available. LAPACK" << endl;
188 cerr << "IFPACK WARNING: allocates memory to store the matrix as" << endl;
189 cerr << "IFPACK WARNING: dense, I hope you have enough memory..." << endl;
190 cerr << "IFPACK WARNING: (file " << __FILE__ << ", line " << __LINE__
191 << ")" << endl;
192 FirstTime = false;
193 }
194 Label_ = "Amesos_Lapack";
195 Solver_ = Teuchos::rcp( Factory.Create("Amesos_Lapack",*Problem_) );
196 }
197 // if empty, give up.
198 if (Solver_ == Teuchos::null)
199 IFPACK_CHK_ERR(-1);
200
201 IFPACK_CHK_ERR(Solver_->SetUseTranspose(UseTranspose_));
202 Solver_->SetParameters(List_);
203 IFPACK_CHK_ERR(Solver_->SymbolicFactorization());
204
205 IsInitialized_ = true;
207 InitializeTime_ += Time_->ElapsedTime();
208 return(0);
209}
210
211//==============================================================================
213{
214
215 if (!IsInitialized())
217
218 if (IsEmpty_) {
219 IsComputed_ = true;
220 ++NumCompute_;
221 return(0);
222 }
223
224 IsComputed_ = false;
225 Time_->ResetStartTime();
226
227 if (Matrix_ == Teuchos::null)
228 IFPACK_CHK_ERR(-1);
229
230 IFPACK_CHK_ERR(Solver_->NumericFactorization());
231
232 IsComputed_ = true;
233 ++NumCompute_;
234 ComputeTime_ += Time_->ElapsedTime();
235 return(0);
236}
237
238//==============================================================================
239int Ifpack_Amesos::SetUseTranspose(bool UseTranspose_in)
240{
241 // store the value in UseTranspose_. If we have the solver, we pass to it
242 // right away, otherwise we wait till when it is created.
243 UseTranspose_ = UseTranspose_in;
244 if (Solver_ != Teuchos::null)
245 IFPACK_CHK_ERR(Solver_->SetUseTranspose(UseTranspose_in));
246
247 return(0);
248}
249
250//==============================================================================
253{
254 // check for maps ???
255 IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
256 return(0);
257}
258
259//==============================================================================
262{
263 if (IsEmpty_) {
265 return(0);
266 }
267
268 if (IsComputed() == false)
269 IFPACK_CHK_ERR(-1);
270
271 if (X.NumVectors() != Y.NumVectors())
272 IFPACK_CHK_ERR(-1); // wrong input
273
274 Time_->ResetStartTime();
275
276 // AztecOO gives X and Y pointing to the same memory location,
277 // need to create an auxiliary vector, Xcopy
278 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
279 if (X.Pointers()[0] == Y.Pointers()[0])
280 Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
281 else
282 Xcopy = Teuchos::rcp( &X, false );
283
284 Problem_->SetLHS(&Y);
285 Problem_->SetRHS((Epetra_MultiVector*)Xcopy.get());
286 IFPACK_CHK_ERR(Solver_->Solve());
287
289 ApplyInverseTime_ += Time_->ElapsedTime();
290
291 return(0);
292}
293
294//==============================================================================
296{
297 return(-1.0);
298}
299
300//==============================================================================
301const char* Ifpack_Amesos::Label() const
302{
303 return((char*)Label_.c_str());
304}
305
306//==============================================================================
308{
309 return(UseTranspose_);
310}
311
312//==============================================================================
314{
315 return(false);
316}
317
318//==============================================================================
320{
321 return(Matrix_->Comm());
322}
323
324//==============================================================================
326{
327 return(Matrix_->OperatorDomainMap());
328}
329
330//==============================================================================
332{
333 return(Matrix_->OperatorRangeMap());
334}
335
336//==============================================================================
338 const int MaxIters, const double Tol,
339 Epetra_RowMatrix* Matrix_in)
340{
341
342 if (!IsComputed()) // cannot compute right now
343 return(-1.0);
344
345 if (Condest_ == -1.0)
346 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
347
348 return(Condest_);
349}
350
351//==============================================================================
352std::ostream& Ifpack_Amesos::Print(std::ostream& os) const
353{
354 using std::endl;
355
356 if (!Comm().MyPID()) {
357 os << endl;
358 os << "================================================================================" << endl;
359 os << "Ifpack_Amesos: " << Label () << endl << endl;
360 os << "Condition number estimate = " << Condest() << endl;
361 os << "Global number of rows = " << Matrix_->NumGlobalRows64() << endl;
362 os << endl;
363 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
364 os << "----- ------- -------------- ------------ --------" << endl;
365 os << "Initialize() " << std::setw(5) << NumInitialize_
366 << " " << std::setw(15) << InitializeTime_
367 << " 0.0 0.0" << endl;
368 os << "Compute() " << std::setw(5) << NumCompute_
369 << " " << std::setw(15) << ComputeTime_
370 << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
371 if (ComputeTime_ != 0.0)
372 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
373 else
374 os << " " << std::setw(15) << 0.0 << endl;
375 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
376 << " " << std::setw(15) << ApplyInverseTime_
377 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
378 if (ApplyInverseTime_ != 0.0)
379 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
380 else
381 os << " " << std::setw(15) << 0.0 << endl;
382 os << "================================================================================" << endl;
383 os << endl;
384 }
385
386 return(os);
387}
static bool FirstTime
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
int NumVectors() const
double ** Pointers() const
Ifpack_Amesos: a class to use Amesos' factorizations as preconditioners.
Teuchos::RefCountPtr< Amesos_BaseSolver > Solver_
Amesos solver, use to apply the inverse of the local matrix.
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object.
double ComputeFlops_
Contains the number of flops for Compute().
int NumInitialize_
Contains the number of successful calls to Initialize().
double Condest_
Contains the estimated condition number.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
Teuchos::ParameterList List_
Contains a copy of the input parameter list.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual const Teuchos::ParameterList & List() const
double ComputeTime_
Contains the time for all successful calls to Compute().
double InitializeTime_
Contains the time for all successful calls to Initialize().
int NumCompute_
Contains the number of successful call to Compute().
virtual const char * Label() const
Returns a character string describing the operator.
virtual int Initialize()
Initializes the preconditioners.
virtual int Compute()
Computes the preconditioners.
std::string Label_
Contains the label of this object.
virtual int SetUseTranspose(bool UseTranspose_in)
If set true, transpose of this operator will be applied (not implemented).
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
bool IsEmpty_
If true, the linear system on this processor is empty, thus the preconditioner is null operation.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Teuchos::RefCountPtr< Epetra_LinearProblem > Problem_
Linear problem required by Solver_.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual double Condest() const
Returns the estimated condition number, never computes it.
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
Ifpack_Amesos(Epetra_RowMatrix *Matrix)
Constructor.
virtual std::ostream & Print(std::ostream &os) const
Prints on ostream basic information about this object.
bool IsComputed_
If true, the preconditioner has been successfully computed.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
bool UseTranspose_
If true, the preconditioner solves for the transpose of the matrix.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
virtual bool IsInitialized() const
Returns true is the preconditioner has been successfully initialized.
#define false