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()
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();
316 if (X.NumVectors() != Y.NumVectors())
322 bool UnitDiagonal =
true;
326 RefCountPtr< const Epetra_MultiVector > Xcopy;
327 if (X.Pointers()[0] == Y.Pointers()[0])
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();
354 if (X.NumVectors() != Y.NumVectors())
360 U_->Multiply(
false, *X1, *Y1);
361 Y1->Update(1.0, *X1, 1.0);
362 Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.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 NumInitialize() const
Returns the number of calls to Initialize().
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
int SetParameters(Teuchos::ParameterList ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
double ElapsedTime(void) const
Ifpack_IC(Epetra_RowMatrix *A)
Ifpack_IC constuctor with variable number of indices per row.
virtual ~Ifpack_IC()
Ifpack_IC Destructor.
const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual double ComputeTime() const
Returns the time spent in Compute().
const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
int Initialize()
Initialize L and U with values from user matrix A.
virtual int MaxNumEntries() const =0
virtual const Epetra_Comm & Comm() const =0
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 NumMyRows() const =0
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized, false otherwise.
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.
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
void ResetStartTime(void)
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.