Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
matrixublas.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-11-13
7 
8  Copyright (C) 2005,2006 EPFL
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __MatrixUBlas_H
30 #define __MatrixUBlas_H 1
31 
32 #include <set>
33 #include <boost/timer.hpp>
34 #include <boost/numeric/ublas/matrix_sparse.hpp>
35 #include <boost/numeric/ublas/matrix_proxy.hpp>
36 //#include <boost/numeric/bindings/traits/traits.hpp>
37 //#include <boost/numeric/bindings/traits/ublas_sparse.hpp>
38 
39 
40 namespace Feel
41 {
42 namespace ublas = boost::numeric::ublas;
43 
51 template<typename T, typename LayoutType>
53 {
54 public:
55 
56 
60 
61  typedef T value_type;
62  typedef ublas::compressed_matrix<value_type, LayoutType> matrix_type;
63 
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;
66 
67  static const bool is_row_major = boost::is_same<ordering_type,
68  boost::numeric::bindings::traits::row_major_t>::value;
69 
70 
71  typedef std::vector<std::set<size_type> > pattern_type;
72 
74 
78 
79  MatrixUBlas()
80  :
81  M_is_initialized( false ),
82  M_mat()
83  {}
84  MatrixUBlas( MatrixUBlas const & m )
85  :
86  M_is_initialized( m.M_is_initialized ),
87  M_mat( m.M_mat )
88  {}
89 
90  ~MatrixUBlas()
91  {
92  }
93 
95 
99 
100  value_type& operator()( size_type i, size_type j )
101  {
102  return M_mat( i, j );
103  }
104 
105  value_type const& operator()( size_type i, size_type j ) const
106  {
107  return M_mat( i, j );
108  }
109 
111 
115 
120  unsigned int size1 () const
121  {
122  return M_mat.size1();
123  }
124 
129  unsigned int size2 () const
130  {
131  return M_mat.size2();
132  }
133 
137  size_type nnz() const
138  {
139  return M_mat.nnz();
140  }
141 
146  unsigned int rowStart () const
147  {
148  return 0;
149  }
150 
155  unsigned int rowStop () const
156  {
157  return 0;
158  }
159 
163  bool isInitialized() const
164  {
165  return M_is_initialized;
166  }
167 
171  matrix_type const& mat () const
172  {
173  return M_mat;
174  }
175 
179  matrix_type & mat ()
180  {
181  return M_mat;
182  }
183 
184 
185 
187 
191 
192 
194 
198 
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 );
213 
220  void clear ()
221  {
222  M_mat.clear();
223  }
224 
229  void zero ()
230  {
231  M_mat = ublas::zero_matrix<value_type>( M_mat.size1(), M_mat.size2() );
232  }
233 
234  void zero ( size_type start1, size_type stop1, size_type start2, size_type stop2 )
235  {
236  ublas::subrange( M_mat, start1, stop1, start2, stop2 ) = ublas::zero_matrix<value_type>( stop1-start1, stop2-start2 );
237  }
238 
247  void add ( const unsigned int i,
248  const unsigned int j,
249  const value_type value )
250  {
251  M_mat( i, j ) += value;
252  }
253 
262  void set ( const unsigned int i,
263  const unsigned int j,
264  const value_type value )
265  {
266  M_mat( i, j ) = value;
267  }
268 
269 
276  void printMatlab( const std::string name="NULL" ) const;
277 
278 
282  void fill( pattern_type const& );
283 
284  void resize( size_type nr, size_type nc, bool preserve = false )
285  {
286  M_mat.resize( nr, nc, preserve );
287  }
288 
292  template<typename VE1, typename VE2>
293  value_type
294  energy( ublas::vector_expression<VE1> const& __v,
295  ublas::vector_expression<VE2> const& __u ) const
296  {
297  return ublas::inner_prod( __v, ublas::prod( M_mat, __u ) );
298  }
299 
303  void diagonalize( size_type );
305 
306 
307 
308 protected:
309 
310 private:
311 
312  bool M_is_initialized;
313 
317  matrix_type M_mat;
318 
319 };
320 template <typename T, typename LayoutType>
321 void MatrixUBlas<T,LayoutType>::init ( const unsigned int m,
322  const unsigned int n,
323  const unsigned int /*m_l*/,
324  const unsigned int /*n_l*/,
325  const unsigned int /*nnz*/,
326  const unsigned int /*noz*/ )
327 {
328  if ( ( m==0 ) || ( n==0 ) )
329  return;
330 
331  M_mat.resize( m,n,false );
332  this->zero ();
333 }
334 
335 template<typename T, typename LayoutType>
336 void
338 {
339  // eliminating row
340  ublas::matrix_row<matrix_type> mr ( M_mat, __dof_index );
341  typedef typename ublas::matrix_row<matrix_type>::iterator r_it;
342 
343  for ( r_it __r = mr.begin(); __r != mr.end(); ++__r )
344  {
345  *__r = 0.0;
346  }
347 
348  // eliminating column
349  ublas::matrix_column<matrix_type> mc ( M_mat, __dof_index );
350 
351  for ( typename ublas::matrix_column<matrix_type>::iterator therow = mc.begin();
352  therow != mc.end(); ++therow )
353  {
354  *therow = 0.0;
355  }
356 
357  // 1 on the diagonal
358  M_mat( __dof_index, __dof_index ) = 1.0;
359 }
360 template<typename T, typename LayoutType>
361 void
362 MatrixUBlas<T, LayoutType>::fill( pattern_type const& __pattern )
363 {
364  namespace bindings = boost::numeric::bindings;
365 
366  boost::timer chrono;
367  typename pattern_type::const_iterator __itl = __pattern.begin();
368  size_type __nnz = 0;
369 #if 0
370  std::for_each( __pattern.begin(),
371  __pattern.end(),
372  ( boost::lambda::var( __nnz ) += boost::lambda::bind( &std::set<size_type>::size,
373  boost::lambda::_1 ) ) );
374 #else
375 
376  for ( size_type __i = 0; __i < __pattern.size(); ++__i )
377  {
378  __nnz += __pattern[__i].size();
379  // std::cout << "line " << __i << "\n";
380  // std::for_each( __pattern[__i].begin(), __pattern[__i].end(), std::cout << lambda::_1 << " " );
381  // std::cout << "\n";
382  }
383 
384 #endif
385 
386  // ublas::unbounded_array<value_type> __val( __nnz );
387  // std::for_each( __val.begin(), __val.end(), boost::lambda::_1 = 0.0 );
388 
389  FEELPP_ASSERT( __nnz >= M_mat.nnz() )( __nnz )( M_mat.nnz() ).error( "incompatible sizes" );
390 
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";
393  // DVLOG(2) << " size val : " << __val.size() << "\n";
394  // save current nonzero entrie of M in the vector val
395  // std::copy( M_mat.value_data().begin(), M_mat.value_data().begin()+M_mat.nnz(), __val.begin() );
396 
397  //std::cout << "values=";
398  // std::for_each( __val.begin(), __val.end(), std::cout << boost::lambda::_1 << "\n" );
399  //std::cout << "\n";
400 
401  matrix_type M_mat_backup( M_mat );
402 
403  DVLOG(2) << "resizing M old : " << M_mat_backup.size1() << "," << M_mat_backup.size2() << " nnz = " << M_mat_backup.nnz() <<"\n";
404 
405  M_mat.reserve( __nnz, false );
406 
407  DVLOG(2) << "resizing M new : " << M_mat.size1() << "," << M_mat.size2() << " nnz = " << M_mat.nnz() <<"\n";
408 
409  std::set<size_type>::const_iterator __it;
410  std::set<size_type>::const_iterator __en;
411 
412  DVLOG(2) << "*** counting Nnz : " << chrono.elapsed() <<"\n";
413  //DVLOG(2) << "*** l size : " << __pattern.size() <<"\n";
414  DVLOG(2) << "*** nnz size : " << __nnz <<"\n";
415  chrono.restart();
416 
417  uint32_type __max_nnz_per_line = 0;
418  uint __row = 0;
419  uint __nnz_entry = 0;
420 
421  size_type __filled1 = 1;
422  size_type __filled2 = 0;
423 
424  size_type thesize;
425 
426  if ( is_row_major )
427  thesize = M_mat.size1();
428 
429  else
430  thesize = M_mat.size2();
431 
432  while ( __row < thesize )
433  {
434  __it = __pattern[__row].begin();
435  __en = __pattern[__row].end();
436 
437  M_mat.index1_data()[__row] = __nnz_entry;
438 
439  ++__filled1;
440  uint32_type __nnz_line = 0;
441 
442  while ( __it != __en )
443  {
444  M_mat.index2_data()[__nnz_entry] = *__it;
445 
446  ++__filled2;
447 
448  namespace bindings = boost::numeric::bindings;
449 
450  typename matrix_type::value_type* __pv = 0;
451 
452  if ( boost::is_same<typename bindings::traits::sparse_matrix_traits<matrix_type>::ordering_type,
453  bindings::traits::row_major_t>::value )
454  {
455  if ( __row < M_mat_backup.size1() && *__it < M_mat_backup.size2() )
456  __pv = M_mat_backup.find_element( __row, *__it );
457  }
458 
459  else if ( boost::is_same<typename bindings::traits::sparse_matrix_traits<matrix_type>::ordering_type,
460  bindings::traits::column_major_t>::value )
461  {
462  if ( __row < M_mat_backup.size2() && *__it < M_mat_backup.size1() )
463  __pv = M_mat_backup.find_element( *__it, __row );
464  }
465 
466  else
467  {
468  std::cout << "ERROR " << __FILE__ << ": " << __LINE__ << "\n";
469  }
470 
471  if ( __pv )
472  M_mat.value_data()[__nnz_entry] = *__pv;
473 
474  else
475  M_mat.value_data()[__nnz_entry] = value_type( 0 );
476 
477  ++__nnz_entry;
478  ++__it;
479  ++__nnz_line;
480  }
481 
482  __max_nnz_per_line = std::max( __nnz_line, __max_nnz_per_line );
483  ++__row;
484  }
485 
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" );
490 
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";
495 }
496 
497 template<typename T, typename LayoutType>
498 void
499 MatrixUBlas<T, LayoutType>::printMatlab( const std::string filename ) const
500 {
501  std::string name = filename;
502  std::string separator = " , ";
503 
504  // check on the file name
505  int i = filename.find( "." );
506 
507  if ( i <= 0 )
508  name = filename + ".m";
509 
510  else
511  {
512  if ( ( unsigned int ) i != filename.size() - 2 ||
513  filename[ i + 1 ] != 'm' )
514  {
515  std::cerr << "Wrong file name ";
516  name = filename + ".m";
517  }
518  }
519 
520  std::ofstream file_out( name.c_str() );
521 
522  FEELPP_ASSERT( file_out )( filename ).error( "[Feel::spy] ERROR: File cannot be opened for writing." );
523 
524  file_out << "S = [ ";
525 
526  for ( typename matrix_type::const_iterator1 i1=M_mat.begin1();
527  i1!=M_mat.end1(); ++i1 )
528  {
529  for ( typename matrix_type::const_iterator2 i2=i1.begin();
530  i2!=i1.end(); ++i2 )
531  file_out << i2.index1() + 1 << separator
532  << i2.index2() + 1 << separator
533  << *i2 << std::endl;
534  }
535 
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;
539 }
540 
541 } // Feel
542 #endif /* __MatrixUBlas_H */

Generated on Fri Oct 25 2013 14:24:18 for Feel++ by doxygen 1.8.4