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__