43 #include "Epetra_ConfigDefs.h" 44 #include "Epetra_Comm.h" 45 #include "Epetra_Map.h" 46 #include "Epetra_CrsGraph.h" 47 #include "Epetra_VbrMatrix.h" 48 #include "Epetra_RowMatrix.h" 49 #include "Epetra_MultiVector.h" 51 #include <Teuchos_ParameterList.hpp> 56 : UserMatrixIsVbr_(
false),
57 UserMatrixIsCrs_(
false),
59 Comm_(Graph_in.Comm()),
63 ValuesInitialized_(
false),
77 : UserMatrixIsVbr_(FactoredMatrix.UserMatrixIsVbr_),
78 UserMatrixIsCrs_(FactoredMatrix.UserMatrixIsCrs_),
79 IsOverlapped_(FactoredMatrix.IsOverlapped_),
80 Graph_(FactoredMatrix.Graph_),
81 IlukRowMap_(FactoredMatrix.IlukRowMap_),
82 IlukDomainMap_(FactoredMatrix.IlukDomainMap_),
83 IlukRangeMap_(FactoredMatrix.IlukRangeMap_),
84 Comm_(FactoredMatrix.Comm_),
85 UseTranspose_(FactoredMatrix.UseTranspose_),
86 NumMyDiagonals_(FactoredMatrix.NumMyDiagonals_),
87 Allocated_(FactoredMatrix.Allocated_),
88 ValuesInitialized_(FactoredMatrix.ValuesInitialized_),
89 Factored_(FactoredMatrix.Factored_),
90 RelaxValue_(FactoredMatrix.RelaxValue_),
91 Athresh_(FactoredMatrix.Athresh_),
92 Rthresh_(FactoredMatrix.Rthresh_),
93 Condest_(FactoredMatrix.Condest_),
94 OverlapMode_(FactoredMatrix.OverlapMode_)
136 if (
Graph().LevelFill()) {
163 bool cerr_warning_if_unused)
198 int MaxNumEntries = OverlapA->MaxNumEntries();
211 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME LONG LONG 240 int MaxNumEntries = OverlapA->MaxNumNonzeros();
255 int NumIn, NumL, NumU;
257 int NumNonzeroDiags = 0;
260 std::vector<int> InI(MaxNumEntries);
261 std::vector<int> LI(MaxNumEntries);
262 std::vector<int> UI(MaxNumEntries);
263 std::vector<double> InV(MaxNumEntries);
264 std::vector<double> LV(MaxNumEntries);
265 std::vector<double> UV(MaxNumEntries);
267 bool ReplaceValues = (
L_->StaticGraph() ||
L_->IndicesAreLocal());
291 for (j=0; j< NumIn; j++) {
315 if (DiagFound) NumNonzeroDiags++;
338 if (!ReplaceValues) {
352 int TotalNonzeroDiags = 0;
355 if (NumNonzeroDiags !=
NumMyRows()) ierr = 1;
373 double MaxDiagonalValue = 1.0/MinDiagonalValue;
377 int * LI=0, * UI = 0;
378 double * LV=0, * UV = 0;
379 int NumIn, NumL, NumU;
382 int MaxNumEntries =
L_->MaxNumEntries() +
U_->MaxNumEntries() + 1;
384 std::vector<int> InI(MaxNumEntries);
385 std::vector<double> InV(MaxNumEntries);
389 ierr =
D_->ExtractView(&DV);
391 #ifdef IFPACK_FLOPCOUNTERS 392 int current_madds = 0;
401 for (j=0; j<
NumMyCols(); j++) colflag[j] = - 1;
407 NumIn = MaxNumEntries;
415 EPETRA_CHK_ERR(
U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));
421 for (j=0; j<NumIn; j++) colflag[InI[j]] = j;
423 double diagmod = 0.0;
425 for (
int jj=0; jj<NumL; jj++) {
427 double multiplier = InV[jj];
434 for (k=0; k<NumUU; k++) {
435 int kk = colflag[UUI[k]];
437 InV[kk] -= multiplier*UUV[k];
438 #ifdef IFPACK_FLOPCOUNTERS 445 for (k=0; k<NumUU; k++) {
446 int kk = colflag[UUI[k]];
447 if (kk>-1) InV[kk] -= multiplier*UUV[k];
448 else diagmod -= multiplier*UUV[k];
449 #ifdef IFPACK_FLOPCOUNTERS 466 if (fabs(DV[i]) > MaxDiagonalValue) {
467 if (DV[i] < 0) DV[i] = - MinDiagonalValue;
468 else DV[i] = MinDiagonalValue;
473 for (j=0; j<NumU; j++) UV[j] *= DV[i];
480 for (j=0; j<NumIn; j++) colflag[InI[j]] = -1;
485 if( !
L_->LowerTriangular() )
487 if( !
U_->UpperTriangular() )
490 #ifdef IFPACK_FLOPCOUNTERS 493 double current_flops = 2 * current_madds;
494 double total_flops = 0;
499 total_flops += (double)
L_->NumGlobalNonzeros();
500 total_flops += (double)
D_->GlobalLength();
501 if (
RelaxValue_!=0.0) total_flops += 2 * (double)
D_->GlobalLength();
520 Teuchos::RefCountPtr<Epetra_MultiVector> X1;
521 Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
526 bool UnitDiagonal =
true;
528 #ifdef IFPACK_FLOPCOUNTERS 531 L_->SetFlopCounter(*counter);
532 Y1->SetFlopCounter(*counter);
533 U_->SetFlopCounter(*counter);
561 Teuchos::RefCountPtr<Epetra_MultiVector> X1;
562 Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
565 #ifdef IFPACK_FLOPCOUNTERS 568 L_->SetFlopCounter(*counter);
569 Y1->SetFlopCounter(*counter);
570 U_->SetFlopCounter(*counter);
626 std::vector<int> tmpIndices(Length);
628 int BlockRow, BlockOffset, NumEntries;
634 for (
int i=0; i<NumMyRows_tmp; i++) {
638 int * ptr = &tmpIndices[0];
647 int jstop =
EPETRA_MIN(NumMyRows_tmp,i+RowDim-BlockOffset);
648 for (
int j= jstart; j< jstop; j++) {*ptr++ = j; NumEntries++;}
651 for (
int j=0; j<NumBlockEntries; j++) {
652 int ColDim = ColElementSizeList[BlockIndices[j]];
653 NumEntries += ColDim;
654 assert(NumEntries<=Length);
655 int Index = ColFirstPointInElementList[BlockIndices[j]];
656 for (
int k=0; k < ColDim; k++) *ptr++ = Index++;
664 for (
int j = jstart; j < jstop; j++) {*ptr++ = j; NumEntries++;}
684 std::vector<int> PtMyGlobalElements_int;
685 std::vector<long long> PtMyGlobalElements_LL;
690 if (PtNumMyElements>0) {
691 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 693 PtMyGlobalElements_int.resize(PtNumMyElements);
694 for (
int i=0; i<NumMyElements; i++) {
695 int StartID = BlockMap.
GID(i)*MaxElementSize;
697 for (
int j=0; j<ElementSize; j++) PtMyGlobalElements_int[curID++] = StartID+j;
702 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 704 PtMyGlobalElements_LL.resize(PtNumMyElements);
705 for (
int i=0; i<NumMyElements; i++) {
706 long long StartID = BlockMap.
GID64(i)*MaxElementSize;
708 for (
int j=0; j<ElementSize; j++) PtMyGlobalElements_LL[curID++] = StartID+j;
713 throw "Ifpack_CrsRiluk::BlockMap2PointMap: GlobalIndices type unknown";
716 assert(curID==PtNumMyElements);
718 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 723 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 728 throw "Ifpack_CrsRiluk::BlockMap2PointMap: GlobalIndices type unknown";
736 Teuchos::RefCountPtr<Epetra_MultiVector>* Xout,
737 Teuchos::RefCountPtr<Epetra_MultiVector>* Yout)
const {
749 if (
VbrX_!=Teuchos::null) {
751 VbrX_ = Teuchos::null;
752 VbrY_ = Teuchos::null;
755 if (
VbrX_==Teuchos::null) {
802 int LevelFill =
A.Graph().LevelFill();
803 int LevelOverlap =
A.Graph().LevelOverlap();
810 os <<
" Level of Fill = "; os << LevelFill;
813 os <<
" Level of Overlap = "; os << LevelOverlap;
817 os <<
" Lower Triangle = ";
823 os <<
" Inverse of Diagonal = ";
829 os <<
" Upper Triangle = ";
const Epetra_Import * Importer() const
Teuchos::RefCountPtr< const Epetra_Map > L_RangeMap_
void UpdateFlops(int Flops_in) const
int NumMyRows() const
Returns the number of local matrix rows.
int BlockMap2PointMap(const Epetra_BlockMap &BlockMap, Teuchos::RefCountPtr< Epetra_Map > *PointMap)
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
void SetFactored(bool Flag)
double double_params[FIRST_INT_PARAM]
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
bool DistributedGlobal() const
Teuchos::RefCountPtr< Epetra_CrsGraph > U_Graph_
Teuchos::RefCountPtr< const Epetra_Map > U_DomainMap_
const Epetra_CrsMatrix & L() const
Returns the address of the L factor associated with this factored matrix.
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapY_
virtual Epetra_CrsGraph & L_Graph()
Returns the graph of lower triangle of the ILU(k) graph as a Epetra_CrsGraph.
Epetra_Flops * GetFlopCounter() const
int PutScalar(double ScalarConstant)
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
int Solve(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_CrsRiluk forward/back solve on a Epetra_MultiVector X in Y (works for ...
#define EPETRA_CHK_ERR(a)
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
int * FirstPointInElementList() const
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
int BlockGraph2PointGraph(const Epetra_CrsGraph &BG, Epetra_CrsGraph &PG, bool Upper)
const Epetra_BlockMap & RowMap() const
int SetAllocated(bool Flag)
int * ElementSizeList() const
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
void SetValuesInitialized(bool Flag)
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
std::ostream & operator<<(std::ostream &os, const Ifpack_CrsRiluk &A)
<< operator will work for Ifpack_CrsRiluk.
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Teuchos::RefCountPtr< Epetra_CrsGraph > L_Graph_
Teuchos::RefCountPtr< Epetra_Map > IlukDomainMap_
Epetra_CombineMode overlap_mode
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
long long GID64(int LID) const
const Epetra_BlockMap & RangeMap() const
const Ifpack_IlukGraph & Graph() const
returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
int FindLocalElementID(int PointID, int &ElementID, int &ElementOffset) const
int InitAllValues(const Epetra_RowMatrix &A, int MaxNumEntries)
int SetParameters(const Teuchos::ParameterList ¶meterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
int Multiply(bool Trans, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of multiplying U, D and L in that order on an Epetra_MultiVector X in Y...
bool GlobalIndicesInt() const
int NumMyElements() const
const Ifpack_IlukGraph & Graph_
int NumMyCols() const
Returns the number of local matrix columns.
Teuchos::RefCountPtr< Epetra_Map > IlukRowMap_
virtual const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapX_
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
const Epetra_Comm & Comm() const
bool IndicesAreLocal() const
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
int MaxMyElementSize() const
const double Epetra_MinDouble
const Epetra_BlockMap & ImportMap() const
long long IndexBase64() const
const Epetra_BlockMap & DomainMap() const
Teuchos::RefCountPtr< Epetra_MultiVector > VbrY_
int MaxNumIndices() const
int MaxElementSize() const
const Epetra_CrsMatrix & U() const
Returns the address of the L factor associated with this factored matrix.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
int GenerateXY(bool Trans, const Epetra_MultiVector &Xin, const Epetra_MultiVector &Yin, Teuchos::RefCountPtr< Epetra_MultiVector > *Xout, Teuchos::RefCountPtr< Epetra_MultiVector > *Yout) const
void set_parameters(const Teuchos::ParameterList ¶meterlist, param_struct ¶ms, bool cerr_warning_if_unused)
Teuchos::RefCountPtr< Epetra_Vector > D_
Teuchos::RefCountPtr< Epetra_MultiVector > VbrX_
virtual ~Ifpack_CrsRiluk()
Ifpack_CrsRiluk Destructor.
bool GlobalIndicesLongLong() const
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
bool PointSameAs(const Epetra_BlockMap &Map) const
Epetra_CombineMode OverlapMode_
Teuchos::RefCountPtr< Epetra_Map > IlukRangeMap_
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
Ifpack_CrsRiluk(const Ifpack_IlukGraph &Graph_in)
Ifpack_CrsRiluk constuctor with variable number of indices per row.
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.