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;
160template<
typename int_type>
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)
185#if !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
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
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]));
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
439#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
441 return TCompute<int>();
445#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
447 return TCompute<long long>();
451 throw "Ifpack_ICT::Compute: GlobalIndices type unknown for A_";
469 Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
473 Xcopy = Teuchos::rcp( &X,
false );
482#ifdef IFPACK_FLOPCOUNTERS
503 const int MaxIters,
const double Tol,
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;
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;
#define EPETRA_CHK_ERR(a)
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
double Ifpack_Condest(const Ifpack_Preconditioner &IFP, const Ifpack_CondestType CT, const int MaxIters, const double Tol, Epetra_RowMatrix *Matrix)
#define IFPACK_CHK_ERR(ifpack_err)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
bool GlobalIndicesInt() const
bool GlobalIndicesLongLong() const
virtual int NumProc() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
double ** Pointers() const
virtual const Epetra_Comm & Comm() const=0
virtual int NumMyRows() const=0
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual long long NumGlobalRows64() const=0
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
void ResetStartTime(void)
double ElapsedTime(void) const
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
double Rthresh_
Relative threshold.
double LevelOfFill() const
Returns the level-of-fill.
const Epetra_RowMatrix & A_
Reference to the matrix to be preconditioned, supposed symmetric.
double ComputeFlops_
Contains the number of flops for Compute().
double DropTolerance_
During factorization, drop all values below this.
double RelaxValue() const
Returns the relaxation value.
double Relax_
Relaxation value.
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.
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
Ifpack_ICT(const Epetra_RowMatrix *A)
Ifpack_ICT constuctor with variable number of indices per row.
double DropTolerance() const
Returns the drop threshold.
double Condest_
Contains the estimate of the condition number, if -1.0 if not computed.
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
int NumInitialize_
Contains the number of successful calls to Initialize().
virtual double ComputeTime() const
Returns the time spent in Compute().
void Destroy()
Destroys all data associated to the preconditioner.
double Athresh_
Absolute threshold.
Teuchos::RefCountPtr< Epetra_SerialComm > SerialComm_
double ComputeTime_
Contains the time for all successful calls to Compute().
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
long long GlobalNonzeros_
Global number of nonzeros in L and U factors.
double RelativeThreshold() const
Returns the relative threshold.
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
int SetParameters(Teuchos::ParameterList ¶meterlis)
Set parameters using a Teuchos::ParameterList object.
bool IsComputed() const
If factor is completed, this query returns true, otherwise it returns false.
double AbsoluteThreshold() const
Returns the absolute threshold.
virtual double ApplyInverseFlops() const
Returns the number of flops in all applications of ApplyInverse().
int Initialize()
Initialize L and U with values from user matrix A.
virtual ~Ifpack_ICT()
Ifpack_ICT Destructor.
bool IsComputed_
If true, the preconditioner has been successfully computed.
Teuchos::RefCountPtr< Epetra_Map > SerialMap_
double LevelOfFill_
Level of fill.
int NumMyRows_
Number of local rows in the matrix.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
bool IsInitialized() const
Returns true is the preconditioner has been successfully initialized.
const char * Label() const
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
Epetra_Time Time_
Used for timing purposes.
int NumCompute_
Contains the number of successful call to Compute().
virtual double ComputeFlops() const
Returns the number of flops in all applications of Compute().
int Compute()
Compute IC factor U using the specified graph, diagonal perturbation thresholds and relaxation parame...
Teuchos::RefCountPtr< Epetra_CrsMatrix > H_
Contains the Cholesky factorization.
std::string Label_
Label of this object.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
const Epetra_RowMatrix & Matrix() const
Returns a reference to the matrix to be preconditioned.
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
double get(const key_type key)
Returns an element from the hash table, or 0.0 if not found.
void arrayify(key_type *key_array, double *val_array)
Converts the contents in array format for both keys and values.
void set(const key_type key, const double value, const bool addToValue=false)
Sets an element in the hash table.
void reset()
Resets the entries of the already allocated memory. This method can be used to clean an array,...
int getNumEntries() const
Returns the number of stored entries.