Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_Polynomial.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_Polynomial.h"
53#include "Ifpack_Utils.h"
54#include "Ifpack_Condest.h"
55#include "Teuchos_LAPACK.hpp"
56#include "Teuchos_SerialDenseMatrix.hpp"
57#include <complex>
58#ifdef HAVE_IFPACK_AZTECOO
60#include "AztecOO.h"
61#endif
62
63#ifdef HAVE_IFPACK_EPETRAEXT
64#include "Epetra_CrsMatrix.h"
65#include "EpetraExt_PointToBlockDiagPermute.h"
66#endif
67
68
69#define ABS(x) ((x)>0?(x):-(x))
70
71// Helper function for normal equations
72inline void Apply_Transpose(Teuchos::RCP<const Epetra_Operator> Operator_,const Epetra_MultiVector &X,Epetra_MultiVector &Y){
73 Epetra_Operator * Operator=const_cast<Epetra_Operator*>(&*Operator_);
74 Operator->SetUseTranspose(true);
75 Operator->Apply(X,Y);
76 Operator->SetUseTranspose(false);
77}
78
79
80
81
82//==============================================================================
83// NOTE: any change to the default values should be committed to the other
84// constructor as well.
86Ifpack_Polynomial(const Epetra_Operator* Operator) :
87 IsInitialized_(false),
88 IsComputed_(false),
89 IsIndefinite_(false),
90 IsComplex_(false),
91 NumInitialize_(0),
92 NumCompute_(0),
93 NumApplyInverse_(0),
94 InitializeTime_(0.0),
95 ComputeTime_(0.0),
96 ApplyInverseTime_(0.0),
97 ComputeFlops_(0.0),
98 ApplyInverseFlops_(0.0),
99 PolyDegree_(3),
100 LSPointsReal_(10),
101 LSPointsImag_(10),
102 UseTranspose_(false),
103 Condest_(-1.0),
104 /* ComputeCondest_(false), (unused; commented out to avoid build warnings) */
105 RealEigRatio_(10.0),
106 ImagEigRatio_(10.0),
107 Label_(),
108 LambdaRealMin_(0.0),
109 LambdaRealMax_(-1.0),
110 LambdaImagMin_(0.0),
111 LambdaImagMax_(0.0),
112 MinDiagonalValue_(0.0),
113 NumMyRows_(0),
114 NumMyNonzeros_(0),
115 NumGlobalRows_(0),
116 NumGlobalNonzeros_(0),
117 Operator_(Teuchos::rcp(Operator,false)),
118 UseBlockMode_(false),
119 SolveNormalEquations_(false),
120 IsRowMatrix_(false),
121 ZeroStartingSolution_(true)
122{
123}
124
125//==============================================================================
126// NOTE: This constructor has been introduced because SWIG does not appear
127// to appreciate dynamic_cast. An instruction of type
128// Matrix_ = dynamic_cast<const Epetra_RowMatrix*> in the
129// other construction does not work in PyTrilinos -- of course
130// it does in any C++ code (for an Epetra_Operator that is also
131// an Epetra_RowMatrix).
132//
133// FIXME: move declarations into a separate method?
135Ifpack_Polynomial(const Epetra_RowMatrix* Operator) :
136 IsInitialized_(false),
137 IsComputed_(false),
138 IsIndefinite_(false),
139 IsComplex_(false),
140 NumInitialize_(0),
141 NumCompute_(0),
142 NumApplyInverse_(0),
143 InitializeTime_(0.0),
144 ComputeTime_(0.0),
145 ApplyInverseTime_(0.0),
146 ComputeFlops_(0.0),
147 ApplyInverseFlops_(0.0),
148 PolyDegree_(3),
149 LSPointsReal_(10),
150 LSPointsImag_(10),
151 UseTranspose_(false),
152 Condest_(-1.0),
153 /* ComputeCondest_(false), (unused; commented out to avoid build warnings) */
154 RealEigRatio_(10.0),
155 ImagEigRatio_(10.0),
156 EigMaxIters_(10),
157 Label_(),
158 LambdaRealMin_(0.0),
159 LambdaRealMax_(-1.0),
160 LambdaImagMin_(0.0),
161 LambdaImagMax_(0.0),
162 MinDiagonalValue_(0.0),
163 NumMyRows_(0),
164 NumMyNonzeros_(0),
165 NumGlobalRows_(0),
166 NumGlobalNonzeros_(0),
167 Operator_(Teuchos::rcp(Operator,false)),
168 Matrix_(Teuchos::rcp(Operator,false)),
169 UseBlockMode_(false),
170 SolveNormalEquations_(false),
171 IsRowMatrix_(true),
172 ZeroStartingSolution_(true)
173{
174}
175
176//==============================================================================
177int Ifpack_Polynomial::SetParameters(Teuchos::ParameterList& List)
178{
179
180 RealEigRatio_ = List.get("polynomial: real eigenvalue ratio", RealEigRatio_);
181 ImagEigRatio_ = List.get("polynomial: imag eigenvalue ratio", ImagEigRatio_);
182 LambdaRealMin_ = List.get("polynomial: min real part", LambdaRealMin_);
183 LambdaRealMax_ = List.get("polynomial: max real part", LambdaRealMax_);
184 LambdaImagMin_ = List.get("polynomial: min imag part", LambdaImagMin_);
185 LambdaImagMax_ = List.get("polynomial: max imag part", LambdaImagMax_);
186 PolyDegree_ = List.get("polynomial: degree",PolyDegree_);
187 LSPointsReal_ = List.get("polynomial: real interp points",LSPointsReal_);
188 LSPointsImag_ = List.get("polynomial: imag interp points",LSPointsImag_);
189 IsIndefinite_ = List.get("polynomial: indefinite",IsIndefinite_);
190 IsComplex_ = List.get("polynomial: complex",IsComplex_);
191 MinDiagonalValue_ = List.get("polynomial: min diagonal value",
193 ZeroStartingSolution_ = List.get("polynomial: zero starting solution",
195
196 Epetra_Vector* ID = List.get("polynomial: operator inv diagonal",
197 (Epetra_Vector*)0);
198 EigMaxIters_ = List.get("polynomial: eigenvalue max iterations",EigMaxIters_);
199
200#ifdef HAVE_IFPACK_EPETRAEXT
201 // This is *always* false if EpetraExt isn't enabled
202 UseBlockMode_ = List.get("polynomial: use block mode",UseBlockMode_);
203 if(!List.isParameter("polynomial: block list")) UseBlockMode_=false;
204 else{
205 BlockList_ = List.get("polynomial: block list",BlockList_);
206
207 // Since we know we're doing a matrix inverse, clobber the block list
208 // w/"invert" if it's set to multiply
209 Teuchos::ParameterList Blist;
210 Blist=BlockList_.get("blockdiagmatrix: list",Blist);
211 std::string dummy("invert");
212 std::string ApplyMode=Blist.get("apply mode",dummy);
213 if(ApplyMode==std::string("multiply")){
214 Blist.set("apply mode","invert");
215 BlockList_.set("blockdiagmatrix: list",Blist);
216 }
217 }
218#endif
219
220 SolveNormalEquations_ = List.get("polynomial: solve normal equations",SolveNormalEquations_);
221
222 if (ID != 0)
223 {
224 InvDiagonal_ = Teuchos::rcp( new Epetra_Vector(*ID) );
225 }
226
227 SetLabel();
228
229 return(0);
230}
231
232//==============================================================================
234{
235 return(Operator_->Comm());
236}
237
238//==============================================================================
240{
241 return(Operator_->OperatorDomainMap());
242}
243
244//==============================================================================
246{
247 return(Operator_->OperatorRangeMap());
248}
249
250//==============================================================================
253{
254 if (IsComputed() == false)
255 IFPACK_CHK_ERR(-3);
256
257 if (X.NumVectors() != Y.NumVectors())
258 IFPACK_CHK_ERR(-2);
259
260 if (IsRowMatrix_)
261 {
262 IFPACK_CHK_ERR(Matrix_->Multiply(UseTranspose(),X,Y));
263 }
264 else
265 {
266 IFPACK_CHK_ERR(Operator_->Apply(X,Y));
267 }
268
269 return(0);
270}
271
272//==============================================================================
274{
275 IsInitialized_ = false;
276
277 if (Operator_ == Teuchos::null)
278 IFPACK_CHK_ERR(-2);
279
280 if (Time_ == Teuchos::null)
281 Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
282
283 if (IsRowMatrix_)
284 {
285 if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
286 IFPACK_CHK_ERR(-2); // only square matrices
287
288 NumMyRows_ = Matrix_->NumMyRows();
289 NumMyNonzeros_ = Matrix_->NumMyNonzeros();
290 NumGlobalRows_ = Matrix_->NumGlobalRows64();
291 NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
292 }
293 else
294 {
295 if (Operator_->OperatorDomainMap().NumGlobalElements64() !=
296 Operator_->OperatorRangeMap().NumGlobalElements64())
297 IFPACK_CHK_ERR(-2); // only square operators
298 }
299
301 InitializeTime_ += Time_->ElapsedTime();
302 IsInitialized_ = true;
303 return(0);
304}
305
306//==============================================================================
308{
309 if (!IsInitialized())
311
312 Time_->ResetStartTime();
313
314 // reset values
315 IsComputed_ = false;
316 Condest_ = -1.0;
317
318 if (PolyDegree_ <= 0)
319 IFPACK_CHK_ERR(-2); // at least one application
320
321#ifdef HAVE_IFPACK_EPETRAEXT
322 // Check to see if we can run in block mode
323 if(IsRowMatrix_ && InvDiagonal_ == Teuchos::null && UseBlockMode_){
324 const Epetra_CrsMatrix *CrsMatrix=dynamic_cast<const Epetra_CrsMatrix*>(&*Matrix_);
325
326 // If we don't have CrsMatrix, we can't use the block preconditioner
327 if(!CrsMatrix) UseBlockMode_=false;
328 else{
329 int ierr;
330 InvBlockDiagonal_=Teuchos::rcp(new EpetraExt_PointToBlockDiagPermute(*CrsMatrix));
331 if(InvBlockDiagonal_==Teuchos::null) IFPACK_CHK_ERR(-6);
332
333 ierr=InvBlockDiagonal_->SetParameters(BlockList_);
334 if(ierr) IFPACK_CHK_ERR(-7);
335
336 ierr=InvBlockDiagonal_->Compute();
337 if(ierr) IFPACK_CHK_ERR(-8);
338 }
339 }
340#endif
341 if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_)
342 {
343 InvDiagonal_ = Teuchos::rcp( new Epetra_Vector(Matrix().Map()) );
344
345 if (InvDiagonal_ == Teuchos::null)
346 IFPACK_CHK_ERR(-5);
347
348 IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*InvDiagonal_));
349
350 // Inverse diagonal elements
351 // Replace zeros with 1.0
352 for (int i = 0 ; i < NumMyRows_ ; ++i) {
353 double diag = (*InvDiagonal_)[i];
354 if (IFPACK_ABS(diag) < MinDiagonalValue_)
355 (*InvDiagonal_)[i] = MinDiagonalValue_;
356 else
357 (*InvDiagonal_)[i] = 1.0 / diag;
358 }
359 }
360
361 // Automatically compute maximum eigenvalue estimate of D^{-1}A if user hasn't provided one
362 double lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max;
363 if (LambdaRealMax_ == -1) {
364 //PowerMethod(Matrix(), *InvDiagonal_, EigMaxIters_, lambda_max);
365 GMRES(Matrix(), *InvDiagonal_, EigMaxIters_, lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max);
366 LambdaRealMin_=lambda_real_min; LambdaImagMin_=lambda_imag_min;
367 LambdaRealMax_=lambda_real_max; LambdaImagMax_=lambda_imag_max;
368 //std::cout<<"LambdaRealMin: "<<LambdaRealMin_<<std::endl;
369 //std::cout<<"LambdaRealMax: "<<LambdaRealMax_<<std::endl;
370 //std::cout<<"LambdaImagMin: "<<LambdaImagMin_<<std::endl;
371 //std::cout<<"LambdaImagMax: "<<LambdaImagMax_<<std::endl;
372 }
373
374 // find least squares polynomial for (LSPointsReal_*LSPointsImag_) zeros
375 // on a rectangle in the complex plane defined as
376 // [LambdaRealMin_,LambdaRealMax_] x [LambdaImagMin_,LambdaImagMax_]
377
378 const std::complex<double> zero(0.0,0.0);
379
380 // Compute points in complex plane
381 double lenx = LambdaRealMax_-LambdaRealMin_;
382 int nx = ceil(lenx*((double) LSPointsReal_));
383 if (nx<2) { nx = 2; }
384 double hx = lenx/((double) nx);
385 std::vector<double> xs;
386 if(abs(lenx)>1.0e-8) {
387 for( int pt=0; pt<=nx; pt++ ) {
388 xs.push_back(hx*pt+LambdaRealMin_);
389 }
390 }
391 else {
392 xs.push_back(LambdaRealMax_);
393 nx=1;
394 }
395 double leny = LambdaImagMax_-LambdaImagMin_;
396 int ny = ceil(leny*((double) LSPointsImag_));
397 if (ny<2) { ny = 2; }
398 double hy = leny/((double) ny);
399 std::vector<double> ys;
400 if(abs(leny)>1.0e-8) {
401 for( int pt=0; pt<=ny; pt++ ) {
402 ys.push_back(hy*pt+LambdaImagMin_);
403 }
404 }
405 else {
406 ys.push_back(LambdaImagMax_);
407 ny=1;
408 }
409 std::vector< std::complex<double> > cpts;
410 for( int jj=0; jj<ny; jj++ ) {
411 for( int ii=0; ii<nx; ii++ ) {
412 std::complex<double> cpt(xs[ii],ys[jj]);
413 cpts.push_back(cpt);
414 }
415 }
416 cpts.push_back(zero);
417
418#ifdef HAVE_TEUCHOS_COMPLEX
419 const std::complex<double> one(1.0,0.0);
420
421 // Construct overdetermined Vandermonde matrix
422 Teuchos::SerialDenseMatrix<int, std::complex<double> > Vmatrix(cpts.size(),PolyDegree_+1);
423 Vmatrix.putScalar(zero);
424 for (int jj = 0; jj <= PolyDegree_; ++jj) {
425 for (int ii = 0; ii < static_cast<int> (cpts.size ()) - 1; ++ii) {
426 if (jj > 0) {
427 Vmatrix(ii,jj) = pow(cpts[ii],jj);
428 }
429 else {
430 Vmatrix(ii,jj) = one;
431 }
432 }
433 }
434 Vmatrix(cpts.size()-1,0)=one;
435
436 // Right hand side: all zero except last entry
437 Teuchos::SerialDenseMatrix< int,std::complex<double> > RHS(cpts.size(),1);
438 RHS.putScalar(zero);
439 RHS(cpts.size()-1,0)=one;
440
441 // Solve least squares problem using LAPACK
442 Teuchos::LAPACK< int, std::complex<double> > lapack;
443 const int N = Vmatrix.numCols();
444 Teuchos::Array<double> singularValues(N);
445 Teuchos::Array<double> rwork(1);
446 rwork.resize (std::max (1, 5 * N));
447 std::complex<double> lworkScalar(1.0,0.0);
448 int info = 0;
449 lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
450 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
451 &lworkScalar, -1, &info);
452 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
453 "_GELSS workspace query returned INFO = "
454 << info << " != 0.");
455 const int lwork = static_cast<int> (real(lworkScalar));
456 TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
457 "_GELSS workspace query returned LWORK = "
458 << lwork << " < 0.");
459 // Allocate workspace. Size > 0 means &work[0] makes sense.
460 Teuchos::Array< std::complex<double> > work (std::max (1, lwork));
461 // Solve the least-squares problem.
462 lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
463 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
464 &work[0], lwork, &info);
465
466 coeff_.resize(PolyDegree_+1);
467 std::complex<double> c0=RHS(0,0);
468 for(int ii=0; ii<=PolyDegree_; ii++) {
469 // test that the imaginary part is nonzero
470 //TEUCHOS_TEST_FOR_EXCEPTION(abs(imag(RHS(ii,0))) > 1e-8, std::logic_error,
471 // "imaginary part of polynomial coefficients is nonzero! coeff = "
472 // << RHS(ii,0));
473 coeff_[ii]=real(RHS(ii,0)/c0);
474 //std::cout<<"coeff["<<ii<<"]="<<coeff_[ii]<<std::endl;
475 }
476
477#else
478
479 // Construct overdetermined Vandermonde matrix
480 Teuchos::SerialDenseMatrix< int, double > Vmatrix(xs.size()+1,PolyDegree_+1);
481 Vmatrix.putScalar(0.0);
482 for( int jj=0; jj<=PolyDegree_; jj++) {
483 for( std::vector<double>::size_type ii=0; ii<xs.size(); ii++) {
484 if(jj>0) {
485 Vmatrix(ii,jj)=pow(xs[ii],jj);
486 }
487 else {
488 Vmatrix(ii,jj)=1.0;
489 }
490 }
491 }
492 Vmatrix(xs.size(),0)=1.0;
493
494 // Right hand side: all zero except last entry
495 Teuchos::SerialDenseMatrix< int, double > RHS(xs.size()+1,1);
496 RHS.putScalar(0.0);
497 RHS(xs.size(),0)=1.0;
498
499 // Solve least squares problem using LAPACK
500 Teuchos::LAPACK< int, double > lapack;
501 const int N = Vmatrix.numCols();
502 Teuchos::Array<double> singularValues(N);
503 Teuchos::Array<double> rwork(1);
504 rwork.resize (std::max (1, 5 * N));
505 double lworkScalar(1.0);
506 int info = 0;
507 lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
508 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
509 &lworkScalar, -1, &info);
510 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
511 "_GELSS workspace query returned INFO = "
512 << info << " != 0.");
513 const int lwork = static_cast<int> (lworkScalar);
514 TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
515 "_GELSS workspace query returned LWORK = "
516 << lwork << " < 0.");
517 // Allocate workspace. Size > 0 means &work[0] makes sense.
518 Teuchos::Array< double > work (std::max (1, lwork));
519 // Solve the least-squares problem.
520 lapack.GELS('N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
521 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
522 &work[0], lwork, &info);
523
524 coeff_.resize(PolyDegree_+1);
525 double c0=RHS(0,0);
526 for(int ii=0; ii<=PolyDegree_; ii++) {
527 // test that the imaginary part is nonzero
528 //TEUCHOS_TEST_FOR_EXCEPTION(abs(imag(RHS(ii,0))) > 1e-8, std::logic_error,
529 // "imaginary part of polynomial coefficients is nonzero! coeff = "
530 // << RHS(ii,0));
531 coeff_[ii]=RHS(ii,0)/c0;
532 }
533
534#endif
535
536#ifdef IFPACK_FLOPCOUNTERS
538#endif
539
540 ++NumCompute_;
541 ComputeTime_ += Time_->ElapsedTime();
542 IsComputed_ = true;
543
544 return(0);
545}
546
547//==============================================================================
548std::ostream& Ifpack_Polynomial::Print(std::ostream & os) const
549{
550 using std::endl;
551
552 double MyMinVal, MyMaxVal;
553 double MinVal, MaxVal;
554
555 if (IsComputed_) {
556 InvDiagonal_->MinValue(&MyMinVal);
557 InvDiagonal_->MaxValue(&MyMaxVal);
558 Comm().MinAll(&MyMinVal,&MinVal,1);
559 Comm().MinAll(&MyMaxVal,&MaxVal,1);
560 }
561
562 if (!Comm().MyPID()) {
563 os << endl;
564 os << "================================================================================" << endl;
565 os << "Ifpack_Polynomial" << endl;
566 os << "Degree of polynomial = " << PolyDegree_ << endl;
567 os << "Condition number estimate = " << Condest() << endl;
568 os << "Global number of rows = " << Operator_->OperatorRangeMap().NumGlobalElements64() << endl;
569 if (IsComputed_) {
570 os << "Minimum value on stored inverse diagonal = " << MinVal << endl;
571 os << "Maximum value on stored inverse diagonal = " << MaxVal << endl;
572 }
574 os << "Using zero starting solution" << endl;
575 else
576 os << "Using input starting solution" << endl;
577 os << endl;
578 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
579 os << "----- ------- -------------- ------------ --------" << endl;
580 os << "Initialize() " << std::setw(5) << NumInitialize_
581 << " " << std::setw(15) << InitializeTime_
582 << " 0.0 0.0" << endl;
583 os << "Compute() " << std::setw(5) << NumCompute_
584 << " " << std::setw(15) << ComputeTime_
585 << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
586 if (ComputeTime_ != 0.0)
587 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
588 else
589 os << " " << std::setw(15) << 0.0 << endl;
590 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
591 << " " << std::setw(15) << ApplyInverseTime_
592 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
593 if (ApplyInverseTime_ != 0.0)
594 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
595 else
596 os << " " << std::setw(15) << 0.0 << endl;
597 os << "================================================================================" << endl;
598 os << endl;
599 }
600
601 return(os);
602}
603
604//==============================================================================
607 const int MaxIters, const double Tol,
608 Epetra_RowMatrix* Matrix_in)
609{
610 if (!IsComputed()) // cannot compute right now
611 return(-1.0);
612
613 // always computes it. Call Condest() with no parameters to get
614 // the previous estimate.
615 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
616
617 return(Condest_);
618}
619
620//==============================================================================
622{
623 Label_ = "IFPACK (Least squares polynomial), degree=" + Ifpack_toString(PolyDegree_);
624}
625
626//==============================================================================
629{
630
631 if (!IsComputed())
632 IFPACK_CHK_ERR(-3);
633
634 if (PolyDegree_ == 0)
635 return 0;
636
637 int nVec = X.NumVectors();
638 if (nVec != Y.NumVectors())
639 IFPACK_CHK_ERR(-2);
640
641 Time_->ResetStartTime();
642
643 Epetra_MultiVector Xcopy(X);
644 if(ZeroStartingSolution_==true) {
645 Y.PutScalar(0.0);
646 }
647
648 // mfh 20 Mar 2014: IBD never gets used, so I'm commenting out the
649 // following lines of code in order to forestall build warnings.
650// #ifdef HAVE_IFPACK_EPETRAEXT
651// EpetraExt_PointToBlockDiagPermute* IBD=0;
652// if (UseBlockMode_) IBD=&*InvBlockDiagonal_;
653// #endif
654
655 Y.Update(-coeff_[1], Xcopy, 1.0);
656 for (int ii = 2; ii < static_cast<int> (coeff_.size ()); ++ii) {
657 const Epetra_MultiVector V(Xcopy);
658 Operator_->Apply(V,Xcopy);
659 Xcopy.Multiply(1.0, *InvDiagonal_, Xcopy, 0.0);
660 // Update Y
661 Y.Update(-coeff_[ii], Xcopy, 1.0);
662 }
663
664 // Flops are updated in each of the following.
666 ApplyInverseTime_ += Time_->ElapsedTime();
667 return(0);
668}
669
670//==============================================================================
672PowerMethod(const Epetra_Operator& Operator,
673 const Epetra_Vector& InvPointDiagonal,
674 const int MaximumIterations,
675 double& lambda_max)
676{
677 // this is a simple power method
678 lambda_max = 0.0;
679 double RQ_top, RQ_bottom, norm;
680 Epetra_Vector x(Operator.OperatorDomainMap());
681 Epetra_Vector y(Operator.OperatorRangeMap());
682 x.Random();
683 x.Norm2(&norm);
684 if (norm == 0.0) IFPACK_CHK_ERR(-1);
685
686 x.Scale(1.0 / norm);
687
688 for (int iter = 0; iter < MaximumIterations; ++iter)
689 {
690 Operator.Apply(x, y);
691 IFPACK_CHK_ERR(y.Multiply(1.0, InvPointDiagonal, y, 0.0));
692 IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
693 IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
694 lambda_max = RQ_top / RQ_bottom;
695 IFPACK_CHK_ERR(y.Norm2(&norm));
696 if (norm == 0.0) IFPACK_CHK_ERR(-1);
697 IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
698 }
699
700 return(0);
701}
702
703//==============================================================================
705CG(const Epetra_Operator& Operator,
706 const Epetra_Vector& InvPointDiagonal,
707 const int MaximumIterations,
708 double& lambda_min, double& lambda_max)
709{
710#ifdef HAVE_IFPACK_AZTECOO
711 Epetra_Vector x(Operator.OperatorDomainMap());
712 Epetra_Vector y(Operator.OperatorRangeMap());
713 x.Random();
714 y.PutScalar(0.0);
715
716 Epetra_LinearProblem LP(const_cast<Epetra_Operator*>(&Operator), &x, &y);
717 AztecOO solver(LP);
718 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
719 solver.SetAztecOption(AZ_output, AZ_none);
720
722 Operator.OperatorRangeMap(),
723 InvPointDiagonal);
724 solver.SetPrecOperator(&diag);
725 solver.Iterate(MaximumIterations, 1e-10);
726
727 const double* status = solver.GetAztecStatus();
728
729 lambda_min = status[AZ_lambda_min];
730 lambda_max = status[AZ_lambda_max];
731
732 return(0);
733#else
734 using std::cout;
735 using std::endl;
736
737 cout << "You need to configure IFPACK with support for AztecOO" << endl;
738 cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
739 cout << "in your configure script." << endl;
740 IFPACK_CHK_ERR(-1);
741#endif
742}
743
744//==============================================================================
745#ifdef HAVE_IFPACK_EPETRAEXT
747PowerMethod(const int MaximumIterations, double& lambda_max)
748{
749
751 // this is a simple power method
752 lambda_max = 0.0;
753 double RQ_top, RQ_bottom, norm;
754 Epetra_Vector x(Operator_->OperatorDomainMap());
755 Epetra_Vector y(Operator_->OperatorRangeMap());
756 Epetra_Vector z(Operator_->OperatorRangeMap());
757 x.Random();
758 x.Norm2(&norm);
759 if (norm == 0.0) IFPACK_CHK_ERR(-1);
760
761 x.Scale(1.0 / norm);
762
763 for (int iter = 0; iter < MaximumIterations; ++iter)
764 {
765 Operator_->Apply(x, z);
766 InvBlockDiagonal_->ApplyInverse(z,y);
768 InvBlockDiagonal_->ApplyInverse(y,z);
770 }
771
772 IFPACK_CHK_ERR(y.Dot(x, &RQ_top));
773 IFPACK_CHK_ERR(x.Dot(x, &RQ_bottom));
774 lambda_max = RQ_top / RQ_bottom;
775 IFPACK_CHK_ERR(y.Norm2(&norm));
776 if (norm == 0.0) IFPACK_CHK_ERR(-1);
777 IFPACK_CHK_ERR(x.Update(1.0 / norm, y, 0.0));
778 }
779
780 return(0);
781}
782#endif
783
784//==============================================================================
785#ifdef HAVE_IFPACK_EPETRAEXT
787CG(const int /* MaximumIterations */,
788 double& /* lambda_min */, double& /* lambda_max */)
789{
790 IFPACK_CHK_ERR(-1);// NTS: This always seems to yield errors in AztecOO, ergo,
791 // I turned it off.
792
793 // mfh 06 Aug 2017: None of the code below here is reachable.
794 // This causes build warnings with some compilers. Thus, I'm
795 // commenting out this code.
796#if 0
797
799
800#ifdef HAVE_IFPACK_AZTECOO
801 Epetra_Vector x(Operator_->OperatorDomainMap());
802 Epetra_Vector y(Operator_->OperatorRangeMap());
803 x.Random();
804 y.PutScalar(0.0);
805 Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
806
807 AztecOO solver(LP);
808 solver.SetAztecOption(AZ_solver, AZ_cg_condnum);
809 solver.SetAztecOption(AZ_output, AZ_none);
810
811 solver.SetPrecOperator(&*InvBlockDiagonal_);
812 solver.Iterate(MaximumIterations, 1e-10);
813
814 const double* status = solver.GetAztecStatus();
815
816 lambda_min = status[AZ_lambda_min];
817 lambda_max = status[AZ_lambda_max];
818
819 return(0);
820#else
821 using std::cout;
822 using std::endl;
823
824 cout << "You need to configure IFPACK with support for AztecOO" << endl;
825 cout << "to use the CG estimator. This may require --enable-aztecoo" << endl;
826 cout << "in your configure script." << endl;
827 IFPACK_CHK_ERR(-1);
828#endif
829
830#endif // 0
831}
832#endif // HAVE_IFPACK_EPETRAEXT
833
834//==============================================================================
836GMRES(const Epetra_Operator& Operator,
837 const Epetra_Vector& InvPointDiagonal,
838 const int MaximumIterations,
839 double& lambda_real_min, double& lambda_real_max,
840 double& lambda_imag_min, double& lambda_imag_max)
841{
842#ifdef HAVE_IFPACK_AZTECOO
843 Epetra_Vector x(Operator_->OperatorDomainMap());
844 Epetra_Vector y(Operator_->OperatorRangeMap());
845 x.Random();
846 y.PutScalar(0.0);
847 Epetra_LinearProblem LP(const_cast<Epetra_RowMatrix*>(&*Matrix_), &x, &y);
848 AztecOO solver(LP);
849 solver.SetAztecOption(AZ_solver, AZ_gmres_condnum);
850 solver.SetAztecOption(AZ_output, AZ_none);
852 Operator.OperatorRangeMap(),
853 InvPointDiagonal);
854 solver.SetPrecOperator(&diag);
855 solver.Iterate(MaximumIterations, 1e-10);
856 const double* status = solver.GetAztecStatus();
857 lambda_real_min = status[AZ_lambda_real_min];
858 lambda_real_max = status[AZ_lambda_real_max];
859 lambda_imag_min = status[AZ_lambda_imag_min];
860 lambda_imag_max = status[AZ_lambda_imag_max];
861 return(0);
862#else
863 using std::cout;
864 using std::endl;
865
866 cout << "You need to configure IFPACK with support for AztecOO" << endl;
867 cout << "to use the GMRES estimator. This may require --enable-aztecoo" << endl;
868 cout << "in your configure script." << endl;
869 IFPACK_CHK_ERR(-1);
870#endif
871}
void Apply_Transpose(Teuchos::RCP< const Epetra_Operator > Operator_, const Epetra_MultiVector &X, Epetra_MultiVector &Y)
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
#define IFPACK_ABS(x)
void Apply_Transpose(Teuchos::RCP< const Epetra_Operator > Operator_, const Epetra_MultiVector &X, Epetra_MultiVector &Y)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
#define RHS(a)
Definition: MatGenFD.c:60
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const=0
int Scale(double ScalarValue)
int NumVectors() const
int Dot(const Epetra_MultiVector &A, double *Result) 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)
int Norm2(double *Result) const
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
Ifpack_DiagPreconditioner: a class for diagonal preconditioning.
bool ZeroStartingSolution_
If true, the starting solution is always the zero vector.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
bool SolveNormalEquations_
Run on the normal equations.
int LSPointsReal_
Contains the number of discretization points of the least squares problem.
Ifpack_Polynomial(const Epetra_Operator *Matrix)
Ifpack_Polynomial constructor with given Epetra_Operator/Epetra_RowMatrix.
bool IsComputed_
If true, the preconditioner has been computed successfully.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int Compute()
Computes the preconditioners.
int NumMyNonzeros_
Number of local nonzeros.
int EigMaxIters_
Max number of iterations to use in eigenvalue estimation (if automatic).
virtual void SetLabel()
Sets the label.
int NumCompute_
Contains the number of successful call to Compute().
int NumMyRows_
Number of local rows.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned as an Epetra_RowMatrix.
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.
int PolyDegree_
Contains the degree of the least squares polynomial.
bool UseBlockMode_
Use Block Preconditioning.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
long long NumGlobalRows_
Number of global rows.
int GMRES(const Epetra_Operator &Operator, const Epetra_Vector &InvPointDiagonal, const int MaximumIterations, double &lambda_real_min, double &lambda_real_max, double &lambda_imag_min, double &lambda_imag_max)
Uses AztecOO's GMRES to estimate the height and width of the spectrum.
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.
double Condest_
Contains the estimated condition number.
Teuchos::RefCountPtr< Epetra_Time > Time_
Time object to track timing.
bool IsIndefinite_
If true, have to compute polynomial for a spectrum with negative eigenvalues.
Teuchos::RefCountPtr< const Epetra_Operator > Operator_
Pointers to the matrix to be preconditioned as an Epetra_Operator.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
double MinDiagonalValue_
Contains the minimum value on the diagonal.
bool IsInitialized_
If true, the preconditioner has been computed successfully.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
long long NumGlobalNonzeros_
Number of global nonzeros.
double ComputeTime_
Contains the time for all successful calls to Compute().
std::string Label_
Contains the label of this object.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
double InitializeTime_
Contains the time for all successful calls to Initialize().
std::vector< double > coeff_
coefficients of the polynomial
double ComputeFlops_
Contains the number of flops for Compute().
double LambdaRealMin_
Bounds on the spectrum.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
bool IsRowMatrix_
If true, the Operator_ is an Epetra_RowMatrix.
Teuchos::RefCountPtr< Epetra_Vector > InvDiagonal_
Contains the inverse of diagonal elements of Matrix.
double RealEigRatio_
Contains the ratio such that the rectangular domain in the complex plane is [-LambdaRealMax_ / EigRat...
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
bool IsComplex_
If true, have to compute polynomial for a spectrum with nonzero imaginary part.
#define true
#define false
const int N