IFPACK  Development
Ifpack_Chebyshev.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_Chebyshev.h"
53 #include "Ifpack_Utils.h"
54 #include "Ifpack_Condest.h"
55 #ifdef HAVE_IFPACK_AZTECOO
56 #include "Ifpack_DiagPreconditioner.h"
57 #include "AztecOO.h"
58 #endif
59 
60 #ifdef HAVE_IFPACK_EPETRAEXT
61 #include "Epetra_CrsMatrix.h"
62 #include "EpetraExt_PointToBlockDiagPermute.h"
63 #endif
64 
65 
66 #define ABS(x) ((x)>0?(x):-(x))
67 
68 // Helper function for normal equations
69 inline void Apply_Transpose(Teuchos::RCP<const Epetra_Operator> Operator_,const Epetra_MultiVector &X,Epetra_MultiVector &Y){
70  Epetra_Operator * Operator=const_cast<Epetra_Operator*>(&*Operator_);
71  Operator->SetUseTranspose(true);
72  Operator->Apply(X,Y);
73  Operator->SetUseTranspose(false);
74 }
75 
76 
77 
78 
79 //==============================================================================
80 // NOTE: any change to the default values should be committed to the other
81 // constructor as well.
83 Ifpack_Chebyshev(const Epetra_Operator* Operator) :
84  IsInitialized_(false),
85  IsComputed_(false),
86  NumInitialize_(0),
87  NumCompute_(0),
88  NumApplyInverse_(0),
89  InitializeTime_(0.0),
90  ComputeTime_(0.0),
91  ApplyInverseTime_(0.0),
92  ComputeFlops_(0.0),
93  ApplyInverseFlops_(0.0),
94  PolyDegree_(1),
95  UseTranspose_(false),
96  Condest_(-1.0),
97  /* ComputeCondest_(false), (Unused; commented out to avoid build warnings) */
98  EigRatio_(30.0),
99  Label_(),
100  LambdaMin_(0.0),
101  LambdaMax_(-1.0),
102  MinDiagonalValue_(0.0),
103  NumMyRows_(0),
104  NumMyNonzeros_(0),
105  NumGlobalRows_(0),
106  NumGlobalNonzeros_(0),
107  Operator_(Teuchos::rcp(Operator,false)),
108  UseBlockMode_(false),
109  SolveNormalEquations_(false),
110  IsRowMatrix_(false),
111  ZeroStartingSolution_(true)
112 {
113 }
114 
115 //==============================================================================
116 // NOTE: This constructor has been introduced because SWIG does not appear
117 // to appreciate dynamic_cast. An instruction of type
118 // Matrix_ = dynamic_cast<const Epetra_RowMatrix*> in the
119 // other construction does not work in PyTrilinos -- of course
120 // it does in any C++ code (for an Epetra_Operator that is also
121 // an Epetra_RowMatrix).
122 //
123 // FIXME: move declarations into a separate method?
125 Ifpack_Chebyshev(const Epetra_RowMatrix* Operator) :
126  IsInitialized_(false),
127  IsComputed_(false),
128  NumInitialize_(0),
129  NumCompute_(0),
130  NumApplyInverse_(0),
131  InitializeTime_(0.0),
132  ComputeTime_(0.0),
133  ApplyInverseTime_(0.0),
134  ComputeFlops_(0.0),
135  ApplyInverseFlops_(0.0),
136  PolyDegree_(1),
137  UseTranspose_(false),
138  Condest_(-1.0),
139  /* ComputeCondest_(false), (Unused; commented out to avoid build warnings) */
140  EigRatio_(30.0),
141  EigMaxIters_(10),
142  Label_(),
143  LambdaMin_(0.0),
144  LambdaMax_(-1.0),
145  MinDiagonalValue_(0.0),
146  NumMyRows_(0),
147  NumMyNonzeros_(0),
148  NumGlobalRows_(0),
149  NumGlobalNonzeros_(0),
150  Operator_(Teuchos::rcp(Operator,false)),
151  Matrix_(Teuchos::rcp(Operator,false)),
152  UseBlockMode_(false),
153  SolveNormalEquations_(false),
154  IsRowMatrix_(true),
155  ZeroStartingSolution_(true)
156 {
157 }
158 
159 //==============================================================================
160 int Ifpack_Chebyshev::SetParameters(Teuchos::ParameterList& List)
161 {
162 
163  EigRatio_ = List.get("chebyshev: ratio eigenvalue", EigRatio_);
164  LambdaMin_ = List.get("chebyshev: min eigenvalue", LambdaMin_);
165  LambdaMax_ = List.get("chebyshev: max eigenvalue", LambdaMax_);
166  PolyDegree_ = List.get("chebyshev: degree",PolyDegree_);
167  MinDiagonalValue_ = List.get("chebyshev: min diagonal value",
168  MinDiagonalValue_);
169  ZeroStartingSolution_ = List.get("chebyshev: zero starting solution",
170  ZeroStartingSolution_);
171 
172  Epetra_Vector* ID = List.get("chebyshev: operator inv diagonal",
173  (Epetra_Vector*)0);
174  EigMaxIters_ = List.get("chebyshev: eigenvalue max iterations",EigMaxIters_);
175 
176 #ifdef HAVE_IFPACK_EPETRAEXT
177  // This is *always* false if EpetraExt isn't enabled
178  UseBlockMode_ = List.get("chebyshev: use block mode",UseBlockMode_);
179  if(!List.isParameter("chebyshev: block list")) UseBlockMode_=false;
180  else{
181  BlockList_ = List.get("chebyshev: block list",BlockList_);
182 
183  // Since we know we're doing a matrix inverse, clobber the block list
184  // w/"invert" if it's set to multiply
185  Teuchos::ParameterList Blist;
186  Blist=BlockList_.get("blockdiagmatrix: list",Blist);
187  std::string dummy("invert");
188  std::string ApplyMode=Blist.get("apply mode",dummy);
189  if(ApplyMode==std::string("multiply")){
190  Blist.set("apply mode","invert");
191  BlockList_.set("blockdiagmatrix: list",Blist);
192  }
193  }
194 #endif
195 
196  SolveNormalEquations_ = List.get("chebyshev: solve normal equations",SolveNormalEquations_);
197 
198  if (ID != 0)
199  {
200  InvDiagonal_ = Teuchos::rcp( new Epetra_Vector(*ID) );
201  }
202 
203  SetLabel();
204 
205  return(0);
206 }
207 
208 //==============================================================================
209 const Epetra_Comm& Ifpack_Chebyshev::Comm() const
210 {
211  return(Operator_->Comm());
212 }
213 
214 //==============================================================================
215 const Epetra_Map& Ifpack_Chebyshev::OperatorDomainMap() const
216 {
217  return(Operator_->OperatorDomainMap());
218 }
219 
220 //==============================================================================
221 const Epetra_Map& Ifpack_Chebyshev::OperatorRangeMap() const
222 {
223  return(Operator_->OperatorRangeMap());
224 }
225 
226 //==============================================================================
228 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
229 {
230  if (IsComputed() == false)
231  IFPACK_CHK_ERR(-3);
232 
233  if (X.NumVectors() != Y.NumVectors())
234  IFPACK_CHK_ERR(-2);
235 
236  if (IsRowMatrix_)
237  {
238  IFPACK_CHK_ERR(Matrix_->Multiply(UseTranspose(),X,Y));
239  }
240  else
241  {
242  IFPACK_CHK_ERR(Operator_->Apply(X,Y));
243  }
244 
245  return(0);
246 }
247 
248 //==============================================================================
250 {
251  IsInitialized_ = false;
252 
253  if (Operator_ == Teuchos::null)
254  IFPACK_CHK_ERR(-2);
255 
256  if (Time_ == Teuchos::null)
257  Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
258 
259  if (IsRowMatrix_)
260  {
261  if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
262  IFPACK_CHK_ERR(-2); // only square matrices
263 
264  NumMyRows_ = Matrix_->NumMyRows();
265  NumMyNonzeros_ = Matrix_->NumMyNonzeros();
266  NumGlobalRows_ = Matrix_->NumGlobalRows64();
267  NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
268  }
269  else
270  {
271  if (Operator_->OperatorDomainMap().NumGlobalElements64() !=
272  Operator_->OperatorRangeMap().NumGlobalElements64())
273  IFPACK_CHK_ERR(-2); // only square operators
274  }
275 
276  ++NumInitialize_;
277  InitializeTime_ += Time_->ElapsedTime();
278  IsInitialized_ = true;
279  return(0);
280 }
281 
282 //==============================================================================
284 {
285  if (!IsInitialized())
286  IFPACK_CHK_ERR(Initialize());
287 
288  Time_->ResetStartTime();
289 
290  // reset values
291  IsComputed_ = false;
292  Condest_ = -1.0;
293 
294  if (PolyDegree_ <= 0)
295  IFPACK_CHK_ERR(-2); // at least one application
296 
297 #ifdef HAVE_IFPACK_EPETRAEXT
298  // Check to see if we can run in block mode
299  if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && UseBlockMode_){
300  const Epetra_CrsMatrix *CrsMatrix = dynamic_cast<const Epetra_CrsMatrix*>(&*Matrix_);
301 
302  // If we don't have CrsMatrix, we can't use the block preconditioner
303  if (!CrsMatrix) {
304  UseBlockMode_ = false;
305 
306  } else {
307  int ierr;
308  InvBlockDiagonal_ = Teuchos::rcp(new EpetraExt_PointToBlockDiagPermute(*CrsMatrix));
309  if (InvBlockDiagonal_ == Teuchos::null)
310  IFPACK_CHK_ERR(-6);
311 
312  ierr = InvBlockDiagonal_->SetParameters(BlockList_);
313  if (ierr)
314  IFPACK_CHK_ERR(-7);
315 
316  ierr = InvBlockDiagonal_->Compute();
317  if (ierr)
318  IFPACK_CHK_ERR(-8);
319  }
320 
321  // Automatically Compute Eigenvalues
322  double lambda_max = 0;
323  if (LambdaMax_ == -1) {
324  PowerMethod(EigMaxIters_, lambda_max);
325  LambdaMax_ = lambda_max;
326  }
327 
328  // Test for Exact Preconditioned case
329  if (ABS(LambdaMax_-1) < 1e-6)
330  LambdaMax_ = LambdaMin_ = 1.0;
331  else
332  LambdaMin_ = LambdaMax_/EigRatio_;
333  }
334 #endif
335 
336  if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_) {
337  InvDiagonal_ = Teuchos::rcp(new Epetra_Vector(Matrix().Map()));
338 
339  if (InvDiagonal_ == Teuchos::null)
340  IFPACK_CHK_ERR(-5);
341 
342  IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*InvDiagonal_));
343 
344  // Inverse diagonal elements
345  // Replace zeros with 1.0
346  for (int i = 0 ; i < NumMyRows_ ; ++i) {
347  double diag = (*InvDiagonal_)[i];
348  if (IFPACK_ABS(diag) < MinDiagonalValue_)
349  (*InvDiagonal_)[i] = MinDiagonalValue_;
350  else
351  (*InvDiagonal_)[i] = 1.0 / diag;
352  }
353  // Automatically compute maximum eigenvalue estimate of D^{-1}A if user hasn't provided one
354  double lambda_max=0;
355  if (LambdaMax_ == -1) {
356  PowerMethod(Matrix(), *InvDiagonal_, EigMaxIters_, lambda_max);
357  LambdaMax_=lambda_max;
358  // Test for Exact Preconditioned case
359  if (ABS(LambdaMax_-1) < 1e-6) LambdaMax_=LambdaMin_=1.0;
360  else LambdaMin_=LambdaMax_/EigRatio_;
361  }
362  // otherwise the inverse of the diagonal has been given by the user
363  }
364 #ifdef IFPACK_FLOPCOUNTERS
365  ComputeFlops_ += NumMyRows_;
366 #endif
367 
368  ++NumCompute_;
369  ComputeTime_ += Time_->ElapsedTime();
370  IsComputed_ = true;
371 
372  SetLabel();
373 
374  return(0);
375 }
376 
377 //==============================================================================
378 std::ostream& Ifpack_Chebyshev::Print(std::ostream & os) const
379 {
380  using std::endl;
381 
382  double MyMinVal, MyMaxVal;
383  double MinVal, MaxVal;
384 
385  if (IsComputed_) {
386  InvDiagonal_->MinValue(&MyMinVal);
387  InvDiagonal_->MaxValue(&MyMaxVal);
388  Comm().MinAll(&MyMinVal,&MinVal,1);
389  Comm().MinAll(&MyMaxVal,&MaxVal,1);
390  }
391 
392  if (!Comm().MyPID()) {
393  os << endl;
394  os << "================================================================================" << endl;
395  os << "Ifpack_Chebyshev" << endl;
396  os << "Degree of polynomial = " << PolyDegree_ << endl;
397  os << "Condition number estimate = " << Condest() << endl;
398  os << "Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
399  if (IsComputed_) {
400  os << "Minimum value on stored inverse diagonal = " << MinVal << endl;
401  os << "Maximum value on stored inverse diagonal = " << MaxVal << endl;
402  }
403  if (ZeroStartingSolution_)
404  os << "Using zero starting solution" << endl;
405  else
406  os << "Using input starting solution" << endl;
407  os << endl;
408  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
409  os << "----- ------- -------------- ------------ --------" << endl;
410  os << "Initialize() " << std::setw(5) << NumInitialize_
411  << " " << std::setw(15) << InitializeTime_
412  << " 0.0 0.0" << endl;
413  os << "Compute() " << std::setw(5) << NumCompute_
414  << " " << std::setw(15) << ComputeTime_
415  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
416  if (ComputeTime_ != 0.0)
417  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
418  else
419  os << " " << std::setw(15) << 0.0 << endl;
420  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
421  << " " << std::setw(15) << ApplyInverseTime_
422  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
423  if (ApplyInverseTime_ != 0.0)
424  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
425  else
426  os << " " << std::setw(15) << 0.0 << endl;
427  os << "================================================================================" << endl;
428  os << endl;
429  }
430 
431  return(os);
432 }
433 
434 //==============================================================================
435 double Ifpack_Chebyshev::
436 Condest(const Ifpack_CondestType CT,
437  const int MaxIters, const double Tol,
438  Epetra_RowMatrix* Matrix_in)
439 {
440  if (!IsComputed()) // cannot compute right now
441  return(-1.0);
442 
443  // always computes it. Call Condest() with no parameters to get
444  // the previous estimate.
445  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
446 
447  return(Condest_);
448 }
449 
450 //==============================================================================
451 void Ifpack_Chebyshev::SetLabel()
452 {
453  std::ostringstream oss;
454  oss << "\"Ifpack Chebyshev polynomial\": {"
455  << "Initialized: " << (IsInitialized() ? "true" : "false")
456  << ", Computed: " << (IsComputed() ? "true" : "false")
457  << ", degree: " << PolyDegree_
458  << ", lambdaMax: " << LambdaMax_
459  << ", alpha: " << EigRatio_
460  << ", lambdaMin: " << LambdaMin_
461  << "}";
462  Label_ = oss.str();
463 }
464 
465 //==============================================================================
467 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
468 {
469  if (!IsComputed())
470  IFPACK_CHK_ERR(-3);
471 
472  if (PolyDegree_ == 0)
473  return 0;
474 
475  int nVec = X.NumVectors();
476  int len = X.MyLength();
477  if (nVec != Y.NumVectors())
478  IFPACK_CHK_ERR(-2);
479 
480  Time_->ResetStartTime();
481 
482  // AztecOO gives X and Y pointing to the same memory location,
483  // need to create an auxiliary vector, Xcopy
484  Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
485  if (X.Pointers()[0] == Y.Pointers()[0])
486  Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
487  else
488  Xcopy = Teuchos::rcp( &X, false );
489 
490  double **xPtr = 0, **yPtr = 0;
491  Xcopy->ExtractView(&xPtr);
492  Y.ExtractView(&yPtr);
493 
494 #ifdef HAVE_IFPACK_EPETRAEXT
495  EpetraExt_PointToBlockDiagPermute* IBD=0;
496  if (UseBlockMode_)
497  IBD = &*InvBlockDiagonal_;
498 #endif
499 
500  //--- Do a quick solve when the matrix is identity
501  double *invDiag = 0;
502  if (!UseBlockMode_)
503  invDiag = InvDiagonal_->Values();
504 
505  if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_)) {
506 #ifdef HAVE_IFPACK_EPETRAEXT
507  if (UseBlockMode_)
508  IBD->ApplyInverse(*Xcopy, Y);
509  else
510 #endif
511  if (nVec == 1) {
512  double *yPointer = yPtr[0], *xPointer = xPtr[0];
513  for (int i = 0; i < len; ++i)
514  yPointer[i] = xPointer[i] * invDiag[i];
515 
516  } else {
517  for (int i = 0; i < len; ++i) {
518  double coeff = invDiag[i];
519  for (int k = 0; k < nVec; ++k)
520  yPtr[k][i] = xPtr[k][i] * coeff;
521  }
522  } // if (nVec == 1)
523 
524  return 0;
525  } // if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_))
526 
527  //--- Initialize coefficients
528  // Note that delta stores the inverse of ML_Cheby::delta
529  double alpha = LambdaMax_ / EigRatio_;
530  double beta = 1.1 * LambdaMax_;
531  double delta = 2.0 / (beta - alpha);
532  double theta = 0.5 * (beta + alpha);
533  double s1 = theta * delta;
534 
535  // Temporary vectors
536  // In ML_Cheby, V corresponds to pAux and W to dk
537  // NOTE: we would like to move the construction to the Compute()
538  // call, but that is not possible as we don't know how many
539  // vectors are in the multivector
540  bool zeroOut = false;
541  Epetra_MultiVector V(X.Map(), X.NumVectors(), zeroOut);
542  Epetra_MultiVector W(X.Map(), X.NumVectors(), zeroOut);
543 #ifdef HAVE_IFPACK_EPETRAEXT
544  Epetra_MultiVector Temp(X.Map(), X.NumVectors(), zeroOut);
545 #endif
546 
547  double *vPointer = V.Values(), *wPointer = W.Values();
548 
549  double oneOverTheta = 1.0/theta;
550 
551  //--- If solving normal equations, multiply RHS by A^T
552  if (SolveNormalEquations_) {
553  Apply_Transpose(Operator_, Y, V);
554  Y = V;
555  }
556 
557  // Do the smoothing when block scaling is turned OFF
558  // --- Treat the initial guess
559  if (ZeroStartingSolution_ == false) {
560  Operator_->Apply(Y, V);
561 
562  // Compute W = invDiag * ( X - V )/ Theta
563 #ifdef HAVE_IFPACK_EPETRAEXT
564  if (UseBlockMode_) {
565  Temp.Update(oneOverTheta, X, -oneOverTheta, V, 0.0);
566  IBD->ApplyInverse(Temp, W);
567 
568  // Perform additional matvecs for normal equations
569  // CMS: Testing this only in block mode FOR NOW
570  if (SolveNormalEquations_){
571  IBD->ApplyInverse(W, Temp);
572  Apply_Transpose(Operator_, Temp, W);
573  }
574  }
575  else
576 #endif
577  if (nVec == 1) {
578  double *xPointer = xPtr[0];
579  for (int i = 0; i < len; ++i)
580  wPointer[i] = invDiag[i] * (xPointer[i] - vPointer[i]) * oneOverTheta;
581 
582  } else {
583  for (int i = 0; i < len; ++i) {
584  double coeff = invDiag[i]*oneOverTheta;
585  double *wi = wPointer + i, *vi = vPointer + i;
586  for (int k = 0; k < nVec; ++k) {
587  *wi = (xPtr[k][i] - (*vi)) * coeff;
588  wi = wi + len; vi = vi + len;
589  }
590  }
591  } // if (nVec == 1)
592  // Update the vector Y
593  Y.Update(1.0, W, 1.0);
594 
595  } else { // if (ZeroStartingSolution_ == false)
596  // Compute W = invDiag * X / Theta
597 #ifdef HAVE_IFPACK_EPETRAEXT
598  if (UseBlockMode_) {
599  IBD->ApplyInverse(X, W);
600 
601  // Perform additional matvecs for normal equations
602  // CMS: Testing this only in block mode FOR NOW
603  if (SolveNormalEquations_) {
604  IBD->ApplyInverse(W, Temp);
605  Apply_Transpose(Operator_, Temp, W);
606  }
607 
608  W.Scale(oneOverTheta);
609  Y.Update(1.0, W, 0.0);
610  }
611  else
612 #endif
613  if (nVec == 1) {
614  double *xPointer = xPtr[0];
615  for (int i = 0; i < len; ++i)
616  wPointer[i] = invDiag[i] * xPointer[i] * oneOverTheta;
617 
618  memcpy(yPtr[0], wPointer, len*sizeof(double));
619 
620  } else {
621  for (int i = 0; i < len; ++i) {
622  double coeff = invDiag[i]*oneOverTheta;
623  double *wi = wPointer + i;
624  for (int k = 0; k < nVec; ++k) {
625  *wi = xPtr[k][i] * coeff;
626  wi = wi + len;
627  }
628  }
629 
630  for (int k = 0; k < nVec; ++k)
631  memcpy(yPtr[k], wPointer + k*len, len*sizeof(double));
632  } // if (nVec == 1)
633  } // if (ZeroStartingSolution_ == false)
634 
635  //--- Apply the polynomial
636  double rhok = 1.0/s1, rhokp1;
637  double dtemp1, dtemp2;
638  int degreeMinusOne = PolyDegree_ - 1;
639  if (nVec == 1) {
640  double *xPointer = xPtr[0];
641  for (int k = 0; k < degreeMinusOne; ++k) {
642  Operator_->Apply(Y, V);
643  rhokp1 = 1.0 / (2.0*s1 - rhok);
644  dtemp1 = rhokp1 * rhok;
645  dtemp2 = 2.0 * rhokp1 * delta;
646  rhok = rhokp1;
647 
648  // Compute W = dtemp1 * W
649  W.Scale(dtemp1);
650 
651  // Compute W = W + dtemp2 * invDiag * ( X - V )
652 #ifdef HAVE_IFPACK_EPETRAEXT
653  if (UseBlockMode_) {
654  //NTS: We can clobber V since it will be reset in the Apply
655  V.Update(dtemp2, X, -dtemp2);
656  IBD->ApplyInverse(V, Temp);
657 
658  // Perform additional matvecs for normal equations
659  // CMS: Testing this only in block mode FOR NOW
660  if (SolveNormalEquations_) {
661  IBD->ApplyInverse(V, Temp);
662  Apply_Transpose(Operator_, Temp, V);
663  }
664 
665  W.Update(1.0, Temp, 1.0);
666  }
667  else
668 #endif
669  for (int i = 0; i < len; ++i)
670  wPointer[i] += dtemp2* invDiag[i] * (xPointer[i] - vPointer[i]);
671 
672  // Update the vector Y
673  Y.Update(1.0, W, 1.0);
674  } // for (k = 0; k < degreeMinusOne; ++k)
675 
676  } else { // if (nVec == 1) {
677  for (int k = 0; k < degreeMinusOne; ++k) {
678  Operator_->Apply(Y, V);
679  rhokp1 = 1.0 / (2.0*s1 - rhok);
680  dtemp1 = rhokp1 * rhok;
681  dtemp2 = 2.0 * rhokp1 * delta;
682  rhok = rhokp1;
683 
684  // Compute W = dtemp1 * W
685  W.Scale(dtemp1);
686 
687  // Compute W = W + dtemp2 * invDiag * ( X - V )
688 #ifdef HAVE_IFPACK_EPETRAEXT
689  if (UseBlockMode_) {
690  // We can clobber V since it will be reset in the Apply
691  V.Update(dtemp2, X, -dtemp2);
692  IBD->ApplyInverse(V, Temp);
693 
694  // Perform additional matvecs for normal equations
695  // CMS: Testing this only in block mode FOR NOW
696  if (SolveNormalEquations_) {
697  IBD->ApplyInverse(V,Temp);
698  Apply_Transpose(Operator_,Temp,V);
699  }
700 
701  W.Update(1.0, Temp, 1.0);
702  }
703  else
704 #endif
705  for (int i = 0; i < len; ++i) {
706  double coeff = invDiag[i]*dtemp2;
707  double *wi = wPointer + i, *vi = vPointer + i;
708  for (int j = 0; j < nVec; ++j) {
709  *wi += (xPtr[j][i] - (*vi)) * coeff;
710  wi = wi + len; vi = vi + len;
711  }
712  }
713 
714  // Update the vector Y
715  Y.Update(1.0, W, 1.0);
716  } // for (k = 0; k < degreeMinusOne; ++k)
717  } // if (nVec == 1)
718 
719 
720  // Flops are updated in each of the following.
721  ++NumApplyInverse_;
722  ApplyInverseTime_ += Time_->ElapsedTime();
723 
724  return(0);
725 }
726 
727 //==============================================================================
729 PowerMethod(const Epetra_Operator& Operator,
730  const Epetra_Vector& InvPointDiagonal,
731  const int MaximumIterations,
732  double& lambda_max,const unsigned int * RngSeed)
733 {
734  // this is a simple power method
735  lambda_max = 0.0;
736  double RQ_top, RQ_bottom, norm;
737  Epetra_Vector x(Operator.OperatorDomainMap());
738  Epetra_Vector y(Operator.OperatorRangeMap());
739  if(RngSeed) x.SetSeed(*RngSeed);
740  x.Random();
741  x.Norm2(&norm);
742  if (norm == 0.0) IFPACK_CHK_ERR(-1);
743 
744  x.Scale(1.0 / norm);
745 
746  for (int iter = 0; iter < MaximumIterations; ++iter)
747  {
748  Operator.Apply(x, y);
749  IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
750  IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
751  IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
752  lambda_max = RQ_top / RQ_bottom;
753  IFPACK_CHK_ERR(y.Norm2(&norm));
754  if (norm == 0.0) IFPACK_CHK_ERR(-1);
755  IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
756  }
757 
758  return(0);
759 }
760 
761 //==============================================================================
763 CG(const Epetra_Operator& Operator,
764  const Epetra_Vector& InvPointDiagonal,
765  const int MaximumIterations,
766  double& lambda_min, double& lambda_max,const unsigned int * RngSeed)
767 {
768 #ifdef HAVE_IFPACK_AZTECOO
769  Epetra_Vector x(Operator.OperatorDomainMap());
770  Epetra_Vector y(Operator.OperatorRangeMap());
771  if(RngSeed) x.SetSeed(*RngSeed);
772  x.Random();
773  y.PutScalar(0.0);
774 
775  Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
776  AztecOO solver(LP);
777  solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
778  solver.SetAztecOption(AZ_output, AZ_none);
779 
780  Ifpack_DiagPreconditioner diag(Operator.OperatorDomainMap(),
781  Operator.OperatorRangeMap(),
782  InvPointDiagonal);
783  solver.SetPrecOperator(&diag);
784  solver.Iterate(MaximumIterations, 1e-10);
785 
786  const double* status = solver.GetAztecStatus();
787 
788  lambda_min = status[AZ_lambda_min];
789  lambda_max = status[AZ_lambda_max];
790 
791  return(0);
792 #else
793  using std::cout;
794  using std::endl;
795 
796  cout << "You need to configure IFPACK with support for AztecOO" << endl;
797  cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
798  cout << "in your configure script." << endl;
799  IFPACK_CHK_ERR(-1);
800 #endif
801 }
802 
803 //==============================================================================
804 #ifdef HAVE_IFPACK_EPETRAEXT
806 PowerMethod(const int MaximumIterations, double& lambda_max,const unsigned int * RngSeed)
807 {
808 
809  if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
810  // this is a simple power method
811  lambda_max = 0.0;
812  double RQ_top, RQ_bottom, norm;
813  Epetra_Vector x(Operator_->OperatorDomainMap());
814  Epetra_Vector y(Operator_->OperatorRangeMap());
815  Epetra_Vector z(Operator_->OperatorRangeMap());
816  if(RngSeed) x.SetSeed(*RngSeed);
817  x.Random();
818  x.Norm2(&norm);
819  if (norm == 0.0) IFPACK_CHK_ERR(-1);
820 
821  x.Scale(1.0 / norm);
822 
823  for (int iter = 0; iter < MaximumIterations; ++iter)
824  {
825  Operator_->Apply(x, z);
826  InvBlockDiagonal_->ApplyInverse(z,y);
827  if(SolveNormalEquations_){
828  InvBlockDiagonal_->ApplyInverse(y,z);
829  Apply_Transpose(Operator_,z, y);
830  }
831 
832  IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
833  IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
834  lambda_max = RQ_top / RQ_bottom;
835  IFPACK_CHK_ERR(y.Norm2(&norm));
836  if (norm == 0.0) IFPACK_CHK_ERR(-1);
837  IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
838  }
839 
840  return(0);
841 }
842 #endif
843 
844 //==============================================================================
845 #ifdef HAVE_IFPACK_EPETRAEXT
847 CG(const int MaximumIterations,
848  double& lambda_min, double& lambda_max,const unsigned int * RngSeed)
849 {
850  IFPACK_CHK_ERR(-1);// NTS: This always seems to yield errors in AztecOO, ergo,
851  // I turned it off.
852 
853  // Prevent build warning for unreachable statements.
854 #if 0
855 
856  if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
857 
858 #ifdef HAVE_IFPACK_AZTECOO
859  Epetra_Vector x(Operator_->OperatorDomainMap());
860  Epetra_Vector y(Operator_->OperatorRangeMap());
861  if(RngSeed) x.SetSeed(*RngSeed);
862  x.Random();
863  y.PutScalar(0.0);
864  Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
865 
866  AztecOO solver(LP);
867  solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
868  solver.SetAztecOption(AZ_output, AZ_none);
869 
870  solver.SetPrecOperator(&*InvBlockDiagonal_);
871  solver.Iterate(MaximumIterations, 1e-10);
872 
873  const double* status = solver.GetAztecStatus();
874 
875  lambda_min = status[AZ_lambda_min];
876  lambda_max = status[AZ_lambda_max];
877 
878  return(0);
879 #else
880  using std::cout;
881  using std::endl;
882 
883  cout << "You need to configure IFPACK with support for AztecOO" << endl;
884  cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
885  cout << "in your configure script." << endl;
886  IFPACK_CHK_ERR(-1);
887 #endif
888 
889 #endif // 0
890 }
891 #endif // HAVE_IFPACK_EPETRAEXT
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 int Compute()
Computes the preconditioners.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
static int CG(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &lambda_min, double &lambda_max, const unsigned int *RngSeed=0)
Uses AztecOO&#39;s CG to estimate lambda_min and lambda_max.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
Ifpack_Chebyshev(const Epetra_Operator *Matrix)
Ifpack_Chebyshev constructor with given Epetra_Operator/Epetra_RowMatrix.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
static int PowerMethod(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &LambdaMax, const unsigned int *RngSeed=0)
Simple power method to compute lambda_max.