34 #ifndef __MatrixEpetra_H
35 #define __MatrixEpetra_H 1
37 #include <feel/feelconfig.h>
45 #if defined( FEELPP_HAS_TRILINOS_EPETRA )
46 #undef PACKAGE_BUGREPORT
49 #undef PACKAGE_TARNAME
50 #undef PACKAGE_VERSION
51 #if defined(FEELPP_HAS_MPI)
52 #include <Epetra_MpiComm.h>
54 #include <Epetra_SerialComm.h>
56 #include <Epetra_Map.h>
57 #include <EpetraExt_MatrixMatrix.h>
58 #include <EpetraExt_RowMatrixOut.h>
59 #include <Epetra_FECrsMatrix.h>
60 #include <Epetra_Vector.h>
61 #undef PACKAGE_BUGREPORT
64 #undef PACKAGE_TARNAME
65 #undef PACKAGE_VERSION
69 template<
typename T>
class VectorEpetra;
84 class MatrixEpetra :
public MatrixSparse<double>
86 typedef MatrixSparse<double> super;
94 static const bool is_row_major =
true;
95 typedef super::value_type value_type;
96 typedef super::real_type real_type;
97 typedef std::vector<std::set<size_type> > pattern_type;
99 typedef super::graph_type graph_type;
100 typedef super::graph_ptrtype graph_ptrtype;
102 typedef Vector<value_type> vector_type;
103 typedef boost::shared_ptr<Vector<value_type> > vector_ptrtype;
105 typedef VectorEpetra<value_type> epetra_vector_type;
106 typedef boost::shared_ptr<epetra_vector_type> epetra_vector_ptrtype;
108 typedef MatrixEpetra epetra_sparse_matrix_type;
125 MatrixEpetra ( Epetra_FECrsMatrix
const& m );
130 MatrixEpetra ( MatrixEpetra
const& T );
135 MatrixEpetra ( Epetra_Map
const& emap,
int nnz = 50 );
137 MatrixEpetra ( Epetra_Vector
const& x );
139 MatrixEpetra ( Epetra_Map
const& rowmap, Epetra_Map
const& colmap );
141 MatrixEpetra( Epetra_Map
const& row_emap, Epetra_Map
const& col_emap,
142 Epetra_Map
const& dom_map, Epetra_Map
const& range_map );
166 MatrixEpetra &
operator = ( MatrixSparse<value_type>
const& M )
181 value_type operator () (
const size_type i,
190 MatrixEpetra&
operator=(
const MatrixEpetra& T )
194 M_emap = T.getRowMap();
195 M_col_emap = T.getColMap();
196 M_dom_map = T.mat().DomainMap();
197 M_range_map = T.mat().RangeMap();
200 this->setInitialized( T.isInitialized() );
240 Epetra_Map getRowMap()
const
248 Epetra_Map getColMap()
const
256 Epetra_Map getDomainMap()
const
258 return M_mat->DomainMap();
264 Epetra_Map getRangeMap()
const
266 return M_mat->RangeMap();
271 std::vector<size_type> bcIndices()
276 boost::shared_ptr<Epetra_FECrsMatrix> matrix()
314 graph_ptrtype
const& graph );
362 void diagonal ( Vector<double>& dest )
const;
364 void setDiagonal( Epetra_Vector
const& x );
391 bool EpetraIndicesAreLocal()
const;
393 bool HaveColMap()
const;
400 void transpose( MatrixSparse<value_type>& Mt )
const;
405 real_type energy ( vector_type
const& v1, vector_type
const& v2,
bool tranpose =
false )
const;
414 real_type l1Norm ()
const;
426 real_type linftyNorm ()
const;
436 const value_type& value );
441 void multiplyMatrix (
const MatrixEpetra& A,
const MatrixEpetra& B );
444 void multiply (
bool trans,
const Vector<T>& v, Vector<T>& r )
const
446 epetra_vector_type
const& ev( dynamic_cast<epetra_vector_type const&>( v ) );
447 epetra_vector_type& er( dynamic_cast<epetra_vector_type&>( r ) );
448 M_mat->Multiply( trans, ev.vec(), er.vec() );
451 void multiply (
bool trans,
const Epetra_MultiVector &X, Epetra_MultiVector &Y )
const
453 M_mat->Multiply( trans, X, Y );
467 const value_type& value );
475 void addMatrix (
double coeff,
const MatrixEpetra& _m )
477 EpetraExt::MatrixMatrix::Add( _m.mat(),
false, coeff, ( *this ).mat(), 1. );
480 void addMatrix (
const ublas::matrix<value_type> &,
481 const std::vector<size_type> &,
482 const std::vector<size_type> & )
494 void addMatrix (
int* rows,
int nrows,
495 int* cols,
int ncols,
502 void addMatrix (
const ublas::matrix<value_type> &dm,
503 const std::vector<size_type> &dof_indices )
505 this->addMatrix ( dm, dof_indices, dof_indices );
520 void addMatrix (
const value_type a, super &X );
526 void scale (
double const scalar );
532 Epetra_FECrsMatrix& mat ()
541 Epetra_FECrsMatrix
const& mat ()
const
553 void printMatlab(
const std::string name=
"NULL" )
const;
554 void printKonsole()
const;
571 void zeroRows( std::vector<int>
const& rows, std::vector<value_type>
const& values, Vector<value_type>& rhs, Context
const& on_context );
577 void updateBlockMat( boost::shared_ptr<MatrixSparse<value_type> > m, std::vector<size_type> start_i, std::vector<size_type> start_j );
583 mpi::communicator M_comm;
607 Epetra_Map M_col_emap;
608 Epetra_Map M_dom_map;
609 Epetra_Map M_range_map;
611 mutable boost::shared_ptr<Epetra_FECrsMatrix> M_mat;
613 std::vector<size_type> M_bc_index;
619 MatrixEpetra::MatrixEpetra()
623 #ifdef FEELPP_HAS_MPI
624 M_emap( Epetra_Map( -1, 0, 0, Epetra_MpiComm( M_comm ) ) ),
625 M_col_emap( Epetra_Map( -1, 0, 0, Epetra_MpiComm( M_comm ) ) ),
626 M_dom_map( Epetra_Map( -1, 0, 0, Epetra_MpiComm( M_comm ) ) ),
627 M_range_map( Epetra_Map( -1, 0, 0, Epetra_MpiComm( M_comm ) ) ),
628 M_mat( new Epetra_FECrsMatrix( Copy, M_emap, 0 ) )
630 M_emap( Epetra_Map( -1, 0, 0, Epetra_SerialComm ) ),
631 M_col_emap( Epetra_Map( -1, 0, 0, Epetra_SerialComm ) ),
632 M_dom_map( Epetra_Map( -1, 0, 0, Epetra_SerialComm ) ),
633 M_range_map( Epetra_Map( -1, 0, 0, Epetra_SerialComm ) ),
634 M_mat( new Epetra_FECrsMatrix( Copy, M_emap, 0 ) )
640 MatrixEpetra::MatrixEpetra( Epetra_Map
const& emap,
int nnz )
647 M_mat( new Epetra_FECrsMatrix( Copy, emap, 0 ) )
654 MatrixEpetra::MatrixEpetra( Epetra_Map
const& row_emap, Epetra_Map
const& col_emap )
658 M_col_emap( col_emap ),
659 M_dom_map( col_emap ),
660 M_range_map( row_emap ),
661 M_mat( new Epetra_FECrsMatrix( Copy, row_emap, col_emap, 0 ) )
668 MatrixEpetra::MatrixEpetra( Epetra_Map
const& row_emap, Epetra_Map
const& col_emap,
669 Epetra_Map
const& dom_map, Epetra_Map
const& range_map )
673 M_col_emap( col_emap ),
674 M_dom_map( dom_map ),
675 M_range_map( range_map ),
676 M_mat( new Epetra_FECrsMatrix( Copy, M_emap, M_col_emap, 0 ) )
684 MatrixEpetra::MatrixEpetra( Epetra_FECrsMatrix
const& M )
687 M_emap( M.RowMap() ),
688 M_col_emap( M.ColMap() ),
689 M_dom_map( M.DomainMap() ),
690 M_range_map( M.RangeMap() ),
691 M_mat( new Epetra_FECrsMatrix( Copy, M_emap, M_col_emap, 0 ) )
697 MatrixEpetra::MatrixEpetra( MatrixEpetra
const& T )
701 M_col_emap( T.M_col_emap ),
702 M_dom_map( T.M_dom_map ),
703 M_range_map( T.M_range_map ),
704 M_mat( new Epetra_FECrsMatrix( T.mat() ) )
711 MatrixEpetra::MatrixEpetra(
const size_type m,
719 M_emap( Epetra_Map( m, m_l, 0, Epetra_MpiComm( M_comm ) ) ),
720 M_col_emap( Epetra_Map( n, n, 0, Epetra_MpiComm( M_comm ) ) ),
722 M_range_map( M_emap ),
723 M_mat( new Epetra_FECrsMatrix( Copy, M_emap, M_col_emap, nnz ) )
729 MatrixEpetra::~MatrixEpetra()
736 void MatrixEpetra::init (
const size_type m,
743 if ( ( m==0 ) || ( n==0 ) )
748 if ( this->isInitialized() )
751 this->setInitialized(
true );
760 void MatrixEpetra::init ()
762 if ( this->isInitialized() )
765 this->setInitialized(
true );
769 void MatrixEpetra::zero ()
771 FEELPP_ASSERT ( this->isInitialized() ).error(
"epetra matrix not properly initialized" ) ;
773 M_bc_index.resize( 0 );
774 M_mat->PutScalar( 0.0 );
780 void MatrixEpetra::clear ()
782 if ( this->isInitialized() )
784 M_mat = boost::shared_ptr<Epetra_FECrsMatrix>(
new Epetra_FECrsMatrix( Copy, M_emap, M_col_emap, 50 ) );
786 this->setInitialized(
false );
794 void MatrixEpetra::close ()
const
798 if ( !this->closed() )
799 ierr = M_mat->GlobalAssemble( M_dom_map, M_range_map,
true );
803 DVLOG(2) <<
"ERRCODE GlobalAssemble: " << ierr <<
"\n";
809 void MatrixEpetra::setDiagonal ( Epetra_Vector
const& x )
811 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );;
813 const Epetra_Map& rowMap( M_mat->RowMatrixRowMap() );
814 const Epetra_Map& colMap( M_mat->RowMatrixColMap() );
818 value_type zero_value =
static_cast<value_type
>( 0.0 );
822 int i_val =
static_cast<int>( rowMap.GID( i ) );
823 int j_val =
static_cast<int>( colMap.GID( i ) );
825 M_mat->InsertGlobalValues( 1, &i_val, 1, &j_val, &zero_value );
830 M_mat->ReplaceDiagonalValues( x );
838 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );;
840 DVLOG(2) <<
"Size in size1(): " << M_mat->NumGlobalRows() <<
"\n";
842 int epetra_m = M_mat->NumGlobalRows();
844 return static_cast<size_type>( epetra_m );
852 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );;
854 int epetra_n = M_mat->NumGlobalCols();
856 return static_cast<size_type>( epetra_n );
862 size_type MatrixEpetra::rowStart ()
const
864 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );;
866 int start = M_emap.MinMyGID();
877 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );;
879 int stop = M_emap.MaxMyGID();
886 void MatrixEpetra::set (
const size_type i,
888 const value_type& value )
890 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );;
892 int i_val =
static_cast<int>( i );
893 int j_val =
static_cast<int>( j );
897 value_type epetra_value =
static_cast<value_type
>( value );
899 ierr = M_mat->ReplaceGlobalValues( 1, &i_val, 1, &j_val, &epetra_value );
903 ierr = M_mat->InsertGlobalValues( 1, &i_val, 1, &j_val, &epetra_value );
907 DVLOG(2) <<
"ERRORCODE InsertGlobalValues: " << ierr <<
" in M(" << i_val <<
"," << j_val <<
") for value "<< epetra_value <<
"." <<
"\n";
918 bool MatrixEpetra::closed()
const
920 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );
923 filled = M_mat->Filled();
925 return static_cast<bool>( filled );
930 bool MatrixEpetra::EpetraIndicesAreLocal()
const
932 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );
935 IsLocal = M_mat->IndicesAreLocal();
937 return static_cast<bool>( IsLocal );
942 bool MatrixEpetra::HaveColMap()
const
944 FEELPP_ASSERT ( this->isInitialized() ).error(
"MatrixEpetra<> not properly initialized" );
947 HaveMap = M_mat->HaveColMap();
949 return static_cast<bool>( HaveMap );
954 MatrixEpetra::addMatrix (
const value_type coeff, super &X )
956 FEELPP_ASSERT ( this->isInitialized() ).error(
"epetra matrix not initialized" );
958 epetra_sparse_matrix_type* A_ptr =
dynamic_cast< epetra_sparse_matrix_type*
>( &X );
960 this->addMatrix( coeff, *A_ptr );
965 MatrixEpetra::scale (
double scalar )
967 FEELPP_ASSERT ( this->isInitialized() ).error(
"epetra matrix not initialized" );
969 M_mat->Scale( scalar );
973 MatrixEpetra::real_type
974 MatrixEpetra::l1Norm()
const
976 FEELPP_ASSERT ( this->isInitialized() ).error(
"epetra matrix not initialized" );
979 real_type NormOne = M_mat->NormOne();
981 return static_cast<real_type
>( NormOne );
985 MatrixEpetra::real_type
986 MatrixEpetra::linftyNorm()
const
988 FEELPP_ASSERT ( this->isInitialized() ).error(
"epetra matrix not initialized" );
991 real_type NormInf = M_mat->NormInf();
993 return static_cast<real_type
>( NormInf );
1000 MatrixEpetra::value_type
1001 MatrixEpetra::operator () (
const size_type i,
1004 FEELPP_ASSERT ( this->isInitialized() ).error(
"epetra matrix not initialized" );
1007 int i_val=
static_cast<int>( i ),
1008 j_val=static_cast<int>( j );
1016 M_mat->ExtractMyRowView( i_val, NumEntries, Values, Indices );
1018 for (
int k = 0; k < NumEntries; ++k )
1020 if ( Indices[k] == j_val )
1021 return static_cast<double> ( Values[k] );
1074 MatrixEpetra::updateBlockMat( boost::shared_ptr<MatrixSparse<value_type> > m, std::vector<size_type> start_i, std::vector<size_type> start_j )
1076 LOG(ERROR) <<
"Invalid call to updateBlockMat, not yet implemented\n";
1082 operator<<( NdebugStream& __os, MatrixEpetra
const& )