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;
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;
#define EPETRA_CHK_ERR(a)
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)
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
void Ifpack_AIJMatrix_dealloc(Ifpack_AIJMatrix *a)
void crout_ict(int n, const Ifpack_AIJMatrix *AL, const double *Adiag, double droptol, int lfil, Ifpack_AIJMatrix *L, double **pdiag)
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.
double LevelOfFill() const
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
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.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
double ComputeFlops_
Contains the number of flops for Compute().
Epetra_Time Time_
Used for timing purposes.
double DropTolerance() const
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.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
double ComputeTime_
Contains the time for all successful calls to Compute().
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
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().
Teuchos::RefCountPtr< Epetra_Vector > D_
const char * Label() const
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
double RelativeThreshold() const
int NumInitialize_
Contains the number of successful calls to Initialize().
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
virtual double ComputeTime() const
Returns the time spent in Compute().
int NumCompute_
Contains the number of successful call to Compute().
double AbsoluteThreshold() const