IFPACK Development
Loading...
Searching...
No Matches
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
57static const int IFPACK_JACOBI = 0;
58static const int IFPACK_GS = 1;
59static 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//==============================================================================
97int 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//==============================================================================
302std::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//==============================================================================
369Condest(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//==============================================================================
384void 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).
456int Ifpack_PointRelaxation::
457ApplyInverseJacobi(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//==============================================================================
498int Ifpack_PointRelaxation::
499ApplyInverseGS(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//==============================================================================
523int Ifpack_PointRelaxation::
524ApplyInverseGS_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//==============================================================================
668int Ifpack_PointRelaxation::
669ApplyInverseGS_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
764int Ifpack_PointRelaxation::
765ApplyInverseGS_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
854int Ifpack_PointRelaxation::
855ApplyInverseGS_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//==============================================================================
945int Ifpack_PointRelaxation::
946ApplyInverseSGS(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//==============================================================================
968int Ifpack_PointRelaxation::
969ApplyInverseSGS_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//==============================================================================
1057int Ifpack_PointRelaxation::
1058ApplyInverseSGS_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()
1145int Ifpack_PointRelaxation::
1146ApplyInverseSGS_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_
1233int Ifpack_PointRelaxation::
1234ApplyInverseSGS_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 int MinAll(double *PartialMins, double *GlobalMins, int Count) const=0
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
bool StorageOptimized() const
const Epetra_BlockMap & Map() const
int NumVectors() const
int ExtractView(double **A, int *MyLDA) const
double ** Pointers() const
int PutScalar(double ScalarConstant)
virtual int MaxNumEntries() const=0
const Epetra_BlockMap & RowMap() const
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual int Compute()
Computes the preconditioners.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual std::ostream & Print(std::ostream &os) const
Prints object to an output stream.
virtual double Condest() const
Returns the condition number estimate, or -1.0 if not computed.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Ifpack_PointRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_PointRelaxation constructor with given Epetra_RowMatrix.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual int Initialize()
Computes all it is necessary to initialize the preconditioner.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.