46 #include "Epetra_ConfigDefs.h"
47 #include "Epetra_Comm.h"
48 #include "Epetra_Map.h"
49 #include "Epetra_RowMatrix.h"
50 #include "Epetra_Vector.h"
51 #include "Epetra_MultiVector.h"
52 #include "Epetra_CrsGraph.h"
53 #include "Epetra_CrsMatrix.h"
54 #include "Teuchos_ParameterList.hpp"
55 #include "Teuchos_RefCountPtr.hpp"
57 using Teuchos::RefCountPtr;
60 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
61 # include "Teuchos_TimeMonitor.hpp"
67 Comm_(Matrix_in->Comm()),
75 IsInitialized_(
false),
82 ApplyInverseTime_(0.0),
84 ApplyInverseFlops_(0.0),
108 sprintf(
Label_,
"IFPACK ILU (fill=%d, relax=%f, athr=%f, rthr=%f)",
118 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
125 if ((
L_.get() == 0) || (
U_.get() == 0) || (
D_.get() == 0))
138 int NumIn, NumL, NumU;
140 int NumNonzeroDiags = 0;
142 std::vector<int> InI(MaxNumEntries);
143 std::vector<int> LI(MaxNumEntries);
144 std::vector<int> UI(MaxNumEntries);
145 std::vector<double> InV(MaxNumEntries);
146 std::vector<double> LV(MaxNumEntries);
147 std::vector<double> UV(MaxNumEntries);
149 bool ReplaceValues = (
L_->StaticGraph() ||
L_->IndicesAreLocal());
172 for (j=0; j< NumIn; j++) {
196 if (DiagFound) NumNonzeroDiags++;
201 (
L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));
210 (
U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));
219 if (!ReplaceValues) {
229 int TotalNonzeroDiags = 0;
232 if (NumNonzeroDiags !=
NumMyRows()) ierr = 1;
242 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
254 if (CrsMatrix == 0) {
259 int size =
A_->MaxNumEntries();
264 std::vector<double> Values(size);
266 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
267 if(
A_->RowMatrixRowMap().GlobalIndicesInt()) {
268 std::vector<int> Indices(size);
271 for (
int i = 0 ; i <
A_->NumMyRows() ; ++i) {
273 int GlobalRow =
A_->RowMatrixRowMap().GID(i);
275 &Values[0], &Indices[0]));
277 for (
int j = 0 ; j < NumEntries ; ++j) {
278 Indices[j] =
A_->RowMatrixColMap().GID(Indices[j]);
286 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
287 if(
A_->RowMatrixRowMap().GlobalIndicesLongLong()) {
288 std::vector<int> Indices_local(size);
289 std::vector<long long> Indices(size);
292 for (
int i = 0 ; i <
A_->NumMyRows() ; ++i) {
294 long long GlobalRow =
A_->RowMatrixRowMap().GID64(i);
296 &Values[0], &Indices_local[0]));
298 for (
int j = 0 ; j < NumEntries ; ++j) {
299 Indices[j] =
A_->RowMatrixColMap().GID64(Indices_local[j]);
307 throw "Ifpack_ILU::Initialize: GlobalIndices type unknown";
310 A_->RowMatrixRowMap()));
337 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
354 double MaxDiagonalValue = 1.0/MinDiagonalValue;
360 int NumIn, NumL, NumU;
363 int MaxNumEntries =
L_->MaxNumEntries() +
U_->MaxNumEntries() + 1;
365 std::vector<int> InI(MaxNumEntries+1);
366 std::vector<double> InV(MaxNumEntries+1);
370 ierr =
D_->ExtractView(&DV);
372 #ifdef IFPACK_FLOPCOUNTERS
373 int current_madds = 0;
384 for (j = 0; j <
NumMyCols(); ++j) colflag[j] = - 1;
390 NumIn = MaxNumEntries;
398 IFPACK_CHK_ERR(
U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));
404 for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
406 double diagmod = 0.0;
408 for (
int jj=0; jj<NumL; jj++) {
410 double multiplier = InV[jj];
417 for (k=0; k<NumUU; k++) {
418 int kk = colflag[UUI[k]];
420 InV[kk] -= multiplier*UUV[k];
421 #ifdef IFPACK_FLOPCOUNTERS
428 for (k=0; k<NumUU; k++) {
429 int kk = colflag[UUI[k]];
430 if (kk>-1) InV[kk] -= multiplier*UUV[k];
431 else diagmod -= multiplier*UUV[k];
432 #ifdef IFPACK_FLOPCOUNTERS
449 if (fabs(DV[i]) > MaxDiagonalValue) {
450 if (DV[i] < 0) DV[i] = - MinDiagonalValue;
451 else DV[i] = MinDiagonalValue;
456 for (j=0; j<NumU; j++) UV[j] *= DV[i];
463 for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
468 if (!
L_->LowerTriangular())
470 if (!
U_->UpperTriangular())
473 #ifdef IFPACK_FLOPCOUNTERS
476 double current_flops = 2 * current_madds;
477 double total_flops = 0;
482 total_flops += (double)
L_->NumGlobalNonzeros();
483 total_flops += (double)
D_->GlobalLength();
484 if (
RelaxValue_!=0.0) total_flops += 2 * (double)
D_->GlobalLength();
504 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
511 bool UnitDiagonal =
true;
539 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
579 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
586 if (X.NumVectors() != Y.NumVectors())
593 Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
594 if (X.Pointers()[0] == Y.Pointers()[0])
602 #ifdef IFPACK_FLOPCOUNTERS
604 (
L_->NumGlobalNonzeros() +
U_->NumGlobalNonzeros());
616 const int MaxIters,
const double Tol,
620 #ifdef IFPACK_TEUCHOS_TIME_MONITOR
638 if (!
Comm().MyPID()) {
640 os <<
"================================================================================" << endl;
641 os <<
"Ifpack_ILU: " <<
Label() << endl << endl;
645 os <<
"Relax value = " <<
RelaxValue() << endl;
646 os <<
"Condition number estimate = " <<
Condest() << endl;
647 os <<
"Global number of rows = " <<
A_->NumGlobalRows64() << endl;
649 os <<
"Number of rows of L, D, U = " <<
L_->NumGlobalRows64() << endl;
651 os <<
"nonzeros / rows = "
655 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
656 os <<
"----- ------- -------------- ------------ --------" << endl;
659 <<
" 0.0 0.0" << endl;
660 os <<
"Compute() " << std::setw(5) <<
NumCompute()
666 os <<
" " << std::setw(15) << 0.0 << endl;
673 os <<
" " << std::setw(15) << 0.0 << endl;
674 os <<
"================================================================================" << endl;
char Label_[160]
Label of this object.
Teuchos::RefCountPtr< Epetra_CrsGraph > CrsGraph_
double InitializeTime_
Contains the time for all successful calls to Initialize().
const Ifpack_IlukGraph & Graph() const
Returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
double Athresh_
absolute threshold
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y...
Teuchos::RefCountPtr< Epetra_Vector > D_
Diagonal of factors.
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
int LevelOfFill() const
Returns the level of fill.
T & get(ParameterList &l, const std::string &name)
double ElapsedTime(void) const
int Multiply(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
bool UseTranspose() const
Returns the current UseTranspose setting.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
virtual const Epetra_Map & OperatorDomainMap() const =0
Teuchos::RefCountPtr< Ifpack_IlukGraph > Graph_
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
const Epetra_Map * U_DomainMap_
bool IsComputed_
If true, the preconditioner has been successfully computed.
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the Epetra_RowMatrix to factorize.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
const Epetra_Map * L_RangeMap_
virtual double ComputeTime() const
Returns the time spent in Compute().
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
virtual int MaxNumEntries() const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int NumCompute_
Contains the number of successful call to Compute().
int Compute()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
int SetParameters(Teuchos::ParameterList ¶meterlist)
Set parameters using a Teuchos::ParameterList object.
int NumMyCols() const
Returns the number of local matrix columns.
Epetra_Time Time_
Used for timing issues.
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
Contains the L factors.
double AbsoluteThreshold() const
Get absolute threshold value.
int Initialize()
Initialize the preconditioner, does not touch matrix values.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
virtual double InitializeTime() const
Returns the time spent in Initialize().
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
Contains the U factors.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
const double Epetra_MinDouble
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
double Condest() const
Returns the computed estimated condition number, or -1.0 if not computed.
virtual std::ostream & Print(std::ostream &os) const
Prints on stream basic information about this object.
Ifpack_ILU(Epetra_RowMatrix *A)
Constructor.
int NumMyRows() const
Returns the number of local matrix rows.
long long NumGlobalNonzeros64() const
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
int NumInitialize_
Contains the number of successful calls to Initialize().
void Destroy()
Destroys all internal data.
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
const Epetra_CrsGraph & Graph() const
double ComputeFlops_
Contains the number of flops for Compute().
const char * Label() const
Returns a character string describing the operator.
#define IFPACK_CHK_ERR(ifpack_err)
double RelativeThreshold() const
Get relative threshold value.
void ResetStartTime(void)
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual int NumInitialize() const
Returns the number of calls to Initialize().
double ComputeTime_
Contains the time for all successful calls to Compute().
int Solve(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_ILU forward/back solve on a Epetra_MultiVector X in Y...
double Rthresh_
relative threshold
double RelaxValue_
Relaxation value.
double RelaxValue() const
Get ILU(k) relaxation parameter.
int LevelOfFill_
Level of fill.
double Condest_
condition number estimate