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 //==============================================================================
62 Ifpack_PointRelaxation::
63 Ifpack_PointRelaxation(const Epetra_RowMatrix* Matrix_in) :
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 //==============================================================================
153 const Epetra_Comm& Ifpack_PointRelaxation::Comm() const
154 {
155  return(Matrix_->Comm());
156 }
157 
158 //==============================================================================
159 const Epetra_Map& Ifpack_PointRelaxation::OperatorDomainMap() const
160 {
161  return(Matrix_->OperatorDomainMap());
162 }
163 
164 //==============================================================================
165 const Epetra_Map& Ifpack_PointRelaxation::OperatorRangeMap() const
166 {
167  return(Matrix_->OperatorRangeMap());
168 }
169 
170 //==============================================================================
171 int Ifpack_PointRelaxation::Initialize()
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 //==============================================================================
201 int Ifpack_PointRelaxation::Compute()
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 //==============================================================================
368 double Ifpack_PointRelaxation::
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 //
413 int Ifpack_PointRelaxation::
414 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
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 }
Epetra_MultiVector
Ifpack_PointRelaxation::Print
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
Definition: Ifpack_PointRelaxation.cpp:302
Epetra_Map
Ifpack_PointRelaxation::Apply
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
Definition: Ifpack_PointRelaxation.h:172
Ifpack_PointRelaxation::Comm
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Definition: Ifpack_PointRelaxation.cpp:153
Epetra_CrsMatrix::StorageOptimized
bool StorageOptimized() const
Ifpack_PointRelaxation::OperatorDomainMap
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Definition: Ifpack_PointRelaxation.cpp:159
Epetra_VbrMatrix
Epetra_Vector
Ifpack_PointRelaxation::Ifpack_PointRelaxation
Ifpack_PointRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_PointRelaxation constructor with given Epetra_RowMatrix.
Definition: Ifpack_PointRelaxation.cpp:63
Epetra_Comm::MinAll
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Import
Ifpack_PointRelaxation::Condest
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
Definition: Ifpack_PointRelaxation.h:264
Ifpack_PointRelaxation::SetParameters
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
Definition: Ifpack_PointRelaxation.cpp:97
Epetra_Time
Epetra_CrsMatrix::ExtractMyRowView
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Epetra_RowMatrix::MaxNumEntries
virtual int MaxNumEntries() const =0
Epetra_Comm
Ifpack_PointRelaxation::IsInitialized
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
Definition: Ifpack_PointRelaxation.h:234
Ifpack_PointRelaxation::Compute
virtual int Compute()
Computes the preconditioners.
Definition: Ifpack_PointRelaxation.cpp:201
Ifpack_PointRelaxation::Initialize
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
Definition: Ifpack_PointRelaxation.cpp:171
Ifpack_PointRelaxation::Matrix
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
Definition: Ifpack_PointRelaxation.h:252
Ifpack_PointRelaxation::IsComputed
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
Definition: Ifpack_PointRelaxation.h:240
Ifpack_Utils.h
Epetra_CrsMatrix
Epetra_RowMatrix
Ifpack_PointRelaxation::OperatorRangeMap
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Definition: Ifpack_PointRelaxation.cpp:165
Ifpack_PointRelaxation::ApplyInverse
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
Definition: Ifpack_PointRelaxation.cpp:414
Ifpack_toString
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
Definition: Ifpack_Utils.cpp:242
Epetra_CrsMatrix::ExtractCrsDataPointers
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
Generated by   doxygen 1.8.5