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;
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;
139 int Ifpack_IC::ComputeSetup()
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++) {
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;
207 Matrix().Comm().MaxAll(&ierr1, &ierr, 1);
208 IFPACK_CHK_ERR(ierr);
219 Time_.ResetStartTime();
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_);
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);
292 ComputeFlops_ += total_flops;
294 ComputeFlops_ += (double) U_->NumGlobalNonzeros();
295 ComputeFlops_ += (double) D_->GlobalLength();
298 ComputeTime_ += Time_.ElapsedTime();
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 337 ApplyInverseFlops_ += 4.0 * U_->NumGlobalNonzeros();
338 ApplyInverseFlops_ += D_->GlobalLength();
342 ApplyInverseTime_ += Time_.ElapsedTime();
350 int Ifpack_IC::Apply(
const Epetra_MultiVector& X,
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)
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;
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
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.
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.
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual int NumInitialize() const
Returns the number of calls to Initialize().
int Initialize()
Initialize L and U with values from user matrix A.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
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...
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
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.
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...
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().