43 #include "Ifpack_ConfigDefs.h" 44 #include "Ifpack_Preconditioner.h" 45 #include "Ifpack_ICT.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" 72 IsInitialized_(false),
81 ApplyInverseTime_(0.0),
83 ApplyInverseFlops_(0.0),
97 void Ifpack_ICT::Destroy()
99 IsInitialized_ =
false;
111 LevelOfFill_ = List.get(
"fact: ict level-of-fill",LevelOfFill_);
112 Athresh_ = List.get(
"fact: absolute threshold", Athresh_);
113 Rthresh_ = List.get(
"fact: relative threshold", Rthresh_);
114 Relax_ = List.get(
"fact: relax value", Relax_);
115 DropTolerance_ = List.get(
"fact: drop tolerance", DropTolerance_);
129 cerr <<
"Caught an exception while parsing the parameter list" << endl;
130 cerr <<
"This typically means that a parameter was set with the" << endl;
131 cerr <<
"wrong type (for example, int instead of double). " << endl;
132 cerr <<
"please check the documentation for the type required by each parameter." << endl;
143 Time_.ResetStartTime();
149 NumMyRows_ =
Matrix().NumMyRows();
152 IsInitialized_ =
true;
154 InitializeTime_ += Time_.ElapsedTime();
160 template<
typename int_type>
161 int Ifpack_ICT::TCompute()
166 Time_.ResetStartTime();
169 NumMyRows_ = A_.NumMyRows();
170 int Length = A_.MaxNumEntries();
171 std::vector<int> RowIndices(Length);
172 std::vector<double> RowValues(Length);
174 bool distributed = (
Comm().NumProc() > 1)?
true:
false;
176 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES) 179 SerialComm_ = Teuchos::rcp(
new Epetra_SerialComm);
180 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) 181 if(A_.RowMatrixRowMap().GlobalIndicesInt())
182 SerialMap_ = Teuchos::rcp(
new Epetra_Map(NumMyRows_, 0, *SerialComm_));
185 #if !defined(EPETRA_NO_64BIT_GLOBAL_INDICES) 186 if(A_.RowMatrixRowMap().GlobalIndicesLongLong())
187 SerialMap_ = Teuchos::rcp(
new Epetra_Map((
long long) NumMyRows_, (
long long) 0, *SerialComm_));
190 throw "Ifpack_ICT::TCompute: Global indices unknown.";
191 assert (SerialComm_.get() != 0);
192 assert (SerialMap_.get() != 0);
195 SerialMap_ = Teuchos::rcp(const_cast<Epetra_Map*>(&A_.RowMatrixRowMap()),
false);
199 #ifdef IFPACK_FLOPCOUNTERS 203 H_ = Teuchos::rcp(
new Epetra_CrsMatrix(Copy,*SerialMap_,0));
208 IFPACK_CHK_ERR(A_.ExtractMyRowCopy(0,Length,RowNnz,
209 &RowValues[0],&RowIndices[0]));
215 for (
int i = 0 ;i < RowNnz ; ++i)
217 if (RowIndices[i] < NumMyRows_){
218 RowIndices[count] = RowIndices[i];
219 RowValues[count] = RowValues[i];
229 double diag_val = 0.0;
230 for (
int i = 0 ;i < RowNnz ; ++i) {
231 if (RowIndices[i] == 0) {
232 double& v = RowValues[i];
239 diag_val = sqrt(diag_val);
240 int_type diag_idx = 0;
241 EPETRA_CHK_ERR(H_->InsertGlobalValues(0,1,&diag_val, &diag_idx));
249 for (
int row_i = 1 ; row_i < NumMyRows_ ; ++row_i) {
252 IFPACK_CHK_ERR(A_.ExtractMyRowCopy(row_i,Length,RowNnz,
253 &RowValues[0],&RowIndices[0]));
259 for (
int i = 0 ;i < RowNnz ; ++i)
261 if (RowIndices[i] < NumMyRows_){
262 RowIndices[count] = RowIndices[i];
263 RowValues[count] = RowValues[i];
275 if (LOF == 0) LOF = 1;
281 for (
int i = 0 ; i < RowNnz ; ++i) {
282 if (RowIndices[i] == row_i) {
283 double& v = RowValues[i];
286 else if (RowIndices[i] < row_i)
288 Hash.set(RowIndices[i], RowValues[i],
true);
295 for (
int col_j = RowIndices[0] ; col_j < row_i ; ++col_j) {
297 double h_ij = 0.0, h_jj = 0.0;
299 h_ij = Hash.get(col_j);
302 int_type* ColIndices;
305 int_type col_j_GID = (int_type) H_->RowMap().GID64(col_j);
306 H_->ExtractGlobalRowView(col_j_GID, ColNnz, ColValues, ColIndices);
308 for (
int k = 0 ; k < ColNnz ; ++k) {
309 int_type col_k = ColIndices[k];
314 double xxx = Hash.get(col_k);
317 h_ij -= ColValues[k] * xxx;
318 #ifdef IFPACK_FLOPCOUNTERS 327 if (IFPACK_ABS(h_ij) > DropTolerance_)
329 Hash.set(col_j, h_ij);
332 #ifdef IFPACK_FLOPCOUNTERS 334 ComputeFlops_ += 2.0 * flops + 1.0;
338 int size = Hash.getNumEntries();
340 std::vector<double> AbsRow(size);
344 std::vector<int_type> keys(size + 1);
345 std::vector<double> values(size + 1);
347 Hash.arrayify(&keys[0], &values[0]);
349 for (
int i = 0 ; i < size ; ++i)
351 AbsRow[i] = IFPACK_ABS(values[i]);
357 nth_element(AbsRow.begin(), AbsRow.begin() + LOF, AbsRow.begin() + count,
359 std::greater<double>());
360 cutoff = AbsRow[LOF];
363 for (
int i = 0 ; i < size ; ++i)
365 h_ii -= values[i] * values[i];
368 if (h_ii < 0.0) h_ii = 1e-12;;
372 #ifdef IFPACK_FLOPCOUNTERS 374 ComputeFlops_ += 2 * size + 1;
377 double DiscardedElements = 0.0;
380 for (
int i = 0 ; i < size ; ++i)
382 if (IFPACK_ABS(values[i]) > cutoff)
384 values[count] = values[i];
385 keys[count] = keys[i];
389 DiscardedElements += values[i];
394 h_ii += DiscardedElements;
397 values[count] = h_ii;
398 keys[count] = (int_type) H_->RowMap().GID64(row_i);
401 H_->InsertGlobalValues((int_type) H_->RowMap().GID64(row_i), count, &(values[0]), (int_type*)&(keys[0]));
404 IFPACK_CHK_ERR(H_->FillComplete());
408 Epetra_Vector LHS(
Matrix().RowMatrixRowMap());
409 Epetra_Vector RHS1(
Matrix().RowMatrixRowMap());
410 Epetra_Vector RHS2(
Matrix().RowMatrixRowMap());
411 Epetra_Vector RHS3(
Matrix().RowMatrixRowMap());
414 Matrix().Multiply(
false,LHS,RHS1);
415 H_->Multiply(
true,LHS,RHS2);
416 H_->Multiply(
false,RHS2,RHS3);
418 RHS1.Update(-1.0, RHS3, 1.0);
422 long long MyNonzeros = H_->NumGlobalNonzeros64();
423 Comm().SumAll(&MyNonzeros, &GlobalNonzeros_, 1);
426 #ifdef IFPACK_FLOPCOUNTERS 428 A_.Comm().SumAll(&flops, &TotalFlops, 1);
429 ComputeFlops_ += TotalFlops;
432 ComputeTime_ += Time_.ElapsedTime();
439 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 440 if(A_.RowMatrixRowMap().GlobalIndicesInt()) {
441 return TCompute<int>();
445 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 446 if(A_.RowMatrixRowMap().GlobalIndicesLongLong()) {
447 return TCompute<long long>();
451 throw "Ifpack_ICT::Compute: GlobalIndices type unknown for A_";
456 Epetra_MultiVector& Y)
const 462 if (X.NumVectors() != Y.NumVectors())
465 Time_.ResetStartTime();
469 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
470 if (X.Pointers()[0] == Y.Pointers()[0])
471 Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X) );
473 Xcopy = Teuchos::rcp( &X,
false );
479 EPETRA_CHK_ERR(H_->Solve(
false,
false,
false,*Xcopy,Y));
480 EPETRA_CHK_ERR(H_->Solve(
false,
true,
false,Y,Y));
482 #ifdef IFPACK_FLOPCOUNTERS 484 ApplyInverseFlops_ += 4.0 * GlobalNonzeros_;
488 ApplyInverseTime_ += Time_.ElapsedTime();
494 int Ifpack_ICT::Apply(
const Epetra_MultiVector& X,
495 Epetra_MultiVector& Y)
const 503 const int MaxIters,
const double Tol,
504 Epetra_RowMatrix* Matrix_in)
510 if (Condest_ == -1.0)
511 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
522 if (!
Comm().MyPID()) {
524 os <<
"================================================================================" << endl;
525 os <<
"Ifpack_ICT: " << Label() << endl << endl;
529 os <<
"Relax value = " <<
RelaxValue() << endl;
530 os <<
"Condition number estimate = " <<
Condest() << endl;
531 os <<
"Global number of rows = " <<
Matrix().NumGlobalRows64() << endl;
533 os <<
"Number of nonzeros of H = " << H_->NumGlobalNonzeros64() << endl;
534 os <<
"nonzeros / rows = " 535 << 1.0 * H_->NumGlobalNonzeros64() / H_->NumGlobalRows64() << endl;
538 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
539 os <<
"----- ------- -------------- ------------ --------" << endl;
542 <<
" 0.0 0.0" << endl;
543 os <<
"Compute() " << std::setw(5) <<
NumCompute()
549 os <<
" " << std::setw(15) << 0.0 << endl;
556 os <<
" " << std::setw(15) << 0.0 << endl;
557 os <<
"================================================================================" << endl;
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual double ComputeFlops() const
Returns the number of flops in all applications of Compute().
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
double AbsoluteThreshold() const
Returns the absolute threshold.
int SetParameters(Teuchos::ParameterList ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
int Initialize()
Initialize L and U with values from user matrix A.
double RelaxValue() const
Returns the relaxation value.
virtual double InitializeTime() const
Returns the time spent in Initialize().
Ifpack_ICT(const Epetra_RowMatrix *A)
Ifpack_ICT constuctor with variable number of indices per row.
double DropTolerance() const
Returns the drop threshold.
virtual int NumCompute() const
Returns the number of calls to Compute().
bool IsInitialized() const
Returns true is the preconditioner has been successfully initialized.
virtual double ApplyInverseFlops() const
Returns the number of flops in all applications of ApplyInverse().
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Ifpack_ICT forward/back solve on a Epetra_MultiVector X in Y...
double RelativeThreshold() const
Returns the relative threshold.
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
virtual ~Ifpack_ICT()
Ifpack_ICT Destructor.
double LevelOfFill() const
Returns the level-of-fill.