Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
44 #include "Ifpack_Preconditioner.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"
57 using Teuchos::RefCountPtr;
58 using 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 {
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 //==========================================================================
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 
133  ++NumInitialize_;
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 //=============================================================================
384 std::ostream&
385 Ifpack_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 }
double ComputeTime_
Contains the time for all successful calls to Compute().
Definition: Ifpack_IC.h:376
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack_IC.h:272
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
Definition: Ifpack_IC.h:164
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Definition: Ifpack_IC.h:248
double Condest_
Definition: Ifpack_IC.h:352
char Label_[160]
Definition: Ifpack_IC.h:361
int SetParameters(Teuchos::ParameterList &parameterlis)
Set parameters using a Teuchos::ParameterList object.
Definition: Ifpack_IC.cpp:110
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
Definition: Ifpack_IC.h:302
double RelativeThreshold() const
Definition: Ifpack_IC.h:336
T & get(ParameterList &l, const std::string &name)
double ElapsedTime(void) const
Epetra_Time Time_
Used for timing purposes.
Definition: Ifpack_IC.h:380
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Definition: Ifpack_IC.h:346
double AbsoluteThreshold() const
Definition: Ifpack_IC.h:331
Ifpack_IC(Epetra_RowMatrix *A)
Ifpack_IC constuctor with variable number of indices per row.
Definition: Ifpack_IC.cpp:61
#define EPETRA_CHK_ERR(a)
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:127
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Definition: Ifpack_IC.h:378
Teuchos::RefCountPtr< Epetra_Vector > D_
Definition: Ifpack_IC.h:349
virtual double ComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack_IC.h:296
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
Definition: Ifpack_IC.h:251
double Athresh_
Definition: Ifpack_IC.h:353
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
bool IsComputed_
Definition: Ifpack_IC.h:364
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
Definition: Ifpack_IC.h:348
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
int Initialize()
Initialize L and U with values from user matrix A.
Definition: Ifpack_IC.cpp:126
virtual int MaxNumEntries() const =0
virtual const Epetra_Comm & Comm() const =0
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
Definition: Ifpack_IC.h:385
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
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
double Droptol_
Definition: Ifpack_IC.h:355
double LevelOfFill() const
Definition: Ifpack_IC.h:326
#define EPETRA_SGN(x)
double ComputeFlops_
Contains the number of flops for Compute().
Definition: Ifpack_IC.h:383
virtual int NumMyRows() const =0
double * Ldiag_
Definition: Ifpack_IC.h:360
int ComputeSetup()
Definition: Ifpack_IC.cpp:139
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
Definition: Ifpack_IC.h:318
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator&lt;&lt;.
Definition: Ifpack_IC.cpp:385
const char * Label() const
Definition: Ifpack_IC.h:257
#define false
int NumCompute_
Contains the number of successful call to Compute().
Definition: Ifpack_IC.h:369
void * Aict_
Definition: Ifpack_IC.h:358
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Definition: Ifpack_IC.cpp:350
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
Definition: Ifpack_IC.h:137
double Rthresh_
Definition: Ifpack_IC.h:354
void * Lict_
Definition: Ifpack_IC.h:359
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
Definition: Ifpack_IC.h:197
double DropTolerance() const
Definition: Ifpack_IC.h:341
virtual double InitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack_IC.h:290
bool IsInitialized_
Definition: Ifpack_IC.h:363
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Definition: Ifpack_IC.h:284
virtual int NumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack_IC.h:278
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
Definition: Ifpack_IC.h:313
double Lfil_
Definition: Ifpack_IC.h:356
#define IFPACK_CHK_ERR(ifpack_err)
void ResetStartTime(void)
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
Definition: Ifpack_IC.cpp:214
int n
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition: Ifpack_IC.h:254
void crout_ict(int n, const Ifpack_AIJMatrix *AL, const double *Adiag, double droptol, int lfil, Ifpack_AIJMatrix *L, double **pdiag)
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)
int NumInitialize_
Contains the number of successful calls to Initialize().
Definition: Ifpack_IC.h:367
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
Definition: Ifpack_IC.h:371