Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Ifpack_SORa.cpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // Ifpack: Object-Oriented Algebraic Preconditioner Package
6 // Copyright (2002) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 #include "Ifpack_SORa.h"
45 #include "Ifpack.h"
46 #include "Ifpack_Utils.h"
47 #include "Teuchos_ParameterList.hpp"
48 #include "Teuchos_RCP.hpp"
49 #include "Epetra_Import.h"
50 #include "Epetra_Export.h"
51 #include "Epetra_IntVector.h"
52 
53 #ifdef HAVE_IFPACK_EPETRAEXT
54 #include "EpetraExt_MatrixMatrix.h"
55 #include "EpetraExt_RowMatrixOut.h"
56 #include "EpetraExt_MultiVectorOut.h"
57 #include "EpetraExt_Transpose_RowMatrix.h"
58 
59 
60 #define ABS(x) ((x)>=0?(x):-(x))
61 #define MIN(x,y) ((x)<(y)?(x):(y))
62 #define MAX(x,y) ((x)>(y)?(x):(y))
63 
64 // Useful functions horked from ML
65 int* FindLocalDiricheltRowsFromOnesAndZeros(const Epetra_CrsMatrix & Matrix, int &numBCRows);
66 void Apply_BCsToMatrixRowsAndColumns(const int *dirichletRows, int numBCRows,const Epetra_IntVector &dirichletColumns,const Epetra_CrsMatrix & Matrix);
67 Epetra_IntVector * FindLocalDirichletColumnsFromRows(const int *dirichletRows, int numBCRows,const Epetra_CrsMatrix & Matrix);
68 
69 Ifpack_SORa::Ifpack_SORa(Epetra_RowMatrix* A):
70  IsInitialized_(false),
71  IsComputed_(false),
72  Label_(),
73  Alpha_(1.5),
74  Gamma_(1.0),
75  NumSweeps_(1),
76  IsParallel_(false),
77  HaveOAZBoundaries_(false),
78  UseInterprocDamping_(false),
79  UseGlobalDamping_(false),
80  LambdaMax_(-1.0),
81  LambdaMaxBoost_(1.0),
82  PowerMethodIters_(10),
83  Time_(A->Comm())
84 {
85  Epetra_CrsMatrix *Acrs=dynamic_cast<Epetra_CrsMatrix*>(A);
86  if (Acrs) {
87  Acrs_= Teuchos::rcp (Acrs, false);
88  }
89  A_ = Teuchos::rcp (A, false);
90 }
91 
92 void Ifpack_SORa::Destroy(){
93 }
94 
95 
96 
97 int Ifpack_SORa::Initialize(){
98  Alpha_ = List_.get("sora: alpha", Alpha_);
99  Gamma_ = List_.get("sora: gamma",Gamma_);
100  NumSweeps_ = List_.get("sora: sweeps",NumSweeps_);
101  HaveOAZBoundaries_ = List_.get("sora: oaz boundaries", HaveOAZBoundaries_);
102  UseInterprocDamping_= List_.get("sora: use interproc damping",UseInterprocDamping_);
103  UseGlobalDamping_ = List_.get("sora: use global damping",UseGlobalDamping_);
104  LambdaMax_ = List_.get("sora: max eigenvalue",LambdaMax_);
105  LambdaMaxBoost_ = List_.get("sora: eigen-analysis: boost",LambdaMaxBoost_);
106  PowerMethodIters_ = List_.get("sora: eigen-analysis: max iters",PowerMethodIters_);
107 
108  if (A_->Comm().NumProc() != 1) IsParallel_ = true;
109  else {
110  IsParallel_ = false;
111  // Don't use interproc damping, for obvious reasons
112  UseInterprocDamping_=false;
113  }
114 
115  // Counters
116  IsInitialized_=true;
117  NumInitialize_++;
118  return 0;
119 }
120 
121 int Ifpack_SORa::SetParameters(Teuchos::ParameterList& parameterlist){
122  List_=parameterlist;
123  return 0;
124 }
125 
126 
127 template<typename int_type>
128 int Ifpack_SORa::TCompute(){
129  using std::cout;
130  using std::endl;
131 
132  Epetra_Map *RowMap=const_cast<Epetra_Map*>(&A_->RowMatrixRowMap());
133  Epetra_Vector Adiag(*RowMap);
134  Epetra_CrsMatrix *Askew2=0, *Aherm2=0,*W=0;
135  int *rowptr_s,*colind_s,*rowptr_h,*colind_h;
136  double *vals_s,*vals_h;
137  bool RowMatrixMode=(Acrs_==Teuchos::null);
138 
139  // Label
140  sprintf(Label_, "IFPACK SORa (alpha=%5.2f gamma=%5.2f)",Alpha_,Gamma_);
141 
142  if(RowMatrixMode){
143  if(!A_->Comm().MyPID()) cout<<"SORa: RowMatrix mode enabled"<<endl;
144  // RowMatrix mode, build Acrs_ the hard way.
145  Epetra_RowMatrix *Arow=&*A_;
146  Epetra_Map *ColMap=const_cast<Epetra_Map*>(&A_->RowMatrixColMap());
147 
148  int Nmax=A_->MaxNumEntries();
149  int length;
150  std::vector<double> values(Nmax);
151  Epetra_CrsMatrix *Acrs=new Epetra_CrsMatrix(Copy,*RowMap,Nmax);
152 
153  std::vector<int> indices_local(Nmax);
154 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
155  if(RowMap->GlobalIndicesInt()) {
156  for(int i=0;i<Arow->NumMyRows();i++) {
157  Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices_local[0]);
158  for(int j=0;j<length;j++)
159  indices_local[j]= ColMap->GID(indices_local[j]);
160  Acrs->InsertGlobalValues(RowMap->GID(i),length,&values[0],&indices_local[0]);
161  }
162  }
163  else
164 #endif
165 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
166  if(RowMap->GlobalIndicesLongLong()) {
167  std::vector<int_type> indices_global(Nmax);
168  for(int i=0;i<Arow->NumMyRows();i++) {
169  Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices_local[0]);
170  for(int j=0;j<length;j++)
171  indices_global[j]= (int_type) ColMap->GID64(indices_local[j]);
172  Acrs->InsertGlobalValues((int_type) RowMap->GID64(i),length,&values[0],&indices_global[0]);
173  }
174  }
175  else
176 #endif
177  throw "Ifpack_SORa::TCompute: GlobalIndices type unknown";
178 
179  Acrs->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
180  Acrs_ = Teuchos::rcp (Acrs, true);
181  }
182 
183  // Create Askew2
184  // Note: Here I let EpetraExt do the thinking for me. Since it gets all the maps correct for the E + F^T stencil.
185  // There are probably more efficient ways to do this but this method has the bonus of allowing code reuse.
186  IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,false,1,*Acrs_,true,-1,Askew2));
187  Askew2->FillComplete();
188 
189  // Create Aherm2
190  IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,false,1,*Acrs_,true,1,Aherm2));
191  Aherm2->FillComplete();
192 
193  int nnz2=Askew2->NumMyNonzeros();
194  int N=Askew2->NumMyRows();
195 
196 
197  // Grab pointers
198  IFPACK_CHK_ERR(Askew2->ExtractCrsDataPointers(rowptr_s,colind_s,vals_s));
199  IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));
200 
201  // Sanity checking: Make sure the sparsity patterns are the same
202 #define SANITY_CHECK
203 #ifdef SANITY_CHECK
204  for(int i=0;i<N;i++)
205  if(rowptr_s[i]!=rowptr_h[i]) IFPACK_CHK_ERR(-2);
206  for(int i=0;i<nnz2;i++)
207  if(colind_s[i]!=colind_h[i]) IFPACK_CHK_ERR(-3);
208 #endif
209 
210  // Dirichlet Detection & Nuking of Aherm2 and Askew2
211  // Note: Relies on Aherm2/Askew2 having identical sparsity patterns (see sanity check above)
212  if(HaveOAZBoundaries_){
213  int numBCRows;
214  int* dirRows=FindLocalDiricheltRowsFromOnesAndZeros(*Acrs_,numBCRows);
215  Epetra_IntVector* dirCols=FindLocalDirichletColumnsFromRows(dirRows,numBCRows,*Aherm2);
216  Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Aherm2);
217  Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Askew2);
218  delete [] dirRows;
219  delete dirCols;
220  }
221 
222  // Grab diagonal of A
223  A_->ExtractDiagonalCopy(Adiag);
224 
225  // Allocate the diagonal for W
226  Epetra_Vector *Wdiag = new Epetra_Vector(*RowMap);
227 
228  // Build the W matrix (lower triangle only)
229  // Note: Relies on EpetraExt giving me identical sparsity patterns for both Askew2 and Aherm2 (see sanity check above)
230  int maxentries=Askew2->MaxNumEntries();
231  int_type* gids=new int_type [maxentries];
232  double* newvals=new double[maxentries];
233  W=new Epetra_CrsMatrix(Copy,*RowMap,0);
234  for(int i=0;i<N;i++){
235  // Build the - (1+alpha)/2 E - (1-alpha)/2 F part of the W matrix
236  int_type rowgid = (int_type) Acrs_->GRID64(i);
237  double c_data=0.0;
238  double ipdamp=0.0;
239  int idx=0;
240 
241  for(int j=rowptr_s[i];j<rowptr_s[i+1];j++){
242  int_type colgid = (int_type) Askew2->GCID64(colind_s[j]);
243  c_data+=fabs(vals_s[j]);
244  if(rowgid>colgid){
245  // Rely on the fact that off-diagonal entries are always numbered last, dropping the entry entirely.
246  if(colind_s[j] < N) {
247  gids[idx]=colgid;
248  newvals[idx]=vals_h[j]/2 + Alpha_ * vals_s[j]/2;
249  idx++;
250  }
251  else{
252  ipdamp+=fabs(vals_h[j]/2 + Alpha_ * vals_s[j]/2);
253  }
254  }
255  }
256  if(idx>0)
257  IFPACK_CHK_ERR(W->InsertGlobalValues(rowgid,idx,newvals,gids));
258 
259  // Do the diagonal
260  double w_val= c_data*Alpha_*Gamma_/4 + Adiag[Acrs_->LRID(rowgid)];
261  if(UseInterprocDamping_) w_val+=ipdamp;
262 
263  W->InsertGlobalValues(rowgid,1,&w_val,&rowgid);
264  IFPACK_CHK_ERR(Wdiag->ReplaceGlobalValues(1,&w_val,&rowgid));
265  }
266  W->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
267  W_= Teuchos::rcp(W);
268  Wdiag_= Teuchos::rcp(Wdiag);
269 
270  // Mark as computed
271  IsComputed_=true;
272 
273  // Global damping, if wanted
274  if(UseGlobalDamping_ && LambdaMax_ == -1) {
275  if(List_.isParameter("sora: eigen-analysis: random seed")) {
276  unsigned int seed=0;
277  seed = List_.get("sora: eigen-analysis: random seed",seed);
278  PowerMethod(PowerMethodIters_,LambdaMax_,&seed);
279  }
280  else {
281  PowerMethod(PowerMethodIters_,LambdaMax_);
282  }
283 
284  LambdaMax_*=LambdaMaxBoost_;
285  if(!A_->Comm().MyPID()) printf("SORa: Global damping parameter = %6.4e (lmax=%6.4e)\n",GetOmega(),LambdaMax_);
286  }
287 
288  // Cleanup
289  delete [] gids;
290  delete [] newvals;
291  delete Aherm2;
292  delete Askew2;
293  if(RowMatrixMode) {
294  Acrs_=Teuchos::null;
295  }
296 
297  // Counters
298  NumCompute_++;
299  ComputeTime_ += Time_.ElapsedTime();
300  return 0;
301 }
302 
303 int Ifpack_SORa::Compute(){
304  if(!IsInitialized_) Initialize();
305  int ret_val = 0;
306 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
307  if(A_->RowMatrixRowMap().GlobalIndicesInt()) {
308  ret_val = TCompute<int>();
309  }
310  else
311 #endif
312 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
313  if(A_->RowMatrixRowMap().GlobalIndicesLongLong()) {
314  ret_val = TCompute<long long>();
315  }
316  else
317 #endif
318  throw "Ifpack_SORa::Compute: GlobalIndices type unknown";
319 
320  return ret_val;
321 }
322 
323 
324 int Ifpack_SORa::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
325  if(!IsComputed_) return -1;
326  Time_.ResetStartTime();
327  bool initial_guess_is_zero=false;
328  const int lclNumRows = W_->NumMyRows();
329  const int NumVectors = X.NumVectors();
330  Epetra_MultiVector Temp(A_->RowMatrixRowMap(),NumVectors);
331 
332  double omega=GetOmega();
333 
334  // need to create an auxiliary vector, Xcopy
336  if (X.Pointers()[0] == Y.Pointers()[0]){
337  Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
338  // Since the user didn't give us anything better, our initial guess is zero.
339  Y.Scale(0.0);
340  initial_guess_is_zero=true;
341  }
342  else
343  Xcopy = Teuchos::rcp( &X, false );
344 
346  // Note: Assuming that the matrix has an importer. Ifpack_PointRelaxation doesn't do this, but given that
347  // I have a CrsMatrix, I'm probably OK.
348  // Note: This is the lazy man's version sacrificing a few extra flops for avoiding if statements to determine
349  // if things are on or off processor.
350  // Note: T2 must be zero'd out
351  if (IsParallel_ && W_->Importer()) T2 = Teuchos::rcp( new Epetra_MultiVector(W_->Importer()->TargetMap(),NumVectors,true));
352  else T2 = Teuchos::rcp( new Epetra_MultiVector(A_->RowMatrixRowMap(),NumVectors,true));
353 
354  // Pointer grabs
355  int* rowptr,*colind;
356  double *values;
357  double **t_ptr,** y_ptr, ** t2_ptr, **x_ptr,*d_ptr;
358  T2->ExtractView(&t2_ptr);
359  Y.ExtractView(&y_ptr);
360  Temp.ExtractView(&t_ptr);
361  Xcopy->ExtractView(&x_ptr);
362  Wdiag_->ExtractView(&d_ptr);
363  IFPACK_CHK_ERR(W_->ExtractCrsDataPointers(rowptr,colind,values));
364 
365 
366  for(int i=0; i<NumSweeps_; i++){
367  // Calculate b-Ax
368  if(!initial_guess_is_zero || i > 0) {
369  A_->Apply(Y,Temp);
370  Temp.Update(1.0,*Xcopy,-1.0);
371  }
372  else
373  Temp.Update(1.0,*Xcopy,0.0);
374 
375  // Note: The off-processor entries of T2 never get touched (they're always zero) and the other entries are updated
376  // in this sweep before they are used, so we don't need to reset T2 to zero here.
377 
378  // Do backsolve & update
379  // x = x + W^{-1} (b - A x)
380  for(int j=0; j<lclNumRows; j++){
381  double diag=d_ptr[j];
382  for (int m=0 ; m<NumVectors; m++) {
383  double dtmp=0.0;
384  // Note: Since the diagonal is in the matrix, we need to zero that entry of T2 here to make sure it doesn't contribute.
385  t2_ptr[m][j]=0.0;
386  for(int k=rowptr[j];k<rowptr[j+1];k++){
387  dtmp+= values[k]*t2_ptr[m][colind[k]];
388  }
389  // Yes, we need to update both of these.
390  t2_ptr[m][j] = (t_ptr[m][j]- dtmp)/diag;
391  y_ptr[m][j] += omega*t2_ptr[m][j];
392  }
393  }
394  }
395 
396  // Counter update
397  NumApplyInverse_++;
398  ApplyInverseTime_ += Time_.ElapsedTime();
399  return 0;
400 }
401 
402 
403 std::ostream& Ifpack_SORa::Print(std::ostream& os) const{
404  using std::endl;
405 
406  os<<"Ifpack_SORa"<<endl;
407  os<<" alpha = "<<Alpha_<<endl;
408  os<<" gamma = "<<Gamma_<<endl;
409  os<<endl;
410  return os;
411 }
412 
413 
414 double Ifpack_SORa::Condest(const Ifpack_CondestType /* CT */,
415  const int /* MaxIters */,
416  const double /* Tol */,
417  Epetra_RowMatrix* /* Matrix_in */){
418  return -1.0;
419 }
420 
421 
422 
423 
424 
425 // ============================================================================
426 inline int* FindLocalDiricheltRowsFromOnesAndZeros(const Epetra_CrsMatrix & Matrix, int &numBCRows){
427  int *dirichletRows = new int[Matrix.NumMyRows()];
428  numBCRows = 0;
429  for (int i=0; i<Matrix.NumMyRows(); i++) {
430  int numEntries, *cols;
431  double *vals;
432  int ierr = Matrix.ExtractMyRowView(i,numEntries,vals,cols);
433  if (ierr == 0) {
434  int nz=0;
435  for (int j=0; j<numEntries; j++) if (vals[j]!=0) nz++;
436  if (nz==1) dirichletRows[numBCRows++] = i;
437  // EXPERIMENTAL: Treat Special Inflow Boundaries as Dirichlet Boundaries
438  if(nz==2) dirichletRows[numBCRows++] = i;
439  }/*end if*/
440  }/*end for*/
441  return dirichletRows;
442 }/*end FindLocalDiricheltRowsFromOnesAndZeros*/
443 
444 
445 // ======================================================================
447 inline Epetra_IntVector * FindLocalDirichletColumnsFromRows(const int *dirichletRows, int numBCRows,const Epetra_CrsMatrix & Matrix){
448  // Zero the columns corresponding to the Dirichlet rows. Completely ignore the matrix maps.
449 
450  // Build rows
451  Epetra_IntVector ZeroRows(Matrix.RowMap());
452  Epetra_IntVector *ZeroCols=new Epetra_IntVector(Matrix.ColMap());
453  ZeroRows.PutValue(0);
454  ZeroCols->PutValue(0);
455 
456  // Easy Case: We're all on one processor
457  if(Matrix.RowMap().SameAs(Matrix.ColMap())){
458  for (int i=0; i < numBCRows; i++)
459  (*ZeroCols)[dirichletRows[i]]=1;
460  return ZeroCols;
461  }
462 
463  // Flag the rows which are zero locally
464  for (int i=0; i < numBCRows; i++)
465  ZeroRows[dirichletRows[i]]=1;
466 
467  // Boundary exchange to move the data
468  if(Matrix.RowMap().SameAs(Matrix.DomainMap())){
469  // Use A's Importer if we have one
470  ZeroCols->Import(ZeroRows,*Matrix.Importer(),Insert);
471  }
472  else{
473  // Use custom importer if we don't
474  Epetra_Import Importer(Matrix.ColMap(),Matrix.RowMap());
475  ZeroCols->Import(ZeroRows,Importer,Insert);
476  }
477  return ZeroCols;
478 }
479 
480 
481 // ======================================================================
482 inline void Apply_BCsToMatrixRowsAndColumns(const int *dirichletRows, int numBCRows,const Epetra_IntVector &dirichletColumns,const Epetra_CrsMatrix & Matrix){
483  /* This function zeros out rows & columns of Matrix.
484  Comments: The graph of Matrix is unchanged.
485  */
486  // Nuke the rows
487  for(int i=0;i<numBCRows;i++){
488  int numEntries, *cols;
489  double *vals;
490  Matrix.ExtractMyRowView(dirichletRows[i],numEntries,vals,cols);
491  for (int j=0; j<numEntries; j++) vals[j]=0.0;
492  }/*end for*/
493 
494  // Nuke the columns
495  for (int i=0; i < Matrix.NumMyRows(); i++) {
496  int numEntries;
497  double *vals;
498  int *cols;
499  Matrix.ExtractMyRowView(i,numEntries,vals,cols);
500  for (int j=0; j < numEntries; j++) {
501  if (dirichletColumns[ cols[j] ] > 0) vals[j] = 0.0;
502  }/*end for*/
503  }/*end for*/
504 }/* end Apply_BCsToMatrixColumns */
505 
506 
507 
508 
509 
510 int Ifpack_SORa::
511 PowerMethod(const int MaximumIterations, double& lambda_max, const unsigned int * RngSeed)
512 {
513  // this is a simple power method
514  lambda_max = 0.0;
515  double RQ_top, RQ_bottom, norm;
516  Epetra_Vector x(A_->OperatorDomainMap());
517  Epetra_Vector y(A_->OperatorRangeMap());
518  Epetra_Vector z(A_->OperatorRangeMap());
519 
520  if(RngSeed) x.SetSeed(*RngSeed);
521  x.Random();
522  x.Norm2(&norm);
523  if (norm == 0.0) IFPACK_CHK_ERR(-1);
524 
525  x.Scale(1.0 / norm);
526 
527  // Only do 1 sweep per PM step, turn off global damping
528  int NumSweepsBackup=NumSweeps_;
529  bool UseGlobalDampingBackup=UseGlobalDamping_;
530  NumSweeps_=1;UseGlobalDamping_=false;
531 
532  for (int iter = 0; iter < MaximumIterations; ++iter)
533  {
534  y.PutScalar(0.0);
535  A_->Apply(x, z);
536  ApplyInverse(z,y);
537  y.Update(1.0,x,-1.0);
538 
539  // Compute the Rayleigh Quotient
540  y.Dot(x, &RQ_top);
541  x.Dot(x, &RQ_bottom);
542  lambda_max = RQ_top / RQ_bottom;
543  y.Norm2(&norm);
544  if (norm == 0.0) IFPACK_CHK_ERR(-1);
545  x.Update(1.0 / norm, y, 0.0);
546 
547  }
548 
549  // Enable if we want to prevent over-relaxation
550  // lambda_max=MAX(lambda_max,1.0);
551 
552  // Reset parameters
553  NumSweeps_=NumSweepsBackup;
554  UseGlobalDamping_=UseGlobalDampingBackup;
555 
556  return(0);
557 }
558 
559 
560 #endif
bool GlobalIndicesLongLong() const
bool SameAs(const Epetra_BlockMap &Map) const
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
const int NumVectors
Definition: performance.cpp:71
int PutValue(int Value)
bool GlobalIndicesInt() const
const Epetra_Map & ColMap() const
int FillComplete(bool OptimizeDataStorage=true)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int ReplaceGlobalValues(int NumEntries, const double *Values, const int *Indices)
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
const Epetra_Map & RowMap() const
long long GID64(int LID) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
const Epetra_Import * Importer() const
int MaxNumEntries() const
int GID(int LID) const
int NumMyRows() const
virtual int NumMyRows() const =0
const int N
#define false
static int Add(const Epetra_CrsMatrix &A, bool transposeA, double scalarA, Epetra_CrsMatrix &B, double scalarB)
TransListIter iter
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the matrix to be preconditioned.
const Epetra_Map & DomainMap() const
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
#define IFPACK_CHK_ERR(ifpack_err)
long long GCID64(int LCID_in) const
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
int NumMyNonzeros() const