42 #include "Ifpack_CrsRiluk.h" 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> 52 #include <ifp_parameters.h> 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_)
99 if (IlukRowMap_!=Teuchos::null) IlukRowMap_ = Teuchos::rcp(
new Epetra_Map(*IlukRowMap_) );
100 if (IlukDomainMap_!=Teuchos::null) IlukDomainMap_ = Teuchos::rcp(
new Epetra_Map(*IlukDomainMap_) );
101 if (IlukRangeMap_!=Teuchos::null) IlukRangeMap_ = Teuchos::rcp(
new Epetra_Map(*IlukRangeMap_) );
107 ValuesInitialized_ =
false;
112 int Ifpack_CrsRiluk::AllocateCrs() {
118 L_Graph_ = Teuchos::null;
119 U_Graph_ = Teuchos::null;
124 int Ifpack_CrsRiluk::AllocateVbr() {
128 EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.
L_Graph().
RowMap(), &IlukRowMap_));
129 EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.
U_Graph().
DomainMap(), &IlukDomainMap_));
130 EPETRA_CHK_ERR(BlockMap2PointMap(Graph_.
L_Graph().
RangeMap(), &IlukRangeMap_));
133 U_DomainMap_ = IlukDomainMap_;
134 L_RangeMap_ = IlukRangeMap_;
136 if (
Graph().LevelFill()) {
137 L_Graph_ = Teuchos::rcp(
new Epetra_CrsGraph(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
138 U_Graph_ = Teuchos::rcp(
new Epetra_CrsGraph(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
139 EPETRA_CHK_ERR(BlockGraph2PointGraph(Graph_.
L_Graph(), *L_Graph_,
false));
140 EPETRA_CHK_ERR(BlockGraph2PointGraph(Graph_.
U_Graph(), *U_Graph_,
true));
142 L_Graph_->FillComplete(*IlukRowMap_, *IlukRangeMap_);
143 U_Graph_->FillComplete(*IlukDomainMap_, *IlukRowMap_);
151 L_ = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
152 U_ = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *IlukRowMap_, *IlukRowMap_, 0) );
154 L_Graph_ = Teuchos::null;
155 U_Graph_ = Teuchos::null;
163 bool cerr_warning_if_unused)
166 params.double_params[Ifpack::relax_value] = RelaxValue_;
167 params.double_params[Ifpack::absolute_threshold] = Athresh_;
168 params.double_params[Ifpack::relative_threshold] = Rthresh_;
169 params.overlap_mode = OverlapMode_;
171 Ifpack::set_parameters(parameterlist, params, cerr_warning_if_unused);
173 RelaxValue_ = params.double_params[Ifpack::relax_value];
174 Athresh_ = params.double_params[Ifpack::absolute_threshold];
175 Rthresh_ = params.double_params[Ifpack::relative_threshold];
176 OverlapMode_ = params.overlap_mode;
184 UserMatrixIsCrs_ =
true;
186 if (!Allocated()) AllocateCrs();
188 Teuchos::RefCountPtr<Epetra_CrsMatrix> OverlapA = Teuchos::rcp( (
Epetra_CrsMatrix *) &A,
false );
194 EPETRA_CHK_ERR(OverlapA->FillComplete());
198 int MaxNumEntries = OverlapA->MaxNumEntries();
201 U_DomainMap_ = Teuchos::rcp( &(A.
DomainMap()),
false );
202 L_RangeMap_ = Teuchos::rcp( &(A.
RangeMap()),
false );
205 EPETRA_CHK_ERR(InitAllValues(*OverlapA, MaxNumEntries));
211 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES // FIXME LONG LONG 214 UserMatrixIsVbr_ =
true;
216 if (!Allocated()) AllocateVbr();
228 Teuchos::RefCountPtr<Epetra_VbrMatrix> OverlapA = Teuchos::rcp( (
Epetra_VbrMatrix *) &A,
false );
234 EPETRA_CHK_ERR(OverlapA->FillComplete());
240 int MaxNumEntries = OverlapA->MaxNumNonzeros();
244 EPETRA_CHK_ERR(InitAllValues(*OverlapA, MaxNumEntries));
251 int Ifpack_CrsRiluk::InitAllValues(
const Epetra_RowMatrix & OverlapA,
int MaxNumEntries) {
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());
276 EPETRA_CHK_ERR(D_->ExtractView(&DV));
283 EPETRA_CHK_ERR(OverlapA.
ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0]));
291 for (j=0; j< NumIn; j++) {
296 DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_;
299 else if (k < 0) {EPETRA_CHK_ERR(-1);}
315 if (DiagFound) NumNonzeroDiags++;
316 else DV[i] = Athresh_;
320 EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));
323 EPETRA_CHK_ERR(L_->InsertMyValues(i, NumL, &LV[0], &LI[0]));
329 EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));
332 EPETRA_CHK_ERR(U_->InsertMyValues(i, NumU, &UV[0], &UI[0]));
338 if (!ReplaceValues) {
343 EPETRA_CHK_ERR(L_->FillComplete(L_->RowMatrixColMap(), *L_RangeMap_));
344 EPETRA_CHK_ERR(U_->FillComplete(*U_DomainMap_, U_->RowMatrixRowMap()));
349 SetValuesInitialized(
true);
352 int TotalNonzeroDiags = 0;
354 NumMyDiagonals_ = NumNonzeroDiags;
355 if (NumNonzeroDiags !=
NumMyRows()) ierr = 1;
367 SetValuesInitialized(
false);
372 double MinDiagonalValue = Epetra_MinDouble;
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;
408 EPETRA_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, &InV[0], &InI[0]));
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];
431 EPETRA_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI));
433 if (RelaxValue_==0.0) {
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 456 EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI));
461 if (RelaxValue_!=0.0) {
462 DV[i] += RelaxValue_*diagmod;
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];
476 EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, UV, UI));
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;
522 EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));
526 bool UnitDiagonal =
true;
528 #ifdef IFPACK_FLOPCOUNTERS 531 L_->SetFlopCounter(*counter);
532 Y1->SetFlopCounter(*counter);
533 U_->SetFlopCounter(*counter);
539 EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *X1, *Y1));
540 EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0));
541 EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *Y1, *Y1));
542 if (IsOverlapped_) {EPETRA_CHK_ERR(Y.
Export(*Y1,*L_->Exporter(), OverlapMode_));}
545 EPETRA_CHK_ERR(U_->Solve(Upper, Trans, UnitDiagonal, *X1, *Y1));
546 EPETRA_CHK_ERR(Y1->Multiply(1.0, *D_, *Y1, 0.0));
547 EPETRA_CHK_ERR(L_->Solve(Lower, Trans, UnitDiagonal, *Y1, *Y1));
548 if (IsOverlapped_) {EPETRA_CHK_ERR(Y.
Export(*Y1,*U_->Importer(), OverlapMode_));}
561 Teuchos::RefCountPtr<Epetra_MultiVector> X1;
562 Teuchos::RefCountPtr<Epetra_MultiVector> Y1;
563 EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));
565 #ifdef IFPACK_FLOPCOUNTERS 568 L_->SetFlopCounter(*counter);
569 Y1->SetFlopCounter(*counter);
570 U_->SetFlopCounter(*counter);
575 EPETRA_CHK_ERR(U_->Multiply(Trans, *X1, *Y1));
576 EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0));
577 EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0));
579 EPETRA_CHK_ERR(L_->Multiply(Trans, Y1temp, *Y1));
580 EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0));
581 if (IsOverlapped_) {EPETRA_CHK_ERR(Y.
Export(*Y1,*L_->Exporter(), OverlapMode_));}
585 EPETRA_CHK_ERR(L_->Multiply(Trans, *X1, *Y1));
586 EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0));
587 EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0));
589 EPETRA_CHK_ERR(U_->Multiply(Trans, Y1temp, *Y1));
590 EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0));
591 if (IsOverlapped_) {EPETRA_CHK_ERR(Y.
Export(*Y1,*L_->Exporter(), OverlapMode_));}
599 ConditionNumberEstimate = Condest_;
607 EPETRA_CHK_ERR(
Solve(Trans, Ones, OnesResult));
608 EPETRA_CHK_ERR(OnesResult.Abs(OnesResult));
609 EPETRA_CHK_ERR(OnesResult.MaxValue(&ConditionNumberEstimate));
610 Condest_ = ConditionNumberEstimate;
626 std::vector<int> tmpIndices(Length);
628 int BlockRow, BlockOffset, NumEntries;
634 for (
int i=0; i<NumMyRows_tmp; i++) {
636 EPETRA_CHK_ERR(BG.
ExtractMyRowView(BlockRow, NumBlockEntries, BlockIndices));
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++;
662 int jstart = EPETRA_MAX(0,i-RowDim+1);
664 for (
int j = jstart; j < jstop; j++) {*ptr++ = j; NumEntries++;}
675 int Ifpack_CrsRiluk::BlockMap2PointMap(
const Epetra_BlockMap & BlockMap, Teuchos::RefCountPtr<Epetra_Map>* PointMap) {
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 720 (*PointMap) = Teuchos::rcp(
new Epetra_Map(-1, PtNumMyElements, &PtMyGlobalElements_int[0], BlockMap.
IndexBase(), BlockMap.
Comm()) );
723 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 725 (*PointMap) = Teuchos::rcp(
new Epetra_Map(-1LL, PtNumMyElements, &PtMyGlobalElements_LL[0], BlockMap.IndexBase64(), BlockMap.
Comm()) );
728 throw "Ifpack_CrsRiluk::BlockMap2PointMap: GlobalIndices type unknown";
730 if (!BlockMap.
PointSameAs(*(*PointMap))) {EPETRA_CHK_ERR(-1);}
734 int Ifpack_CrsRiluk::GenerateXY(
bool Trans,
736 Teuchos::RefCountPtr<Epetra_MultiVector>* Xout,
737 Teuchos::RefCountPtr<Epetra_MultiVector>* Yout)
const {
746 if (!IsOverlapped_ && UserMatrixIsCrs_)
return(0);
748 if (UserMatrixIsVbr_) {
749 if (VbrX_!=Teuchos::null) {
751 VbrX_ = Teuchos::null;
752 VbrY_ = Teuchos::null;
755 if (VbrX_==Teuchos::null) {
756 VbrX_ = Teuchos::rcp(
new Epetra_MultiVector(View, *U_DomainMap_, (*Xout)->Pointers(), (*Xout)->NumVectors()) );
757 VbrY_ = Teuchos::rcp(
new Epetra_MultiVector(View, *L_RangeMap_, (*Yout)->Pointers(), (*Yout)->NumVectors()) );
760 EPETRA_CHK_ERR(VbrX_->ResetView((*Xout)->Pointers()));
761 EPETRA_CHK_ERR(VbrY_->ResetView((*Yout)->Pointers()));
769 if (OverlapX_!=Teuchos::null) {
770 if (OverlapX_->NumVectors()!=Xin.
NumVectors()) {
771 OverlapX_ = Teuchos::null;
772 OverlapY_ = Teuchos::null;
775 if (OverlapX_==Teuchos::null) {
776 OverlapX_ = Teuchos::rcp(
new Epetra_MultiVector(U_->RowMatrixColMap(), (*Xout)->NumVectors()) );
777 OverlapY_ = Teuchos::rcp(
new Epetra_MultiVector(L_->RowMatrixRowMap(), (*Yout)->NumVectors()) );
780 EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*U_->Importer(),
Insert));
783 EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*L_->Exporter(),
Insert));
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
void UpdateFlops(int Flops_in) const
int NumMyRows() const
Returns the number of local matrix rows.
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
int * ElementSizeList() const
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
bool DistributedGlobal() const
const Epetra_CrsMatrix & L() const
Returns the address of the L factor associated with this factored matrix.
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 ...
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...
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
const Epetra_BlockMap & RowMap() const
const Epetra_Map & RangeMap() const
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
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)
const Epetra_BlockMap & RangeMap() const
const Ifpack_IlukGraph & Graph() const
returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
int FindLocalElementID(int PointID, int &ElementID, int &ElementOffset) const
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
int NumMyCols() const
Returns the number of local matrix columns.
virtual const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
const Epetra_Comm & Comm() const
bool IndicesAreLocal() const
int MaxMyElementSize() const
int * FirstPointInElementList() const
const Epetra_BlockMap & ImportMap() const
const Epetra_BlockMap & DomainMap() const
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.
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
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
const Epetra_Map & DomainMap() const
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.