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();
503 UpdateFlops(total_flops);
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 {
741 if (Xin.NumVectors()!=Yin.NumVectors()) EPETRA_CHK_ERR(-1);
746 if (!IsOverlapped_ && UserMatrixIsCrs_)
return(0);
748 if (UserMatrixIsVbr_) {
749 if (VbrX_!=Teuchos::null) {
750 if (VbrX_->NumVectors()!=Xin.NumVectors()) {
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 = ";
bool PointSameAs(const Epetra_BlockMap &Map) const
const Epetra_BlockMap & RangeMap() const
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 DistributedGlobal() const
const Epetra_Map & RangeMap() const
Ifpack_CrsRiluk: A class for constructing and using an incomplete lower/upper (ILU) factorization of ...
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
bool GlobalIndicesLongLong() 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.
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
const Epetra_Import * Importer() const
Ifpack_IlukGraph: A class for constructing level filled graphs for use with ILU(k) class precondition...
const Epetra_BlockMap & ImportMap() const
bool GlobalIndicesInt() const
const Epetra_BlockMap & DomainMap() const
int * FirstPointInElementList() const
virtual const Epetra_BlockMap & DomainMap() const
Returns the Epetra_BlockMap object associated with the domain of this matrix operator.
virtual int LevelFill() const
Returns the level of fill used to construct this graph.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
int SetParameters(const Teuchos::ParameterList ¶meterlist, bool cerr_warning_if_unused=false)
Set parameters using a Teuchos::ParameterList object.
int NumMyCols() const
Returns the number of local matrix columns.
int InitValues(const Epetra_CrsMatrix &A)
Initialize L and U with values from user matrix A.
const Epetra_CrsMatrix & U() const
Returns the address of the L factor associated with this factored matrix.
int * ElementSizeList() const
const Epetra_BlockMap & RowMap() const
int FindLocalElementID(int PointID, int &ElementID, int &ElementOffset) const
virtual int LevelOverlap() const
Returns the level of overlap used to construct this graph.
const Epetra_Comm & Comm() const
virtual Epetra_Import * OverlapImporter() const
Returns the importer used to create the overlapped graph.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
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 ...
bool Factored() const
If factor is completed, this query returns true, otherwise it returns false.
int NumMyRows() const
Returns the number of local matrix rows.
bool ValuesInitialized() const
If values have been initialized, this query returns true, otherwise it returns false.
const Epetra_Vector & D() const
Returns the address of the D factor associated with this factored matrix.
int MaxNumIndices() const
const Epetra_Map & DomainMap() const
int MaxElementSize() const
int MaxMyElementSize() const
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
virtual ~Ifpack_CrsRiluk()
Ifpack_CrsRiluk Destructor.
int Factor()
Compute ILU factors L and U using the specified graph, diagonal perturbation thresholds and relaxatio...
const Ifpack_IlukGraph & Graph() const
returns the address of the Ifpack_IlukGraph associated with this factored matrix. ...
Ifpack_CrsRiluk(const Ifpack_IlukGraph &Graph_in)
Ifpack_CrsRiluk constuctor with variable number of indices per row.
bool IndicesAreLocal() const