IFPACK Development
Loading...
Searching...
No Matches
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
69inline 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.
83Ifpack_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?
125Ifpack_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//==============================================================================
160int 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//==============================================================================
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())
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//==============================================================================
378std::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//==============================================================================
436Condest(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//==============================================================================
451void 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//==============================================================================
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//==============================================================================
729PowerMethod(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//==============================================================================
763CG(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
806PowerMethod(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
847CG(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 int MinAll(double *PartialMins, double *GlobalMins, int Count) const=0
const Epetra_BlockMap & Map() const
int Scale(double ScalarValue)
int NumVectors() const
int Dot(const Epetra_MultiVector &A, double *Result) const
int MyLength() const
int ExtractView(double **A, int *MyLDA) const
double * Values() const
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
double ** Pointers() const
int Norm2(double *Result) const
int SetSeed(unsigned int Seed_in)
int PutScalar(double ScalarConstant)
virtual int SetUseTranspose(bool UseTranspose)=0
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
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 SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
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.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Ifpack_Chebyshev(const Epetra_Operator *Matrix)
Ifpack_Chebyshev constructor with given Epetra_Operator/Epetra_RowMatrix.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual int Compute()
Computes the preconditioners.
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's CG to estimate lambda_min and lambda_max.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.