43#include "Ifpack_ConfigDefs.h"
44#include "Ifpack_Preconditioner.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;
73 IsInitialized_(false),
83 ApplyInverseFlops_(0.0)
85 Teuchos::ParameterList List;
103 if (Ldiag_ != 0)
delete [] Ldiag_;
105 IsInitialized_ =
false;
114 Lfil_ = List.get(
"fact: ict level-of-fill", Lfil_);
115 Athresh_ = List.get(
"fact: absolute threshold", Athresh_);
116 Rthresh_ = List.get(
"fact: relative threshold", Rthresh_);
117 Droptol_ = List.get(
"fact: drop tolerance", Droptol_);
120 sprintf(Label_,
"IFPACK IC (fill=%f, drop=%f)",
129 IsInitialized_ =
false;
134 IsInitialized_ =
true;
139int Ifpack_IC::ComputeSetup()
143 Matrix().RowMatrixRowMap(), 0));
146 if (U_.get() == 0 || D_.get() == 0)
153 int NumNonzeroDiags = 0;
158 std::vector<int> InI(MaxNumEntries);
159 std::vector<int> UI(MaxNumEntries);
160 std::vector<double> InV(MaxNumEntries);
161 std::vector<double> UV(MaxNumEntries);
164 ierr = D_->ExtractView(&DV);
170 for (i=0; i< NumRows; i++) {
178 for (j=0; j< NumIn; j++) {
184 DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_;
187 else if (k < 0)
return(-1);
188 else if (i<k && k<NumRows) {
197 if (DiagFound) NumNonzeroDiags++;
198 if (NumU) U_->InsertMyValues(i, NumU, &UV[0], &UI[0]);
206 if (NumNonzeroDiags<U_->NumMyRows()) ierr1 = 1;
208 IFPACK_CHK_ERR(ierr);
223 IFPACK_CHK_ERR(ComputeSetup());
227 int m, n, nz, Nrhs, ldrhs, ldlhs;
229 double * val, * rhs, * lhs;
231 int ierr = Epetra_Util_ExtractHbData(U_.get(), 0, 0, m, n, nz, ptr, ind,
232 val, Nrhs, rhs, ldrhs, lhs, ldlhs);
234 IFPACK_CHK_ERR(ierr);
239 Aict_ = (
void *) Aict;
245 Lict_ = (
void *) Lict;
249 Ifpack_AIJMatrix_dealloc( Lict );
251 if (Ldiag_ != 0)
delete [] Ldiag_;
256 EPETRA_CHK_ERR(D_->ExtractView(&DV));
261 int lfil = (Lfil_ * nz)/n + .5;
263 crout_ict(m, Aict, DV, Droptol_, lfil, Lict, &Ldiag_);
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);
283 U_->FillComplete(A_->OperatorDomainMap(), A_->OperatorRangeMap());
286#ifdef IFPACK_FLOPCOUNTERS
287 double current_flops = 2 * nz;
288 double total_flops = 0;
290 A_->Comm().SumAll(¤t_flops, &total_flops, 1);
292 ComputeFlops_ += total_flops;
294 ComputeFlops_ += (double) U_->NumGlobalNonzeros();
295 ComputeFlops_ += (double) D_->GlobalLength();
322 bool UnitDiagonal =
true;
326 RefCountPtr< const Epetra_MultiVector > Xcopy;
330 Xcopy = rcp( &X,
false );
332 U_->Solve(Upper,
true, UnitDiagonal, *Xcopy, Y);
334 U_->Solve(Upper,
false, UnitDiagonal, Y, Y);
336#ifdef IFPACK_FLOPCOUNTERS
337 ApplyInverseFlops_ += 4.0 * U_->NumGlobalNonzeros();
338 ApplyInverseFlops_ += D_->GlobalLength();
360 U_->Multiply(
false, *X1, *Y1);
361 Y1->
Update(1.0, *X1, 1.0);
364 U_->Multiply(
true, Y1temp, *Y1);
365 Y1->
Update(1.0, Y1temp, 1.0);
371 const int MaxIters,
const double Tol,
377 if (Condest_ == -1.0)
378 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
389 if (!
Comm().MyPID()) {
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;
400 os <<
"Number of nonzeros of H = " << U_->NumGlobalNonzeros64() << endl;
401 os <<
"nonzeros / rows = "
402 << 1.0 * U_->NumGlobalNonzeros64() / U_->NumGlobalRows64() << endl;
405 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
406 os <<
"----- ------- -------------- ------------ --------" << endl;
409 <<
" 0.0 0.0" << endl;
410 os <<
"Compute() " << std::setw(5) <<
NumCompute()
416 os <<
" " << std::setw(15) << 0.0 << endl;
423 os <<
" " << std::setw(15) << 0.0 << endl;
424 os <<
"================================================================================" << endl;
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const=0
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().
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
int Initialize()
Initialize L and U with values from user matrix A.
int SetParameters(Teuchos::ParameterList ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
virtual int NumCompute() const
Returns the number of calls to Compute().
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
Ifpack_IC(Epetra_RowMatrix *A)
Ifpack_IC constuctor with variable number of indices per row.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
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.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
virtual ~Ifpack_IC()
Ifpack_IC Destructor.
const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
virtual double ComputeTime() const
Returns the time spent in Compute().