IFPACK Development
Loading...
Searching...
No Matches
Ifpack_ICT.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 "Ifpack_Preconditioner.h"
45#include "Ifpack_ICT.h"
46#include "Ifpack_Condest.h"
47#include "Ifpack_Utils.h"
48#include "Ifpack_HashTable.h"
49#include "Epetra_SerialComm.h"
50#include "Epetra_Comm.h"
51#include "Epetra_Map.h"
52#include "Epetra_RowMatrix.h"
53#include "Epetra_CrsMatrix.h"
54#include "Epetra_Vector.h"
55#include "Epetra_MultiVector.h"
56#include "Epetra_Util.h"
57#include "Teuchos_ParameterList.hpp"
58#include "Teuchos_RefCountPtr.hpp"
59#include <functional>
60
61//==============================================================================
62// FIXME: allocate Comm_ and Time_ the first Initialize() call
64 A_(*A),
65 Comm_(A_.Comm()),
66 Condest_(-1.0),
67 Athresh_(0.0),
68 Rthresh_(1.0),
69 LevelOfFill_(1.0),
70 DropTolerance_(0.0),
71 Relax_(0.0),
72 IsInitialized_(false),
73 IsComputed_(false),
74 UseTranspose_(false),
75 NumMyRows_(0),
76 NumInitialize_(0),
77 NumCompute_(0),
78 NumApplyInverse_(0),
79 InitializeTime_(0.0),
80 ComputeTime_(0.0),
81 ApplyInverseTime_(0.0),
82 ComputeFlops_(0.0),
83 ApplyInverseFlops_(0.0),
84 Time_(Comm()),
85 GlobalNonzeros_(0)
86{
87 // do nothing here
88}
89
90//==============================================================================
92{
93 Destroy();
94}
95
96//==============================================================================
97void Ifpack_ICT::Destroy()
98{
99 IsInitialized_ = false;
100 IsComputed_ = false;
101}
102
103//==========================================================================
104int Ifpack_ICT::SetParameters(Teuchos::ParameterList& List)
105{
106 using std::cerr;
107 using std::endl;
108
109 try
110 {
111 LevelOfFill_ = List.get("fact: ict level-of-fill",LevelOfFill_);
112 Athresh_ = List.get("fact: absolute threshold", Athresh_);
113 Rthresh_ = List.get("fact: relative threshold", Rthresh_);
114 Relax_ = List.get("fact: relax value", Relax_);
115 DropTolerance_ = List.get("fact: drop tolerance", DropTolerance_);
116
117 // set label
118 Label_ = "ICT (fill=" + Ifpack_toString(LevelOfFill())
119 + ", athr=" + Ifpack_toString(AbsoluteThreshold())
120 + ", rthr=" + Ifpack_toString(RelativeThreshold())
121 + ", relax=" + Ifpack_toString(RelaxValue())
122 + ", droptol=" + Ifpack_toString(DropTolerance())
123 + ")";
124
125 return(0);
126 }
127 catch (...)
128 {
129 cerr << "Caught an exception while parsing the parameter list" << endl;
130 cerr << "This typically means that a parameter was set with the" << endl;
131 cerr << "wrong type (for example, int instead of double). " << endl;
132 cerr << "please check the documentation for the type required by each parameter." << endl;
133 IFPACK_CHK_ERR(-1);
134 }
135}
136
137//==========================================================================
139{
140 // clean data if present
141 Destroy();
142
143 Time_.ResetStartTime();
144
145 // matrix must be square. Check only on one processor
146 if (Comm().NumProc() == 1 && Matrix().NumMyRows() != Matrix().NumMyCols())
147 IFPACK_CHK_ERR(-2);
148
149 NumMyRows_ = Matrix().NumMyRows();
150
151 // nothing else to do here
152 IsInitialized_ = true;
153 ++NumInitialize_;
154 InitializeTime_ += Time_.ElapsedTime();
155
156 return(0);
157}
158
159//==========================================================================
160template<typename int_type>
161int Ifpack_ICT::TCompute()
162{
163 if (!IsInitialized())
164 IFPACK_CHK_ERR(Initialize());
165
166 Time_.ResetStartTime();
167 IsComputed_ = false;
168
169 NumMyRows_ = A_.NumMyRows();
170 int Length = A_.MaxNumEntries();
171 std::vector<int> RowIndices(Length);
172 std::vector<double> RowValues(Length);
173
174 bool distributed = (Comm().NumProc() > 1)?true:false;
175
176#if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
177 if (distributed)
178 {
179 SerialComm_ = Teuchos::rcp(new Epetra_SerialComm);
180#if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES)
182 SerialMap_ = Teuchos::rcp(new Epetra_Map(NumMyRows_, 0, *SerialComm_));
183 else
184#endif
185#if !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
187 SerialMap_ = Teuchos::rcp(new Epetra_Map((long long) NumMyRows_, (long long) 0, *SerialComm_));
188 else
189#endif
190 throw "Ifpack_ICT::TCompute: Global indices unknown.";
191 assert (SerialComm_.get() != 0);
192 assert (SerialMap_.get() != 0);
193 }
194 else
195 SerialMap_ = Teuchos::rcp(const_cast<Epetra_Map*>(&A_.RowMatrixRowMap()), false);
196#endif
197
198 int RowNnz;
199#ifdef IFPACK_FLOPCOUNTERS
200 double flops = 0.0;
201#endif
202
203 H_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy,*SerialMap_,0));
204 if (H_.get() == 0)
205 IFPACK_CHK_ERR(-5); // memory allocation error
206
207 // get A(0,0) element and insert it (after sqrt)
208 IFPACK_CHK_ERR(A_.ExtractMyRowCopy(0,Length,RowNnz,
209 &RowValues[0],&RowIndices[0]));
210
211 // skip off-processor elements
212 if (distributed)
213 {
214 int count = 0;
215 for (int i = 0 ;i < RowNnz ; ++i)
216 {
217 if (RowIndices[i] < NumMyRows_){
218 RowIndices[count] = RowIndices[i];
219 RowValues[count] = RowValues[i];
220 ++count;
221 }
222 else
223 continue;
224 }
225 RowNnz = count;
226 }
227
228 // modify diagonal
229 double diag_val = 0.0;
230 for (int i = 0 ;i < RowNnz ; ++i) {
231 if (RowIndices[i] == 0) {
232 double& v = RowValues[i];
233 diag_val = AbsoluteThreshold() * EPETRA_SGN(v) +
234 RelativeThreshold() * v;
235 break;
236 }
237 }
238
239 diag_val = sqrt(diag_val);
240 int_type diag_idx = 0;
241 EPETRA_CHK_ERR(H_->InsertGlobalValues(0,1,&diag_val, &diag_idx));
242
243 // The 10 is just a small constant to limit collisons as the actual keys
244 // we store are the indices and not integers
245 // [0..A_.MaxNumEntries()*LevelofFill()].
247
248 // start factorization for line 1
249 for (int row_i = 1 ; row_i < NumMyRows_ ; ++row_i) {
250
251 // get row `row_i' of the matrix
252 IFPACK_CHK_ERR(A_.ExtractMyRowCopy(row_i,Length,RowNnz,
253 &RowValues[0],&RowIndices[0]));
254
255 // skip off-processor elements
256 if (distributed)
257 {
258 int count = 0;
259 for (int i = 0 ;i < RowNnz ; ++i)
260 {
261 if (RowIndices[i] < NumMyRows_){
262 RowIndices[count] = RowIndices[i];
263 RowValues[count] = RowValues[i];
264 ++count;
265 }
266 else
267 continue;
268 }
269 RowNnz = count;
270 }
271
272 // number of nonzeros in this row are defined as the nonzeros
273 // of the matrix, plus the level of fill
274 int LOF = (int)(LevelOfFill() * RowNnz);
275 if (LOF == 0) LOF = 1;
276
277 // convert line `row_i' into hash for fast access
278 Hash.reset();
279
280 double h_ii = 0.0;
281 for (int i = 0 ; i < RowNnz ; ++i) {
282 if (RowIndices[i] == row_i) {
283 double& v = RowValues[i];
284 h_ii = AbsoluteThreshold() * EPETRA_SGN(v) + RelativeThreshold() * v;
285 }
286 else if (RowIndices[i] < row_i)
287 {
288 Hash.set(RowIndices[i], RowValues[i], true);
289 }
290 }
291
292 // form element (row_i, col_j)
293 // I start from the first row that has a nonzero column
294 // index in row_i.
295 for (int col_j = RowIndices[0] ; col_j < row_i ; ++col_j) {
296
297 double h_ij = 0.0, h_jj = 0.0;
298 // note: get() returns 0.0 if col_j is not found
299 h_ij = Hash.get(col_j);
300
301 // get pointers to row `col_j'
302 int_type* ColIndices;
303 double* ColValues;
304 int ColNnz;
305 int_type col_j_GID = (int_type) H_->RowMap().GID64(col_j);
306 H_->ExtractGlobalRowView(col_j_GID, ColNnz, ColValues, ColIndices);
307
308 for (int k = 0 ; k < ColNnz ; ++k) {
309 int_type col_k = ColIndices[k];
310
311 if (col_k == col_j)
312 h_jj = ColValues[k];
313 else {
314 double xxx = Hash.get(col_k);
315 if (xxx != 0.0)
316 {
317 h_ij -= ColValues[k] * xxx;
318#ifdef IFPACK_FLOPCOUNTERS
319 flops += 2.0;
320#endif
321 }
322 }
323 }
324
325 h_ij /= h_jj;
326
327 if (IFPACK_ABS(h_ij) > DropTolerance_)
328 {
329 Hash.set(col_j, h_ij);
330 }
331
332#ifdef IFPACK_FLOPCOUNTERS
333 // only approx
334 ComputeFlops_ += 2.0 * flops + 1.0;
335#endif
336 }
337
338 int size = Hash.getNumEntries();
339
340 std::vector<double> AbsRow(size);
341 int count = 0;
342
343 // +1 because I use the extra position for diagonal in insert
344 std::vector<int_type> keys(size + 1);
345 std::vector<double> values(size + 1);
346
347 Hash.arrayify(&keys[0], &values[0]);
348
349 for (int i = 0 ; i < size ; ++i)
350 {
351 AbsRow[i] = IFPACK_ABS(values[i]);
352 }
353 count = size;
354
355 double cutoff = 0.0;
356 if (count > LOF) {
357 nth_element(AbsRow.begin(), AbsRow.begin() + LOF, AbsRow.begin() + count,
358
359 std::greater<double>());
360 cutoff = AbsRow[LOF];
361 }
362
363 for (int i = 0 ; i < size ; ++i)
364 {
365 h_ii -= values[i] * values[i];
366 }
367
368 if (h_ii < 0.0) h_ii = 1e-12;;
369
370 h_ii = sqrt(h_ii);
371
372#ifdef IFPACK_FLOPCOUNTERS
373 // only approx, + 1 == sqrt
374 ComputeFlops_ += 2 * size + 1;
375#endif
376
377 double DiscardedElements = 0.0;
378
379 count = 0;
380 for (int i = 0 ; i < size ; ++i)
381 {
382 if (IFPACK_ABS(values[i]) > cutoff)
383 {
384 values[count] = values[i];
385 keys[count] = keys[i];
386 ++count;
387 }
388 else
389 DiscardedElements += values[i];
390 }
391
392 if (RelaxValue() != 0.0) {
393 DiscardedElements *= RelaxValue();
394 h_ii += DiscardedElements;
395 }
396
397 values[count] = h_ii;
398 keys[count] = (int_type) H_->RowMap().GID64(row_i);
399 ++count;
400
401 H_->InsertGlobalValues((int_type) H_->RowMap().GID64(row_i), count, &(values[0]), (int_type*)&(keys[0]));
402 }
403
404 IFPACK_CHK_ERR(H_->FillComplete());
405
406#if 0
407 // to check the complete factorization
408 Epetra_Vector LHS(Matrix().RowMatrixRowMap());
409 Epetra_Vector RHS1(Matrix().RowMatrixRowMap());
410 Epetra_Vector RHS2(Matrix().RowMatrixRowMap());
411 Epetra_Vector RHS3(Matrix().RowMatrixRowMap());
412 LHS.Random();
413
414 Matrix().Multiply(false,LHS,RHS1);
415 H_->Multiply(true,LHS,RHS2);
416 H_->Multiply(false,RHS2,RHS3);
417
418 RHS1.Update(-1.0, RHS3, 1.0);
419 std::cout << endl;
420 std::cout << RHS1;
421#endif
422 long long MyNonzeros = H_->NumGlobalNonzeros64();
423 Comm().SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
424
425 IsComputed_ = true;
426#ifdef IFPACK_FLOPCOUNTERS
427 double TotalFlops; // sum across all the processors
428 A_.Comm().SumAll(&flops, &TotalFlops, 1);
429 ComputeFlops_ += TotalFlops;
430#endif
431 ++NumCompute_;
432 ComputeTime_ += Time_.ElapsedTime();
433
434 return(0);
435
436}
437
439#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
441 return TCompute<int>();
442 }
443 else
444#endif
445#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
447 return TCompute<long long>();
448 }
449 else
450#endif
451 throw "Ifpack_ICT::Compute: GlobalIndices type unknown for A_";
452}
453
454//=============================================================================
456 Epetra_MultiVector& Y) const
457{
458
459 if (!IsComputed())
460 IFPACK_CHK_ERR(-3); // compute preconditioner first
461
462 if (X.NumVectors() != Y.NumVectors())
463 IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size
464
465 Time_.ResetStartTime();
466
467 // AztecOO gives X and Y pointing to the same memory location,
468 // need to create an auxiliary vector, Xcopy
469 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
470 if (X.Pointers()[0] == Y.Pointers()[0])
471 Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
472 else
473 Xcopy = Teuchos::rcp( &X, false );
474
475 // NOTE: H_ is based on SerialMap_, while Xcopy is based
476 // on A.Map()... which are in general different. However, Solve()
477 // does not seem to care... which is fine with me.
478 //
479 EPETRA_CHK_ERR(H_->Solve(false,false,false,*Xcopy,Y));
480 EPETRA_CHK_ERR(H_->Solve(false,true,false,Y,Y));
481
482#ifdef IFPACK_FLOPCOUNTERS
483 // these are global flop count
484 ApplyInverseFlops_ += 4.0 * GlobalNonzeros_;
485#endif
486
487 ++NumApplyInverse_;
488 ApplyInverseTime_ += Time_.ElapsedTime();
489
490 return(0);
491}
492//=============================================================================
493// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
494int Ifpack_ICT::Apply(const Epetra_MultiVector& /* X */,
495 Epetra_MultiVector& /* Y */) const
496{
497
498 IFPACK_CHK_ERR(-98);
499}
500
501//=============================================================================
502double Ifpack_ICT::Condest(const Ifpack_CondestType CT,
503 const int MaxIters, const double Tol,
504 Epetra_RowMatrix* Matrix_in)
505{
506 if (!IsComputed()) // cannot compute right now
507 return(-1.0);
508
509 // NOTE: this is computing the *local* condest
510 if (Condest_ == -1.0)
511 Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
512
513 return(Condest_);
514}
515
516//=============================================================================
517std::ostream&
518Ifpack_ICT::Print(std::ostream& os) const
519{
520 using std::endl;
521
522 if (!Comm().MyPID()) {
523 os << endl;
524 os << "================================================================================" << endl;
525 os << "Ifpack_ICT: " << Label() << endl << endl;
526 os << "Level-of-fill = " << LevelOfFill() << endl;
527 os << "Absolute threshold = " << AbsoluteThreshold() << endl;
528 os << "Relative threshold = " << RelativeThreshold() << endl;
529 os << "Relax value = " << RelaxValue() << endl;
530 os << "Condition number estimate = " << Condest() << endl;
531 os << "Global number of rows = " << Matrix().NumGlobalRows64() << endl;
532 if (IsComputed_) {
533 os << "Number of nonzeros of H = " << H_->NumGlobalNonzeros64() << endl;
534 os << "nonzeros / rows = "
535 << 1.0 * H_->NumGlobalNonzeros64() / H_->NumGlobalRows64() << endl;
536 }
537 os << endl;
538 os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
539 os << "----- ------- -------------- ------------ --------" << endl;
540 os << "Initialize() " << std::setw(5) << NumInitialize()
541 << " " << std::setw(15) << InitializeTime()
542 << " 0.0 0.0" << endl;
543 os << "Compute() " << std::setw(5) << NumCompute()
544 << " " << std::setw(15) << ComputeTime()
545 << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
546 if (ComputeTime() != 0.0)
547 os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
548 else
549 os << " " << std::setw(15) << 0.0 << endl;
550 os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
551 << " " << std::setw(15) << ApplyInverseTime()
552 << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
553 if (ApplyInverseTime() != 0.0)
554 os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
555 else
556 os << " " << std::setw(15) << 0.0 << endl;
557 os << "================================================================================" << endl;
558 os << endl;
559 }
560
561
562 return(os);
563}
bool GlobalIndicesInt() const
bool GlobalIndicesLongLong() const
virtual int NumProc() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
int NumVectors() const
double ** Pointers() const
virtual const Epetra_Comm & Comm() const=0
virtual int NumMyRows() const=0
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Map & RowMatrixRowMap() 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
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition Ifpack_ICT.h:224
double LevelOfFill() const
Returns the level-of-fill.
Definition Ifpack_ICT.h:299
double RelaxValue() const
Returns the relaxation value.
Definition Ifpack_ICT.h:317
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_ICT forward/back solve on a Epetra_MultiVector X in Y.
virtual int NumCompute() const
Returns the number of calls to Compute().
Definition Ifpack_ICT.h:248
virtual double InitializeTime() const
Returns the time spent in Initialize().
Definition Ifpack_ICT.h:260
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
Ifpack_ICT(const Epetra_RowMatrix *A)
Ifpack_ICT constuctor with variable number of indices per row.
double DropTolerance() const
Returns the drop threshold.
Definition Ifpack_ICT.h:323
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
Definition Ifpack_ICT.h:272
double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
Definition Ifpack_ICT.h:175
virtual double ComputeTime() const
Returns the time spent in Compute().
Definition Ifpack_ICT.h:266
double RelativeThreshold() const
Returns the relative threshold.
Definition Ifpack_ICT.h:311
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
Definition Ifpack_ICT.h:254
int SetParameters(Teuchos::ParameterList &parameterlis)
Set parameters using a Teuchos::ParameterList object.
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
Definition Ifpack_ICT.h:142
double AbsoluteThreshold() const
Returns the absolute threshold.
Definition Ifpack_ICT.h:305
virtual double ApplyInverseFlops() const
Returns the number of flops in all applications of ApplyInverse().
Definition Ifpack_ICT.h:290
int Initialize()
Initialize L and U with values from user matrix A.
virtual ~Ifpack_ICT()
Ifpack_ICT Destructor.
bool IsInitialized() const
Returns true is the preconditioner has been successfully initialized.
Definition Ifpack_ICT.h:116
virtual double ComputeFlops() const
Returns the number of flops in all applications of Compute().
Definition Ifpack_ICT.h:284
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
virtual int NumInitialize() const
Returns the number of calls to Initialize().
Definition Ifpack_ICT.h:242
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
Definition Ifpack_ICT.h:110