43#ifndef IFPACK_SUPPORTGRAPH_H
44#define IFPACK_SUPPORTGRAPH_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>
68using Teuchos::RefCountPtr;
70typedef std::pair<int, int>
E;
73typedef adjacency_list < vecS, vecS, undirectedS,
74 no_property, property < edge_weight_t, double > >
Graph;
75typedef graph_traits < Graph >::edge_descriptor
Edge;
76typedef graph_traits < Graph >::vertex_descriptor
Vertex;
146 virtual const char *
Label()
const;
217 const int MaxIters = 1550,
218 const double Tol = 1e-9,
234 virtual std::ostream&
Print(std::ostream&)
const;
299 Teuchos::RefCountPtr<const Epetra_RowMatrix>
Matrix_;
350 Teuchos::RefCountPtr<Epetra_Time>
Time_;
377Matrix_(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;
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()));
663 IsInitialized_ =
true;
665 InitializeTime_ += Time_->ElapsedTime();
674 if (IsInitialized() ==
false)
677 Time_->ResetStartTime();
685 ComputeTime_ += Time_->ElapsedTime();
696 UseTranspose_ = UseTranspose_in;
699 if (Inverse_!=Teuchos::null)
716 return(Label_.c_str());
726 Time_->ResetStartTime();
729 Inverse_->ApplyInverse(X,Y);
732 ApplyInverseTime_ += Time_->ElapsedTime();
739Print(std::ostream& os)
const
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;
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
@ Ifpack_Cheap
cheap 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)
graph_traits< Graph >::edge_descriptor Edge
graph_traits< Graph >::vertex_descriptor Vertex
adjacency_list< vecS, vecS, undirectedS, no_property, property< edge_weight_t, double > > Graph
int PutScalar(double ScalarConstant)
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
double DiagPertRel_
Relative diagonal pertubation.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the preconditioner to X, returns the result in Y.
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual double Condest() const
Returns the computed condition number.
int NumApplyInverse_
Contains the number of successful call to ApplyInverse().
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully initialized.
virtual int SetUseTranspose(bool UseTranspose_in)
If set true, transpose of this operator will be applied (not implemented).
bool IsInitialized_
If true, the preconditioner has been successfully initialized.
Ifpack_SupportGraph(Epetra_RowMatrix *Matrix_in)
Constructor.
int NumForests_
Contains the number of forests in the support graph.
double Condest_
Contains the estimated condition number.
virtual const char * Label() const
Returns a character string describing the operator.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
Teuchos::RefCountPtr< Epetra_CrsMatrix > Support_
Pointers to the matrix of the support graph.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
int NumInitialize_
Contains the number of successful calls to Initialize().
int NumCompute_
Contains the number of successful call to Compute().
double KeepDiag_
Contains the option to keep the diagonal of original matrix, or weighted average.
std::string Label_
Contains the label of this object.
Teuchos::RefCountPtr< T > Inverse_
Pointer to the local solver.
virtual int Initialize()
Initialize the preconditioner.
virtual double ComputeFlops() const
Returns the total number of flops to compute the preconditioner.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
double InitializeFlops_
Contains the number of flops for Initialize().
Teuchos::RefCountPtr< Epetra_Time > Time_
Object used for timing purposes.
virtual double ApplyInverseFlops() const
Returns the total number of flops to apply the preconditioner.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual bool HasNormInf() const
Returns true if this object can provide an approximate Inf-norm, false otherwise.
virtual const Epetra_RowMatrix & Matrix() const
Returns a const reference to the internally stored matrix.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
int Randomize_
Option to add random pertubation to edge weights, to get random spanning trees.
int FindSupport()
Compute the support graph.
double ComputeFlops_
Contains the number of flops for Compute().
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
virtual std::ostream & Print(std::ostream &) const
Prints on ostream basic information about this object.
virtual int Compute()
Computes the preconditioners.
double DiagPertAbs_
Absolute diagonal pertubation.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
bool IsComputed_
If true, the preconditioner has been successfully computed.
double ApplyInverseFlops_
Contain sthe number of flops for ApplyInverse().
double InitializeTime_
Contains the time for all successful calls to Initialize().
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
Teuchos::ParameterList List_
Stores a copy of the list given in SetParameters()
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
double ApplyInverseTime_
Contains the time for all successful calls to ApplyInverse().
virtual double ApplyInverseTime() const
Returns the total time spent in ApplyInverse().
virtual double InitializeFlops() const
Returns the number of flops in the initialization phase.
double ComputeTime_
Contains the time for all successful calls to Compute().
virtual double ComputeTime() const
Returns the total time spent in Compute().
bool UseTranspose_
If true, solve with the transpose (not supported by all solvers).
Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_
Pointers to the matrix to be preconditioned.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
virtual double InitializeTime() const
Returns the total time spent in Initialize().