43 #ifndef IFPACK_SUPPORTGRAPH_H
44 #define IFPACK_SUPPORTGRAPH_H
46 #include "Ifpack_ConfigDefs.h"
47 #include "Ifpack_Condest.h"
48 #include "Ifpack_Preconditioner.h"
49 #include "Ifpack_Amesos.h"
50 #include "Ifpack_Condest.h"
51 #include "Epetra_Map.h"
52 #include "Epetra_Comm.h"
53 #include "Epetra_Time.h"
54 #include "Epetra_Vector.h"
55 #include "Epetra_MultiVector.h"
56 #include "Epetra_LinearProblem.h"
57 #include "Epetra_RowMatrix.h"
58 #include "Epetra_CrsMatrix.h"
60 #include "Teuchos_ParameterList.hpp"
61 #include "Teuchos_RefCountPtr.hpp"
63 #include <boost/graph/adjacency_list.hpp>
64 #include <boost/graph/kruskal_min_spanning_tree.hpp>
65 #include <boost/graph/prim_minimum_spanning_tree.hpp>
66 #include <boost/config.hpp>
68 using Teuchos::RefCountPtr;
70 typedef std::pair<int, int> E;
71 using namespace boost;
73 typedef adjacency_list < vecS, vecS, undirectedS,
74 no_property, property < edge_weight_t, double > > Graph;
75 typedef graph_traits < Graph >::edge_descriptor Edge;
76 typedef graph_traits < Graph >::vertex_descriptor Vertex;
105 virtual int SetUseTranspose(
bool UseTranspose_in);
146 virtual const char * Label()
const;
171 return(IsInitialized_);
193 virtual int SetParameters(Teuchos::ParameterList& List);
199 virtual int Initialize();
204 virtual int Compute();
216 virtual double Condest(
const Ifpack_CondestType CT = Ifpack_Cheap,
217 const int MaxIters = 1550,
218 const double Tol = 1e-9,
234 virtual std::ostream& Print(std::ostream&)
const;
239 return(NumInitialize_);
251 return(NumApplyInverse_);
257 return(InitializeTime_);
263 return(ComputeTime_);
269 return(ApplyInverseTime_);
275 return(InitializeFlops_);
281 return(ComputeFlops_);
287 return(ApplyInverseFlops_);
299 Teuchos::RefCountPtr<const Epetra_RowMatrix>
Matrix_;
350 Teuchos::RefCountPtr<Epetra_Time>
Time_;
377 Matrix_(rcp(Matrix_in,false)),
378 IsInitialized_(false),
380 UseTranspose_(false),
385 InitializeTime_(0.0),
387 ApplyInverseTime_(0.0),
388 InitializeFlops_(0.0),
390 ApplyInverseFlops_(0.0),
398 Teuchos::ParameterList List_in;
407 long long rows = (*Matrix_).NumGlobalRows64();
408 long long cols = (*Matrix_).NumGlobalCols64();
409 int num_edges = ((*Matrix_).NumMyNonzeros() - (*Matrix_).NumMyDiagonals())/2;
410 std::cout <<
"global num rows " << rows << std::endl;
413 IFPACK_CHK_ERR((rows == cols));
415 if(rows > std::numeric_limits<int>::max())
417 std::cerr <<
"Ifpack_SupportGraph<T>::FindSupport: global num rows won't fit an int. " << rows << std::endl;
422 int num_verts = (int) rows;
425 E *edge_array =
new E[num_edges];
426 double *weights =
new double[num_edges];
429 int max_num_entries = (*Matrix_).MaxNumEntries();
430 double *values =
new double[max_num_entries];
431 int *indices =
new int[max_num_entries];
433 double * diagonal =
new double[num_verts];
436 for(
int i = 0; i < max_num_entries; i++)
444 for(
int i = 0; i < num_verts; i++)
446 (*Matrix_).ExtractMyRowCopy(i,max_num_entries,num_entries,values,indices);
448 for(
int j = 0; j < num_entries; j++)
453 diagonal[i] = values[j];
456 diagonal[i] *= DiagPertRel_;
458 diagonal[i] += DiagPertAbs_;
463 edge_array[k] = E(i,indices[j]);
464 weights[k] = values[j];
468 weights[k] *= (1.0 + 1e-8 * drand48());
477 Graph g(edge_array, edge_array + num_edges, weights, num_verts);
480 property_map < Graph, edge_weight_t >::type weight =
get(edge_weight, g);
482 std::vector < Edge > spanning_tree;
485 kruskal_minimum_spanning_tree(g, std::back_inserter(spanning_tree));
488 std::vector<int> NumNz(num_verts,1);
491 for (std::vector < Edge >::iterator ei = spanning_tree.begin();
492 ei != spanning_tree.end(); ++ei)
494 NumNz[source(*ei,g)] = NumNz[source(*ei,g)] + 1;
495 NumNz[target(*ei,g)] = NumNz[target(*ei,g)] + 1;
500 std::vector< std::vector< int > > Indices(num_verts);
505 std::vector< std::vector< double > > Values(num_verts);
507 for(
int i = 0; i < num_verts; i++)
509 std::vector<int> temp(NumNz[i],0);
510 std::vector<double> temp2(NumNz[i],0);
515 int *l =
new int[num_verts];
516 for(
int i = 0; i < num_verts; i++)
524 for(
int i = 0; i < NumForests_; i++)
528 spanning_tree.clear();
529 kruskal_minimum_spanning_tree(g,std::back_inserter(spanning_tree));
530 for(std::vector < Edge >::iterator ei = spanning_tree.begin();
531 ei != spanning_tree.end(); ++ei)
533 NumNz[source(*ei,g)] = NumNz[source(*ei,g)] + 1;
534 NumNz[target(*ei,g)] = NumNz[target(*ei,g)] + 1;
536 for(
int i = 0; i < num_verts; i++)
538 Indices[i].resize(NumNz[i]);
539 Values[i].resize(NumNz[i]);
543 for (std::vector < Edge >::iterator ei = spanning_tree.begin();
544 ei != spanning_tree.end(); ++ei)
548 Indices[source(*ei,g)][0] = source(*ei,g);
549 Values[source(*ei,g)][0] = Values[source(*ei,g)][0] - weight[*ei];
550 Indices[target(*ei,g)][0] = target(*ei,g);
551 Values[target(*ei,g)][0] = Values[target(*ei,g)][0] - weight[*ei];
553 Indices[source(*ei,g)][l[source(*ei,g)]] = target(*ei,g);
554 Values[source(*ei,g)][l[source(*ei,g)]] = weight[*ei];
555 l[source(*ei,g)] = l[source(*ei,g)] + 1;
557 Indices[target(*ei,g)][l[target(*ei,g)]] = source(*ei,g);
558 Values[target(*ei,g)][l[target(*ei,g)]] = weight[*ei];
559 l[target(*ei,g)] = l[target(*ei,g)] + 1;
576 Matrix_->Multiply(
false, ones, surplus);
578 for(
int i = 0; i < num_verts; i++)
580 Values[i][0] += surplus[i];
581 Values[i][0] = KeepDiag_*diagonal[i] +
582 (1.-KeepDiag_) * Values[i][0];
588 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
589 if((*Matrix_).RowMatrixRowMap().GlobalIndicesLongLong())
592 for(
int i = 0; i < num_verts; i++)
594 std::vector<long long> IndicesLL(l[i]);
595 for(
int k = 0; k < l[i]; ++k)
596 IndicesLL[k] = Indices[i][k];
598 (*Support_).InsertGlobalValues(i,l[i],&Values[i][0],&IndicesLL[0]);
603 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
604 if((*Matrix_).RowMatrixRowMap().GlobalIndicesInt())
607 for(
int i = 0; i < num_verts; i++)
609 (*Support_).InsertGlobalValues(i,l[i],&Values[i][0],&Indices[i][0]);
614 throw "Ifpack_SupportGraph::FindSupport: GlobalIndices unknown.";;
616 (*Support_).FillComplete();
632 NumForests_ = List_in.get(
"MST: forest number", NumForests_);
633 KeepDiag_ = List_in.get(
"MST: keep diagonal", KeepDiag_);
634 Randomize_ = List_in.get(
"MST: randomize", Randomize_);
636 DiagPertRel_ = List_in.get(
"fact: relative threshold", DiagPertRel_);
637 DiagPertAbs_ = List_in.get(
"fact: absolute threshold", DiagPertAbs_);
645 IsInitialized_ =
false;
649 if (Time_ == Teuchos::null)
655 Time_->ResetStartTime();
659 Inverse_ = Teuchos::rcp(
new T(Support_.get()));
661 IFPACK_CHK_ERR(Inverse_->Initialize());
663 IsInitialized_ =
true;
665 InitializeTime_ += Time_->ElapsedTime();
674 if (IsInitialized() ==
false)
675 IFPACK_CHK_ERR(Initialize());
677 Time_->ResetStartTime();
681 IFPACK_CHK_ERR(Inverse_->Compute());
685 ComputeTime_ += Time_->ElapsedTime();
696 UseTranspose_ = UseTranspose_in;
699 if (Inverse_!=Teuchos::null)
700 IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose_in));
709 IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
716 return(Label_.c_str());
726 Time_->ResetStartTime();
729 Inverse_->ApplyInverse(X,Y);
732 ApplyInverseTime_ += Time_->ElapsedTime();
741 os <<
"================================================================================" << std::endl;
742 os <<
"Ifpack_SupportGraph: " << Label () << std::endl << std::endl;
743 os <<
"Condition number estimate = " << Condest() << std::endl;
744 os <<
"Global number of rows = " << Matrix_->NumGlobalRows64() << std::endl;
745 os <<
"Number of edges in support graph = " << (Support_->NumGlobalNonzeros64()-Support_->NumGlobalDiagonals64())/2 << std::endl;
746 os <<
"Fraction of off diagonals of support graph/off diagonals of original = "
747 << ((double)Support_->NumGlobalNonzeros64()-Support_->NumGlobalDiagonals64())/(Matrix_->NumGlobalNonzeros64()-Matrix_->NumGlobalDiagonals64());
749 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << std::endl;
750 os <<
"----- ------- -------------- ------------ --------" << std::endl;
751 os <<
"Initialize() " << std::setw(10) << NumInitialize_
752 <<
" " << std::setw(15) << InitializeTime_
753 <<
" 0.0 0.0" << std::endl;
754 os <<
"Compute() " << std::setw(10) << NumCompute_
755 <<
" " << std::setw(22) << ComputeTime_
756 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
757 if (ComputeTime_ != 0.0)
758 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << std::endl;
760 os <<
" " << std::setw(15) << 0.0 << std::endl;
761 os <<
"ApplyInverse() " << std::setw(10) << NumApplyInverse_
762 <<
" " << std::setw(22) << ApplyInverseTime_
763 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
764 if (ApplyInverseTime_ != 0.0)
765 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << std::endl;
767 os <<
" " << std::setw(15) << 0.0 << std::endl;
769 os << std::endl << std::endl;
770 os <<
"Now calling the underlying preconditioner's print()" << std::endl;
779 Condest(
const Ifpack_CondestType CT,
const int MaxIters,
787 Condest_ = Ifpack_Condest(*
this, CT, MaxIters, Tol, Matrix_in);
792 #endif // IFPACK_SUPPORTGRAPH_H
double ComputeTime_
Contains the time for all successful calls to Compute().
virtual double Condest() const
Returns the computed condition number.
Teuchos::RefCountPtr< Epetra_Time > Time_
Object used for timing purposes.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
bool IsComputed_
If true, the preconditioner has been successfully computed.
virtual const char * Label() const
Returns a character string describing the operator.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual double ApplyInverseFlops() const
Returns the total number of flops to apply the preconditioner.
int NumInitialize_
Contains the number of successful calls to Initialize().
Teuchos::ParameterList List_
Stores a copy of the list given in SetParameters()
virtual double ApplyInverseTime() const
Returns the total time spent in ApplyInverse().
std::string Label_
Contains the label of this object.
int Randomize_
Option to add random pertubation to edge weights, to get random spanning trees.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
Ifpack_SupportGraph(Epetra_RowMatrix *Matrix_in)
Constructor.
virtual double InitializeTime() const
Returns the total time spent in Initialize().
Teuchos::RefCountPtr< Epetra_CrsMatrix > Support_
Pointers to the matrix of the support graph.
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual double ComputeFlops() const
Returns the total number of flops to compute the preconditioner.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
double Condest_
Contains the estimated condition number.
virtual bool HasNormInf() const
Returns true if this object can provide an approximate Inf-norm, false otherwise. ...
virtual int SetUseTranspose(bool UseTranspose_in)
If set true, transpose of this operator will be applied (not implemented).
virtual const Epetra_RowMatrix & Matrix() const
Returns a const reference to the internally stored matrix.
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
int NumCompute_
Contains the number of successful call to Compute().
int NumForests_
Contains the number of forests in the support graph.
Teuchos::RefCountPtr< T > Inverse_
Pointer to the local solver.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
double InitializeFlops_
Contains the number of flops for Initialize().
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual std::ostream & Print(std::ostream &) const
Prints on ostream basic information about this object.
double ComputeFlops_
Contains the number of flops for Compute().
double KeepDiag_
Contains the option to keep the diagonal of original matrix, or weighted average. ...
double InitializeTime_
Contains the time for all successful calls to Initialize().
bool UseTranspose_
If true, solve with the transpose (not supported by all solvers).
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual int Initialize()
Initialize the preconditioner.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
virtual int Compute()
Computes the preconditioners.
double DiagPertAbs_
Absolute diagonal pertubation.
virtual double ComputeTime() const
Returns the total time spent in Compute().
int FindSupport()
Compute the support graph.
double DiagPertRel_
Relative diagonal pertubation.
virtual double InitializeFlops() const
Returns the number of flops in the initialization phase.