IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_Amesos.cpp
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 {
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 //==============================================================================
117 int 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;
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())
216  IFPACK_CHK_ERR(Initialize());
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_) {
264  ++NumApplyInverse_;
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 
288  ++NumApplyInverse_;
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 //==============================================================================
337 double Ifpack_Amesos::Condest(const Ifpack_CondestType CT,
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.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
Ifpack_Amesos(Epetra_RowMatrix *Matrix)
Constructor.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
Ifpack_Amesos: a class to use Amesos&#39; factorizations as preconditioners.
Definition: Ifpack_Amesos.h:81
virtual int Compute()
Computes the preconditioners.
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.
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.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual bool IsInitialized() const
Returns true is the preconditioner has been successfully initialized.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
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.