IFPACK  Development
 All Classes Namespaces 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.
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), (unused; commented out to avoid build warnings) */
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?
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), (unused; commented out to avoid build warnings) */
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  std::string dummy("invert");
212  std::string ApplyMode=Blist.get("apply mode",dummy);
213  if(ApplyMode==std::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 //==============================================================================
234 {
235  return(Operator_->Comm());
236 }
237 
238 //==============================================================================
240 {
241  return(Operator_->OperatorDomainMap());
242 }
243 
244 //==============================================================================
246 {
247  return(Operator_->OperatorRangeMap());
248 }
249 
250 //==============================================================================
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  const std::complex<double> zero(0.0,0.0);
379 
380  // Compute points in complex plane
381  double lenx = LambdaRealMax_-LambdaRealMin_;
382  int nx = ceil(lenx*((double) LSPointsReal_));
383  if (nx<2) { nx = 2; }
384  double hx = lenx/((double) nx);
385  std::vector<double> xs;
386  if(abs(lenx)>1.0e-8) {
387  for( int pt=0; pt<=nx; pt++ ) {
388  xs.push_back(hx*pt+LambdaRealMin_);
389  }
390  }
391  else {
392  xs.push_back(LambdaRealMax_);
393  nx=1;
394  }
395  double leny = LambdaImagMax_-LambdaImagMin_;
396  int ny = ceil(leny*((double) LSPointsImag_));
397  if (ny<2) { ny = 2; }
398  double hy = leny/((double) ny);
399  std::vector<double> ys;
400  if(abs(leny)>1.0e-8) {
401  for( int pt=0; pt<=ny; pt++ ) {
402  ys.push_back(hy*pt+LambdaImagMin_);
403  }
404  }
405  else {
406  ys.push_back(LambdaImagMax_);
407  ny=1;
408  }
409  std::vector< std::complex<double> > cpts;
410  for( int jj=0; jj<ny; jj++ ) {
411  for( int ii=0; ii<nx; ii++ ) {
412  std::complex<double> cpt(xs[ii],ys[jj]);
413  cpts.push_back(cpt);
414  }
415  }
416  cpts.push_back(zero);
417 
418 #ifdef HAVE_TEUCHOS_COMPLEX
419  const std::complex<double> one(1.0,0.0);
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( std::vector<double>::size_type 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 std::ostream& Ifpack_Polynomial::Print(std::ostream & os) const
549 {
550  using std::endl;
551 
552  double MyMinVal, MyMaxVal;
553  double MinVal, MaxVal;
554 
555  if (IsComputed_) {
556  InvDiagonal_->MinValue(&MyMinVal);
557  InvDiagonal_->MaxValue(&MyMaxVal);
558  Comm().MinAll(&MyMinVal,&MinVal,1);
559  Comm().MinAll(&MyMaxVal,&MaxVal,1);
560  }
561 
562  if (!Comm().MyPID()) {
563  os << endl;
564  os << "================================================================================" << endl;
565  os << "Ifpack_Polynomial" << endl;
566  os << "Degree of polynomial = " << PolyDegree_ << endl;
567  os << "Condition number estimate = " << Condest() << endl;
568  os << "Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
569  if (IsComputed_) {
570  os << "Minimum value on stored inverse diagonal = " << MinVal << endl;
571  os << "Maximum value on stored inverse diagonal = " << MaxVal << endl;
572  }
573  if (ZeroStartingSolution_)
574  os << "Using zero starting solution" << endl;
575  else
576  os << "Using input starting solution" << endl;
577  os << endl;
578  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
579  os << "----- ------- -------------- ------------ --------" << endl;
580  os << "Initialize() " << std::setw(5) << NumInitialize_
581  << " " << std::setw(15) << InitializeTime_
582  << " 0.0 0.0" << endl;
583  os << "Compute() " << std::setw(5) << NumCompute_
584  << " " << std::setw(15) << ComputeTime_
585  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
586  if (ComputeTime_ != 0.0)
587  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
588  else
589  os << " " << std::setw(15) << 0.0 << endl;
590  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
591  << " " << std::setw(15) << ApplyInverseTime_
592  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
593  if (ApplyInverseTime_ != 0.0)
594  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
595  else
596  os << " " << std::setw(15) << 0.0 << endl;
597  os << "================================================================================" << endl;
598  os << endl;
599  }
600 
601  return(os);
602 }
603 
604 //==============================================================================
605 double Ifpack_Polynomial::
606 Condest(const Ifpack_CondestType CT,
607  const int MaxIters, const double Tol,
608  Epetra_RowMatrix* Matrix_in)
609 {
610  if (!IsComputed()) // cannot compute right now
611  return(-1.0);
612 
613  // always computes it. Call Condest() with no parameters to get
614  // the previous estimate.
615  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
616 
617  return(Condest_);
618 }
619 
620 //==============================================================================
621 void Ifpack_Polynomial::SetLabel()
622 {
623  Label_ = "IFPACK (Least squares polynomial), degree=" + Ifpack_toString(PolyDegree_);
624 }
625 
626 //==============================================================================
629 {
630 
631  if (!IsComputed())
632  IFPACK_CHK_ERR(-3);
633 
634  if (PolyDegree_ == 0)
635  return 0;
636 
637  int nVec = X.NumVectors();
638  if (nVec != Y.NumVectors())
639  IFPACK_CHK_ERR(-2);
640 
641  Time_->ResetStartTime();
642 
643  Epetra_MultiVector Xcopy(X);
644  if(ZeroStartingSolution_==true) {
645  Y.PutScalar(0.0);
646  }
647 
648  // mfh 20 Mar 2014: IBD never gets used, so I'm commenting out the
649  // following lines of code in order to forestall build warnings.
650 // #ifdef HAVE_IFPACK_EPETRAEXT
651 // EpetraExt_PointToBlockDiagPermute* IBD=0;
652 // if (UseBlockMode_) IBD=&*InvBlockDiagonal_;
653 // #endif
654 
655  Y.Update(-coeff_[1], Xcopy, 1.0);
656  for (int ii = 2; ii < static_cast<int> (coeff_.size ()); ++ii) {
657  const Epetra_MultiVector V(Xcopy);
658  Operator_->Apply(V,Xcopy);
659  Xcopy.Multiply(1.0, *InvDiagonal_, Xcopy, 0.0);
660  // Update Y
661  Y.Update(-coeff_[ii], Xcopy, 1.0);
662  }
663 
664  // Flops are updated in each of the following.
665  ++NumApplyInverse_;
666  ApplyInverseTime_ += Time_->ElapsedTime();
667  return(0);
668 }
669 
670 //==============================================================================
672 PowerMethod(const Epetra_Operator& Operator,
673  const Epetra_Vector& InvPointDiagonal,
674  const int MaximumIterations,
675  double& lambda_max)
676 {
677  // this is a simple power method
678  lambda_max = 0.0;
679  double RQ_top, RQ_bottom, norm;
680  Epetra_Vector x(Operator.OperatorDomainMap());
681  Epetra_Vector y(Operator.OperatorRangeMap());
682  x.Random();
683  x.Norm2(&norm);
684  if (norm == 0.0) IFPACK_CHK_ERR(-1);
685 
686  x.Scale(1.0 / norm);
687 
688  for (int iter = 0; iter < MaximumIterations; ++iter)
689  {
690  Operator.Apply(x, y);
691  IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
692  IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
693  IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
694  lambda_max = RQ_top / RQ_bottom;
695  IFPACK_CHK_ERR(y.Norm2(&norm));
696  if (norm == 0.0) IFPACK_CHK_ERR(-1);
697  IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
698  }
699 
700  return(0);
701 }
702 
703 //==============================================================================
705 CG(const Epetra_Operator& Operator,
706  const Epetra_Vector& InvPointDiagonal,
707  const int MaximumIterations,
708  double& lambda_min, double& lambda_max)
709 {
710 #ifdef HAVE_IFPACK_AZTECOO
711  Epetra_Vector x(Operator.OperatorDomainMap());
712  Epetra_Vector y(Operator.OperatorRangeMap());
713  x.Random();
714  y.PutScalar(0.0);
715 
716  Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
717  AztecOO solver(LP);
718  solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
719  solver.SetAztecOption(AZ_output, AZ_none);
720 
722  Operator.OperatorRangeMap(),
723  InvPointDiagonal);
724  solver.SetPrecOperator(&diag);
725  solver.Iterate(MaximumIterations, 1e-10);
726 
727  const double* status = solver.GetAztecStatus();
728 
729  lambda_min = status[AZ_lambda_min];
730  lambda_max = status[AZ_lambda_max];
731 
732  return(0);
733 #else
734  using std::cout;
735  using std::endl;
736 
737  cout << "You need to configure IFPACK with support for AztecOO" << endl;
738  cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
739  cout << "in your configure script." << endl;
740  IFPACK_CHK_ERR(-1);
741 #endif
742 }
743 
744 //==============================================================================
745 #ifdef HAVE_IFPACK_EPETRAEXT
747 PowerMethod(const int MaximumIterations, double& lambda_max)
748 {
749 
750  if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
751  // this is a simple power method
752  lambda_max = 0.0;
753  double RQ_top, RQ_bottom, norm;
754  Epetra_Vector x(Operator_->OperatorDomainMap());
755  Epetra_Vector y(Operator_->OperatorRangeMap());
756  Epetra_Vector z(Operator_->OperatorRangeMap());
757  x.Random();
758  x.Norm2(&norm);
759  if (norm == 0.0) IFPACK_CHK_ERR(-1);
760 
761  x.Scale(1.0 / norm);
762 
763  for (int iter = 0; iter < MaximumIterations; ++iter)
764  {
765  Operator_->Apply(x, z);
766  InvBlockDiagonal_->ApplyInverse(z,y);
767  if(SolveNormalEquations_){
768  InvBlockDiagonal_->ApplyInverse(y,z);
769  Apply_Transpose(Operator_,z, y);
770  }
771 
772  IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
773  IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
774  lambda_max = RQ_top / RQ_bottom;
775  IFPACK_CHK_ERR(y.Norm2(&norm));
776  if (norm == 0.0) IFPACK_CHK_ERR(-1);
777  IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
778  }
779 
780  return(0);
781 }
782 #endif
783 
784 //==============================================================================
785 #ifdef HAVE_IFPACK_EPETRAEXT
787 CG(const int /* MaximumIterations */,
788  double& /* lambda_min */, double& /* lambda_max */)
789 {
790  IFPACK_CHK_ERR(-1);// NTS: This always seems to yield errors in AztecOO, ergo,
791  // I turned it off.
792 
793  // mfh 06 Aug 2017: None of the code below here is reachable.
794  // This causes build warnings with some compilers. Thus, I'm
795  // commenting out this code.
796 #if 0
797 
798  if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
799 
800 #ifdef HAVE_IFPACK_AZTECOO
801  Epetra_Vector x(Operator_->OperatorDomainMap());
802  Epetra_Vector y(Operator_->OperatorRangeMap());
803  x.Random();
804  y.PutScalar(0.0);
805  Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
806 
807  AztecOO solver(LP);
808  solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
809  solver.SetAztecOption(AZ_output, AZ_none);
810 
811  solver.SetPrecOperator(&*InvBlockDiagonal_);
812  solver.Iterate(MaximumIterations, 1e-10);
813 
814  const double* status = solver.GetAztecStatus();
815 
816  lambda_min = status[AZ_lambda_min];
817  lambda_max = status[AZ_lambda_max];
818 
819  return(0);
820 #else
821  using std::cout;
822  using std::endl;
823 
824  cout << "You need to configure IFPACK with support for AztecOO" << endl;
825  cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
826  cout << "in your configure script." << endl;
827  IFPACK_CHK_ERR(-1);
828 #endif
829 
830 #endif // 0
831 }
832 #endif // HAVE_IFPACK_EPETRAEXT
833 
834 //==============================================================================
836 GMRES(const Epetra_Operator& Operator,
837  const Epetra_Vector& InvPointDiagonal,
838  const int MaximumIterations,
839  double& lambda_real_min, double& lambda_real_max,
840  double& lambda_imag_min, double& lambda_imag_max)
841 {
842 #ifdef HAVE_IFPACK_AZTECOO
843  Epetra_Vector x(Operator_->OperatorDomainMap());
844  Epetra_Vector y(Operator_->OperatorRangeMap());
845  x.Random();
846  y.PutScalar(0.0);
847  Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
848  AztecOO solver(LP);
849  solver.SetAztecOption(AZ_solver, AZ_gmres_condnum);
850  solver.SetAztecOption(AZ_output, AZ_none);
852  Operator.OperatorRangeMap(),
853  InvPointDiagonal);
854  solver.SetPrecOperator(&diag);
855  solver.Iterate(MaximumIterations, 1e-10);
856  const double* status = solver.GetAztecStatus();
857  lambda_real_min = status[AZ_lambda_real_min];
858  lambda_real_max = status[AZ_lambda_real_max];
859  lambda_imag_min = status[AZ_lambda_imag_min];
860  lambda_imag_max = status[AZ_lambda_imag_max];
861  return(0);
862 #else
863  using std::cout;
864  using std::endl;
865 
866  cout << "You need to configure IFPACK with support for AztecOO" << endl;
867  cout << "to use the GMRES estimator. This may require --enable-aztecoo" << endl;
868  cout << "in your configure script." << endl;
869  IFPACK_CHK_ERR(-1);
870 #endif
871 }
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&#39;s GMRES to estimate the height and width of the spectrum.
virtual int SetUseTranspose(bool UseTranspose)=0
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
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.
virtual const Epetra_Map & OperatorDomainMap() const =0
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
static int CG(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &lambda_min, double &lambda_max)
Uses AztecOO&#39;s CG to estimate lambda_min and lambda_max.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
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 std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
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.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
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.