IFPACK  Development
 All Classes 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),
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),
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 //==============================================================================
159 const Epetra_Comm& Ifpack_Krylov::Comm() const
160 {
161  return(Operator_->Comm());
162 }
163 
164 //==============================================================================
165 const Epetra_Map& Ifpack_Krylov::OperatorDomainMap() const
166 {
167  return(Operator_->OperatorDomainMap());
168 }
169 
170 //==============================================================================
171 const Epetra_Map& Ifpack_Krylov::OperatorRangeMap() const
172 {
173  return(Operator_->OperatorRangeMap());
174 }
175 
176 //==============================================================================
177 int Ifpack_Krylov::
178 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
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  cout << "You need to configure IFPACK with support for AztecOO" << endl;
291  cout << "to use this preconditioner. This may require --enable-aztecoo" << endl;
292  cout << "in your configure script." << endl;
293  IFPACK_CHK_ERR(-1);
294 #endif
295 
296  // reset values
297  IsComputed_ = false;
298  Condest_ = -1.0;
299  ++NumCompute_;
300  ComputeTime_ += Time_->ElapsedTime();
301  IsComputed_ = true;
302 
303  return(0);
304 }
305 
306 //==============================================================================
307 ostream& Ifpack_Krylov::Print(ostream & os) const
308 {
309 
310  if (!Comm().MyPID()) {
311  os << endl;
312  os << "================================================================================" << endl;
313  os << "Ifpack_Krylov" << endl;
314  os << "Number of iterations = " << Iterations_ << endl;
315  os << "Residual Tolerance = " << Tolerance_ << endl;
316  os << "Solver type (O for CG, 1 for GMRES) = " << SolverType_ << endl;
317  os << "Preconditioner type = " << PreconditionerType_ << endl;
318  os << "(0 for none, 1 for Jacobi, 2 for GS, 3 for SGS )" << endl;
319  os << "Condition number estimate = " << Condest() << endl;
320  os << "Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
321  if (ZeroStartingSolution_)
322  os << "Using zero starting solution" << endl;
323  else
324  os << "Using input starting solution" << endl;
325  os << endl;
326  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
327  os << "----- ------- -------------- ------------ --------" << endl;
328  os << "Initialize() " << std::setw(5) << NumInitialize_
329  << " " << std::setw(15) << InitializeTime_
330  << " 0.0 0.0" << endl;
331  os << "Compute() " << std::setw(5) << NumCompute_
332  << " " << std::setw(15) << ComputeTime_
333  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
334  if (ComputeTime_ != 0.0)
335  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
336  else
337  os << " " << std::setw(15) << 0.0 << endl;
338  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
339  << " " << std::setw(15) << ApplyInverseTime_
340  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
341  if (ApplyInverseTime_ != 0.0)
342  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
343  else
344  os << " " << std::setw(15) << 0.0 << endl;
345  os << "================================================================================" << endl;
346  os << endl;
347  }
348 
349  return(os);
350 }
351 
352 //==============================================================================
353 double Ifpack_Krylov::
354 Condest(const Ifpack_CondestType CT,
355  const int MaxIters, const double Tol,
356  Epetra_RowMatrix* Matrix_in)
357 {
358  if (!IsComputed()) // cannot compute right now
359  return(-1.0);
360 
361  // always computes it. Call Condest() with no parameters to get
362  // the previous estimate.
363  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
364 
365  return(Condest_);
366 }
367 
368 //==============================================================================
369 void Ifpack_Krylov::SetLabel()
370 {
371  Label_ = "IFPACK (Krylov smoother), iterations=" + Ifpack_toString(Iterations_);
372 }
373 
374 //==============================================================================
375 int Ifpack_Krylov::
376 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
377 {
378 
379  if (!IsComputed())
380  IFPACK_CHK_ERR(-3);
381 
382  if (Iterations_ == 0)
383  return 0;
384 
385  int nVec = X.NumVectors();
386  if (nVec != Y.NumVectors())
387  IFPACK_CHK_ERR(-2);
388 
389  Time_->ResetStartTime();
390 
391  // AztecOO gives X and Y pointing to the same memory location,
392  // need to create an auxiliary vector, Xcopy
393  Teuchos::RCP<Epetra_MultiVector> Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
394  if(ZeroStartingSolution_==true) {
395  Y.PutScalar(0.0);
396  }
397 
398 #ifdef HAVE_IFPACK_AZTECOO
399  AztecSolver_ -> SetLHS(&Y);
400  AztecSolver_ -> SetRHS(&*Xcopy);
401  AztecSolver_ -> Iterate(Iterations_,Tolerance_);
402 #else
403  cout << "You need to configure IFPACK with support for AztecOO" << endl;
404  cout << "to use this preconditioner. This may require --enable-aztecoo" << endl;
405  cout << "in your configure script." << endl;
406  IFPACK_CHK_ERR(-1);
407 #endif
408 
409  // Flops are updated in each of the following.
410  ++NumApplyInverse_;
411  ApplyInverseTime_ += Time_->ElapsedTime();
412  return(0);
413 }
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix'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 ostream & Print(ostream &os) const
Prints object to an output stream.
string Ifpack_toString(const int &x)
Convertes an integer to string.
virtual int Compute()
Computes the preconditioners.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
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'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.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.