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;
152 IsInitialized_ =
true;
160 template<
typename int_type>
161 int Ifpack_ICT::TCompute()
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)
182 SerialMap_ = Teuchos::rcp(
new Epetra_Map(NumMyRows_, 0, *SerialComm_));
185 #if !defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
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
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) {
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());
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
429 ComputeFlops_ += TotalFlops;
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 );
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_;
503 const int MaxIters,
const double Tol,
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;