Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
double Condest_
Contains the estimated condition number.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
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)
Ifpack_Amesos(Epetra_RowMatrix *Matrix)
Constructor.
Teuchos::ParameterList List_
Contains a copy of the input parameter list.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
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().
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual const Teuchos::ParameterList & List() const
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:87
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 int SetUseTranspose(bool UseTranspose_in)
If set true, transpose of this operator will be applied (not implemented).
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual const char * Label() const
Returns a character string describing the operator.
virtual double Condest() const
Returns the estimated condition number, never computes it.
Teuchos::RefCountPtr< Epetra_LinearProblem > Problem_
Linear problem required by Solver_.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
#define false
bool IsComputed_
If true, the preconditioner has been successfully computed.
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
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().
virtual bool IsInitialized() const
Returns true is the preconditioner has been successfully initialized.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
bool IsEmpty_
If true, the linear system on this processor is empty, thus the preconditioner is null operation...
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
#define IFPACK_CHK_ERR(ifpack_err)
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual std::ostream & Print(std::ostream &os) const
Prints on ostream basic information about this object.
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object.