IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_Krylov.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 <iomanip>
45 #include "Epetra_Operator.h"
46 #include "Epetra_RowMatrix.h"
47 #include "Epetra_Comm.h"
48 #include "Epetra_Map.h"
49 #include "Epetra_MultiVector.h"
50 #include "Epetra_Vector.h"
51 #include "Epetra_Time.h"
52 #include "Ifpack_Krylov.h"
53 #include "Ifpack_Utils.h"
54 #include "Ifpack_Condest.h"
55 #ifdef HAVE_IFPACK_AZTECOO
56 #include "AztecOO.h"
57 #endif
58 
59 #ifdef HAVE_IFPACK_EPETRAEXT
60 #include "Epetra_CrsMatrix.h"
61 #include "EpetraExt_PointToBlockDiagPermute.h"
62 #endif
63 
64 //==============================================================================
65 // NOTE: any change to the default values should be committed to the other
66 // constructor as well.
67 Ifpack_Krylov::
68 Ifpack_Krylov(Epetra_Operator* Operator) :
69  IsInitialized_(false),
70  IsComputed_(false),
71  NumInitialize_(0),
72  NumCompute_(0),
73  NumApplyInverse_(0),
74  InitializeTime_(0.0),
75  ComputeTime_(0.0),
76  ApplyInverseTime_(0.0),
77  ComputeFlops_(0.0),
78  ApplyInverseFlops_(0.0),
79  Iterations_(5),
80  Tolerance_(1e-6),
81  SolverType_(1),
82  PreconditionerType_(1),
83  NumSweeps_(0),
84  BlockSize_(1),
85  DampingParameter_(1.0),
86  UseTranspose_(false),
87  Condest_(-1.0),
88  /* ComputeCondest_(false), (unused; commented out to avoid build warnings) */
89  Label_(),
90  NumMyRows_(0),
91  NumMyNonzeros_(0),
92  NumGlobalRows_(0),
93  NumGlobalNonzeros_(0),
94  Operator_(Teuchos::rcp(Operator,false)),
95  IsRowMatrix_(false),
96  ZeroStartingSolution_(true)
97 {
98 }
99 
100 //==============================================================================
101 // NOTE: This constructor has been introduced because SWIG does not appear
102 // to appreciate dynamic_cast. An instruction of type
103 // Matrix_ = dynamic_cast<const Epetra_RowMatrix*> in the
104 // other construction does not work in PyTrilinos -- of course
105 // it does in any C++ code (for an Epetra_Operator that is also
106 // an Epetra_RowMatrix).
107 //
108 // FIXME: move declarations into a separate method?
109 Ifpack_Krylov::
110 Ifpack_Krylov(Epetra_RowMatrix* Operator) :
111  IsInitialized_(false),
112  IsComputed_(false),
113  NumInitialize_(0),
114  NumCompute_(0),
115  NumApplyInverse_(0),
116  InitializeTime_(0.0),
117  ComputeTime_(0.0),
118  ApplyInverseTime_(0.0),
119  ComputeFlops_(0.0),
120  ApplyInverseFlops_(0.0),
121  Iterations_(5),
122  Tolerance_(1e-6),
123  SolverType_(1),
124  PreconditionerType_(1),
125  NumSweeps_(0),
126  BlockSize_(1),
127  DampingParameter_(1.0),
128  UseTranspose_(false),
129  Condest_(-1.0),
130  /* ComputeCondest_(false), (unused; commented out to avoid build warnings) */
131  Label_(),
132  NumMyRows_(0),
133  NumMyNonzeros_(0),
134  NumGlobalRows_(0),
135  NumGlobalNonzeros_(0),
136  Operator_(Teuchos::rcp(Operator,false)),
137  Matrix_(Teuchos::rcp(Operator,false)),
138  IsRowMatrix_(true),
139  ZeroStartingSolution_(true)
140 {
141 }
142 
143 //==============================================================================
144 int Ifpack_Krylov::SetParameters(Teuchos::ParameterList& List)
145 {
146  Iterations_ = List.get("krylov: iterations",Iterations_);
147  Tolerance_ = List.get("krylov: tolerance",Tolerance_);
148  SolverType_ = List.get("krylov: solver",SolverType_);
149  PreconditionerType_ = List.get("krylov: preconditioner",PreconditionerType_);
150  NumSweeps_ = List.get("krylov: number of sweeps",NumSweeps_);
151  BlockSize_ = List.get("krylov: block size",BlockSize_);
152  DampingParameter_ = List.get("krylov: damping parameter",DampingParameter_);
153  ZeroStartingSolution_ = List.get("krylov: zero starting solution",ZeroStartingSolution_);
154  SetLabel();
155  return(0);
156 }
157 
158 //==============================================================================
160 {
161  return(Operator_->Comm());
162 }
163 
164 //==============================================================================
166 {
167  return(Operator_->OperatorDomainMap());
168 }
169 
170 //==============================================================================
172 {
173  return(Operator_->OperatorRangeMap());
174 }
175 
176 //==============================================================================
177 int Ifpack_Krylov::
179 {
180  if (IsComputed() == false)
181  IFPACK_CHK_ERR(-3);
182 
183  if (X.NumVectors() != Y.NumVectors())
184  IFPACK_CHK_ERR(-2);
185 
186  if (IsRowMatrix_)
187  {
188  IFPACK_CHK_ERR(Matrix_->Multiply(UseTranspose(),X,Y));
189  }
190  else
191  {
192  IFPACK_CHK_ERR(Operator_->Apply(X,Y));
193  }
194 
195  return(0);
196 }
197 
198 //==============================================================================
200 {
201  IsInitialized_ = false;
202 
203  if (Operator_ == Teuchos::null)
204  IFPACK_CHK_ERR(-2);
205 
206  if (Time_ == Teuchos::null)
207  Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
208 
209  if (IsRowMatrix_)
210  {
211  if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
212  IFPACK_CHK_ERR(-2); // only square matrices
213 
214  NumMyRows_ = Matrix_->NumMyRows();
215  NumMyNonzeros_ = Matrix_->NumMyNonzeros();
216  NumGlobalRows_ = Matrix_->NumGlobalRows64();
217  NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
218  }
219  else
220  {
221  if (Operator_->OperatorDomainMap().NumGlobalElements64() !=
222  Operator_->OperatorRangeMap().NumGlobalElements64())
223  IFPACK_CHK_ERR(-2); // only square operators
224  }
225 
226  ++NumInitialize_;
227  InitializeTime_ += Time_->ElapsedTime();
228  IsInitialized_ = true;
229  return(0);
230 }
231 
232 //==============================================================================
234 {
235  if (!IsInitialized())
236  IFPACK_CHK_ERR(Initialize());
237 
238  Time_->ResetStartTime();
239 
240 #ifdef HAVE_IFPACK_AZTECOO
241  // setup Aztec solver
242  AztecSolver_ = Teuchos::rcp( new AztecOO );
243  if(IsRowMatrix_==true) {
244  AztecSolver_ -> SetUserMatrix(&*Matrix_);
245  }
246  else {
247  AztecSolver_ -> SetUserOperator(&*Operator_);
248  }
249  if(SolverType_==0) {
250  AztecSolver_ -> SetAztecOption(AZ_solver, AZ_cg);
251  }
252  else {
253  AztecSolver_ -> SetAztecOption(AZ_solver, AZ_gmres);
254  }
255  AztecSolver_ -> SetAztecOption(AZ_output, AZ_none);
256  // setup preconditioner
257  Teuchos::ParameterList List;
258  List.set("relaxation: damping factor", DampingParameter_);
259  List.set("relaxation: sweeps",NumSweeps_);
260  if(PreconditionerType_==0) { }
261  else if(PreconditionerType_==1) { List.set("relaxation: type", "Jacobi" ); }
262  else if(PreconditionerType_==2) { List.set("relaxation: type", "Gauss-Seidel" ); }
263  else if(PreconditionerType_==3) { List.set("relaxation: type", "symmetric Gauss-Seidel"); }
264  if(BlockSize_==1) {
265  IfpackPrec_ = Teuchos::rcp( new Ifpack_PointRelaxation(&*Matrix_) );
266  }
267  else {
268  IfpackPrec_ = Teuchos::rcp( new Ifpack_BlockRelaxation< Ifpack_DenseContainer > (&*Matrix_) );
269  int NumRows;
270  if(IsRowMatrix_==true) {
271  NumRows = Matrix_->NumMyRows();
272  }
273  else {
274  long long NumRows_LL = Operator_->OperatorDomainMap().NumGlobalElements64();
275  if(NumRows_LL > std::numeric_limits<int>::max())
276  throw "Ifpack_Krylov::Compute: NumGlobalElements don't fit an int";
277  else
278  NumRows = static_cast<int>(NumRows_LL);
279  }
280  List.set("partitioner: type", "linear");
281  List.set("partitioner: local parts", NumRows/BlockSize_);
282  }
283  if(PreconditionerType_>0) {
284  IfpackPrec_ -> SetParameters(List);
285  IfpackPrec_ -> Initialize();
286  IfpackPrec_ -> Compute();
287  AztecSolver_ -> SetPrecOperator(&*IfpackPrec_);
288  }
289 #else
290  using std::cout;
291  using std::endl;
292 
293  cout << "You need to configure IFPACK with support for AztecOO" << endl;
294  cout << "to use this preconditioner. This may require --enable-aztecoo" << endl;
295  cout << "in your configure script." << endl;
296  IFPACK_CHK_ERR(-1);
297 #endif
298 
299  // reset values
300  IsComputed_ = false;
301  Condest_ = -1.0;
302  ++NumCompute_;
303  ComputeTime_ += Time_->ElapsedTime();
304  IsComputed_ = true;
305 
306  return(0);
307 }
308 
309 //==============================================================================
310 std::ostream& Ifpack_Krylov::Print(std::ostream & os) const
311 {
312  using std::endl;
313 
314  if (!Comm().MyPID()) {
315  os << endl;
316  os << "================================================================================" << endl;
317  os << "Ifpack_Krylov" << endl;
318  os << "Number of iterations = " << Iterations_ << endl;
319  os << "Residual Tolerance = " << Tolerance_ << endl;
320  os << "Solver type (O for CG, 1 for GMRES) = " << SolverType_ << endl;
321  os << "Preconditioner type = " << PreconditionerType_ << endl;
322  os << "(0 for none, 1 for Jacobi, 2 for GS, 3 for SGS )" << endl;
323  os << "Condition number estimate = " << Condest() << endl;
324  os << "Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
325  if (ZeroStartingSolution_)
326  os << "Using zero starting solution" << endl;
327  else
328  os << "Using input starting solution" << endl;
329  os << endl;
330  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
331  os << "----- ------- -------------- ------------ --------" << endl;
332  os << "Initialize() " << std::setw(5) << NumInitialize_
333  << " " << std::setw(15) << InitializeTime_
334  << " 0.0 0.0" << endl;
335  os << "Compute() " << std::setw(5) << NumCompute_
336  << " " << std::setw(15) << ComputeTime_
337  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
338  if (ComputeTime_ != 0.0)
339  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
340  else
341  os << " " << std::setw(15) << 0.0 << endl;
342  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
343  << " " << std::setw(15) << ApplyInverseTime_
344  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
345  if (ApplyInverseTime_ != 0.0)
346  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
347  else
348  os << " " << std::setw(15) << 0.0 << endl;
349  os << "================================================================================" << endl;
350  os << endl;
351  }
352 
353  return(os);
354 }
355 
356 //==============================================================================
357 double Ifpack_Krylov::
358 Condest(const Ifpack_CondestType CT,
359  const int MaxIters, const double Tol,
360  Epetra_RowMatrix* Matrix_in)
361 {
362  if (!IsComputed()) // cannot compute right now
363  return(-1.0);
364 
365  // always computes it. Call Condest() with no parameters to get
366  // the previous estimate.
367  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
368 
369  return(Condest_);
370 }
371 
372 //==============================================================================
373 void Ifpack_Krylov::SetLabel()
374 {
375  Label_ = "IFPACK (Krylov smoother), iterations=" + Ifpack_toString(Iterations_);
376 }
377 
378 //==============================================================================
379 int Ifpack_Krylov::
381 {
382 
383  if (!IsComputed())
384  IFPACK_CHK_ERR(-3);
385 
386  if (Iterations_ == 0)
387  return 0;
388 
389  int nVec = X.NumVectors();
390  if (nVec != Y.NumVectors())
391  IFPACK_CHK_ERR(-2);
392 
393  Time_->ResetStartTime();
394 
395  // AztecOO gives X and Y pointing to the same memory location,
396  // need to create an auxiliary vector, Xcopy
397  Teuchos::RCP<Epetra_MultiVector> Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
398  if(ZeroStartingSolution_==true) {
399  Y.PutScalar(0.0);
400  }
401 
402 #ifdef HAVE_IFPACK_AZTECOO
403  AztecSolver_ -> SetLHS(&Y);
404  AztecSolver_ -> SetRHS(&*Xcopy);
405  AztecSolver_ -> Iterate(Iterations_,Tolerance_);
406 #else
407  using std::cout;
408  using std::endl;
409 
410  cout << "You need to configure IFPACK with support for AztecOO" << endl;
411  cout << "to use this preconditioner. This may require --enable-aztecoo" << endl;
412  cout << "in your configure script." << endl;
413  IFPACK_CHK_ERR(-1);
414 #endif
415 
416  // Flops are updated in each of the following.
417  ++NumApplyInverse_;
418  ApplyInverseTime_ += Time_->ElapsedTime();
419  return(0);
420 }
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix&#39;s.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual int Compute()
Computes the preconditioners.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Ifpack_PointRelaxation: a class to define point relaxation preconditioners of for Epetra_RowMatrix&#39;s...
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.