IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_PointRelaxation.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 <cmath>
46 #include "Epetra_Operator.h"
47 #include "Epetra_CrsMatrix.h"
48 #include "Epetra_VbrMatrix.h"
49 #include "Epetra_Comm.h"
50 #include "Epetra_Map.h"
51 #include "Epetra_MultiVector.h"
52 #include "Ifpack_Preconditioner.h"
53 #include "Ifpack_PointRelaxation.h"
54 #include "Ifpack_Utils.h"
55 #include "Ifpack_Condest.h"
56 
57 static const int IFPACK_JACOBI = 0;
58 static const int IFPACK_GS = 1;
59 static const int IFPACK_SGS = 2;
60 
61 //==============================================================================
64  IsInitialized_(false),
65  IsComputed_(false),
66  NumInitialize_(0),
67  NumCompute_(0),
68  NumApplyInverse_(0),
69  InitializeTime_(0.0),
70  ComputeTime_(0.0),
71  ApplyInverseTime_(0.0),
72  ComputeFlops_(0.0),
73  ApplyInverseFlops_(0.0),
74  NumSweeps_(1),
75  DampingFactor_(1.0),
76  UseTranspose_(false),
77  Condest_(-1.0),
78  /* ComputeCondest_(false), (unused; commented out to avoid build warnings) */
79  PrecType_(IFPACK_JACOBI),
80  MinDiagonalValue_(0.0),
81  NumMyRows_(0),
82  NumMyNonzeros_(0),
83  NumGlobalRows_(0),
84  NumGlobalNonzeros_(0),
85  Matrix_(Teuchos::rcp(Matrix_in,false)),
86  IsParallel_(false),
87  ZeroStartingSolution_(true),
88  DoBackwardGS_(false),
89  DoL1Method_(false),
90  L1Eta_(1.5),
91  NumLocalSmoothingIndices_(Matrix_in->NumMyRows()),
92  LocalSmoothingIndices_(0)
93 {
94 }
95 
96 //==============================================================================
97 int Ifpack_PointRelaxation::SetParameters(Teuchos::ParameterList& List)
98 {
99  using std::cout;
100  using std::endl;
101 
102  std::string PT;
103  if (PrecType_ == IFPACK_JACOBI)
104  PT = "Jacobi";
105  else if (PrecType_ == IFPACK_GS)
106  PT = "Gauss-Seidel";
107  else if (PrecType_ == IFPACK_SGS)
108  PT = "symmetric Gauss-Seidel";
109 
110  PT = List.get("relaxation: type", PT);
111 
112  if (PT == "Jacobi")
113  PrecType_ = IFPACK_JACOBI;
114  else if (PT == "Gauss-Seidel")
115  PrecType_ = IFPACK_GS;
116  else if (PT == "symmetric Gauss-Seidel")
117  PrecType_ = IFPACK_SGS;
118  else {
119  IFPACK_CHK_ERR(-2);
120  }
121 
122  NumSweeps_ = List.get("relaxation: sweeps",NumSweeps_);
123  DampingFactor_ = List.get("relaxation: damping factor",
124  DampingFactor_);
125  MinDiagonalValue_ = List.get("relaxation: min diagonal value",
126  MinDiagonalValue_);
127  ZeroStartingSolution_ = List.get("relaxation: zero starting solution",
128  ZeroStartingSolution_);
129 
130  DoBackwardGS_ = List.get("relaxation: backward mode",DoBackwardGS_);
131 
132  DoL1Method_ = List.get("relaxation: use l1",DoL1Method_);
133 
134  L1Eta_ = List.get("relaxation: l1 eta",L1Eta_);
135 
136 
137  // This (partial) reordering would require a very different implementation of Jacobi than we have no, so we're not
138  // going to do it.
139  NumLocalSmoothingIndices_= List.get("relaxation: number of local smoothing indices",NumLocalSmoothingIndices_);
140  LocalSmoothingIndices_ = List.get("relaxation: local smoothing indices",LocalSmoothingIndices_);
141  if(PrecType_ == IFPACK_JACOBI && LocalSmoothingIndices_) {
142  NumLocalSmoothingIndices_=NumMyRows_;
143  LocalSmoothingIndices_=0;
144  if(!Comm().MyPID()) cout<<"Ifpack_PointRelaxation: WARNING: Reordered/Local Smoothing not implemented for Jacobi. Defaulting to regular Jacobi"<<endl;
145  }
146 
147  SetLabel();
148 
149  return(0);
150 }
151 
152 //==============================================================================
154 {
155  return(Matrix_->Comm());
156 }
157 
158 //==============================================================================
160 {
161  return(Matrix_->OperatorDomainMap());
162 }
163 
164 //==============================================================================
166 {
167  return(Matrix_->OperatorRangeMap());
168 }
169 
170 //==============================================================================
172 {
173  IsInitialized_ = false;
174 
175  if (Matrix_ == Teuchos::null)
176  IFPACK_CHK_ERR(-2);
177 
178  if (Time_ == Teuchos::null)
179  Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
180 
181  if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
182  IFPACK_CHK_ERR(-2); // only square matrices
183 
184  NumMyRows_ = Matrix_->NumMyRows();
185  NumMyNonzeros_ = Matrix_->NumMyNonzeros();
186  NumGlobalRows_ = Matrix_->NumGlobalRows64();
187  NumGlobalNonzeros_ = Matrix_->NumGlobalNonzeros64();
188 
189  if (Comm().NumProc() != 1)
190  IsParallel_ = true;
191  else
192  IsParallel_ = false;
193 
194  ++NumInitialize_;
195  InitializeTime_ += Time_->ElapsedTime();
196  IsInitialized_ = true;
197  return(0);
198 }
199 
200 //==============================================================================
202 {
203  int ierr = 0;
204  if (!IsInitialized())
205  IFPACK_CHK_ERR(Initialize());
206 
207  Time_->ResetStartTime();
208 
209  // reset values
210  IsComputed_ = false;
211  Condest_ = -1.0;
212 
213  if (NumSweeps_ == 0) ierr = 1; // Warning: no sweeps performed.
214  if (NumSweeps_ < 0)
215  IFPACK_CHK_ERR(-2); // at least one application
216 
217  // NOTE: RowMatrixRowMap doesn't work correctly for Epetra_VbrMatrix
218 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
219  const Epetra_VbrMatrix * VbrMat = dynamic_cast<const Epetra_VbrMatrix*>(&Matrix());
220  if(VbrMat) Diagonal_ = Teuchos::rcp( new Epetra_Vector(VbrMat->RowMap()) );
221  else
222 #endif
223  Diagonal_ = Teuchos::rcp( new Epetra_Vector(Matrix().RowMatrixRowMap()) );
224 
225  if (Diagonal_ == Teuchos::null)
226  IFPACK_CHK_ERR(-5);
227 
228  IFPACK_CHK_ERR(Matrix().ExtractDiagonalCopy(*Diagonal_));
229 
230  // Setup for L1 Methods.
231  // Here we add half the value of the off-processor entries in the row,
232  // but only if diagonal isn't sufficiently large.
233  //
234  // Note: This is only done in the slower-but-more-general "RowMatrix" mode.
235  //
236  // This follows from Equation (6.5) in:
237  // Baker, Falgout, Kolev and Yang. Multigrid Smoothers for Ultraparallel Computing.
238  // SIAM J. Sci. Comput., Vol. 33, No. 5. (2011), pp. 2864--2887.
239  if(DoL1Method_ && IsParallel_) {
240  int maxLength = Matrix().MaxNumEntries();
241  std::vector<int> Indices(maxLength);
242  std::vector<double> Values(maxLength);
243  int NumEntries;
244 
245  for (int i = 0 ; i < NumMyRows_ ; ++i) {
246  IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, maxLength,NumEntries,
247  &Values[0], &Indices[0]));
248  double diagonal_boost=0.0;
249  for (int k = 0 ; k < NumEntries ; ++k)
250  if(Indices[k] >= NumMyRows_)
251  diagonal_boost+=std::abs(Values[k]/2.0);
252  if ((*Diagonal_)[i] < L1Eta_*diagonal_boost)
253  (*Diagonal_)[i]+=diagonal_boost;
254  }
255  }
256 
257  // check diagonal elements, store the inverses, and verify that
258  // no zeros are around. If an element is zero, then by default
259  // its inverse is zero as well (that is, the row is ignored).
260  for (int i = 0 ; i < NumMyRows_ ; ++i) {
261  double& diag = (*Diagonal_)[i];
262  if (IFPACK_ABS(diag) < MinDiagonalValue_)
263  diag = MinDiagonalValue_;
264  if (diag != 0.0)
265  diag = 1.0 / diag;
266  }
267 #ifdef IFPACK_FLOPCOUNTERS
268  ComputeFlops_ += NumMyRows_;
269 #endif
270 
271 #if 0
272  // some methods require the inverse of the diagonal, compute it
273  // now and store it in `Diagonal_'
274  if ((PrecType_ == IFPACK_JACOBI) || (PrecType_ == IFPACK_GS)) {
275  Diagonal_->Reciprocal(*Diagonal_);
276  // update flops
277  ComputeFlops_ += NumMyRows_;
278  }
279 #endif
280 
281  // We need to import data from external processors. Here I create an
282  // Epetra_Import object because I cannot assume that Matrix_ has one.
283  // This is a bit of waste of resources (but the code is more robust).
284  // Note that I am doing some strange stuff to set the components of Y
285  // from Y2 (to save some time).
286  //
287  if (IsParallel_ && ((PrecType_ == IFPACK_GS) || (PrecType_ == IFPACK_SGS))) {
288  Importer_ = Teuchos::rcp( new Epetra_Import(Matrix().RowMatrixColMap(),
289  Matrix().RowMatrixRowMap()) );
290  if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
291  }
292 
293  ++NumCompute_;
294  ComputeTime_ += Time_->ElapsedTime();
295  IsComputed_ = true;
296 
297  IFPACK_CHK_ERR(ierr);
298  return(0);
299 }
300 
301 //==============================================================================
302 std::ostream& Ifpack_PointRelaxation::Print(std::ostream & os) const
303 {
304  using std::endl;
305 
306  double MyMinVal, MyMaxVal;
307  double MinVal, MaxVal;
308 
309  if (IsComputed_) {
310  Diagonal_->MinValue(&MyMinVal);
311  Diagonal_->MaxValue(&MyMaxVal);
312  Comm().MinAll(&MyMinVal,&MinVal,1);
313  Comm().MinAll(&MyMaxVal,&MaxVal,1);
314  }
315 
316  if (!Comm().MyPID()) {
317  os << endl;
318  os << "================================================================================" << endl;
319  os << "Ifpack_PointRelaxation" << endl;
320  os << "Sweeps = " << NumSweeps_ << endl;
321  os << "damping factor = " << DampingFactor_ << endl;
322  if (PrecType_ == IFPACK_JACOBI)
323  os << "Type = Jacobi" << endl;
324  else if (PrecType_ == IFPACK_GS)
325  os << "Type = Gauss-Seidel" << endl;
326  else if (PrecType_ == IFPACK_SGS)
327  os << "Type = symmetric Gauss-Seidel" << endl;
328  if (DoBackwardGS_)
329  os << "Using backward mode (GS only)" << endl;
330  if (ZeroStartingSolution_)
331  os << "Using zero starting solution" << endl;
332  else
333  os << "Using input starting solution" << endl;
334  os << "Condition number estimate = " << Condest() << endl;
335  os << "Global number of rows = " << Matrix_->NumGlobalRows64() << endl;
336  if (IsComputed_) {
337  os << "Minimum value on stored diagonal = " << MinVal << endl;
338  os << "Maximum value on stored diagonal = " << MaxVal << endl;
339  }
340  os << endl;
341  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
342  os << "----- ------- -------------- ------------ --------" << endl;
343  os << "Initialize() " << std::setw(5) << NumInitialize_
344  << " " << std::setw(15) << InitializeTime_
345  << " 0.0 0.0" << endl;
346  os << "Compute() " << std::setw(5) << NumCompute_
347  << " " << std::setw(15) << ComputeTime_
348  << " " << std::setw(15) << 1.0e-6 * ComputeFlops_;
349  if (ComputeTime_ != 0.0)
350  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
351  else
352  os << " " << std::setw(15) << 0.0 << endl;
353  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse_
354  << " " << std::setw(15) << ApplyInverseTime_
355  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
356  if (ApplyInverseTime_ != 0.0)
357  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
358  else
359  os << " " << std::setw(15) << 0.0 << endl;
360  os << "================================================================================" << endl;
361  os << endl;
362  }
363 
364  return(os);
365 }
366 
367 //==============================================================================
369 Condest(const Ifpack_CondestType CT,
370  const int MaxIters, const double Tol,
371  Epetra_RowMatrix* Matrix_in)
372 {
373  if (!IsComputed()) // cannot compute right now
374  return(-1.0);
375 
376  // always computes it. Call Condest() with no parameters to get
377  // the previous estimate.
378  Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
379 
380  return(Condest_);
381 }
382 
383 //==============================================================================
384 void Ifpack_PointRelaxation::SetLabel()
385 {
386  std::string PT;
387  if (PrecType_ == IFPACK_JACOBI)
388  PT = "Jacobi";
389  else if (PrecType_ == IFPACK_GS){
390  PT = "GS";
391  if(DoBackwardGS_)
392  PT = "Backward " + PT;
393  }
394  else if (PrecType_ == IFPACK_SGS)
395  PT = "SGS";
396 
397  if(NumLocalSmoothingIndices_!=NumMyRows_) PT="Local " + PT;
398  else if(LocalSmoothingIndices_) PT="Reordered " + PT;
399 
400  Label_ = "IFPACK (" + PT + ", sweeps=" + Ifpack_toString(NumSweeps_)
401  + ", damping=" + Ifpack_toString(DampingFactor_) + ")";
402 }
403 
404 //==============================================================================
405 // Note that Ifpack_PointRelaxation and Jacobi is much faster than
406 // Ifpack_AdditiveSchwarz<Ifpack_PointRelaxation> (because of the
407 // way the matrix-vector produce is performed).
408 //
409 // Another ML-related observation is that the starting solution (in Y)
410 // is NOT supposed to be zero. This may slow down the application of just
411 // one sweep of Jacobi.
412 //
415 {
416  if (!IsComputed())
417  IFPACK_CHK_ERR(-3);
418 
419  if (X.NumVectors() != Y.NumVectors())
420  IFPACK_CHK_ERR(-2);
421 
422  Time_->ResetStartTime();
423 
424  // AztecOO gives X and Y pointing to the same memory location,
425  // need to create an auxiliary vector, Xcopy
426  Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
427  if (X.Pointers()[0] == Y.Pointers()[0])
428  Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
429  else
430  Xcopy = Teuchos::rcp( &X, false );
431 
432  // Flops are updated in each of the following.
433  switch (PrecType_) {
434  case IFPACK_JACOBI:
435  IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
436  break;
437  case IFPACK_GS:
438  IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
439  break;
440  case IFPACK_SGS:
441  IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
442  break;
443  default:
444  IFPACK_CHK_ERR(-1); // something wrong
445  }
446 
447  ++NumApplyInverse_;
448  ApplyInverseTime_ += Time_->ElapsedTime();
449  return(0);
450 }
451 
452 //==============================================================================
453 // This preconditioner can be much slower than AztecOO and ML versions
454 // if the matrix-vector product requires all ExtractMyRowCopy()
455 // (as done through Ifpack_AdditiveSchwarz).
456 int Ifpack_PointRelaxation::
457 ApplyInverseJacobi(const Epetra_MultiVector& RHS, Epetra_MultiVector& LHS) const
458 {
459  int NumVectors = LHS.NumVectors();
460 
461  int startIter = 0;
462  if (NumSweeps_ > 0 && ZeroStartingSolution_) {
463  // When we have a zero initial guess, we can skip the first
464  // matrix apply call and zero initialization
465  for (int v = 0; v < NumVectors; v++)
466  IFPACK_CHK_ERR(LHS(v)->Multiply(DampingFactor_, *(RHS(v)), *Diagonal_, 0.0));
467 
468  startIter = 1;
469  }
470 
471  bool zeroOut = false;
472  Epetra_MultiVector A_times_LHS(LHS.Map(), NumVectors, zeroOut);
473  for (int j = startIter; j < NumSweeps_ ; j++) {
474  IFPACK_CHK_ERR(Apply(LHS, A_times_LHS));
475  IFPACK_CHK_ERR(A_times_LHS.Update(1.0,RHS,-1.0));
476  for (int v = 0 ; v < NumVectors ; ++v)
477  IFPACK_CHK_ERR(LHS(v)->Multiply(DampingFactor_, *(A_times_LHS(v)),
478  *Diagonal_, 1.0));
479 
480  }
481 
482  // Flops:
483  // - matrix vector (2 * NumGlobalNonzeros_)
484  // - update (2 * NumGlobalRows_)
485  // - Multiply:
486  // - DampingFactor (NumGlobalRows_)
487  // - Diagonal (NumGlobalRows_)
488  // - A + B (NumGlobalRows_)
489  // - 1.0 (NumGlobalRows_)
490 #ifdef IFPACK_FLOPCOUNTERS
491  ApplyInverseFlops_ += NumVectors * (6 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
492 #endif
493 
494  return(0);
495 }
496 
497 //==============================================================================
498 int Ifpack_PointRelaxation::
499 ApplyInverseGS(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
500 {
501  if (ZeroStartingSolution_)
502  Y.PutScalar(0.0);
503 
504  const Epetra_CrsMatrix* CrsMatrix = dynamic_cast<const Epetra_CrsMatrix*>(&*Matrix_);
505  // try to pick the best option; performances may be improved
506  // if several sweeps are used.
507  if (CrsMatrix != 0)
508  {
509  if (CrsMatrix->StorageOptimized() && LocalSmoothingIndices_)
510  return(ApplyInverseGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
511  else if (CrsMatrix->StorageOptimized())
512  return(ApplyInverseGS_FastCrsMatrix(CrsMatrix, X, Y));
513  else
514  return(ApplyInverseGS_CrsMatrix(CrsMatrix, X, Y));
515  }
516  else
517  return(ApplyInverseGS_RowMatrix(X, Y));
518 } //ApplyInverseGS()
519 
520 
521 
522 //==============================================================================
523 int Ifpack_PointRelaxation::
524 ApplyInverseGS_RowMatrix(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
525 {
526  int NumVectors = X.NumVectors();
527 
528  int Length = Matrix().MaxNumEntries();
529  std::vector<int> Indices(Length);
530  std::vector<double> Values(Length);
531 
532  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
533  if (IsParallel_)
534  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
535  else
536  Y2 = Teuchos::rcp( &Y, false );
537 
538  // extract views (for nicer and faster code)
539  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
540  X.ExtractView(&x_ptr);
541  Y.ExtractView(&y_ptr);
542  Y2->ExtractView(&y2_ptr);
543  Diagonal_->ExtractView(&d_ptr);
544 
545  for (int j = 0; j < NumSweeps_ ; j++) {
546 
547  // data exchange is here, once per sweep
548  if (IsParallel_)
549  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
550 
551  // FIXME: do I really need this code below?
552  if (NumVectors == 1) {
553 
554  double* y0_ptr = y_ptr[0];
555  double* y20_ptr = y2_ptr[0];
556  double* x0_ptr = x_ptr[0];
557 
558  if(!DoBackwardGS_){
559  /* Forward Mode */
560  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
561  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
562 
563  int NumEntries;
564  int col;
565  IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
566  &Values[0], &Indices[0]));
567 
568  double dtemp = 0.0;
569  for (int k = 0 ; k < NumEntries ; ++k) {
570 
571  col = Indices[k];
572  dtemp += Values[k] * y20_ptr[col];
573  }
574 
575  y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
576  }
577  }
578  else {
579  /* Backward Mode */
580  for (int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
581  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
582 
583  int NumEntries;
584  // int col; // unused
585  // (void) col; // Forestall compiler warning for unused variable.
586  IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
587  &Values[0], &Indices[0]));
588  double dtemp = 0.0;
589  for (int k = 0 ; k < NumEntries ; ++k) {
590  //col = Indices[k]; // unused
591  dtemp += Values[k] * y20_ptr[i];
592  }
593 
594  y20_ptr[i] += DampingFactor_ * d_ptr[i] * (x0_ptr[i] - dtemp);
595  }
596  }
597 
598  // using Export() sounded quite expensive
599  if (IsParallel_)
600  for (int i = 0 ; i < NumMyRows_ ; ++i)
601  y0_ptr[i] = y20_ptr[i];
602 
603  }
604  else {
605  if(!DoBackwardGS_){
606  /* Forward Mode */
607  for (int i = 0 ; i < NumMyRows_ ; ++i) {
608 
609  int NumEntries;
610  int col;
611  IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
612  &Values[0], &Indices[0]));
613 
614  for (int m = 0 ; m < NumVectors ; ++m) {
615 
616  double dtemp = 0.0;
617  for (int k = 0 ; k < NumEntries ; ++k) {
618 
619  col = Indices[k];
620  dtemp += Values[k] * y2_ptr[m][col];
621  }
622 
623  y2_ptr[m][i] += DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
624  }
625  }
626  }
627  else {
628  /* Backward Mode */
629  for (int i = NumMyRows_ - 1 ; i > -1 ; --i) {
630  int NumEntries;
631  int col;
632  IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
633  &Values[0], &Indices[0]));
634 
635  for (int m = 0 ; m < NumVectors ; ++m) {
636 
637  double dtemp = 0.0;
638  for (int k = 0 ; k < NumEntries ; ++k) {
639 
640  col = Indices[k];
641  dtemp += Values[k] * y2_ptr[m][col];
642  }
643 
644  y2_ptr[m][i] += DampingFactor_ * d_ptr[i] * (x_ptr[m][i] - dtemp);
645 
646  }
647  }
648  }
649 
650  // using Export() sounded quite expensive
651  if (IsParallel_)
652  for (int m = 0 ; m < NumVectors ; ++m)
653  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
654  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
655  y_ptr[m][i] = y2_ptr[m][i];
656  }
657  }
658  }
659 
660 #ifdef IFPACK_FLOPCOUNTERS
661  ApplyInverseFlops_ += NumVectors * (4 * NumGlobalRows_ + 2 * NumGlobalNonzeros_);
662 #endif
663 
664  return(0);
665 } //ApplyInverseGS_RowMatrix()
666 
667 //==============================================================================
668 int Ifpack_PointRelaxation::
669 ApplyInverseGS_CrsMatrix(const Epetra_CrsMatrix* A, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
670 {
671  int NumVectors = X.NumVectors();
672 
673  int* Indices;
674  double* Values;
675 
676  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
677  if (IsParallel_) {
678  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
679  }
680  else
681  Y2 = Teuchos::rcp( &Y, false );
682 
683  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
684  X.ExtractView(&x_ptr);
685  Y.ExtractView(&y_ptr);
686  Y2->ExtractView(&y2_ptr);
687  Diagonal_->ExtractView(&d_ptr);
688 
689  for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
690 
691  // only one data exchange per sweep
692  if (IsParallel_)
693  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
694 
695  if(!DoBackwardGS_){
696  /* Forward Mode */
697  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
698  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
699 
700  int NumEntries;
701  int col;
702  double diag = d_ptr[i];
703 
704  IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
705 
706  for (int m = 0 ; m < NumVectors ; ++m) {
707 
708  double dtemp = 0.0;
709 
710  for (int k = 0; k < NumEntries; ++k) {
711 
712  col = Indices[k];
713  dtemp += Values[k] * y2_ptr[m][col];
714  }
715 
716  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
717  }
718  }
719  }
720  else {
721  /* Backward Mode */
722  for (int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
723  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
724 
725  int NumEntries;
726  int col;
727  double diag = d_ptr[i];
728 
729  IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
730 
731  for (int m = 0 ; m < NumVectors ; ++m) {
732 
733  double dtemp = 0.0;
734  for (int k = 0; k < NumEntries; ++k) {
735 
736  col = Indices[k];
737  dtemp += Values[k] * y2_ptr[m][col];
738  }
739 
740  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
741 
742  }
743  }
744  }
745 
746  if (IsParallel_)
747  for (int m = 0 ; m < NumVectors ; ++m)
748  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
749  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
750  y_ptr[m][i] = y2_ptr[m][i];
751  }
752  }
753 
754 #ifdef IFPACK_FLOPCOUNTERS
755  ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
756 #endif
757  return(0);
758 } //ApplyInverseGS_CrsMatrix()
759 
760 //
761 //==============================================================================
762 // ApplyInverseGS_FastCrsMatrix() requires Epetra_CrsMatrix + OptimizeStorage()
763 
764 int Ifpack_PointRelaxation::
765 ApplyInverseGS_FastCrsMatrix(const Epetra_CrsMatrix* A, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
766 {
767  int* IndexOffset;
768  int* Indices;
769  double* Values;
770  IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
771 
772  int NumVectors = X.NumVectors();
773 
774  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
775  if (IsParallel_) {
776  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
777  }
778  else
779  Y2 = Teuchos::rcp( &Y, false );
780 
781  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
782  X.ExtractView(&x_ptr);
783  Y.ExtractView(&y_ptr);
784  Y2->ExtractView(&y2_ptr);
785  Diagonal_->ExtractView(&d_ptr);
786 
787  for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
788 
789  // only one data exchange per sweep
790  if (IsParallel_)
791  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
792 
793  if(!DoBackwardGS_){
794  /* Forward Mode */
795  for (int i = 0 ; i < NumMyRows_ ; ++i) {
796 
797  int col;
798  double diag = d_ptr[i];
799 
800  for (int m = 0 ; m < NumVectors ; ++m) {
801 
802  double dtemp = 0.0;
803 
804  for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
805 
806  col = Indices[k];
807  dtemp += Values[k] * y2_ptr[m][col];
808  }
809 
810  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
811  }
812  }
813  }
814  else {
815  /* Backward Mode */
816  for (int i = NumMyRows_ - 1 ; i > -1 ; --i) {
817 
818  int col;
819  double diag = d_ptr[i];
820 
821  for (int m = 0 ; m < NumVectors ; ++m) {
822 
823  double dtemp = 0.0;
824  for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
825 
826  col = Indices[k];
827  dtemp += Values[k] * y2_ptr[m][col];
828  }
829 
830  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
831 
832  }
833  }
834  }
835 
836 
837  if (IsParallel_)
838  for (int m = 0 ; m < NumVectors ; ++m)
839  for (int i = 0 ; i < NumMyRows_ ; ++i)
840  y_ptr[m][i] = y2_ptr[m][i];
841  }
842 
843 #ifdef IFPACK_FLOPCOUNTERS
844  ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
845 #endif
846  return(0);
847 } //ApplyInverseGS_FastCrsMatrix()
848 
849 
850 //
851 //==============================================================================
852 // ApplyInverseGS_LocalFastCrsMatrix() requires Epetra_CrsMatrix + OptimizeStorage() + LocalSmoothingIndices
853 
854 int Ifpack_PointRelaxation::
855 ApplyInverseGS_LocalFastCrsMatrix(const Epetra_CrsMatrix* A, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
856 {
857  int* IndexOffset;
858  int* Indices;
859  double* Values;
860  IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
861 
862  int NumVectors = X.NumVectors();
863 
864  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
865  if (IsParallel_) {
866  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
867  }
868  else
869  Y2 = Teuchos::rcp( &Y, false );
870 
871  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
872  X.ExtractView(&x_ptr);
873  Y.ExtractView(&y_ptr);
874  Y2->ExtractView(&y2_ptr);
875  Diagonal_->ExtractView(&d_ptr);
876 
877  for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
878 
879  // only one data exchange per sweep
880  if (IsParallel_)
881  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
882 
883  if(!DoBackwardGS_){
884  /* Forward Mode */
885  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
886  int i=LocalSmoothingIndices_[ii];
887 
888  int col;
889  double diag = d_ptr[i];
890 
891  for (int m = 0 ; m < NumVectors ; ++m) {
892 
893  double dtemp = 0.0;
894 
895  for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
896 
897  col = Indices[k];
898  dtemp += Values[k] * y2_ptr[m][col];
899  }
900 
901  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
902  }
903  }
904  }
905  else {
906  /* Backward Mode */
907  for (int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
908  int i=LocalSmoothingIndices_[ii];
909 
910  int col;
911  double diag = d_ptr[i];
912 
913  for (int m = 0 ; m < NumVectors ; ++m) {
914 
915  double dtemp = 0.0;
916  for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
917 
918  col = Indices[k];
919  dtemp += Values[k] * y2_ptr[m][col];
920  }
921 
922  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
923 
924  }
925  }
926  }
927 
928 
929  if (IsParallel_)
930  for (int m = 0 ; m < NumVectors ; ++m)
931  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
932  int i=LocalSmoothingIndices_[ii];
933  y_ptr[m][i] = y2_ptr[m][i];
934  }
935  }
936 
937 #ifdef IFPACK_FLOPCOUNTERS
938  ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
939 #endif
940  return(0);
941 } //ApplyInverseGS_LocalFastCrsMatrix()
942 
943 
944 //==============================================================================
945 int Ifpack_PointRelaxation::
946 ApplyInverseSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
947 {
948  if (ZeroStartingSolution_)
949  Y.PutScalar(0.0);
950 
951  const Epetra_CrsMatrix* CrsMatrix = dynamic_cast<const Epetra_CrsMatrix*>(&*Matrix_);
952  // try to pick the best option; performances may be improved
953  // if several sweeps are used.
954  if (CrsMatrix != 0)
955  {
956  if (CrsMatrix->StorageOptimized() && LocalSmoothingIndices_)
957  return(ApplyInverseSGS_LocalFastCrsMatrix(CrsMatrix, X, Y));
958  else if (CrsMatrix->StorageOptimized())
959  return(ApplyInverseSGS_FastCrsMatrix(CrsMatrix, X, Y));
960  else
961  return(ApplyInverseSGS_CrsMatrix(CrsMatrix, X, Y));
962  }
963  else
964  return(ApplyInverseSGS_RowMatrix(X, Y));
965 }
966 
967 //==============================================================================
968 int Ifpack_PointRelaxation::
969 ApplyInverseSGS_RowMatrix(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
970 {
971  int NumVectors = X.NumVectors();
972  int Length = Matrix().MaxNumEntries();
973  std::vector<int> Indices(Length);
974  std::vector<double> Values(Length);
975 
976  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
977  if (IsParallel_) {
978  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
979  }
980  else
981  Y2 = Teuchos::rcp( &Y, false );
982 
983  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
984  X.ExtractView(&x_ptr);
985  Y.ExtractView(&y_ptr);
986  Y2->ExtractView(&y2_ptr);
987  Diagonal_->ExtractView(&d_ptr);
988 
989  for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
990 
991  // only one data exchange per sweep
992  if (IsParallel_)
993  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
994 
995  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
996  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
997 
998  int NumEntries;
999  int col;
1000  double diag = d_ptr[i];
1001 
1002  IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
1003  &Values[0], &Indices[0]));
1004 
1005  for (int m = 0 ; m < NumVectors ; ++m) {
1006 
1007  double dtemp = 0.0;
1008 
1009  for (int k = 0 ; k < NumEntries ; ++k) {
1010 
1011  col = Indices[k];
1012  dtemp += Values[k] * y2_ptr[m][col];
1013  }
1014 
1015  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1016  }
1017  }
1018  for (int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1019  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1020 
1021  int NumEntries;
1022  int col;
1023  double diag = d_ptr[i];
1024 
1025  IFPACK_CHK_ERR(Matrix_->ExtractMyRowCopy(i, Length,NumEntries,
1026  &Values[0], &Indices[0]));
1027 
1028  for (int m = 0 ; m < NumVectors ; ++m) {
1029 
1030  double dtemp = 0.0;
1031  for (int k = 0 ; k < NumEntries ; ++k) {
1032 
1033  col = Indices[k];
1034  dtemp += Values[k] * y2_ptr[m][col];
1035  }
1036 
1037  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1038 
1039  }
1040  }
1041 
1042  if (IsParallel_)
1043  for (int m = 0 ; m < NumVectors ; ++m)
1044  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1045  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1046  y_ptr[m][i] = y2_ptr[m][i];
1047  }
1048  }
1049 
1050 #ifdef IFPACK_FLOPCOUNTERS
1051  ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1052 #endif
1053  return(0);
1054 }
1055 
1056 //==============================================================================
1057 int Ifpack_PointRelaxation::
1058 ApplyInverseSGS_CrsMatrix(const Epetra_CrsMatrix* A, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
1059 {
1060  int NumVectors = X.NumVectors();
1061 
1062  int* Indices;
1063  double* Values;
1064 
1065  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1066  if (IsParallel_) {
1067  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1068  }
1069  else
1070  Y2 = Teuchos::rcp( &Y, false );
1071 
1072  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1073  X.ExtractView(&x_ptr);
1074  Y.ExtractView(&y_ptr);
1075  Y2->ExtractView(&y2_ptr);
1076  Diagonal_->ExtractView(&d_ptr);
1077 
1078  for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1079 
1080  // only one data exchange per sweep
1081  if (IsParallel_)
1082  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1083 
1084  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1085  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1086 
1087  int NumEntries;
1088  int col;
1089  double diag = d_ptr[i];
1090 
1091  IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
1092 
1093  for (int m = 0 ; m < NumVectors ; ++m) {
1094 
1095  double dtemp = 0.0;
1096 
1097  for (int k = 0; k < NumEntries; ++k) {
1098 
1099  col = Indices[k];
1100  dtemp += Values[k] * y2_ptr[m][col];
1101  }
1102 
1103  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1104  }
1105  }
1106  for (int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1107  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1108 
1109  int NumEntries;
1110  int col;
1111  double diag = d_ptr[i];
1112 
1113  IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
1114 
1115  for (int m = 0 ; m < NumVectors ; ++m) {
1116 
1117  double dtemp = 0.0;
1118  for (int k = 0; k < NumEntries; ++k) {
1119 
1120  col = Indices[k];
1121  dtemp += Values[k] * y2_ptr[m][col];
1122  }
1123 
1124  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1125 
1126  }
1127  }
1128 
1129  if (IsParallel_)
1130  for (int m = 0 ; m < NumVectors ; ++m)
1131  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1132  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1133  y_ptr[m][i] = y2_ptr[m][i];
1134  }
1135  }
1136 
1137 #ifdef IFPACK_FLOPCOUNTERS
1138  ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1139 #endif
1140  return(0);
1141 }
1142 
1143 //==============================================================================
1144 // this requires Epetra_CrsMatrix + OptimizeStorage()
1145 int Ifpack_PointRelaxation::
1146 ApplyInverseSGS_FastCrsMatrix(const Epetra_CrsMatrix* A, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
1147 {
1148  int* IndexOffset;
1149  int* Indices;
1150  double* Values;
1151  IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1152 
1153  int NumVectors = X.NumVectors();
1154 
1155  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1156  if (IsParallel_) {
1157  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1158  }
1159  else
1160  Y2 = Teuchos::rcp( &Y, false );
1161 
1162  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1163  X.ExtractView(&x_ptr);
1164  Y.ExtractView(&y_ptr);
1165  Y2->ExtractView(&y2_ptr);
1166  Diagonal_->ExtractView(&d_ptr);
1167 
1168  for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1169 
1170  // only one data exchange per sweep
1171  if (IsParallel_)
1172  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1173 
1174  for (int i = 0 ; i < NumMyRows_ ; ++i) {
1175 
1176  int col;
1177  double diag = d_ptr[i];
1178 
1179  // no need to extract row i
1180  //IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
1181 
1182  for (int m = 0 ; m < NumVectors ; ++m) {
1183 
1184  double dtemp = 0.0;
1185 
1186  for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1187 
1188  col = Indices[k];
1189  dtemp += Values[k] * y2_ptr[m][col];
1190  }
1191 
1192  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1193  }
1194  }
1195 
1196  for (int i = NumMyRows_ - 1 ; i > -1 ; --i) {
1197 
1198  int col;
1199  double diag = d_ptr[i];
1200 
1201  // no need to extract row i
1202  //IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
1203 
1204  for (int m = 0 ; m < NumVectors ; ++m) {
1205 
1206  double dtemp = 0.0;
1207  for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1208 
1209  col = Indices[k];
1210  dtemp += Values[k] * y2_ptr[m][col];
1211  }
1212 
1213  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1214 
1215  }
1216  }
1217 
1218  if (IsParallel_)
1219  for (int m = 0 ; m < NumVectors ; ++m)
1220  for (int i = 0 ; i < NumMyRows_ ; ++i)
1221  y_ptr[m][i] = y2_ptr[m][i];
1222  }
1223 
1224 #ifdef IFPACK_FLOPCOUNTERS
1225  ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1226 #endif
1227  return(0);
1228 }
1229 
1230 
1231 //==============================================================================
1232 // this requires Epetra_CrsMatrix + OptimizeStorage() + LocalSmoothingIndices_
1233 int Ifpack_PointRelaxation::
1234 ApplyInverseSGS_LocalFastCrsMatrix(const Epetra_CrsMatrix* A, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
1235 {
1236  int* IndexOffset;
1237  int* Indices;
1238  double* Values;
1239  IFPACK_CHK_ERR(A->ExtractCrsDataPointers(IndexOffset, Indices, Values));
1240 
1241  int NumVectors = X.NumVectors();
1242 
1243  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
1244  if (IsParallel_) {
1245  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
1246  }
1247  else
1248  Y2 = Teuchos::rcp( &Y, false );
1249 
1250  double** y_ptr, ** y2_ptr, ** x_ptr, *d_ptr;
1251  X.ExtractView(&x_ptr);
1252  Y.ExtractView(&y_ptr);
1253  Y2->ExtractView(&y2_ptr);
1254  Diagonal_->ExtractView(&d_ptr);
1255 
1256  for (int iter = 0 ; iter < NumSweeps_ ; ++iter) {
1257 
1258  // only one data exchange per sweep
1259  if (IsParallel_)
1260  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
1261 
1262  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1263  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1264 
1265  int col;
1266  double diag = d_ptr[i];
1267 
1268  // no need to extract row i
1269  //IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
1270 
1271  for (int m = 0 ; m < NumVectors ; ++m) {
1272 
1273  double dtemp = 0.0;
1274 
1275  for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1276 
1277  col = Indices[k];
1278  dtemp += Values[k] * y2_ptr[m][col];
1279  }
1280 
1281  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1282  }
1283  }
1284 
1285  for (int ii = NumLocalSmoothingIndices_ - 1 ; ii > -1 ; --ii) {
1286  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1287 
1288  int col;
1289  double diag = d_ptr[i];
1290 
1291  // no need to extract row i
1292  //IFPACK_CHK_ERR(A->ExtractMyRowView(i, NumEntries, Values, Indices));
1293 
1294  for (int m = 0 ; m < NumVectors ; ++m) {
1295 
1296  double dtemp = 0.0;
1297  for (int k = IndexOffset[i] ; k < IndexOffset[i + 1] ; ++k) {
1298 
1299  col = Indices[k];
1300  dtemp += Values[k] * y2_ptr[m][col];
1301  }
1302 
1303  y2_ptr[m][i] += DampingFactor_ * (x_ptr[m][i] - dtemp) * diag;
1304 
1305  }
1306  }
1307 
1308  if (IsParallel_)
1309  for (int m = 0 ; m < NumVectors ; ++m)
1310  for (int ii = 0 ; ii < NumLocalSmoothingIndices_ ; ++ii) {
1311  int i = (!LocalSmoothingIndices_)? ii : LocalSmoothingIndices_[ii];
1312  y_ptr[m][i] = y2_ptr[m][i];
1313  }
1314  }
1315 
1316 #ifdef IFPACK_FLOPCOUNTERS
1317  ApplyInverseFlops_ += NumVectors * (8 * NumGlobalRows_ + 4 * NumGlobalNonzeros_);
1318 #endif
1319  return(0);
1320 }
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
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.
bool StorageOptimized() const
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Ifpack_PointRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_PointRelaxation constructor with given Epetra_RowMatrix.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual int MaxNumEntries() const =0
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual int Compute()
Computes the preconditioners.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const