Ifpack Package Browser (Single Doxygen Collection)  Development
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"
44 #include "Ifpack_Preconditioner.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"
51 #include "Epetra_LinearProblem.h"
52 #include "Epetra_RowMatrix.h"
53 #include "Epetra_Time.h"
54 #include "Teuchos_ParameterList.hpp"
55 
56 static 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 {
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 
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 //==============================================================================
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;
159  ++NumInitialize_;
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;
206  ++NumInitialize_;
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 //==============================================================================
239 int 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 //==============================================================================
251 int Ifpack_Amesos::
253 {
254  // check for maps ???
255  IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
256  return(0);
257 }
258 
259 //==============================================================================
260 int Ifpack_Amesos::
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 //==============================================================================
301 const 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 //==============================================================================
352 std::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 }
virtual std::ostream & Print(std::ostream &os) const
Prints on ostream basic information about this object.
virtual const char * Label() const
Returns a character string describing the operator.
double Condest_
Contains the estimated condition number.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual bool IsInitialized() const
Returns true is the preconditioner has been successfully initialized.
bool UseTranspose_
If true, the preconditioner solves for the transpose of the matrix.
double ComputeTime_
Contains the time for all successful calls to Compute().
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
std::string Label_
Contains the label of this object.
int NumInitialize_
Contains the number of successful calls to Initialize().
T & get(ParameterList &l, const std::string &name)
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
Ifpack_Amesos(Epetra_RowMatrix *Matrix)
Constructor.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Teuchos::ParameterList List_
Contains a copy of the input parameter list.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
Teuchos::RefCountPtr< Amesos_BaseSolver > Solver_
Amesos solver, use to apply the inverse of the local matrix.
double InitializeTime_
Contains the time for all successful calls to Initialize().
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
Ifpack_Amesos: a class to use Amesos&#39; factorizations as preconditioners.
Definition: Ifpack_Amesos.h:81
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual int Compute()
Computes the preconditioners.
double ComputeFlops_
Contains the number of flops for Compute().
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static bool FirstTime
virtual int Initialize()
Initializes the preconditioners.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
virtual int SetUseTranspose(bool UseTranspose_in)
If set true, transpose of this operator will be applied (not implemented).
int NumVectors() const
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Teuchos::RefCountPtr< Epetra_LinearProblem > Problem_
Linear problem required by Solver_.
#define false
bool IsComputed_
If true, the preconditioner has been successfully computed.
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
int NumCompute_
Contains the number of successful call to Compute().
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
virtual double Condest() const
Returns the estimated condition number, never computes it.
bool IsEmpty_
If true, the linear system on this processor is empty, thus the preconditioner is null operation...
#define IFPACK_CHK_ERR(ifpack_err)
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
double ** Pointers() const
virtual const Teuchos::ParameterList & List() const
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object.