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)",
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++) {
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;
227 int m, n, nz, Nrhs, ldrhs, ldlhs;
229 double * val, * rhs, * lhs;
232 val, Nrhs, rhs, ldrhs, lhs, ldlhs);
239 Aict_ = (
void *) Aict;
245 Lict_ = (
void *) Lict;
261 int lfil = (
Lfil_ * nz)/n + .5;
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);
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 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,
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.
#define EPETRA_CHK_ERR(a)
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().
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const=0
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 MaxNumEntries() const=0
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual const Epetra_Comm & Comm() const=0
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().
virtual int NumMyRows() const=0
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.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
#define IFPACK_CHK_ERR(ifpack_err)
void ResetStartTime(void)
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
double ** Pointers() const
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)
double ElapsedTime(void) const
int Epetra_Util_ExtractHbData(Epetra_CrsMatrix *A, Epetra_MultiVector *LHS, Epetra_MultiVector *RHS, int &M, int &N, int &nz, int *&ptr, int *&ind, double *&val, int &Nrhs, double *&rhs, int &ldrhs, double *&lhs, int &ldlhs)
int NumInitialize_
Contains the number of successful calls to Initialize().
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
int ReciprocalMultiply(double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)