49 #ifndef __sparse_matrix_h__
50 #define __sparse_matrix_h__
56 #include <boost/numeric/ublas/matrix.hpp>
57 #include <boost/fusion/include/fold.hpp>
63 #include <feel/feelcore/feel.hpp>
64 #include <feel/feelcore/traits.hpp>
65 #include <feel/feelcore/context.hpp>
67 #include <feel/feelalg/enums.hpp>
69 #include <feel/feelalg/vector.hpp>
73 namespace ublas = boost::numeric::ublas;
77 template <
typename T>
inline std::ostream& operator << ( std::ostream& os, const MatrixSparse<T>& m );
97 typedef typename type_traits<T>::real_type real_type;
99 typedef boost::shared_ptr<graph_type> graph_ptrtype;
101 typedef boost::shared_ptr<Vector<T> > vector_ptrtype;
104 typedef boost::shared_ptr<datamap_type> datamap_ptrtype;
120 MatrixSparse( datamap_ptrtype
const& dmRow, datamap_ptrtype
const& dmCol, WorldComm
const& worldComm=Environment::worldComm() );
160 void setMapRow( datamap_ptrtype
const& d )
164 void setMapCol( datamap_ptrtype
const& d )
206 graph_ptrtype
const&
graph ) = 0;
212 template<
typename DomainSpace,
typename ImageSpace>
213 void initIndexSplit( DomainSpace
const& dm, ImageSpace
const& im )
215 auto nSpace = DomainSpace::element_type::nSpaces;
220 std::vector < std::vector<int> > is( nSpace );
223 auto result = boost::make_tuple( cptSpaces,start );
225 std::vector < std::vector<int> > indexSplit( nSpace );
228 boost::fusion::fold( dm->functionSpaces(), result, cndof );
230 this->setIndexSplit( indexSplit );
238 virtual void setIndexSplit( std::vector< std::vector<size_type> >
const &_indexSplit )
240 M_IndexSplit=_indexSplit;
246 std::vector< std::vector<size_type> > indexSplit()
const
328 bool haveConsistentProperties()
const
332 return ( p1 ==
false ) && ( p2 == false );
336 return M_mprop.test(
DENSE );
338 void checkProperties()
const
340 if ( !haveConsistentProperties() )
342 std::ostringstream ostr;
343 ostr <<
"Invalid matrix properties:\n"
348 <<
" DENSE: " << isDense() <<
"\n";
349 throw std::logic_error( ostr.str() );
365 virtual void clear () = 0;
370 virtual void zero () = 0;
382 virtual void close ()
const = 0;
416 const value_type& value ) = 0;
428 const value_type& value ) = 0;
436 virtual void addMatrix (
const ublas::matrix<value_type> &dm,
437 const std::vector<size_type> &rows,
438 const std::vector<size_type> &cols ) = 0;
446 virtual void addMatrix (
int* rows,
int nrows,
447 int* cols,
int ncols,
448 value_type* data ) = 0;
454 virtual void addMatrix (
const ublas::matrix<value_type> &dm,
455 const std::vector<size_type> &dof_indices ) = 0;
474 virtual void scale (
const T ) = 0;
488 boost::shared_ptr<
Vector<T> >& dest )
const
523 virtual void diagonal ( Vector<T>& dest )
const = 0;
545 boost::shared_ptr<MatrixSparse<T> > Mt;
594 real_type
energy ( vector_ptrtype
const& v,
595 vector_ptrtype
const& u,
596 bool _transpose =
false )
const
598 return this->
energy( *v, *u, _transpose );
608 virtual real_type
l1Norm ()
const = 0;
626 virtual bool closed()
const = 0;
633 void print( std::ostream& os=std::cout )
const;
639 template <
typename U>
640 friend std::ostream& operator << ( std::ostream& os, const MatrixSparse<U>& m );
648 std::cerr <<
"ERROR: Not Implemented in base class yet!" << std::endl;
649 FEELPP_ASSERT( 0 ).error(
"invalid call" );
660 std::cerr <<
"ERROR: Not Implemented in base class yet!" << std::endl;
661 std::cerr <<
"ERROR writing MATLAB file " << name << std::endl;
662 FEELPP_ASSERT( 0 ).error(
"invalid call" );
671 const std::vector<size_type>& rows,
672 const std::vector<size_type>& cols )
const
687 const std::vector<size_type>& rows,
688 const std::vector<size_type>& cols )
const
719 template<
typename DomainSpace,
typename ImageSpace>
720 void updateIndexSplit( DomainSpace
const& dm, ImageSpace
const& im )
722 auto nSpace = DomainSpace::element_type::nSpaces;
727 std::vector < std::vector<int> > is( nSpace );
730 auto result = boost::make_tuple( cptSpaces,start );
732 std::vector < std::vector<int> > indexSplit( nSpace );
734 detail::computeNDofForEachSpace cndof( indexSplit );
735 boost::fusion::fold( dm->functionSpaces(), result, cndof );
737 this->setIndexSplit( indexSplit );
748 const std::vector<size_type>& ,
749 const std::vector<size_type>& ,
752 std::cerr <<
"Error! This function is not yet implemented in the base class!"
754 FEELPP_ASSERT( 0 ).error(
"invalid call" );
767 graph_ptrtype M_graph;
771 std::vector < std::vector<size_type> > M_IndexSplit;
777 datamap_ptrtype M_mapCol;
783 typedef boost::shared_ptr<d_sparse_matrix_type> d_sparse_matrix_ptrtype;
784 typedef boost::shared_ptr<d_sparse_matrix_type> sparse_matrix_ptrtype;
789 template <
typename T>
792 M_is_initialized( false ),
796 template <
typename T>
799 M_worldComm( worldComm ),
800 M_is_initialized( false ),
806 template <
typename T>
813 template <
typename T>
817 assert ( this->isInitialized() );
819 for (
size_type i=0; i<this->size1(); i++ )
821 for (
size_type j=0; j<this->size2(); j++ )
822 os << std::setw( 8 ) << ( *this )( i,j ) <<
" ";
828 template <
typename T>
833 this->multAddVector( arg,dest );
838 template <
typename T>
856 std::cout <<
"Real part:" << std::endl;
858 for (
size_type i=0; i<this->size1(); i++ )
860 for (
size_type j=0; j<this->size2(); j++ )
861 os << std::setw( 8 ) << ( *this )( i,j ).real() <<
" ";
866 os << std::endl <<
"Imaginary part:" << std::endl;
868 for (
size_type i=0; i<this->size1(); i++ )
870 for (
size_type j=0; j<this->size2(); j++ )
871 os << std::setw( 8 ) << ( *this )( i,j ).imag() <<
" ";
881 template <
typename T>
883 std::ostream& operator << ( std::ostream& os, const MatrixSparse<T>& m )
891 #endif // #ifndef __sparse_matrix_h__
virtual void zeroRows(std::vector< int > const &rows, Vector< value_type > const &values, Vector< value_type > &rhs, Context const &on_context)=0
bool isHermitianPositiveDefinite() const
Definition: matrixsparse.hpp:304
Definition: solverlinear.hpp:33
bool isSingular() const
Definition: matrixsparse.hpp:313
virtual void addVector(int *i, int n, value_type *v)=0
virtual T operator()(const size_type i, const size_type j) const =0
datamap_ptrtype const & mapColPtr() const
Definition: matrixsparse.hpp:155
boost::shared_ptr< MatrixSparse< T > > transpose() const
Definition: matrixsparse.hpp:543
WorldComm const & comm() const
Definition: matrixsparse.hpp:356
void setGraph(graph_ptrtype const &graph)
Definition: matrixsparse.hpp:270
void transpose(boost::shared_ptr< MatrixSparse< value_type > > &Mt) const
Definition: matrixsparse.hpp:556
virtual void init(const size_type m, const size_type n, const size_type m_l, const size_type n_l, const size_type nnz=30, const size_type noz=10)=0
virtual void reinitSubmatrix(MatrixSparse< T > &submatrix, const std::vector< size_type > &rows, const std::vector< size_type > &cols) const
Definition: matrixsparse.hpp:686
void multVector(const boost::shared_ptr< Vector< T > > &arg, boost::shared_ptr< Vector< T > > &dest) const
Definition: matrixsparse.hpp:487
virtual ~MatrixSparse()
Definition: matrixsparse.hpp:808
virtual real_type l1Norm() const =0
bool isHermitian() const
Definition: matrixsparse.hpp:286
MatrixSparse()
Definition: matrixsparse.hpp:791
Definition: feelalg/enums.hpp:68
virtual real_type linftyNorm() const =0
virtual void close() const =0
bool isNonHermitian() const
Definition: matrixsparse.hpp:295
virtual void addMatrix(const ublas::matrix< value_type > &dm, const std::vector< size_type > &rows, const std::vector< size_type > &cols)=0
datamap_ptrtype M_mapRow
Definition: matrixsparse.hpp:776
virtual size_type size1() const =0
virtual real_type energy(vector_type const &v, vector_type const &u, bool transpose=false) const =0
void symmetricPart(boost::shared_ptr< MatrixSparse< value_type > > &Ms) const
Definition: matrixsparse.hpp:569
graph_ptrtype const & graph() const
Definition: matrixsparse.hpp:262
Definition: feelalg/enums.hpp:69
Definition: feelalg/enums.hpp:70
WorldComm M_worldComm
mpi communicator
Definition: matrixsparse.hpp:759
bool hasGraph() const
Definition: matrixsparse.hpp:254
datamap_type const & mapRow() const
Definition: matrixsparse.hpp:131
virtual void diagonal(Vector< T > &dest) const =0
Definition: feelalg/enums.hpp:67
datamap_ptrtype const & mapRowPtr() const
Definition: matrixsparse.hpp:147
virtual void set(const size_type i, const size_type j, const value_type &value)=0
data layout in a multi-processor environnement
Definition: datamap.hpp:46
datamap_type const & mapCol() const
Definition: matrixsparse.hpp:139
virtual void updateSparsityPattern(const std::vector< std::vector< size_type > > &)
Definition: matrixsparse.hpp:183
virtual void printMatlab(const std::string name="NULL") const
Definition: matrixsparse.hpp:658
real_type energy(vector_ptrtype const &v, vector_ptrtype const &u, bool _transpose=false) const
Definition: matrixsparse.hpp:594
virtual size_type rowStop() const =0
virtual void createSubmatrix(MatrixSparse< T > &submatrix, const std::vector< size_type > &rows, const std::vector< size_type > &cols) const
Definition: matrixsparse.hpp:670
Definition: feelalg/enums.hpp:66
bool isPositiveDefinite() const
Definition: matrixsparse.hpp:322
void diagonal(boost::shared_ptr< Vector< T > > &dest) const
Definition: matrixsparse.hpp:528
virtual void updateBlockMat(boost::shared_ptr< MatrixSparse< T > > m, std::vector< size_type > start_i, std::vector< size_type > start_j)=0
virtual bool closed() const =0
void addMatrix(const T &s, boost::shared_ptr< MatrixSparse< T > > &m)
Definition: matrixsparse.hpp:469
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
virtual void symmetricPart(MatrixSparse< value_type > &Ms) const
Definition: matrixsparse.hpp:564
Definition: matrixsparse.hpp:76
virtual bool isInitialized() const
Definition: matrixsparse.hpp:173
void multAddVector(const Vector< T > &arg, Vector< T > &dest) const
Definition: matrixsparse.hpp:839
virtual void _get_submatrix(MatrixSparse< T > &, const std::vector< size_type > &, const std::vector< size_type > &, const bool) const
Definition: matrixsparse.hpp:747
virtual size_type size2() const =0
void multVector(const Vector< T > &arg, Vector< T > &dest) const
Definition: matrixsparse.hpp:829
bool M_is_initialized
Definition: matrixsparse.hpp:765
void setInitialized(bool _init)
Definition: matrixsparse.hpp:714
void setMatrixProperties(size_type p)
Definition: matrixsparse.hpp:278
Context class.
Definition: feelcore/context.hpp:63
void print(std::ostream &os=std::cout) const
Definition: matrixsparse.hpp:815
Graph representation of the Compressed Sparse Row format.
Definition: graphcsr.hpp:57
virtual size_type rowStart() const =0
virtual void printPersonal(std::ostream &=std::cout) const
Definition: matrixsparse.hpp:646
virtual void add(const size_type i, const size_type j, const value_type &value)=0