IFPACK  Development
 All Classes Files Functions Variables Enumerations Friends Pages
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),
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),
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  string dummy("invert");
188  string ApplyMode=Blist.get("apply mode",dummy);
189  if(ApplyMode==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  PowerMethod(EigMaxIters_, lambda_max);
324  LambdaMax_ = lambda_max;
325 
326  // Test for Exact Preconditioned case
327  if (ABS(LambdaMax_-1) < 1e-6)
328  LambdaMax_ = LambdaMin_ = 1.0;
329  else
330  LambdaMin_ = LambdaMax_/EigRatio_;
331  }
332 #endif
333 
334  if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_) {
335  InvDiagonal_ = Teuchos::rcp(new Epetra_Vector(Matrix().Map()));
336 
337  if (InvDiagonal_ == Teuchos::null)
338  IFPACK_CHK_ERR(-5);
339 
340  IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*InvDiagonal_));
341 
342  // Inverse diagonal elements
343  // Replace zeros with 1.0
344  for (int i = 0 ; i < NumMyRows_ ; ++i) {
345  double diag = (*InvDiagonal_)[i];
346  if (IFPACK_ABS(diag) < MinDiagonalValue_)
347  (*InvDiagonal_)[i] = MinDiagonalValue_;
348  else
349  (*InvDiagonal_)[i] = 1.0 / diag;
350  }
351  // Automatically compute maximum eigenvalue estimate of D^{-1}A if user hasn't provided one
352  double lambda_max=0;
353  if (LambdaMax_ == -1) {
354  PowerMethod(Matrix(), *InvDiagonal_, EigMaxIters_, lambda_max);
355  LambdaMax_=lambda_max;
356  // Test for Exact Preconditioned case
357  if (ABS(LambdaMax_-1) < 1e-6) LambdaMax_=LambdaMin_=1.0;
358  else LambdaMin_=LambdaMax_/EigRatio_;
359  }
360  // otherwise the inverse of the diagonal has been given by the user
361  }
362 #ifdef IFPACK_FLOPCOUNTERS
363  ComputeFlops_ += NumMyRows_;
364 #endif
365 
366  ++NumCompute_;
367  ComputeTime_ += Time_->ElapsedTime();
368  IsComputed_ = true;
369 
370  return(0);
371 }
372 
373 //==============================================================================
374 ostream& Ifpack_Chebyshev::Print(ostream & os) const
375 {
376 
377  double MyMinVal, MyMaxVal;
378  double MinVal, MaxVal;
379 
380  if (IsComputed_) {
381  InvDiagonal_->MinValue(&MyMinVal);
382  InvDiagonal_->MaxValue(&MyMaxVal);
383  Comm().MinAll(&MyMinVal,&MinVal,1);
384  Comm().MinAll(&MyMaxVal,&MaxVal,1);
385  }
386 
387  if (!Comm().MyPID()) {
388  os << endl;
389  os << "================================================================================" << endl;
390  os << "Ifpack_Chebyshev" << endl;
391  os << "Degree of polynomial = " << PolyDegree_ << endl;
392  os << "Condition number estimate = " << Condest() << endl;
393  os << "Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
394  if (IsComputed_) {
395  os << "Minimum value on stored inverse diagonal = " << MinVal << endl;
396  os << "Maximum value on stored inverse diagonal = " << MaxVal << endl;
397  }
398  if (ZeroStartingSolution_)
399  os << "Using zero starting solution" << endl;
400  else
401  os << "Using input starting solution" << endl;
402  os << endl;
403  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
404  os << "----- ------- -------------- ------------ --------" << endl;
405  os << "Initialize() " << std::setw(5) << NumInitialize_
406  << " " << std::setw(15) << InitializeTime_
407  << " 0.0 0.0" << endl;
408  os << "Compute() " << std::setw(5) << NumCompute_
409  << " " << std::setw(15) << ComputeTime_
410  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
411  if (ComputeTime_ != 0.0)
412  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
413  else
414  os << " " << std::setw(15) << 0.0 << endl;
415  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
416  << " " << std::setw(15) << ApplyInverseTime_
417  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
418  if (ApplyInverseTime_ != 0.0)
419  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
420  else
421  os << " " << std::setw(15) << 0.0 << endl;
422  os << "================================================================================" << endl;
423  os << endl;
424  }
425 
426  return(os);
427 }
428 
429 //==============================================================================
430 double Ifpack_Chebyshev::
431 Condest(const Ifpack_CondestType CT,
432  const int MaxIters, const double Tol,
433  Epetra_RowMatrix* Matrix_in)
434 {
435  if (!IsComputed()) // cannot compute right now
436  return(-1.0);
437 
438  // always computes it. Call Condest() with no parameters to get
439  // the previous estimate.
440  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
441 
442  return(Condest_);
443 }
444 
445 //==============================================================================
446 void Ifpack_Chebyshev::SetLabel()
447 {
448  Label_ = "IFPACK (Chebyshev polynomial), degree=" + Ifpack_toString(PolyDegree_);
449 }
450 
451 //==============================================================================
453 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
454 {
455  if (!IsComputed())
456  IFPACK_CHK_ERR(-3);
457 
458  if (PolyDegree_ == 0)
459  return 0;
460 
461  int nVec = X.NumVectors();
462  int len = X.MyLength();
463  if (nVec != Y.NumVectors())
464  IFPACK_CHK_ERR(-2);
465 
466  Time_->ResetStartTime();
467 
468  // AztecOO gives X and Y pointing to the same memory location,
469  // need to create an auxiliary vector, Xcopy
470  Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
471  if (X.Pointers()[0] == Y.Pointers()[0])
472  Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
473  else
474  Xcopy = Teuchos::rcp( &X, false );
475 
476  double **xPtr = 0, **yPtr = 0;
477  Xcopy->ExtractView(&xPtr);
478  Y.ExtractView(&yPtr);
479 
480 #ifdef HAVE_IFPACK_EPETRAEXT
481  EpetraExt_PointToBlockDiagPermute* IBD=0;
482  if (UseBlockMode_)
483  IBD = &*InvBlockDiagonal_;
484 #endif
485 
486  //--- Do a quick solve when the matrix is identity
487  double *invDiag = 0;
488  if (!UseBlockMode_)
489  invDiag = InvDiagonal_->Values();
490 
491  if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_)) {
492 #ifdef HAVE_IFPACK_EPETRAEXT
493  if (UseBlockMode_)
494  IBD->ApplyInverse(*Xcopy, Y);
495  else
496 #endif
497  if (nVec == 1) {
498  double *yPointer = yPtr[0], *xPointer = xPtr[0];
499  for (int i = 0; i < len; ++i)
500  yPointer[i] = xPointer[i] * invDiag[i];
501 
502  } else {
503  for (int i = 0; i < len; ++i) {
504  double coeff = invDiag[i];
505  for (int k = 0; k < nVec; ++k)
506  yPtr[k][i] = xPtr[k][i] * coeff;
507  }
508  } // if (nVec == 1)
509 
510  return 0;
511  } // if ((LambdaMin_ == 1.0) && (LambdaMax_ == LambdaMin_))
512 
513  //--- Initialize coefficients
514  // Note that delta stores the inverse of ML_Cheby::delta
515  double alpha = LambdaMax_ / EigRatio_;
516  double beta = 1.1 * LambdaMax_;
517  double delta = 2.0 / (beta - alpha);
518  double theta = 0.5 * (beta + alpha);
519  double s1 = theta * delta;
520 
521  // Temporary vectors
522  // In ML_Cheby, V corresponds to pAux and W to dk
523  // NOTE: we would like to move the construction to the Compute()
524  // call, but that is not possible as we don't know how many
525  // vectors are in the multivector
526  bool zeroOut = false;
527  Epetra_MultiVector V(X.Map(), X.NumVectors(), zeroOut);
528  Epetra_MultiVector W(X.Map(), X.NumVectors(), zeroOut);
529 #ifdef HAVE_IFPACK_EPETRAEXT
530  Epetra_MultiVector Temp(X.Map(), X.NumVectors(), zeroOut);
531 #endif
532 
533  double *vPointer = V.Values(), *wPointer = W.Values();
534 
535  double oneOverTheta = 1.0/theta;
536 
537  //--- If solving normal equations, multiply RHS by A^T
538  if (SolveNormalEquations_) {
539  Apply_Transpose(Operator_, Y, V);
540  Y = V;
541  }
542 
543  // Do the smoothing when block scaling is turned OFF
544  // --- Treat the initial guess
545  if (ZeroStartingSolution_ == false) {
546  Operator_->Apply(Y, V);
547 
548  // Compute W = invDiag * ( X - V )/ Theta
549 #ifdef HAVE_IFPACK_EPETRAEXT
550  if (UseBlockMode_) {
551  Temp.Update(oneOverTheta, X, -oneOverTheta, V, 0.0);
552  IBD->ApplyInverse(Temp, W);
553 
554  // Perform additional matvecs for normal equations
555  // CMS: Testing this only in block mode FOR NOW
556  if (SolveNormalEquations_){
557  IBD->ApplyInverse(W, Temp);
558  Apply_Transpose(Operator_, Temp, W);
559  }
560  }
561  else
562 #endif
563  if (nVec == 1) {
564  double *xPointer = xPtr[0];
565  for (int i = 0; i < len; ++i)
566  wPointer[i] = invDiag[i] * (xPointer[i] - vPointer[i]) * oneOverTheta;
567 
568  } else {
569  for (int i = 0; i < len; ++i) {
570  double coeff = invDiag[i]*oneOverTheta;
571  double *wi = wPointer + i, *vi = vPointer + i;
572  for (int k = 0; k < nVec; ++k) {
573  *wi = (xPtr[k][i] - (*vi)) * coeff;
574  wi = wi + len; vi = vi + len;
575  }
576  }
577  } // if (nVec == 1)
578  // Update the vector Y
579  Y.Update(1.0, W, 1.0);
580 
581  } else { // if (ZeroStartingSolution_ == false)
582  // Compute W = invDiag * X / Theta
583 #ifdef HAVE_IFPACK_EPETRAEXT
584  if (UseBlockMode_) {
585  IBD->ApplyInverse(X, W);
586 
587  // Perform additional matvecs for normal equations
588  // CMS: Testing this only in block mode FOR NOW
589  if (SolveNormalEquations_) {
590  IBD->ApplyInverse(W, Temp);
591  Apply_Transpose(Operator_, Temp, W);
592  }
593 
594  W.Scale(oneOverTheta);
595  Y.Update(1.0, W, 0.0);
596  }
597  else
598 #endif
599  if (nVec == 1) {
600  double *xPointer = xPtr[0];
601  for (int i = 0; i < len; ++i)
602  wPointer[i] = invDiag[i] * xPointer[i] * oneOverTheta;
603 
604  memcpy(yPtr[0], wPointer, len*sizeof(double));
605 
606  } else {
607  for (int i = 0; i < len; ++i) {
608  double coeff = invDiag[i]*oneOverTheta;
609  double *wi = wPointer + i;
610  for (int k = 0; k < nVec; ++k) {
611  *wi = xPtr[k][i] * coeff;
612  wi = wi + len;
613  }
614  }
615 
616  for (int k = 0; k < nVec; ++k)
617  memcpy(yPtr[k], wPointer + k*len, len*sizeof(double));
618  } // if (nVec == 1)
619  } // if (ZeroStartingSolution_ == false)
620 
621  //--- Apply the polynomial
622  double rhok = 1.0/s1, rhokp1;
623  double dtemp1, dtemp2;
624  int degreeMinusOne = PolyDegree_ - 1;
625  if (nVec == 1) {
626  double *xPointer = xPtr[0];
627  for (int k = 0; k < degreeMinusOne; ++k) {
628  Operator_->Apply(Y, V);
629  rhokp1 = 1.0 / (2.0*s1 - rhok);
630  dtemp1 = rhokp1 * rhok;
631  dtemp2 = 2.0 * rhokp1 * delta;
632  rhok = rhokp1;
633 
634  // Compute W = dtemp1 * W
635  W.Scale(dtemp1);
636 
637  // Compute W = W + dtemp2 * invDiag * ( X - V )
638 #ifdef HAVE_IFPACK_EPETRAEXT
639  if (UseBlockMode_) {
640  //NTS: We can clobber V since it will be reset in the Apply
641  V.Update(dtemp2, X, -dtemp2);
642  IBD->ApplyInverse(V, Temp);
643 
644  // Perform additional matvecs for normal equations
645  // CMS: Testing this only in block mode FOR NOW
646  if (SolveNormalEquations_) {
647  IBD->ApplyInverse(V, Temp);
648  Apply_Transpose(Operator_, Temp, V);
649  }
650 
651  W.Update(1.0, Temp, 1.0);
652  }
653  else
654 #endif
655  for (int i = 0; i < len; ++i)
656  wPointer[i] += dtemp2* invDiag[i] * (xPointer[i] - vPointer[i]);
657 
658  // Update the vector Y
659  Y.Update(1.0, W, 1.0);
660  } // for (k = 0; k < degreeMinusOne; ++k)
661 
662  } else { // if (nVec == 1) {
663  for (int k = 0; k < degreeMinusOne; ++k) {
664  Operator_->Apply(Y, V);
665  rhokp1 = 1.0 / (2.0*s1 - rhok);
666  dtemp1 = rhokp1 * rhok;
667  dtemp2 = 2.0 * rhokp1 * delta;
668  rhok = rhokp1;
669 
670  // Compute W = dtemp1 * W
671  W.Scale(dtemp1);
672 
673  // Compute W = W + dtemp2 * invDiag * ( X - V )
674 #ifdef HAVE_IFPACK_EPETRAEXT
675  if (UseBlockMode_) {
676  // We can clobber V since it will be reset in the Apply
677  V.Update(dtemp2, X, -dtemp2);
678  IBD->ApplyInverse(V, Temp);
679 
680  // Perform additional matvecs for normal equations
681  // CMS: Testing this only in block mode FOR NOW
682  if (SolveNormalEquations_) {
683  IBD->ApplyInverse(V,Temp);
684  Apply_Transpose(Operator_,Temp,V);
685  }
686 
687  W.Update(1.0, Temp, 1.0);
688  }
689  else
690 #endif
691  for (int i = 0; i < len; ++i) {
692  double coeff = invDiag[i]*dtemp2;
693  double *wi = wPointer + i, *vi = vPointer + i;
694  for (int j = 0; j < nVec; ++j) {
695  *wi += (xPtr[j][i] - (*vi)) * coeff;
696  wi = wi + len; vi = vi + len;
697  }
698  }
699 
700  // Update the vector Y
701  Y.Update(1.0, W, 1.0);
702  } // for (k = 0; k < degreeMinusOne; ++k)
703  } // if (nVec == 1)
704 
705 
706  // Flops are updated in each of the following.
707  ++NumApplyInverse_;
708  ApplyInverseTime_ += Time_->ElapsedTime();
709 
710  return(0);
711 }
712 
713 //==============================================================================
715 PowerMethod(const Epetra_Operator& Operator,
716  const Epetra_Vector& InvPointDiagonal,
717  const int MaximumIterations,
718  double& lambda_max)
719 {
720  // this is a simple power method
721  lambda_max = 0.0;
722  double RQ_top, RQ_bottom, norm;
723  Epetra_Vector x(Operator.OperatorDomainMap());
724  Epetra_Vector y(Operator.OperatorRangeMap());
725  x.Random();
726  x.Norm2(&norm);
727  if (norm == 0.0) IFPACK_CHK_ERR(-1);
728 
729  x.Scale(1.0 / norm);
730 
731  for (int iter = 0; iter < MaximumIterations; ++iter)
732  {
733  Operator.Apply(x, y);
734  IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
735  IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
736  IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
737  lambda_max = RQ_top / RQ_bottom;
738  IFPACK_CHK_ERR(y.Norm2(&norm));
739  if (norm == 0.0) IFPACK_CHK_ERR(-1);
740  IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
741  }
742 
743  return(0);
744 }
745 
746 //==============================================================================
748 CG(const Epetra_Operator& Operator,
749  const Epetra_Vector& InvPointDiagonal,
750  const int MaximumIterations,
751  double& lambda_min, double& lambda_max)
752 {
753 #ifdef HAVE_IFPACK_AZTECOO
754  Epetra_Vector x(Operator.OperatorDomainMap());
755  Epetra_Vector y(Operator.OperatorRangeMap());
756  x.Random();
757  y.PutScalar(0.0);
758 
759  Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
760  AztecOO solver(LP);
761  solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
762  solver.SetAztecOption(AZ_output, AZ_none);
763 
764  Ifpack_DiagPreconditioner diag(Operator.OperatorDomainMap(),
765  Operator.OperatorRangeMap(),
766  InvPointDiagonal);
767  solver.SetPrecOperator(&diag);
768  solver.Iterate(MaximumIterations, 1e-10);
769 
770  const double* status = solver.GetAztecStatus();
771 
772  lambda_min = status[AZ_lambda_min];
773  lambda_max = status[AZ_lambda_max];
774 
775  return(0);
776 #else
777  cout << "You need to configure IFPACK with support for AztecOO" << endl;
778  cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
779  cout << "in your configure script." << endl;
780  IFPACK_CHK_ERR(-1);
781 #endif
782 }
783 
784 //==============================================================================
785 #ifdef HAVE_IFPACK_EPETRAEXT
787 PowerMethod(const int MaximumIterations, double& lambda_max)
788 {
789 
790  if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
791  // this is a simple power method
792  lambda_max = 0.0;
793  double RQ_top, RQ_bottom, norm;
794  Epetra_Vector x(Operator_->OperatorDomainMap());
795  Epetra_Vector y(Operator_->OperatorRangeMap());
796  Epetra_Vector z(Operator_->OperatorRangeMap());
797  x.Random();
798  x.Norm2(&norm);
799  if (norm == 0.0) IFPACK_CHK_ERR(-1);
800 
801  x.Scale(1.0 / norm);
802 
803  for (int iter = 0; iter < MaximumIterations; ++iter)
804  {
805  Operator_->Apply(x, z);
806  InvBlockDiagonal_->ApplyInverse(z,y);
807  if(SolveNormalEquations_){
808  InvBlockDiagonal_->ApplyInverse(y,z);
809  Apply_Transpose(Operator_,z, y);
810  }
811 
812  IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
813  IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
814  lambda_max = RQ_top / RQ_bottom;
815  IFPACK_CHK_ERR(y.Norm2(&norm));
816  if (norm == 0.0) IFPACK_CHK_ERR(-1);
817  IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
818  }
819 
820  return(0);
821 }
822 #endif
823 
824 //==============================================================================
825 #ifdef HAVE_IFPACK_EPETRAEXT
827 CG(const int MaximumIterations,
828  double& lambda_min, double& lambda_max)
829 {
830  IFPACK_CHK_ERR(-1);// NTS: This always seems to yield errors in AztecOO, ergo,
831  // I turned it off.
832 
833  if(!UseBlockMode_) IFPACK_CHK_ERR(-1);
834 
835 #ifdef HAVE_IFPACK_AZTECOO
836  Epetra_Vector x(Operator_->OperatorDomainMap());
837  Epetra_Vector y(Operator_->OperatorRangeMap());
838  x.Random();
839  y.PutScalar(0.0);
840  Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
841 
842  AztecOO solver(LP);
843  solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
844  solver.SetAztecOption(AZ_output, AZ_none);
845 
846  solver.SetPrecOperator(&*InvBlockDiagonal_);
847  solver.Iterate(MaximumIterations, 1e-10);
848 
849  const double* status = solver.GetAztecStatus();
850 
851  lambda_min = status[AZ_lambda_min];
852  lambda_max = status[AZ_lambda_max];
853 
854  return(0);
855 #else
856  cout << "You need to configure IFPACK with support for AztecOO" << endl;
857  cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
858  cout << "in your configure script." << endl;
859  IFPACK_CHK_ERR(-1);
860 #endif
861 }
862 #endif
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 ostream & Print(ostream &os) const
Prints object to an output stream.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
string Ifpack_toString(const int &x)
Convertes an integer to string.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
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_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
Ifpack_Chebyshev(const Epetra_Operator *Matrix)
Ifpack_Chebyshev constructor with given Epetra_Operator/Epetra_RowMatrix.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
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 int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.