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;
73 IsInitialized_(
false),
83 ApplyInverseFlops_(0.0)
85 Teuchos::ParameterList List;
114 Lfil_ = List.get(
"fact: ict level-of-fill",
Lfil_);
120 sprintf(
Label_,
"IFPACK IC (fill=%f, drop=%f)",
142 U_ = rcp(
new Epetra_CrsMatrix(Copy,
Matrix().RowMatrixRowMap(),
143 Matrix().RowMatrixRowMap(), 0));
144 D_ = rcp(
new Epetra_Vector(
Matrix().RowMatrixRowMap()));
146 if (
U_.get() == 0 ||
D_.get() == 0)
153 int NumNonzeroDiags = 0;
156 int MaxNumEntries =
Matrix().MaxNumEntries();
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);
168 int NumRows =
Matrix().NumMyRows();
170 for (i=0; i< NumRows; i++) {
172 Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0]);
178 for (j=0; j< NumIn; j++) {
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;
207 Matrix().Comm().MaxAll(&ierr1, &ierr, 1);
219 Time_.ResetStartTime();
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);
239 Aict_ = (
void *) Aict;
245 Lict_ = (
void *) Lict;
256 EPETRA_CHK_ERR(
D_->ExtractView(&DV));
261 int lfil = (
Lfil_ * nz)/n + .5;
269 U_ = rcp(
new Epetra_CrsMatrix(View,
A_->RowMatrixRowMap(),
A_->RowMatrixRowMap(),0));
270 D_ = rcp(
new Epetra_Vector(View,
A_->RowMatrixRowMap(),
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);
310 Epetra_MultiVector& Y)
const 316 if (X.NumVectors() != Y.NumVectors())
319 Time_.ResetStartTime();
322 bool UnitDiagonal =
true;
326 RefCountPtr< const Epetra_MultiVector > Xcopy;
327 if (X.Pointers()[0] == Y.Pointers()[0])
328 Xcopy = rcp(
new Epetra_MultiVector(X) );
330 Xcopy = rcp( &X,
false );
332 U_->Solve(Upper,
true, UnitDiagonal, *Xcopy, Y);
333 Y.Multiply(1.0, *
D_, Y, 0.0);
334 U_->Solve(Upper,
false, UnitDiagonal, Y, Y);
336 #ifdef IFPACK_FLOPCOUNTERS 351 Epetra_MultiVector& Y)
const 354 if (X.NumVectors() != Y.NumVectors())
357 Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
358 Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;
360 U_->Multiply(
false, *X1, *Y1);
361 Y1->Update(1.0, *X1, 1.0);
362 Y1->ReciprocalMultiply(1.0, *
D_, *Y1, 0.0);
363 Epetra_MultiVector Y1temp(*Y1);
364 U_->Multiply(
true, Y1temp, *Y1);
365 Y1->Update(1.0, Y1temp, 1.0);
371 const int MaxIters,
const double Tol,
372 Epetra_RowMatrix* Matrix_in)
389 if (!
Comm().MyPID()) {
391 os <<
"================================================================================" << endl;
392 os <<
"Ifpack_IC: " <<
Label() << endl << 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;
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
double ComputeTime_
Contains the time for all successful calls to Compute().
double DropTolerance() const
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
int SetParameters(Teuchos::ParameterList ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Epetra_Time Time_
Used for timing purposes.
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
virtual double ComputeTime() const
Returns the time spent in Compute().
const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
Ifpack_IC(Epetra_RowMatrix *A)
Ifpack_IC constuctor with variable number of indices per row.
virtual ~Ifpack_IC()
Ifpack_IC Destructor.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
Teuchos::RefCountPtr< Epetra_Vector > D_
const char * Label() const
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual int NumInitialize() const
Returns the number of calls to Initialize().
void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
double AbsoluteThreshold() const
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.
double RelativeThreshold() const
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
double LevelOfFill() const
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
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...
double ComputeFlops_
Contains the number of flops for Compute().
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
int NumCompute_
Contains the number of successful call to Compute().
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
virtual int NumCompute() const
Returns the number of calls to Compute().
#define IFPACK_CHK_ERR(ifpack_err)
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
void crout_ict(int n, const Ifpack_AIJMatrix *AL, const double *Adiag, double droptol, int lfil, Ifpack_AIJMatrix *L, double **pdiag)
int NumInitialize_
Contains the number of successful calls to Initialize().
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().