Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_Chebyshev.cpp
Go to the documentation of this file.
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
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
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.
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?
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 //==============================================================================
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",
169  ZeroStartingSolution_ = List.get("chebyshev: zero starting solution",
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
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  {
201  }
202 
203  SetLabel();
204 
205  return(0);
206 }
207 
208 //==============================================================================
210 {
211  return(Operator_->Comm());
212 }
213 
214 //==============================================================================
216 {
217  return(Operator_->OperatorDomainMap());
218 }
219 
220 //==============================================================================
222 {
223  return(Operator_->OperatorRangeMap());
224 }
225 
226 //==============================================================================
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())
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
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
333  }
334 #endif
335 
338 
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;
361  }
362  // otherwise the inverse of the diagonal has been given by the user
363  }
364 #ifdef IFPACK_FLOPCOUNTERS
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  }
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::
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 //==============================================================================
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 //==============================================================================
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
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
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.
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 
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);
828  InvBlockDiagonal_->ApplyInverse(y,z);
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
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
bool UseBlockMode_
Use Block Preconditioning.
virtual int Compute()
Computes the preconditioners.
int NumMyNonzeros_
Number of local nonzeros.
int NumCompute_
Contains the number of successful call to Compute().
double EigRatio_
Contains the ratio such that [LambdaMax_ / EigRatio_, LambdaMax_] is the interval of interest for the...
virtual int SetUseTranspose(bool UseTranspose)=0
bool ZeroStartingSolution_
If true, the starting solution is always the zero vector.
std::string Label_
Contains the label of this object.
void Apply_Transpose(Teuchos::RCP< const Epetra_Operator > Operator_, const Epetra_MultiVector &X, Epetra_MultiVector &Y)
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
T & get(ParameterList &l, const std::string &name)
int EigMaxIters_
Max number of iterations to use in eigenvalue estimation (if automatic).
long long NumGlobalNonzeros_
Number of global nonzeros.
long long NumGlobalRows_
Number of global rows.
double LambdaMin_
Contains an approximation to the smallest eigenvalue.
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.
double ComputeFlops_
Contains the number of flops for Compute().
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
virtual const Epetra_Map & OperatorDomainMap() const =0
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
int PolyDegree_
Contains the degree of Chebyshev polynomial.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
bool IsRowMatrix_
If true, the Operator_ is an Epetra_RowMatrix.
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
Teuchos::RefCountPtr< const Epetra_Operator > Operator_
Pointers to the matrix to be preconditioned as an Epetra_Operator.
bool isParameter(const std::string &name) const
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define ABS(x)
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual const Epetra_BlockMap & Map() const =0
double InitializeTime_
Contains the time for all successful calls to Initialize().
double Condest_
Contains the estimated condition number.
#define false
bool SolveNormalEquations_
Run on the normal equations.
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.
int NumMyRows_
Number of local rows.
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 UseTranspose() const
Returns the current UseTranspose setting.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
TransListIter iter
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual void SetLabel()
Sets the label.
bool IsInitialized_
If true, the preconditioner has been computed successfully.
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object to track timing.
#define IFPACK_ABS(x)
double LambdaMax_
Contains an approximation to the largest eigenvalue.
Ifpack_Chebyshev(const Epetra_Operator *Matrix)
Ifpack_Chebyshev constructor with given Epetra_Operator/Epetra_RowMatrix.
Teuchos::RefCountPtr< Epetra_Vector > InvDiagonal_
Contains the inverse of diagonal elements of Matrix.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
double ComputeTime_
Contains the time for all successful calls to Compute().
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.
#define IFPACK_CHK_ERR(ifpack_err)
bool IsComputed_
If true, the preconditioner has been computed successfully.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
#define true
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.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned as an Epetra_RowMatrix.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
double MinDiagonalValue_
Contains the minimum value on the diagonal.