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),
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();
160 template<
typename int_type>
166 Time_.ResetStartTime();
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) 180 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) 181 if(
A_.RowMatrixRowMap().GlobalIndicesInt())
185 #if !defined(EPETRA_NO_64BIT_GLOBAL_INDICES) 186 if(
A_.RowMatrixRowMap().GlobalIndicesLongLong())
190 throw "Ifpack_ICT::TCompute: Global indices unknown.";
199 #ifdef IFPACK_FLOPCOUNTERS 209 &RowValues[0],&RowIndices[0]));
215 for (
int i = 0 ;i < RowNnz ; ++i)
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;
249 for (
int row_i = 1 ; row_i <
NumMyRows_ ; ++row_i) {
253 &RowValues[0],&RowIndices[0]));
259 for (
int i = 0 ;i < RowNnz ; ++i)
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 329 Hash.set(col_j, h_ij);
332 #ifdef IFPACK_FLOPCOUNTERS 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)
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 377 double DiscardedElements = 0.0;
380 for (
int i = 0 ; i < size ; ++i)
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]));
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();
426 #ifdef IFPACK_FLOPCOUNTERS 428 A_.Comm().SumAll(&flops, &TotalFlops, 1);
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])
482 #ifdef IFPACK_FLOPCOUNTERS 495 Epetra_MultiVector& Y)
const 503 const int MaxIters,
const double Tol,
504 Epetra_RowMatrix* 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().
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
virtual double ComputeTime() const
Returns the time spent in Compute().
int NumCompute_
Contains the number of successful call to Compute().
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
double ComputeTime_
Contains the time for all successful calls to Compute().
T & get(ParameterList &l, const std::string &name)
double LevelOfFill_
Level of fill.
Epetra_Time Time_
Used for timing purposes.
double AbsoluteThreshold() const
Returns the absolute threshold.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
int SetParameters(Teuchos::ParameterList ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
double ComputeFlops_
Contains the number of flops for Compute().
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.
Teuchos::RefCountPtr< Epetra_CrsMatrix > H_
Contains the Cholesky factorization.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
int Initialize()
Initialize L and U with values from user matrix A.
double DropTolerance_
During factorization, drop all values below this.
Teuchos::RefCountPtr< Epetra_Map > SerialMap_
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double RelaxValue() const
Returns the relaxation value.
long long GlobalNonzeros_
Global number of nonzeros in L and U factors.
virtual double InitializeTime() const
Returns the time spent in Initialize().
std::string Label_
Label of this object.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ifpack_ICT(const Epetra_RowMatrix *A)
Ifpack_ICT constuctor with variable number of indices per row.
double DropTolerance() const
Returns the drop threshold.
double Rthresh_
Relative threshold.
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual int NumCompute() const
Returns the number of calls to Compute().
bool IsInitialized() const
Returns true is the preconditioner has been successfully initialized.
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual double ApplyInverseFlops() const
Returns the number of flops in all applications of ApplyInverse().
#define EPETRA_CHK_ERR(xxx)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
void Destroy()
Destroys all data associated to the preconditioner.
double Relax_
Relaxation value.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
const char * Label() const
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.
Teuchos::RefCountPtr< Epetra_SerialComm > SerialComm_
double Condest_
Contains the estimate of the condition number, if -1.0 if not computed.
bool IsComputed_
If true, the preconditioner has been successfully computed.
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
int NumMyRows_
Number of local rows in the matrix.
#define IFPACK_CHK_ERR(ifpack_err)
const Epetra_RowMatrix & A_
Reference to the matrix to be preconditioned, supposed symmetric.
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.
double Athresh_
Absolute threshold.