Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Ifpack_IC.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack: Object-Oriented Algebraic Preconditioner Package
5// Copyright (2002) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42
43#include "Ifpack_ConfigDefs.h"
45#include "Ifpack_IC.h"
46#include "Ifpack_IC_Utils.h"
47#include "Ifpack_Condest.h"
48#include "Epetra_Comm.h"
49#include "Epetra_Map.h"
50#include "Epetra_RowMatrix.h"
51#include "Epetra_CrsMatrix.h"
52#include "Epetra_Vector.h"
53#include "Epetra_MultiVector.h"
54#include "Epetra_Util.h"
55#include "Teuchos_ParameterList.hpp"
56#include "Teuchos_RefCountPtr.hpp"
57using Teuchos::RefCountPtr;
58using Teuchos::rcp;
59
60//==============================================================================
62 A_(rcp(A,false)),
63 Comm_(A->Comm()),
64 UseTranspose_(false),
65 Condest_(-1.0),
66 Athresh_(0.0),
67 Rthresh_(1.0),
68 Droptol_(0.0),
69 Lfil_(1.0),
70 Aict_(0),
71 Lict_(0),
72 Ldiag_(0),
73 IsInitialized_(false),
74 IsComputed_(false),
75 NumInitialize_(0),
76 NumCompute_(0),
77 NumApplyInverse_(0),
78 InitializeTime_(0.0),
79 ComputeTime_(0),
80 ApplyInverseTime_(0),
81 Time_(Comm()),
82 ComputeFlops_(0.0),
83 ApplyInverseFlops_(0.0)
84{
85 Teuchos::ParameterList List;
86 SetParameters(List);
87
88}
89//==============================================================================
91{
92 if (Lict_ != 0) {
94 delete [] Lict->ptr;
95 delete [] Lict->col;
96 delete [] Lict->val;
97 delete Lict;
98 }
99 if (Aict_ != 0) {
101 delete Aict;
102 }
103 if (Ldiag_ != 0) delete [] Ldiag_;
104
105 IsInitialized_ = false;
106 IsComputed_ = false;
107}
108
109//==========================================================================
110int Ifpack_IC::SetParameters(Teuchos::ParameterList& List)
111{
112
113 // Lfil_ = List.get("fact: level-of-fill",Lfil_); // Confusing parameter since Ifpack_IC can only do ICT not IC(k)
114 Lfil_ = List.get("fact: ict level-of-fill", Lfil_); // Lfil_ is now the fill ratio as in ICT (1.0 means same #nonzeros as A)
115 Athresh_ = List.get("fact: absolute threshold", Athresh_);
116 Rthresh_ = List.get("fact: relative threshold", Rthresh_);
117 Droptol_ = List.get("fact: drop tolerance", Droptol_);
118
119 // set label
120 sprintf(Label_, "IFPACK IC (fill=%f, drop=%f)",
121 Lfil_, Droptol_);
122 return(0);
123}
124
125//==========================================================================
127{
128
129 IsInitialized_ = false;
130
131 // FIXME: construction of ILUK graph must be here
132
134 IsInitialized_ = true;
135 return(0);
136}
137
138//==========================================================================
140{
141 // (re)allocate memory for ICT factors
142 U_ = rcp(new Epetra_CrsMatrix(Copy, Matrix().RowMatrixRowMap(),
143 Matrix().RowMatrixRowMap(), 0));
144 D_ = rcp(new Epetra_Vector(Matrix().RowMatrixRowMap()));
145
146 if (U_.get() == 0 || D_.get() == 0)
147 IFPACK_CHK_ERR(-5); // memory allocation error
148
149 int ierr = 0;
150 int i, j;
151 int NumIn, NumU;
152 bool DiagFound;
153 int NumNonzeroDiags = 0;
154
155 // Get Maximun Row length
156 int MaxNumEntries = Matrix().MaxNumEntries();
157
158 std::vector<int> InI(MaxNumEntries); // Allocate temp space
159 std::vector<int> UI(MaxNumEntries);
160 std::vector<double> InV(MaxNumEntries);
161 std::vector<double> UV(MaxNumEntries);
162
163 double *DV;
164 ierr = D_->ExtractView(&DV); // Get view of diagonal
165
166 // First we copy the user's matrix into diagonal vector and U, regardless of fill level
167
168 int NumRows = Matrix().NumMyRows();
169
170 for (i=0; i< NumRows; i++) {
171
172 Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0]); // Get Values and Indices
173
174 // Split into L and U (we don't assume that indices are ordered).
175 NumU = 0;
176 DiagFound = false;
177
178 for (j=0; j< NumIn; j++) {
179 int k = InI[j];
180
181 if (k==i) {
182 DiagFound = true;
183 // Store perturbed diagonal in Epetra_Vector D_
184 DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_;
185 }
186
187 else if (k < 0) return(-1); // Out of range
188 else if (i<k && k<NumRows) {
189 UI[NumU] = k;
190 UV[NumU] = InV[j];
191 NumU++;
192 }
193 }
194
195 // Check in things for this row of L and U
196
197 if (DiagFound) NumNonzeroDiags++;
198 if (NumU) U_->InsertMyValues(i, NumU, &UV[0], &UI[0]);
199
200 }
201
202 U_->FillComplete(Matrix().OperatorDomainMap(),
204
205 int ierr1 = 0;
206 if (NumNonzeroDiags<U_->NumMyRows()) ierr1 = 1;
207 Matrix().Comm().MaxAll(&ierr1, &ierr, 1);
208 IFPACK_CHK_ERR(ierr);
209
210 return(0);
211}
212
213//==========================================================================
215
216 if (!IsInitialized())
218
220 IsComputed_ = false;
221
222 // copy matrix into L and U.
224
225 int i;
226
227 int m, n, nz, Nrhs, ldrhs, ldlhs;
228 int * ptr=0, * ind;
229 double * val, * rhs, * lhs;
230
231 int ierr = Epetra_Util_ExtractHbData(U_.get(), 0, 0, m, n, nz, ptr, ind,
232 val, Nrhs, rhs, ldrhs, lhs, ldlhs);
233 if (ierr < 0)
234 IFPACK_CHK_ERR(ierr);
235
236 Ifpack_AIJMatrix * Aict;
237 if (Aict_==0) {
238 Aict = new Ifpack_AIJMatrix;
239 Aict_ = (void *) Aict;
240 }
241 else Aict = (Ifpack_AIJMatrix *) Aict_;
242 Ifpack_AIJMatrix * Lict;
243 if (Lict_==0) {
244 Lict = new Ifpack_AIJMatrix;
245 Lict_ = (void *) Lict;
246 }
247 else {
248 Lict = (Ifpack_AIJMatrix *) Lict_;
249 Ifpack_AIJMatrix_dealloc( Lict ); // Delete storage, crout_ict will allocate it again.
250 }
251 if (Ldiag_ != 0) delete [] Ldiag_; // Delete storage, crout_ict will allocate it again.
252 Aict->val = val;
253 Aict->col = ind;
254 Aict->ptr = ptr;
255 double *DV;
256 EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
257
258 // lfil is average number of nonzeros per row times fill ratio.
259 // Currently crout_ict keeps a constant number of nonzeros per row.
260 // TODO: Pass Lfil_ and make crout_ict keep variable #nonzeros per row.
261 int lfil = (Lfil_ * nz)/n + .5;
262
263 crout_ict(m, Aict, DV, Droptol_, lfil, Lict, &Ldiag_);
264
265 // Get rid of unnecessary data
266 delete [] ptr;
267
268 // Create Epetra View of L from crout_ict
269 U_ = rcp(new Epetra_CrsMatrix(View, A_->RowMatrixRowMap(), A_->RowMatrixRowMap(),0));
270 D_ = rcp(new Epetra_Vector(View, A_->RowMatrixRowMap(), Ldiag_));
271
272 ptr = Lict->ptr;
273 ind = Lict->col;
274 val = Lict->val;
275
276 for (i=0; i< m; i++) {
277 int NumEntries = ptr[i+1]-ptr[i];
278 int * Indices = ind+ptr[i];
279 double * Values = val+ptr[i];
280 U_->InsertMyValues(i, NumEntries, Values, Indices);
281 }
282
283 U_->FillComplete(A_->OperatorDomainMap(), A_->OperatorRangeMap());
284 D_->Reciprocal(*D_); // Put reciprocal of diagonal in this vector
285
286#ifdef IFPACK_FLOPCOUNTERS
287 double current_flops = 2 * nz; // Just an estimate
288 double total_flops = 0;
289
290 A_->Comm().SumAll(&current_flops, &total_flops, 1); // Get total madds across all PEs
291
292 ComputeFlops_ += total_flops;
293 // Now count the rest. NOTE: those flops are *global*
294 ComputeFlops_ += (double) U_->NumGlobalNonzeros(); // Accounts for multiplier above
295 ComputeFlops_ += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal
296#endif
297 ++NumCompute_;
299
300
301 IsComputed_ = true;
302
303 return(0);
304
305}
306
307//=============================================================================
308// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
310 Epetra_MultiVector& Y) const
311{
312
313 if (!IsComputed())
314 IFPACK_CHK_ERR(-3); // compute preconditioner first
315
316 if (X.NumVectors() != Y.NumVectors())
317 IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
318
320
321 bool Upper = true;
322 bool UnitDiagonal = true;
323
324 // AztecOO gives X and Y pointing to the same memory location,
325 // need to create an auxiliary vector, Xcopy
326 RefCountPtr< const Epetra_MultiVector > Xcopy;
327 if (X.Pointers()[0] == Y.Pointers()[0])
328 Xcopy = rcp( new Epetra_MultiVector(X) );
329 else
330 Xcopy = rcp( &X, false );
331
332 U_->Solve(Upper, true, UnitDiagonal, *Xcopy, Y);
333 Y.Multiply(1.0, *D_, Y, 0.0); // y = D*y (D_ has inverse of diagonal)
334 U_->Solve(Upper, false, UnitDiagonal, Y, Y); // Solve Uy = y
335
336#ifdef IFPACK_FLOPCOUNTERS
337 ApplyInverseFlops_ += 4.0 * U_->NumGlobalNonzeros();
338 ApplyInverseFlops_ += D_->GlobalLength();
339#endif
340
343
344 return(0);
345
346}
347
348//=============================================================================
349// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
351 Epetra_MultiVector& Y) const
352{
353
354 if (X.NumVectors() != Y.NumVectors())
355 IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
356
359
360 U_->Multiply(false, *X1, *Y1);
361 Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
362 Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
363 Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
364 U_->Multiply(true, Y1temp, *Y1);
365 Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
366 return(0);
367}
368
369//=============================================================================
371 const int MaxIters, const double Tol,
372 Epetra_RowMatrix* Matrix_in)
373{
374 if (!IsComputed()) // cannot compute right now
375 return(-1.0);
376
377 if (Condest_ == -1.0)
378 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
379
380 return(Condest_);
381}
382
383//=============================================================================
384std::ostream&
385Ifpack_IC::Print(std::ostream& os) const
386{
387 using std::endl;
388
389 if (!Comm().MyPID()) {
390 os << endl;
391 os << "================================================================================" << endl;
392 os << "Ifpack_IC: " << Label() << endl << endl;
393 os << "Level-of-fill = " << LevelOfFill() << endl;
394 os << "Absolute threshold = " << AbsoluteThreshold() << endl;
395 os << "Relative threshold = " << RelativeThreshold() << endl;
396 os << "Drop tolerance = " << DropTolerance() << endl;
397 os << "Condition number estimate = " << Condest() << endl;
398 os << "Global number of rows = " << A_->NumGlobalRows64() << endl;
399 if (IsComputed_) {
400 os << "Number of nonzeros of H = " << U_->NumGlobalNonzeros64() << endl;
401 os << "nonzeros / rows = "
402 << 1.0 * U_->NumGlobalNonzeros64() / U_->NumGlobalRows64() << endl;
403 }
404 os << endl;
405 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
406 os << "----- ------- -------------- ------------ --------" << endl;
407 os << "Initialize() " << std::setw(5) << NumInitialize()
408 << " " << std::setw(15) << InitializeTime()
409 << " 0.0 0.0" << endl;
410 os << "Compute() " << std::setw(5) << NumCompute()
411 << " " << std::setw(15) << ComputeTime()
412 << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
413 if (ComputeTime() != 0.0)
414 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
415 else
416 os << " " << std::setw(15) << 0.0 << endl;
417 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
418 << " " << std::setw(15) << ApplyInverseTime()
419 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
420 if (ApplyInverseTime() != 0.0)
421 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
422 else
423 os << " " << std::setw(15) << 0.0 << endl;
424 os << "================================================================================" << endl;
425 os << endl;
426 }
427
428
429 return(os);
430}
#define EPETRA_SGN(x)
#define EPETRA_CHK_ERR(a)
View
Copy
int Epetra_Util_ExtractHbData(Epetra_CrsMatrix *A, Epetra_MultiVector *LHS, Epetra_MultiVector *RHS, int &M, int &N, int &nz, int *&ptr, int *&ind, double *&val, int &Nrhs, double *&rhs, int &ldrhs, double *&lhs, int &ldlhs)
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
void crout_ict(int n, const Ifpack_AIJMatrix *AL, const double *Adiag, double droptol, int lfil, Ifpack_AIJMatrix *L, double **pdiag)
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const=0
int NumVectors() const
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
double ** Pointers() const
int ReciprocalMultiply(double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
virtual const Epetra_Comm & Comm() const=0
virtual int NumMyRows() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
void ResetStartTime(void)
double ElapsedTime(void) const
virtual double InitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack_IC.h:284
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Definition: Ifpack_IC.h:278
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
Definition: Ifpack_IC.h:131
int Initialize()
Initialize L and U with values from user matrix A.
Definition: Ifpack_IC.cpp:126
void * Lict_
Definition: Ifpack_IC.h:353
int SetParameters(Teuchos::ParameterList &parameterlis)
Set parameters using a Teuchos::ParameterList object.
Definition: Ifpack_IC.cpp:110
int ComputeSetup()
Definition: Ifpack_IC.cpp:139
double LevelOfFill() const
Definition: Ifpack_IC.h:320
bool IsInitialized_
Definition: Ifpack_IC.h:357
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Definition: Ifpack_IC.cpp:350
virtual int NumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack_IC.h:272
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
Definition: Ifpack_IC.cpp:214
double Droptol_
Definition: Ifpack_IC.h:349
Ifpack_IC(Epetra_RowMatrix *A)
Ifpack_IC constuctor with variable number of indices per row.
Definition: Ifpack_IC.cpp:61
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Definition: Ifpack_IC.h:312
char Label_[160]
Definition: Ifpack_IC.h:355
double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
Definition: Ifpack_IC.h:191
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_IC forward/back solve on a Epetra_MultiVector X in Y.
Definition: Ifpack_IC.cpp:309
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
Definition: Ifpack_IC.h:365
double ComputeFlops_
Contains the number of flops for Compute().
Definition: Ifpack_IC.h:377
Epetra_Time Time_
Used for timing purposes.
Definition: Ifpack_IC.h:374
double DropTolerance() const
Definition: Ifpack_IC.h:335
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack_IC.h:266
double Rthresh_
Definition: Ifpack_IC.h:348
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Definition: Ifpack_IC.h:245
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Definition: Ifpack_IC.h:242
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
Definition: Ifpack_IC.h:379
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
Definition: Ifpack_IC.h:158
double ComputeTime_
Contains the time for all successful calls to Compute().
Definition: Ifpack_IC.h:370
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Definition: Ifpack_IC.h:340
virtual ~Ifpack_IC()
Ifpack_IC Destructor.
Definition: Ifpack_IC.cpp:90
const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
Definition: Ifpack_IC.h:121
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
Definition: Ifpack_IC.h:307
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition: Ifpack_IC.h:248
double Condest_
Definition: Ifpack_IC.h:346
double Lfil_
Definition: Ifpack_IC.h:350
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
Definition: Ifpack_IC.h:296
Teuchos::RefCountPtr< Epetra_Vector > D_
Definition: Ifpack_IC.h:343
const char * Label() const
Definition: Ifpack_IC.h:251
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Definition: Ifpack_IC.h:372
void * Aict_
Definition: Ifpack_IC.h:352
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
Definition: Ifpack_IC.cpp:385
double * Ldiag_
Definition: Ifpack_IC.h:354
double Athresh_
Definition: Ifpack_IC.h:347
double RelativeThreshold() const
Definition: Ifpack_IC.h:330
int NumInitialize_
Contains the number of successful calls to Initialize().
Definition: Ifpack_IC.h:361
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
Definition: Ifpack_IC.h:342
virtual double ComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack_IC.h:290
int NumCompute_
Contains the number of successful call to Compute().
Definition: Ifpack_IC.h:363
double AbsoluteThreshold() const
Definition: Ifpack_IC.h:325
bool IsComputed_
Definition: Ifpack_IC.h:358
#define false