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();
503 UpdateFlops(total_flops);
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 {
750 if (
VbrX_->NumVectors()!=Xin.NumVectors()) {
770 if (
OverlapX_->NumVectors()!=Xin.NumVectors()) {
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...
int GenerateXY(bool Trans, const Epetra_MultiVector &Xin, const Epetra_MultiVector &Yin, Teuchos::RefCountPtr< Epetra_MultiVector > *Xout, Teuchos::RefCountPtr< Epetra_MultiVector > *Yout) const
int * FirstPointInElementList() const
bool DistributedGlobal() const
const Epetra_Map & RangeMap() const
Teuchos::RefCountPtr< const Epetra_Map > L_RangeMap_
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 ...
int Condest(bool Trans, double &ConditionNumberEstimate) const
Returns the maximum over all the condition number estimate for each local ILU set of factors...
void SetFactored(bool Flag)
double double_params[FIRST_INT_PARAM]
bool GlobalIndicesLongLong() const
Teuchos::RefCountPtr< Epetra_CrsGraph > U_Graph_
long long IndexBase64() const
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.
virtual Epetra_CrsGraph * OverlapGraph() const
Returns the the overlapped graph.
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
const Epetra_Import * Importer() const
#define EPETRA_CHK_ERR(a)
int * ElementSizeList() 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 BlockGraph2PointGraph(const Epetra_CrsGraph &BG, Epetra_CrsGraph &PG, bool Upper)
int SetAllocated(bool Flag)
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.
void SetValuesInitialized(bool Flag)
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
int NumMyElements() const
long long GID64(int LID) const
virtual Epetra_CrsGraph & U_Graph()
Returns the graph of upper triangle of the ILU(k) graph as a Epetra_CrsGraph.
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)
Teuchos::RefCountPtr< Epetra_CrsMatrix > U_
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.
const Ifpack_IlukGraph & Graph_
Teuchos::RefCountPtr< Epetra_Map > IlukRowMap_
int NumMyCols() const
Returns the number of local matrix columns.
Teuchos::RefCountPtr< Epetra_MultiVector > OverlapX_
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.
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.
std::ostream & operator<<(std::ostream &os, const Ifpack_Container &obj)
Teuchos::RefCountPtr< Epetra_CrsMatrix > L_
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 ...
const double Epetra_MinDouble
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
Teuchos::RefCountPtr< Epetra_MultiVector > VbrY_
int MaxElementSize() const
int MaxMyElementSize() const
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
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.
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...
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