IFPACK  Development
 All Classes Files Functions Variables Enumerations Friends Pages
Ifpack_Polynomial.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_Polynomial.h"
53 #include "Ifpack_Utils.h"
54 #include "Ifpack_Condest.h"
55 #include "Teuchos_LAPACK.hpp"
56 #include "Teuchos_SerialDenseMatrix.hpp"
57 #include <complex>
58 #ifdef HAVE_IFPACK_AZTECOO
59 #include "Ifpack_DiagPreconditioner.h"
60 #include "AztecOO.h"
61 #endif
62 
63 #ifdef HAVE_IFPACK_EPETRAEXT
64 #include "Epetra_CrsMatrix.h"
65 #include "EpetraExt_PointToBlockDiagPermute.h"
66 #endif
67 
68 
69 #define ABS(x) ((x)>0?(x):-(x))
70 
71 // Helper function for normal equations
72 inline void Apply_Transpose(Teuchos::RCP<const Epetra_Operator> Operator_,const Epetra_MultiVector &X,Epetra_MultiVector &Y){
73  Epetra_Operator * Operator=const_cast<Epetra_Operator*>(&*Operator_);
74  Operator->SetUseTranspose(true);
75  Operator->Apply(X,Y);
76  Operator->SetUseTranspose(false);
77 }
78 
79 
80 
81 
82 //==============================================================================
83 // NOTE: any change to the default values should be committed to the other
84 // constructor as well.
86 Ifpack_Polynomial(const Epetra_Operator* Operator) :
87  IsInitialized_(false),
88  IsComputed_(false),
89  IsIndefinite_(false),
90  IsComplex_(false),
91  NumInitialize_(0),
92  NumCompute_(0),
93  NumApplyInverse_(0),
94  InitializeTime_(0.0),
95  ComputeTime_(0.0),
96  ApplyInverseTime_(0.0),
97  ComputeFlops_(0.0),
98  ApplyInverseFlops_(0.0),
99  PolyDegree_(3),
100  LSPointsReal_(10),
101  LSPointsImag_(10),
102  UseTranspose_(false),
103  Condest_(-1.0),
104  ComputeCondest_(false),
105  RealEigRatio_(10.0),
106  ImagEigRatio_(10.0),
107  Label_(),
108  LambdaRealMin_(0.0),
109  LambdaRealMax_(-1.0),
110  LambdaImagMin_(0.0),
111  LambdaImagMax_(0.0),
112  MinDiagonalValue_(0.0),
113  NumMyRows_(0),
114  NumMyNonzeros_(0),
115  NumGlobalRows_(0),
116  NumGlobalNonzeros_(0),
117  Operator_(Teuchos::rcp(Operator,false)),
118  UseBlockMode_(false),
119  SolveNormalEquations_(false),
120  IsRowMatrix_(false),
121  ZeroStartingSolution_(true)
122 {
123 }
124 
125 //==============================================================================
126 // NOTE: This constructor has been introduced because SWIG does not appear
127 // to appreciate dynamic_cast. An instruction of type
128 // Matrix_ = dynamic_cast<const Epetra_RowMatrix*> in the
129 // other construction does not work in PyTrilinos -- of course
130 // it does in any C++ code (for an Epetra_Operator that is also
131 // an Epetra_RowMatrix).
132 //
133 // FIXME: move declarations into a separate method?
135 Ifpack_Polynomial(const Epetra_RowMatrix* Operator) :
136  IsInitialized_(false),
137  IsComputed_(false),
138  IsIndefinite_(false),
139  IsComplex_(false),
140  NumInitialize_(0),
141  NumCompute_(0),
142  NumApplyInverse_(0),
143  InitializeTime_(0.0),
144  ComputeTime_(0.0),
145  ApplyInverseTime_(0.0),
146  ComputeFlops_(0.0),
147  ApplyInverseFlops_(0.0),
148  PolyDegree_(3),
149  LSPointsReal_(10),
150  LSPointsImag_(10),
151  UseTranspose_(false),
152  Condest_(-1.0),
153  ComputeCondest_(false),
154  RealEigRatio_(10.0),
155  ImagEigRatio_(10.0),
156  EigMaxIters_(10),
157  Label_(),
158  LambdaRealMin_(0.0),
159  LambdaRealMax_(-1.0),
160  LambdaImagMin_(0.0),
161  LambdaImagMax_(0.0),
162  MinDiagonalValue_(0.0),
163  NumMyRows_(0),
164  NumMyNonzeros_(0),
165  NumGlobalRows_(0),
166  NumGlobalNonzeros_(0),
167  Operator_(Teuchos::rcp(Operator,false)),
168  Matrix_(Teuchos::rcp(Operator,false)),
169  UseBlockMode_(false),
170  SolveNormalEquations_(false),
171  IsRowMatrix_(true),
172  ZeroStartingSolution_(true)
173 {
174 }
175 
176 //==============================================================================
177 int Ifpack_Polynomial::SetParameters(Teuchos::ParameterList& List)
178 {
179 
180  RealEigRatio_ = List.get("polynomial: real eigenvalue ratio", RealEigRatio_);
181  ImagEigRatio_ = List.get("polynomial: imag eigenvalue ratio", ImagEigRatio_);
182  LambdaRealMin_ = List.get("polynomial: min real part", LambdaRealMin_);
183  LambdaRealMax_ = List.get("polynomial: max real part", LambdaRealMax_);
184  LambdaImagMin_ = List.get("polynomial: min imag part", LambdaImagMin_);
185  LambdaImagMax_ = List.get("polynomial: max imag part", LambdaImagMax_);
186  PolyDegree_ = List.get("polynomial: degree",PolyDegree_);
187  LSPointsReal_ = List.get("polynomial: real interp points",LSPointsReal_);
188  LSPointsImag_ = List.get("polynomial: imag interp points",LSPointsImag_);
189  IsIndefinite_ = List.get("polynomial: indefinite",IsIndefinite_);
190  IsComplex_ = List.get("polynomial: complex",IsComplex_);
191  MinDiagonalValue_ = List.get("polynomial: min diagonal value",
192  MinDiagonalValue_);
193  ZeroStartingSolution_ = List.get("polynomial: zero starting solution",
194  ZeroStartingSolution_);
195 
196  Epetra_Vector* ID = List.get("polynomial: operator inv diagonal",
197  (Epetra_Vector*)0);
198  EigMaxIters_ = List.get("polynomial: eigenvalue max iterations",EigMaxIters_);
199 
200 #ifdef HAVE_IFPACK_EPETRAEXT
201  // This is *always* false if EpetraExt isn't enabled
202  UseBlockMode_ = List.get("polynomial: use block mode",UseBlockMode_);
203  if(!List.isParameter("polynomial: block list")) UseBlockMode_=false;
204  else{
205  BlockList_ = List.get("polynomial: block list",BlockList_);
206 
207  // Since we know we're doing a matrix inverse, clobber the block list
208  // w/"invert" if it's set to multiply
209  Teuchos::ParameterList Blist;
210  Blist=BlockList_.get("blockdiagmatrix: list",Blist);
211  string dummy("invert");
212  string ApplyMode=Blist.get("apply mode",dummy);
213  if(ApplyMode==string("multiply")){
214  Blist.set("apply mode","invert");
215  BlockList_.set("blockdiagmatrix: list",Blist);
216  }
217  }
218 #endif
219 
220  SolveNormalEquations_ = List.get("polynomial: solve normal equations",SolveNormalEquations_);
221 
222  if (ID != 0)
223  {
224  InvDiagonal_ = Teuchos::rcp( new Epetra_Vector(*ID) );
225  }
226 
227  SetLabel();
228 
229  return(0);
230 }
231 
232 //==============================================================================
233 const Epetra_Comm& Ifpack_Polynomial::Comm() const
234 {
235  return(Operator_->Comm());
236 }
237 
238 //==============================================================================
239 const Epetra_Map& Ifpack_Polynomial::OperatorDomainMap() const
240 {
241  return(Operator_->OperatorDomainMap());
242 }
243 
244 //==============================================================================
245 const Epetra_Map& Ifpack_Polynomial::OperatorRangeMap() const
246 {
247  return(Operator_->OperatorRangeMap());
248 }
249 
250 //==============================================================================
252 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
253 {
254  if (IsComputed() == false)
255  IFPACK_CHK_ERR(-3);
256 
257  if (X.NumVectors() != Y.NumVectors())
258  IFPACK_CHK_ERR(-2);
259 
260  if (IsRowMatrix_)
261  {
262  IFPACK_CHK_ERR(Matrix_->Multiply(UseTranspose(),X,Y));
263  }
264  else
265  {
266  IFPACK_CHK_ERR(Operator_->Apply(X,Y));
267  }
268 
269  return(0);
270 }
271 
272 //==============================================================================
274 {
275  IsInitialized_ = false;
276 
277  if (Operator_ == Teuchos::null)
278  IFPACK_CHK_ERR(-2);
279 
280  if (Time_ == Teuchos::null)
281  Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
282 
283  if (IsRowMatrix_)
284  {
285  if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
286  IFPACK_CHK_ERR(-2); // only square matrices
287 
288  NumMyRows_ = Matrix_->NumMyRows();
289  NumMyNonzeros_ = Matrix_->NumMyNonzeros();
290  NumGlobalRows_ = Matrix_->NumGlobalRows64();
291  NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
292  }
293  else
294  {
295  if (Operator_->OperatorDomainMap().NumGlobalElements64() !=
296  Operator_->OperatorRangeMap().NumGlobalElements64())
297  IFPACK_CHK_ERR(-2); // only square operators
298  }
299 
300  ++NumInitialize_;
301  InitializeTime_ += Time_->ElapsedTime();
302  IsInitialized_ = true;
303  return(0);
304 }
305 
306 //==============================================================================
308 {
309  if (!IsInitialized())
310  IFPACK_CHK_ERR(Initialize());
311 
312  Time_->ResetStartTime();
313 
314  // reset values
315  IsComputed_ = false;
316  Condest_ = -1.0;
317 
318  if (PolyDegree_ <= 0)
319  IFPACK_CHK_ERR(-2); // at least one application
320 
321 #ifdef HAVE_IFPACK_EPETRAEXT
322  // Check to see if we can run in block mode
323  if(IsRowMatrix_ && InvDiagonal_ == Teuchos::null && UseBlockMode_){
324  const Epetra_CrsMatrix *CrsMatrix=dynamic_cast<const Epetra_CrsMatrix*>(&*Matrix_);
325 
326  // If we don't have CrsMatrix, we can't use the block preconditioner
327  if(!CrsMatrix) UseBlockMode_=false;
328  else{
329  int ierr;
330  InvBlockDiagonal_=Teuchos::rcp(new EpetraExt_PointToBlockDiagPermute(*CrsMatrix));
331  if(InvBlockDiagonal_==Teuchos::null) IFPACK_CHK_ERR(-6);
332 
333  ierr=InvBlockDiagonal_->SetParameters(BlockList_);
334  if(ierr) IFPACK_CHK_ERR(-7);
335 
336  ierr=InvBlockDiagonal_->Compute();
337  if(ierr) IFPACK_CHK_ERR(-8);
338  }
339  }
340 #endif
341  if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_)
342  {
343  InvDiagonal_ = Teuchos::rcp( new Epetra_Vector(Matrix().Map()) );
344 
345  if (InvDiagonal_ == Teuchos::null)
346  IFPACK_CHK_ERR(-5);
347 
348  IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*InvDiagonal_));
349 
350  // Inverse diagonal elements
351  // Replace zeros with 1.0
352  for (int i = 0 ; i < NumMyRows_ ; ++i) {
353  double diag = (*InvDiagonal_)[i];
354  if (IFPACK_ABS(diag) < MinDiagonalValue_)
355  (*InvDiagonal_)[i] = MinDiagonalValue_;
356  else
357  (*InvDiagonal_)[i] = 1.0 / diag;
358  }
359  }
360 
361  // Automatically compute maximum eigenvalue estimate of D^{-1}A if user hasn't provided one
362  double lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max;
363  if (LambdaRealMax_ == -1) {
364  //PowerMethod(Matrix(), *InvDiagonal_, EigMaxIters_, lambda_max);
365  GMRES(Matrix(), *InvDiagonal_, EigMaxIters_, lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max);
366  LambdaRealMin_=lambda_real_min; LambdaImagMin_=lambda_imag_min;
367  LambdaRealMax_=lambda_real_max; LambdaImagMax_=lambda_imag_max;
368  //std::cout<<"LambdaRealMin: "<<LambdaRealMin_<<std::endl;
369  //std::cout<<"LambdaRealMax: "<<LambdaRealMax_<<std::endl;
370  //std::cout<<"LambdaImagMin: "<<LambdaImagMin_<<std::endl;
371  //std::cout<<"LambdaImagMax: "<<LambdaImagMax_<<std::endl;
372  }
373 
374  // find least squares polynomial for (LSPointsReal_*LSPointsImag_) zeros
375  // on a rectangle in the complex plane defined as
376  // [LambdaRealMin_,LambdaRealMax_] x [LambdaImagMin_,LambdaImagMax_]
377 
378  std::complex<double> zero(0.0,0.0);
379  std::complex<double> one(1.0,0.0);
380 
381  // Compute points in complex plane
382  double lenx = LambdaRealMax_-LambdaRealMin_;
383  int nx = ceil(lenx*((double) LSPointsReal_));
384  if (nx<2) { nx = 2; }
385  double hx = lenx/((double) nx);
386  std::vector<double> xs;
387  if(abs(lenx)>1.0e-8) {
388  for( int pt=0; pt<=nx; pt++ ) {
389  xs.push_back(hx*pt+LambdaRealMin_);
390  }
391  }
392  else {
393  xs.push_back(LambdaRealMax_);
394  nx=1;
395  }
396  double leny = LambdaImagMax_-LambdaImagMin_;
397  int ny = ceil(leny*((double) LSPointsImag_));
398  if (ny<2) { ny = 2; }
399  double hy = leny/((double) ny);
400  std::vector<double> ys;
401  if(abs(leny)>1.0e-8) {
402  for( int pt=0; pt<=ny; pt++ ) {
403  ys.push_back(hy*pt+LambdaImagMin_);
404  }
405  }
406  else {
407  ys.push_back(LambdaImagMax_);
408  ny=1;
409  }
410  std::vector< std::complex<double> > cpts;
411  for( int jj=0; jj<ny; jj++ ) {
412  for( int ii=0; ii<nx; ii++ ) {
413  std::complex<double> cpt(xs[ii],ys[jj]);
414  cpts.push_back(cpt);
415  }
416  }
417  cpts.push_back(zero);
418 
419 #ifdef HAVE_TEUCHOS_COMPLEX
420 
421  // Construct overdetermined Vandermonde matrix
422  Teuchos::SerialDenseMatrix<int, std::complex<double> > Vmatrix(cpts.size(),PolyDegree_+1);
423  Vmatrix.putScalar(zero);
424  for (int jj = 0; jj <= PolyDegree_; ++jj) {
425  for (int ii = 0; ii < static_cast<int> (cpts.size ()) - 1; ++ii) {
426  if (jj > 0) {
427  Vmatrix(ii,jj) = pow(cpts[ii],jj);
428  }
429  else {
430  Vmatrix(ii,jj) = one;
431  }
432  }
433  }
434  Vmatrix(cpts.size()-1,0)=one;
435 
436  // Right hand side: all zero except last entry
437  Teuchos::SerialDenseMatrix< int,std::complex<double> > RHS(cpts.size(),1);
438  RHS.putScalar(zero);
439  RHS(cpts.size()-1,0)=one;
440 
441  // Solve least squares problem using LAPACK
442  Teuchos::LAPACK< int, std::complex<double> > lapack;
443  const int N = Vmatrix.numCols();
444  Teuchos::Array<double> singularValues(N);
445  Teuchos::Array<double> rwork(1);
446  rwork.resize (std::max (1, 5 * N));
447  std::complex<double> lworkScalar(1.0,0.0);
448  int info = 0;
449  lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
450  Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
451  &lworkScalar, -1, &info);
452  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
453  "_GELSS workspace query returned INFO = "
454  << info << " != 0.");
455  const int lwork = static_cast<int> (real(lworkScalar));
456  TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
457  "_GELSS workspace query returned LWORK = "
458  << lwork << " < 0.");
459  // Allocate workspace. Size > 0 means &work[0] makes sense.
460  Teuchos::Array< std::complex<double> > work (std::max (1, lwork));
461  // Solve the least-squares problem.
462  lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
463  Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
464  &work[0], lwork, &info);
465 
466  coeff_.resize(PolyDegree_+1);
467  std::complex<double> c0=RHS(0,0);
468  for(int ii=0; ii<=PolyDegree_; ii++) {
469  // test that the imaginary part is nonzero
470  //TEUCHOS_TEST_FOR_EXCEPTION(abs(imag(RHS(ii,0))) > 1e-8, std::logic_error,
471  // "imaginary part of polynomial coefficients is nonzero! coeff = "
472  // << RHS(ii,0));
473  coeff_[ii]=real(RHS(ii,0)/c0);
474  //std::cout<<"coeff["<<ii<<"]="<<coeff_[ii]<<std::endl;
475  }
476 
477 #else
478 
479  // Construct overdetermined Vandermonde matrix
480  Teuchos::SerialDenseMatrix< int, double > Vmatrix(xs.size()+1,PolyDegree_+1);
481  Vmatrix.putScalar(0.0);
482  for( int jj=0; jj<=PolyDegree_; jj++) {
483  for( int ii=0; ii<xs.size(); ii++) {
484  if(jj>0) {
485  Vmatrix(ii,jj)=pow(xs[ii],jj);
486  }
487  else {
488  Vmatrix(ii,jj)=1.0;
489  }
490  }
491  }
492  Vmatrix(xs.size(),0)=1.0;
493 
494  // Right hand side: all zero except last entry
495  Teuchos::SerialDenseMatrix< int, double > RHS(xs.size()+1,1);
496  RHS.putScalar(0.0);
497  RHS(xs.size(),0)=1.0;
498 
499  // Solve least squares problem using LAPACK
500  Teuchos::LAPACK< int, double > lapack;
501  const int N = Vmatrix.numCols();
502  Teuchos::Array<double> singularValues(N);
503  Teuchos::Array<double> rwork(1);
504  rwork.resize (std::max (1, 5 * N));
505  double lworkScalar(1.0);
506  int info = 0;
507  lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
508  Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
509  &lworkScalar, -1, &info);
510  TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
511  "_GELSS workspace query returned INFO = "
512  << info << " != 0.");
513  const int lwork = static_cast<int> (lworkScalar);
514  TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
515  "_GELSS workspace query returned LWORK = "
516  << lwork << " < 0.");
517  // Allocate workspace. Size > 0 means &work[0] makes sense.
518  Teuchos::Array< double > work (std::max (1, lwork));
519  // Solve the least-squares problem.
520  lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
521  Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
522  &work[0], lwork, &info);
523 
524  coeff_.resize(PolyDegree_+1);
525  double c0=RHS(0,0);
526  for(int ii=0; ii<=PolyDegree_; ii++) {
527  // test that the imaginary part is nonzero
528  //TEUCHOS_TEST_FOR_EXCEPTION(abs(imag(RHS(ii,0))) > 1e-8, std::logic_error,
529  // "imaginary part of polynomial coefficients is nonzero! coeff = "
530  // << RHS(ii,0));
531  coeff_[ii]=RHS(ii,0)/c0;
532  }
533 
534 #endif
535 
536 #ifdef IFPACK_FLOPCOUNTERS
537  ComputeFlops_ += NumMyRows_;
538 #endif
539 
540  ++NumCompute_;
541  ComputeTime_ += Time_->ElapsedTime();
542  IsComputed_ = true;
543 
544  return(0);
545 }
546 
547 //==============================================================================
548 ostream& Ifpack_Polynomial::Print(ostream & os) const
549 {
550 
551  double MyMinVal, MyMaxVal;
552  double MinVal, MaxVal;
553 
554  if (IsComputed_) {
555  InvDiagonal_->MinValue(&MyMinVal);
556  InvDiagonal_->MaxValue(&MyMaxVal);
557  Comm().MinAll(&MyMinVal,&MinVal,1);
558  Comm().MinAll(&MyMaxVal,&MaxVal,1);
559  }
560 
561  if (!Comm().MyPID()) {
562  os << endl;
563  os << "================================================================================" << endl;
564  os << "Ifpack_Polynomial" << endl;
565  os << "Degree of polynomial = " << PolyDegree_ << endl;
566  os << "Condition number estimate = " << Condest() << endl;
567  os << "Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
568  if (IsComputed_) {
569  os << "Minimum value on stored inverse diagonal = " << MinVal << endl;
570  os << "Maximum value on stored inverse diagonal = " << MaxVal << endl;
571  }
572  if (ZeroStartingSolution_)
573  os << "Using zero starting solution" << endl;
574  else
575  os << "Using input starting solution" << endl;
576  os << endl;
577  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
578  os << "----- ------- -------------- ------------ --------" << endl;
579  os << "Initialize() " << std::setw(5) << NumInitialize_
580  << " " << std::setw(15) << InitializeTime_
581  << " 0.0 0.0" << endl;
582  os << "Compute() " << std::setw(5) << NumCompute_
583  << " " << std::setw(15) << ComputeTime_
584  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
585  if (ComputeTime_ != 0.0)
586  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
587  else
588  os << " " << std::setw(15) << 0.0 << endl;
589  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
590  << " " << std::setw(15) << ApplyInverseTime_
591  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
592  if (ApplyInverseTime_ != 0.0)
593  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
594  else
595  os << " " << std::setw(15) << 0.0 << endl;
596  os << "================================================================================" << endl;
597  os << endl;
598  }
599 
600  return(os);
601 }
602 
603 //==============================================================================
604 double Ifpack_Polynomial::
605 Condest(const Ifpack_CondestType CT,
606  const int MaxIters, const double Tol,
607  Epetra_RowMatrix* Matrix_in)
608 {
609  if (!IsComputed()) // cannot compute right now
610  return(-1.0);
611 
612  // always computes it. Call Condest() with no parameters to get
613  // the previous estimate.
614  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
615 
616  return(Condest_);
617 }
618 
619 //==============================================================================
620 void Ifpack_Polynomial::SetLabel()
621 {
622  Label_ = "IFPACK (Least squares polynomial), degree=" + Ifpack_toString(PolyDegree_);
623 }
624 
625 //==============================================================================
627 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
628 {
629 
630  if (!IsComputed())
631  IFPACK_CHK_ERR(-3);
632 
633  if (PolyDegree_ == 0)
634  return 0;
635 
636  int nVec = X.NumVectors();
637  if (nVec != Y.NumVectors())
638  IFPACK_CHK_ERR(-2);
639 
640  Time_->ResetStartTime();
641 
642  Epetra_MultiVector Xcopy(X);
643  if(ZeroStartingSolution_==true) {
644  Y.PutScalar(0.0);
645  }
646 
647  // mfh 20 Mar 2014: IBD never gets used, so I'm commenting out the
648  // following lines of code in order to forestall build warnings.
649 // #ifdef HAVE_IFPACK_EPETRAEXT
650 // EpetraExt_PointToBlockDiagPermute* IBD=0;
651 // if (UseBlockMode_) IBD=&*InvBlockDiagonal_;
652 // #endif
653 
654  Y.Update(-coeff_[1], Xcopy, 1.0);
655  for (int ii = 2; ii < static_cast<int> (coeff_.size ()); ++ii) {
656  const Epetra_MultiVector V(Xcopy);
657  Operator_->Apply(V,Xcopy);
658  Xcopy.Multiply(1.0, *InvDiagonal_, Xcopy, 0.0);
659  // Update Y
660  Y.Update(-coeff_[ii], Xcopy, 1.0);
661  }
662 
663  // Flops are updated in each of the following.
664  ++NumApplyInverse_;
665  ApplyInverseTime_ += Time_->ElapsedTime();
666  return(0);
667 }
668 
669 //==============================================================================
671 PowerMethod(const Epetra_Operator& Operator,
672  const Epetra_Vector& InvPointDiagonal,
673  const int MaximumIterations,
674  double& lambda_max)
675 {
676  // this is a simple power method
677  lambda_max = 0.0;
678  double RQ_top, RQ_bottom, norm;
679  Epetra_Vector x(Operator.OperatorDomainMap());
680  Epetra_Vector y(Operator.OperatorRangeMap());
681  x.Random();
682  x.Norm2(&norm);
683  if (norm == 0.0) IFPACK_CHK_ERR(-1);
684 
685  x.Scale(1.0 / norm);
686 
687  for (int iter = 0; iter < MaximumIterations; ++iter)
688  {
689  Operator.Apply(x, y);
690  IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
691  IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
692  IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
693  lambda_max = RQ_top / RQ_bottom;
694  IFPACK_CHK_ERR(y.Norm2(&norm));
695  if (norm == 0.0) IFPACK_CHK_ERR(-1);
696  IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
697  }
698 
699  return(0);
700 }
701 
702 //==============================================================================
704 CG(const Epetra_Operator& Operator,
705  const Epetra_Vector& InvPointDiagonal,
706  const int MaximumIterations,
707  double& lambda_min, double& lambda_max)
708 {
709 #ifdef HAVE_IFPACK_AZTECOO
710  Epetra_Vector x(Operator.OperatorDomainMap());
711  Epetra_Vector y(Operator.OperatorRangeMap());
712  x.Random();
713  y.PutScalar(0.0);
714 
715  Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
716  AztecOO solver(LP);
717  solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
718  solver.SetAztecOption(AZ_output, AZ_none);
719 
720  Ifpack_DiagPreconditioner diag(Operator.OperatorDomainMap(),
721  Operator.OperatorRangeMap(),
722  InvPointDiagonal);
723  solver.SetPrecOperator(&diag);
724  solver.Iterate(MaximumIterations, 1e-10);
725 
726  const double* status = solver.GetAztecStatus();
727 
728  lambda_min = status[AZ_lambda_min];
729  lambda_max = status[AZ_lambda_max];
730 
731  return(0);
732 #else
733  cout << "You need to configure IFPACK with support for AztecOO" << endl;
734  cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
735  cout << "in your configure script." << endl;
736  IFPACK_CHK_ERR(-1);
737 #endif
738 }
739 
740 //==============================================================================
741 #ifdef HAVE_IFPACK_EPETRAEXT
743 PowerMethod(const int MaximumIterations, double& lambda_max)
744 {
745 
746  if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
747  // this is a simple power method
748  lambda_max = 0.0;
749  double RQ_top, RQ_bottom, norm;
750  Epetra_Vector x(Operator_->OperatorDomainMap());
751  Epetra_Vector y(Operator_->OperatorRangeMap());
752  Epetra_Vector z(Operator_->OperatorRangeMap());
753  x.Random();
754  x.Norm2(&norm);
755  if (norm == 0.0) IFPACK_CHK_ERR(-1);
756 
757  x.Scale(1.0 / norm);
758 
759  for (int iter = 0; iter < MaximumIterations; ++iter)
760  {
761  Operator_->Apply(x, z);
762  InvBlockDiagonal_->ApplyInverse(z,y);
763  if(SolveNormalEquations_){
764  InvBlockDiagonal_->ApplyInverse(y,z);
765  Apply_Transpose(Operator_,z, y);
766  }
767 
768  IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
769  IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
770  lambda_max = RQ_top / RQ_bottom;
771  IFPACK_CHK_ERR(y.Norm2(&norm));
772  if (norm == 0.0) IFPACK_CHK_ERR(-1);
773  IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
774  }
775 
776  return(0);
777 }
778 #endif
779 
780 //==============================================================================
781 #ifdef HAVE_IFPACK_EPETRAEXT
783 CG(const int MaximumIterations,
784  double& lambda_min, double& lambda_max)
785 {
786  IFPACK_CHK_ERR(-1);// NTS: This always seems to yield errors in AztecOO, ergo,
787  // I turned it off.
788 
789  if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
790 
791 #ifdef HAVE_IFPACK_AZTECOO
792  Epetra_Vector x(Operator_->OperatorDomainMap());
793  Epetra_Vector y(Operator_->OperatorRangeMap());
794  x.Random();
795  y.PutScalar(0.0);
796  Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
797 
798  AztecOO solver(LP);
799  solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
800  solver.SetAztecOption(AZ_output, AZ_none);
801 
802  solver.SetPrecOperator(&*InvBlockDiagonal_);
803  solver.Iterate(MaximumIterations, 1e-10);
804 
805  const double* status = solver.GetAztecStatus();
806 
807  lambda_min = status[AZ_lambda_min];
808  lambda_max = status[AZ_lambda_max];
809 
810  return(0);
811 #else
812  cout << "You need to configure IFPACK with support for AztecOO" << endl;
813  cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
814  cout << "in your configure script." << endl;
815  IFPACK_CHK_ERR(-1);
816 #endif
817 }
818 #endif
819 
820 //==============================================================================
822 GMRES(const Epetra_Operator& Operator,
823  const Epetra_Vector& InvPointDiagonal,
824  const int MaximumIterations,
825  double& lambda_real_min, double& lambda_real_max,
826  double& lambda_imag_min, double& lambda_imag_max)
827 {
828 #ifdef HAVE_IFPACK_AZTECOO
829  Epetra_Vector x(Operator_->OperatorDomainMap());
830  Epetra_Vector y(Operator_->OperatorRangeMap());
831  x.Random();
832  y.PutScalar(0.0);
833  Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
834  AztecOO solver(LP);
835  solver.SetAztecOption(AZ_solver, AZ_gmres_condnum);
836  solver.SetAztecOption(AZ_output, AZ_none);
837  Ifpack_DiagPreconditioner diag(Operator.OperatorDomainMap(),
838  Operator.OperatorRangeMap(),
839  InvPointDiagonal);
840  solver.SetPrecOperator(&diag);
841  solver.Iterate(MaximumIterations, 1e-10);
842  const double* status = solver.GetAztecStatus();
843  lambda_real_min = status[AZ_lambda_real_min];
844  lambda_real_max = status[AZ_lambda_real_max];
845  lambda_imag_min = status[AZ_lambda_imag_min];
846  lambda_imag_max = status[AZ_lambda_imag_max];
847  return(0);
848 #else
849  cout << "You need to configure IFPACK with support for AztecOO" << endl;
850  cout << "to use the GMRES estimator. This may require --enable-aztecoo" << endl;
851  cout << "in your configure script." << endl;
852  IFPACK_CHK_ERR(-1);
853 #endif
854 }
int GMRES(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &lambda_real_min, double &lambda_real_max, double &lambda_imag_min, double &lambda_imag_max)
Uses AztecOO's GMRES to estimate the height and width of the spectrum.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual ostream & Print(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.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
string Ifpack_toString(const int &x)
Convertes an integer to string.
static int CG(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &lambda_min, double &lambda_max)
Uses AztecOO's CG to estimate lambda_min and lambda_max.
virtual int Compute()
Computes the preconditioners.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
Ifpack_Polynomial(const Epetra_Operator *Matrix)
Ifpack_Polynomial constructor with given Epetra_Operator/Epetra_RowMatrix.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
static int PowerMethod(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &LambdaMax)
Simple power method to compute lambda_max.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.