43 #include "Ifpack_ConfigDefs.h"
44 #include "Ifpack_Preconditioner.h"
45 #include "Ifpack_IKLU.h"
46 #include "Ifpack_Condest.h"
48 #include "Ifpack_HashTable.h"
49 #include "Epetra_SerialComm.h"
50 #include "Epetra_Comm.h"
51 #include "Epetra_Map.h"
52 #include "Epetra_RowMatrix.h"
53 #include "Epetra_CrsMatrix.h"
54 #include "Epetra_Vector.h"
55 #include "Epetra_MultiVector.h"
56 #include "Epetra_Util.h"
57 #include "Teuchos_ParameterList.hpp"
58 #include "Teuchos_RefCountPtr.hpp"
62 using namespace Teuchos;
74 DropTolerance_(1e-12),
75 IsInitialized_(false),
84 ApplyInverseTime_(0.0),
86 ApplyInverseFlops_(0.0),
103 void Ifpack_IKLU::Destroy()
105 IsInitialized_ =
false;
123 LevelOfFill_ = List.get<
double>(
"fact: ilut level-of-fill", LevelOfFill());
124 if (LevelOfFill_ <= 0.0)
127 Athresh_ = List.get<
double>(
"fact: absolute threshold", Athresh_);
128 Rthresh_ = List.get<
double>(
"fact: relative threshold", Rthresh_);
129 Relax_ = List.get<
double>(
"fact: relax value", Relax_);
130 DropTolerance_ = List.get<
double>(
"fact: drop tolerance", DropTolerance_);
142 cerr <<
"Caught an exception while parsing the parameter list" << endl;
143 cerr <<
"This typically means that a parameter was set with the" << endl;
144 cerr <<
"wrong type (for example, int instead of double). " << endl;
145 cerr <<
"please check the documentation for the type required by each parameer." << endl;
165 cout <<
" There are too many processors !!! " << endl;
166 cerr <<
"Ifpack_IKLU can be used with Comm().NumProc() == 1" << endl;
167 cerr <<
"only. This class is a subdomain solver for Ifpack_AdditiveSchwarz," << endl;
168 cerr <<
"and it is currently not meant to be used otherwise." << endl;
180 std::vector<int> RowIndices(Length);
181 std::vector<double> RowValues(Length);
185 csrA_ = csr_spalloc( NumMyRows_, NumMyRows_, NumMyNonzeros_, 1, 0 );
190 for (
int i = 0; i < NumMyRows_; ++i ) {
193 &RowValues[0],&RowIndices[0]));
194 for (
int j = 0 ; j < RowNnz ; ++j) {
195 csrA_->j[count++] = RowIndices[j];
198 csrA_->p[i+1] = csrA_->p[i] + RowNnz;
203 cssS_ = csr_sqr( order, csrA_ );
204 for (
int i = 0; i < NumMyRows_; ++i ) {
205 cout <<
"AMD Perm (from inside KLU) [" << i <<
"] = " << cssS_->q[i] << endl;
209 IsInitialized_ =
true;
220 inline bool operator()(
const double& x,
const double& y)
222 return(IFPACK_ABS(x) > IFPACK_ABS(y));
242 bool distributed = (
Comm().
NumProc() > 1)?
true:
false;
243 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
247 SerialMap_ = rcp(
new Epetra_Map(NumMyRows_, 0, *SerialComm_));
248 assert (SerialComm_.get() != 0);
249 assert (SerialMap_.get() != 0);
252 SerialMap_ = rcp(const_cast<Epetra_Map*>(&A_.
RowMatrixRowMap()),
false);
256 std::vector<int> RowIndices(Length);
257 std::vector<double> RowValues(Length);
261 for (
int i = 0; i < NumMyRows_; ++i ) {
264 &RowValues[0],&RowIndices[0]));
266 if (RowNnz != (csrA_->p[i+1]-csrA_->p[i])) {
267 cout <<
"The number of nonzeros for this row does not math the expected number of nonzeros!!!" << endl;
269 for (
int j = 0 ; j < RowNnz ; ++j) {
271 csrA_->x[count++] = RowValues[j];
278 csrnN_ = csr_lu( &*csrA_, &*cssS_, tol );
281 csr* L_tmp = csrnN_->L;
282 csr* U_tmp = csrnN_->U;
283 std::vector<int> numEntriesL( NumMyRows_ ), numEntriesU( NumMyRows_ );
284 for (
int i=0; i < NumMyRows_; ++i) {
285 numEntriesL[i] = ( L_tmp->p[i+1] - L_tmp->p[i] );
286 numEntriesU[i] = ( U_tmp->p[i+1] - U_tmp->p[i] );
292 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
293 if(SerialMap_->GlobalIndicesInt()) {
294 for (
int i=0; i < NumMyRows_; ++i) {
295 L_->InsertGlobalValues( i, numEntriesL[i], &(L_tmp->x[L_tmp->p[i]]), &(L_tmp->j[L_tmp->p[i]]) );
296 U_->InsertGlobalValues( i, numEntriesU[i], &(U_tmp->x[U_tmp->p[i]]), &(U_tmp->j[U_tmp->p[i]]) );
301 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
302 if(SerialMap_->GlobalIndicesLongLong()) {
304 const int MaxNumEntries_L_U = std::max(L_->MaxNumEntries(), U_->MaxNumEntries());
305 std::vector<long long> entries(MaxNumEntries_L_U);
307 for (
int i=0; i < NumMyRows_; ++i) {
308 std::copy(&(L_tmp->j[L_tmp->p[i]]), &(L_tmp->j[L_tmp->p[i]]) + numEntriesL[i], entries.begin());
309 L_->InsertGlobalValues( i, numEntriesL[i], &(L_tmp->x[L_tmp->p[i]]), &(entries[0]) );
311 std::copy(&(U_tmp->j[U_tmp->p[i]]), &(U_tmp->j[U_tmp->p[i]]) + numEntriesU[i], entries.begin());
312 U_->InsertGlobalValues( i, numEntriesU[i], &(U_tmp->x[U_tmp->p[i]]), &(entries[0]) );
317 throw "Ifpack_IKLU::Compute: GlobalIndices type unknown for SerialMap_";
319 IFPACK_CHK_ERR(L_->FillComplete());
320 IFPACK_CHK_ERR(U_->FillComplete());
322 long long MyNonzeros = L_->NumGlobalNonzeros64() + U_->NumGlobalNonzeros64();
323 Comm().
SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
341 if (X.NumVectors() != Y.NumVectors())
352 std::vector<int> invq( NumMyRows_ );
354 for (
int i=0; i<NumMyRows_; ++i ) {
355 csrnN_->perm[ csrnN_->pinv[i] ] = i;
356 invq[ cssS_->q[i] ] = i;
359 Teuchos::RefCountPtr<Epetra_MultiVector> Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X.Map(),X.NumVectors()),
false );
360 Teuchos::RefCountPtr<Epetra_MultiVector> Ytemp = Teuchos::rcp(
new Epetra_MultiVector(Y.Map(),Y.NumVectors()) );
362 for (
int i=0; i<NumMyRows_; ++i) {
363 for (
int j=0; j<X.NumVectors(); ++j) {
364 Xcopy->ReplaceMyValue( invq[i], j, (*X(j))[i] );
371 IFPACK_CHK_ERR(L_->Solve(
false,
false,
false,*Xcopy,*Ytemp));
372 IFPACK_CHK_ERR(U_->Solve(
true,
false,
false,*Ytemp,*Ytemp));
377 IFPACK_CHK_ERR(U_->Solve(
true,
true,
false,*Xcopy,*Ytemp));
378 IFPACK_CHK_ERR(L_->Solve(
false,
true,
false,*Ytemp,*Ytemp));
382 for (
int i=0; i<NumMyRows_; ++i) {
383 for (
int j=0; j<Y.NumVectors(); ++j) {
384 Y.ReplaceMyValue( csrnN_->perm[i], j, (*(*Ytemp)(j))[i] );
389 #ifdef IFPACK_FLOPCOUNTERS
390 ApplyInverseFlops_ += X.NumVectors() * 2 * GlobalNonzeros_;
408 const int MaxIters,
const double Tol,
415 if (Condest_ == -1.0)
416 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
427 if (!
Comm().MyPID()) {
429 os <<
"================================================================================" << endl;
430 os <<
"Ifpack_IKLU: " <<
Label() << endl << endl;
431 os <<
"Level-of-fill = " << LevelOfFill() << endl;
434 os <<
"Relax value = " <<
RelaxValue() << endl;
435 os <<
"Condition number estimate = " <<
Condest() << endl;
436 os <<
"Global number of rows = " << A_.NumGlobalRows64() << endl;
438 os <<
"Number of nonzeros in A = " << A_.NumGlobalNonzeros64() << endl;
439 os <<
"Number of nonzeros in L + U = " << NumGlobalNonzeros64()
440 <<
" ( = " << 100.0 * NumGlobalNonzeros64() / A_.NumGlobalNonzeros64()
441 <<
" % of A)" << endl;
442 os <<
"nonzeros / rows = "
443 << 1.0 * NumGlobalNonzeros64() / U_->NumGlobalRows64() << endl;
446 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
447 os <<
"----- ------- -------------- ------------ --------" << endl;
450 <<
" 0.0 0.0" << endl;
451 os <<
"Compute() " << std::setw(5) <<
NumCompute()
457 os <<
" " << std::setw(15) << 0.0 << endl;
464 os <<
" " << std::setw(15) << 0.0 << endl;
465 os <<
"================================================================================" << endl;
Ifpack_IKLU(const Epetra_RowMatrix *A)
Ifpack_IKLU constuctor with variable number of indices per row.
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
virtual const Epetra_Map & RowMatrixRowMap() const =0
double ElapsedTime(void) const
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
int SetParameters(Teuchos::ParameterList ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
virtual int MaxNumEntries() const =0
virtual const Epetra_Comm & Comm() const =0
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
const char * Label() const
Returns the label of this object.
bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
virtual int NumMyRows() const =0
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual ~Ifpack_IKLU()
Ifpack_IKLU Destructor.
double DropTolerance() const
Gets the dropping tolerance.
virtual int NumProc() const =0
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_IKLU forward/back solve on a Epetra_MultiVector X in Y...
double RelativeThreshold() const
Get relative threshold value.
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
double RelaxValue() const
Set relative threshold value.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
void ResetStartTime(void)
double AbsoluteThreshold() const
Get absolute threshold value.
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
int Initialize()
Initialize L and U with values from user matrix A.
virtual int NumMyNonzeros() const =0
double Condest() const
Returns the computed estimated condition number, or -1.0 if no computed.
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.