44 #include "Epetra_Comm.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_CrsGraph.h"
47 #include "Epetra_CrsMatrix.h"
48 #include "Epetra_Vector.h"
49 #include "Epetra_MultiVector.h"
51 #include <Teuchos_ParameterList.hpp>
60 ValuesInitialized_(
false),
75 : A_(FactoredMatrix.A_),
76 Graph_(FactoredMatrix.Graph_),
77 UseTranspose_(FactoredMatrix.UseTranspose_),
78 Allocated_(FactoredMatrix.Allocated_),
79 ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
80 Factored_(FactoredMatrix.Factored_),
81 RelaxValue_(FactoredMatrix.RelaxValue_),
82 Condest_(FactoredMatrix.Condest_),
83 Athresh_(FactoredMatrix.Athresh_),
84 Rthresh_(FactoredMatrix.Rthresh_),
87 OverlapMode_(FactoredMatrix.OverlapMode_)
122 bool cerr_warning_if_unused)
147 int * InI=0, * LI=0, * UI = 0;
148 double * InV=0, * LV=0, * UV = 0;
149 int NumIn, NumL, NumU;
151 int NumNonzeroDiags = 0;
164 InI =
new int[MaxNumEntries];
165 UI =
new int[MaxNumEntries];
166 InV =
new double[MaxNumEntries];
167 UV =
new double[MaxNumEntries];
185 for (j=0; j< NumIn; j++) {
193 else if (k < 0)
return(-1);
203 if (DiagFound) NumNonzeroDiags++;
223 int TotalNonzeroDiags = 0;
246 double MaxDiagonalValue = 1.0/MinDiagonalValue;
257 int * InI =
new int[MaxNumEntries];
258 double * InV =
new double[MaxNumEntries];
264 #ifdef IFPACK_FLOPCOUNTERS
265 int current_madds = 0;
274 for (j=0; j<
NumMyCols(); j++) colflag[j] = - 1;
280 NumIn = MaxNumEntries;
281 IFPACK_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, InV, InI)==0);
294 for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
296 double diagmod = 0.0;
298 for (
int jj=0; jj<NumL; jj++) {
300 double multiplier = InV[jj];
307 for (k=0; k<NumUU; k++) {
308 int kk = colflag[UUI[k]];
310 InV[kk] -= multiplier*UUV[k];
311 #ifdef IFPACK_FLOPCOUNTERS
319 for (k=0; k<NumUU; k++) {
320 int kk = colflag[UUI[k]];
321 if (kk>-1) InV[kk] -= multiplier*UUV[k];
322 else diagmod -= multiplier*UUV[k];
323 #ifdef IFPACK_FLOPCOUNTERS
339 if (fabs(DV[i]) > MaxDiagonalValue) {
340 if (DV[i] < 0) DV[i] = - MinDiagonalValue;
341 else DV[i] = MinDiagonalValue;
346 for (j=0; j<NumU; j++) UV[j] *= DV[i];
353 for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
357 #ifdef IFPACK_FLOPCOUNTERS
360 double current_flops = 2 * current_madds;
361 double total_flops = 0;
366 total_flops += (double) L_->NumGlobalNonzeros();
367 total_flops += (double)
D_->GlobalLength();
368 if (
RelaxValue_!=0.0) total_flops += 2 * (double)
D_->GlobalLength();
370 UpdateFlops(total_flops);
392 bool UnitDiagonal =
true;
407 #ifdef IFPACK_FLOPCOUNTERS
410 L_->SetFlopCounter(*counter);
411 Y1->SetFlopCounter(*counter);
412 U_->SetFlopCounter(*counter);
418 L_->
Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
419 Y1->Multiply(1.0, *
D_, *Y1, 0.0);
420 U_->
Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1);
424 U_->
Solve(Upper, Trans, UnitDiagonal, *X1, *Y1);
425 Y1->Multiply(1.0, *
D_, *Y1, 0.0);
426 L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
448 bool UnitDiagonal =
true;
456 if (
OverlapX_->NumVectors()!=X.NumVectors()) {
470 #ifdef IFPACK_FLOPCOUNTERS
473 L_->SetFlopCounter(*counter);
474 Y1->SetFlopCounter(*counter);
475 U_->SetFlopCounter(*counter);
481 L_->
Solve(Lower, Trans, UnitDiagonal, *X1, *Y1);
482 Y1->Multiply(1.0, *
D_, *Y1, 0.0);
483 U_->
Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1);
487 U_->
Solve(Upper, Trans, UnitDiagonal, *X1, *Y1);
488 Y1->Multiply(1.0, *
D_, *Y1, 0.0);
489 L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1);
509 bool UnitDiagonal =
true;
517 if (
OverlapX_->NumVectors()!=X.NumVectors()) {
531 #ifdef IFPACK_FLOPCOUNTERS
534 L_->SetFlopCounter(*counter);
535 Y1->SetFlopCounter(*counter);
536 U_->SetFlopCounter(*counter);
543 Y1->Update(1.0, *X1, 1.0);
544 Y1->ReciprocalMultiply(1.0, *
D_, *Y1, 0.0);
547 Y1->Update(1.0, Y1temp, 1.0);
551 Y1->Update(1.0, *X1, 1.0);
552 Y1->ReciprocalMultiply(1.0, *
D_, *Y1, 0.0);
554 L_->Multiply(Trans, Y1temp, *Y1);
555 Y1->Update(1.0, Y1temp, 1.0);
597 os <<
" Level of Fill = "; os << LevelFill;
600 os <<
" Level of Overlap = "; os << LevelOverlap;
604 os <<
" Lower Triangle = ";
610 os <<
" Inverse of Diagonal = ";
616 os <<
" Upper Triangle = ";
bool DistributedGlobal() const
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
const Epetra_CrsMatrix & U() const
Returns the address of the U factor associated with this factored matrix.
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
double double_params[FIRST_INT_PARAM]
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Epetra_MultiVector * OverlapY_
int SetAllocated(bool Flag)
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
Epetra_MultiVector * OverlapX_
#define EPETRA_CHK_ERR(a)
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int NumMyRows() const
Returns the number of local matrix rows.
const Epetra_BlockMap & DomainMap() const
const Ifpack_IlukGraph & Graph() const
Returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
int InitValues()
Initialize L and U with values from user matrix A.
int FillComplete(bool OptimizeDataStorage=true)
const Epetra_CrsMatrix & A_
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
const Epetra_Map & RowMap() const
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
const Ifpack_IlukGraph & Graph_
int SetParameters(const Teuchos::ParameterList ¶meterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
int NumGlobalRows() const
Returns the number of global matrix rows.
Epetra_CombineMode overlap_mode
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
int MaxNumEntries() const
virtual int NumGlobalDiagonals() const
Returns the number of diagonal entries found in the global input graph.
const Epetra_BlockMap & RowMap() const
adjacency_list< vecS, vecS, undirectedS, no_property, property< edge_weight_t, double > > Graph
Ifpack_CrsRick: A class for constructing and using an incomplete Cholesky (IC) factorization of a giv...
Epetra_CombineMode OverlapMode_
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
std::ostream & operator<<(std::ostream &os, const Ifpack_Container &obj)
const Epetra_Comm & Comm() const
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.
Ifpack_CrsRick(const Epetra_CrsMatrix &A, const Ifpack_IlukGraph &Graph)
Ifpack_CrsRick constuctor with variable number of indices per row.
const double Epetra_MinDouble
virtual int NumMyDiagonals() const
Returns the number of diagonal entries found in the local input graph.
virtual ~Ifpack_CrsRick()
Ifpack_CrsRick Destructor.
void SetFactored(bool Flag)
int NumMyCols() const
Returns the number of local matrix columns.
int Solve(bool Trans, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Ifpack_CrsRick forward/back solve on a Epetra_Vector x in y...
void SetValuesInitialized(bool Flag)
int Factor()
Compute IC factor L using the specified graph, diagonal perturbation thresholds and relaxation parame...
int Multiply(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of multiplying U, D and U^T in that order on an Epetra_MultiVector X in Y...
void set_parameters(const Teuchos::ParameterList ¶meterlist, param_struct ¶ms, bool cerr_warning_if_unused)
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
#define IFPACK_CHK_ERR(ifpack_err)
int ExtractView(double **V) const
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local IC set of factors...
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.