29 #ifndef __MatrixUBlas_H
30 #define __MatrixUBlas_H 1
33 #include <boost/timer.hpp>
34 #include <boost/numeric/ublas/matrix_sparse.hpp>
35 #include <boost/numeric/ublas/matrix_proxy.hpp>
42 namespace ublas = boost::numeric::ublas;
51 template<
typename T,
typename LayoutType>
62 typedef ublas::compressed_matrix<value_type, LayoutType> matrix_type;
64 typedef typename boost::numeric::bindings::traits::sparse_matrix_traits<matrix_type>::ordering_type ordering_type;
65 typedef typename boost::numeric::bindings::traits::sparse_matrix_traits<matrix_type>::layout_type layout_type;
67 static const bool is_row_major = boost::is_same<ordering_type,
68 boost::numeric::bindings::traits::row_major_t>::value;
71 typedef std::vector<std::set<size_type> > pattern_type;
81 M_is_initialized(
false ),
86 M_is_initialized( m.M_is_initialized ),
102 return M_mat( i, j );
107 return M_mat( i, j );
122 return M_mat.size1();
131 return M_mat.size2();
165 return M_is_initialized;
171 matrix_type
const&
mat ()
const
207 void init (
const unsigned int m,
208 const unsigned int n,
209 const unsigned int m_l,
210 const unsigned int n_l,
211 const unsigned int nnz=30,
212 const unsigned int noz=10 );
231 M_mat = ublas::zero_matrix<value_type>( M_mat.size1(), M_mat.size2() );
236 ublas::subrange( M_mat, start1, stop1, start2, stop2 ) = ublas::zero_matrix<value_type>( stop1-start1, stop2-start2 );
247 void add (
const unsigned int i,
248 const unsigned int j,
249 const value_type value )
251 M_mat( i, j ) += value;
262 void set (
const unsigned int i,
263 const unsigned int j,
264 const value_type value )
266 M_mat( i, j ) = value;
276 void printMatlab(
const std::string name=
"NULL" )
const;
282 void fill( pattern_type
const& );
286 M_mat.resize( nr, nc, preserve );
292 template<
typename VE1,
typename VE2>
294 energy( ublas::vector_expression<VE1>
const& __v,
295 ublas::vector_expression<VE2>
const& __u )
const
297 return ublas::inner_prod( __v, ublas::prod( M_mat, __u ) );
312 bool M_is_initialized;
320 template <
typename T,
typename LayoutType>
322 const unsigned int n,
328 if ( ( m==0 ) || ( n==0 ) )
331 M_mat.resize( m,n,
false );
335 template<
typename T,
typename LayoutType>
340 ublas::matrix_row<matrix_type> mr ( M_mat, __dof_index );
341 typedef typename ublas::matrix_row<matrix_type>::iterator r_it;
343 for ( r_it __r = mr.begin(); __r != mr.end(); ++__r )
349 ublas::matrix_column<matrix_type> mc ( M_mat, __dof_index );
351 for (
typename ublas::matrix_column<matrix_type>::iterator therow = mc.begin();
352 therow != mc.end(); ++therow )
358 M_mat( __dof_index, __dof_index ) = 1.0;
360 template<
typename T,
typename LayoutType>
364 namespace bindings = boost::numeric::bindings;
367 typename pattern_type::const_iterator __itl = __pattern.begin();
370 std::for_each( __pattern.begin(),
372 ( boost::lambda::var( __nnz ) += boost::lambda::bind( &std::set<size_type>::size,
373 boost::lambda::_1 ) ) );
376 for (
size_type __i = 0; __i < __pattern.size(); ++__i )
378 __nnz += __pattern[__i].size();
389 FEELPP_ASSERT( __nnz >= M_mat.nnz() )( __nnz )( M_mat.nnz() ).error(
"incompatible sizes" );
391 DVLOG(2) <<
"number of nnz in old M : " << M_mat.nnz() <<
", " << M_mat.nnz_capacity() <<
"\n";
392 DVLOG(2) <<
"size M.value_data() : " << M_mat.value_data().size() <<
"\n";
401 matrix_type M_mat_backup( M_mat );
403 DVLOG(2) <<
"resizing M old : " << M_mat_backup.size1() <<
"," << M_mat_backup.size2() <<
" nnz = " << M_mat_backup.nnz() <<
"\n";
405 M_mat.reserve( __nnz,
false );
407 DVLOG(2) <<
"resizing M new : " << M_mat.size1() <<
"," << M_mat.size2() <<
" nnz = " << M_mat.nnz() <<
"\n";
409 std::set<size_type>::const_iterator __it;
410 std::set<size_type>::const_iterator __en;
412 DVLOG(2) <<
"*** counting Nnz : " << chrono.elapsed() <<
"\n";
414 DVLOG(2) <<
"*** nnz size : " << __nnz <<
"\n";
417 uint32_type __max_nnz_per_line = 0;
419 uint __nnz_entry = 0;
427 thesize = M_mat.size1();
430 thesize = M_mat.size2();
432 while ( __row < thesize )
434 __it = __pattern[__row].begin();
435 __en = __pattern[__row].end();
437 M_mat.index1_data()[__row] = __nnz_entry;
440 uint32_type __nnz_line = 0;
442 while ( __it != __en )
444 M_mat.index2_data()[__nnz_entry] = *__it;
448 namespace bindings = boost::numeric::bindings;
450 typename matrix_type::value_type* __pv = 0;
452 if ( boost::is_same<
typename bindings::traits::sparse_matrix_traits<matrix_type>::ordering_type,
453 bindings::traits::row_major_t>::value )
455 if ( __row < M_mat_backup.size1() && *__it < M_mat_backup.size2() )
456 __pv = M_mat_backup.find_element( __row, *__it );
459 else if ( boost::is_same<
typename bindings::traits::sparse_matrix_traits<matrix_type>::ordering_type,
460 bindings::traits::column_major_t>::value )
462 if ( __row < M_mat_backup.size2() && *__it < M_mat_backup.size1() )
463 __pv = M_mat_backup.find_element( *__it, __row );
468 std::cout <<
"ERROR " << __FILE__ <<
": " << __LINE__ <<
"\n";
472 M_mat.value_data()[__nnz_entry] = *__pv;
475 M_mat.value_data()[__nnz_entry] = value_type( 0 );
482 __max_nnz_per_line = std::max( __nnz_line, __max_nnz_per_line );
486 M_mat.index1_data()[thesize] = __filled2;
487 FEELPP_ASSERT( thesize+1 == __filled1 )( thesize )( __filled1 ).error(
"invalid matrix storage" );
488 M_mat.set_filled( __filled1, __filled2 );
489 FEELPP_ASSERT( M_mat.nnz() == __filled2 )( M_mat.nnz() )( __filled2 ).error(
"inconsistent matrix storage" );
491 DVLOG(2) <<
"*** value data size : " << M_mat.value_data().size() <<
"\n";
492 DVLOG(2) <<
"*** nnz : " << M_mat.nnz() <<
"\n";
493 DVLOG(2) <<
"*** max nnz per line : " << __max_nnz_per_line <<
"\n";
494 DVLOG(2) <<
"*** fillMatrixFromPattern() done in " << chrono.elapsed() <<
"s\n";
497 template<
typename T,
typename LayoutType>
501 std::string name = filename;
502 std::string separator =
" , ";
505 int i = filename.find(
"." );
508 name = filename +
".m";
512 if ( (
unsigned int ) i != filename.size() - 2 ||
513 filename[ i + 1 ] !=
'm' )
515 std::cerr <<
"Wrong file name ";
516 name = filename +
".m";
520 std::ofstream file_out( name.c_str() );
522 FEELPP_ASSERT( file_out )( filename ).error(
"[Feel::spy] ERROR: File cannot be opened for writing." );
524 file_out <<
"S = [ ";
526 for (
typename matrix_type::const_iterator1 i1=M_mat.begin1();
527 i1!=M_mat.end1(); ++i1 )
529 for (
typename matrix_type::const_iterator2 i2=i1.begin();
531 file_out << i2.index1() + 1 << separator
532 << i2.index2() + 1 << separator
536 file_out <<
"];" << std::endl;
537 file_out <<
"I=S(:,1); J=S(:,2); S=S(:,3);" << std::endl;
538 file_out <<
"A=sparse(I,J,S); spy(A);" << std::endl;