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
matrixsparse.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-27
7 
8  Copyright (C) 2005,2006 EPFL
9  Copyright (C) 2007-2011 Université Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 3.0 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 // $Id: sparse_matrix.h,v 1.12 2005/05/24 12:54:57 jwpeterson Exp $
31 
32 // The libMesh Finite Element Library.
33 // Copyright (C) 2002-2005 Benjamin S. Kirk, John W. Peterson
34 
35 // This library is free software; you can redistribute it and/or
36 // modify it under the terms of the GNU Lesser General Public
37 // License as published by the Free Software Foundation; either
38 // version 3.0 of the License, or (at your option) any later version.
39 
40 // This library is distributed in the hope that it will be useful,
41 // but WITHOUT ANY WARRANTY; without even the implied warranty of
42 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
43 // Lesser General Public License for more details.
44 
45 // You should have received a copy of the GNU Lesser General Public
46 // License along with this library; if not, write to the Free Software
47 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
48 
49 #ifndef __sparse_matrix_h__
50 #define __sparse_matrix_h__
51 
52 // C++ includes
53 #include <iomanip>
54 #include <vector>
55 
56 #include <boost/numeric/ublas/matrix.hpp>
57 #include <boost/fusion/include/fold.hpp>
58 
59 #include <Eigen/Core>
60 
61 
62 
63 #include <feel/feelcore/feel.hpp>
64 #include <feel/feelcore/traits.hpp>
65 #include <feel/feelcore/context.hpp>
66 
67 #include <feel/feelalg/enums.hpp>
69 #include <feel/feelalg/vector.hpp>
70 
71 namespace Feel
72 {
73 namespace ublas = boost::numeric::ublas;
74 
75 // forward declarations
76 template <typename T> class MatrixSparse;
77 template <typename T> inline std::ostream& operator << ( std::ostream& os, const MatrixSparse<T>& m );
78 
79 
91 template <typename T>
92 class MatrixSparse
93 {
94 public:
95 
96  typedef T value_type;
97  typedef typename type_traits<T>::real_type real_type;
98  typedef GraphCSR graph_type;
99  typedef boost::shared_ptr<graph_type> graph_ptrtype;
100  typedef Vector<T> vector_type;
101  typedef boost::shared_ptr<Vector<T> > vector_ptrtype;
102 
103  typedef DataMap datamap_type;
104  typedef boost::shared_ptr<datamap_type> datamap_ptrtype;
105 
118  MatrixSparse ();
119 
120  MatrixSparse( datamap_ptrtype const& dmRow, datamap_ptrtype const& dmCol, WorldComm const& worldComm=Environment::worldComm() );
121 
126  virtual ~MatrixSparse ();
127 
131  datamap_type const& mapRow() const
132  {
133  return *M_mapRow;
134  }
135 
139  datamap_type const& mapCol() const
140  {
141  return *M_mapCol;
142  }
143 
147  datamap_ptrtype const& mapRowPtr() const
148  {
149  return M_mapRow;
150  }
151 
155  datamap_ptrtype const& mapColPtr() const
156  {
157  return M_mapCol;
158  }
159 
160  void setMapRow( datamap_ptrtype const& d )
161  {
162  M_mapRow=d;
163  }
164  void setMapCol( datamap_ptrtype const& d )
165  {
166  M_mapCol=d;
167  }
168 
173  virtual bool isInitialized() const
174  {
175  return M_is_initialized;
176  }
177 
183  virtual void updateSparsityPattern ( const std::vector<std::vector<size_type> >& ) {}
184 
192  virtual void init ( const size_type m,
193  const size_type n,
194  const size_type m_l,
195  const size_type n_l,
196  const size_type nnz=30,
197  const size_type noz=10 ) = 0;
198 
202  virtual void init ( const size_type m,
203  const size_type n,
204  const size_type m_l,
205  const size_type n_l,
206  graph_ptrtype const& graph ) = 0;
207 
208 #if 0
209 
212  template<typename DomainSpace, typename ImageSpace>
213  void initIndexSplit( DomainSpace const& dm, ImageSpace const& im )
214  {
215  auto nSpace = DomainSpace::element_type::nSpaces;
216 
217  if ( nSpace>1 )
218  {
219  //std::cout << "\n Debug : nSpace " << nSpace << "\n";
220  std::vector < std::vector<int> > is( nSpace );
221  uint cptSpaces=0;
222  uint start=0;
223  auto result = boost::make_tuple( cptSpaces,start );
224 
225  std::vector < std::vector<int> > indexSplit( nSpace );
226  //detail::computeNDofForEachSpace cndof(nSpace);
227  //detail::computeNDofForEachSpace cndof(indexSplit);
228  boost::fusion::fold( dm->functionSpaces(), result, cndof );
229 
230  this->setIndexSplit( indexSplit );
231  }
232 
233  }
234 #endif
235 
238  virtual void setIndexSplit( std::vector< std::vector<size_type> > const &_indexSplit )
239  {
240  M_IndexSplit=_indexSplit;
241  }
242 
246  std::vector< std::vector<size_type> > indexSplit() const
247  {
248  return M_IndexSplit;
249  }
250 
254  bool hasGraph() const
255  {
256  return M_graph != 0;
257  }
258 
262  graph_ptrtype const& graph() const
263  {
264  return M_graph;
265  }
266 
270  void setGraph( graph_ptrtype const& graph )
271  {
272  M_graph = graph;
273  }
274 
279  {
280  M_mprop = ( size_type )p;
281  checkProperties();
282  }
286  bool isHermitian() const
287  {
288  checkProperties();
289  return M_mprop.test( HERMITIAN );
290  }
291 
295  bool isNonHermitian() const
296  {
297  checkProperties();
298  return M_mprop.test( NON_HERMITIAN );
299  }
300 
305  {
306  checkProperties();
307  return M_mprop.test( HERMITIAN | POSITIVE_DEFINITE );
308  }
309 
313  bool isSingular() const
314  {
315  checkProperties();
316  return M_mprop.test( SINGULAR );
317  }
318 
322  bool isPositiveDefinite() const
323  {
324  checkProperties();
325  return M_mprop.test( POSITIVE_DEFINITE );
326  }
327 
328  bool haveConsistentProperties() const
329  {
330  bool p1 = M_mprop.test( SINGULAR ) && M_mprop.test( POSITIVE_DEFINITE );
331  bool p2 = M_mprop.test( HERMITIAN ) && M_mprop.test( NON_HERMITIAN );
332  return ( p1 == false ) && ( p2 == false );
333  }
334  bool isDense() const
335  {
336  return M_mprop.test( DENSE );
337  }
338  void checkProperties() const
339  {
340  if ( !haveConsistentProperties() )
341  {
342  std::ostringstream ostr;
343  ostr << "Invalid matrix properties:\n"
344  << " HERMITIAN: " << isHermitian() << "\n"
345  << " NON_HERMITIAN: " << isNonHermitian() << "\n"
346  << " SINGULAR: " << isSingular() << "\n"
347  << " POSITIVE_DEFINITE: " << isPositiveDefinite() << "\n"
348  << " DENSE: " << isDense() << "\n";
349  throw std::logic_error( ostr.str() );
350  }
351  }
355  //mpi::communicator const& comm() const { return M_comm; }
356  WorldComm const& comm() const
357  {
358  return M_worldComm;
359  }
360 
365  virtual void clear () = 0;
366 
370  virtual void zero () = 0;
371 
375  virtual void zero ( size_type start1, size_type size1,
376  size_type start2, size_type size2 ) = 0;
377 
382  virtual void close () const = 0;
383 
388  virtual size_type size1 () const = 0;
389 
394  virtual size_type size2 () const = 0;
395 
400  virtual size_type rowStart () const = 0;
401 
406  virtual size_type rowStop () const = 0;
407 
414  virtual void set ( const size_type i,
415  const size_type j,
416  const value_type& value ) = 0;
417 
426  virtual void add ( const size_type i,
427  const size_type j,
428  const value_type& value ) = 0;
429 
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;
439 
446  virtual void addMatrix ( int* rows, int nrows,
447  int* cols, int ncols,
448  value_type* data ) = 0;
449 
454  virtual void addMatrix ( const ublas::matrix<value_type> &dm,
455  const std::vector<size_type> &dof_indices ) = 0;
456 
462  virtual void addMatrix ( const T, MatrixSparse<T> & ) = 0;
463 
469  void addMatrix ( const T& s, boost::shared_ptr<MatrixSparse<T> > & m )
470  {
471  this->addMatrix( s, *m );
472  }
473 
474  virtual void scale ( const T ) = 0;
475 
480  void multVector ( const Vector<T>& arg,
481  Vector<T>& dest ) const;
482 
487  void multVector ( const boost::shared_ptr<Vector<T> >& arg,
488  boost::shared_ptr<Vector<T> >& dest ) const
489  {
490  this->multVector( *arg, *dest );
491  }
492 
496  void multAddVector ( const Vector<T>& arg,
497  Vector<T>& dest ) const;
498 
499 
511  virtual T operator () ( const size_type i,
512  const size_type j ) const = 0;
513 
514  virtual MatrixSparse<T> & operator = ( MatrixSparse<value_type> const& M ) = 0;
515  MatrixSparse<T> & operator = ( boost::shared_ptr<MatrixSparse<value_type> > const& M )
516  {
517  *this = *M;
518  return *this;
519  }
523  virtual void diagonal ( Vector<T>& dest ) const = 0;
524 
528  void diagonal ( boost::shared_ptr<Vector<T> >& dest ) const
529  {
530  diagonal( *dest );
531  }
532 
538  virtual void transpose( MatrixSparse<value_type>& Mt ) const = 0;
539 
543  boost::shared_ptr<MatrixSparse<T> > transpose() const
544  {
545  boost::shared_ptr<MatrixSparse<T> > Mt;
546  transpose( *Mt );
547  return Mt;
548  }
549 
556  void transpose( boost::shared_ptr<MatrixSparse<value_type> >& Mt ) const
557  {
558  this->transpose( *Mt );
559  }
560 
564  virtual void symmetricPart( MatrixSparse<value_type>& Ms ) const {}
565 
569  void symmetricPart( boost::shared_ptr<MatrixSparse<value_type> >& Ms ) const
570  {
571  symmetricPart( *Ms );
572  }
573 
582  virtual real_type energy ( vector_type const& v,
583  vector_type const& u,
584  bool transpose = false ) const = 0;
585 
594  real_type energy ( vector_ptrtype const& v,
595  vector_ptrtype const& u,
596  bool _transpose = false ) const
597  {
598  return this->energy( *v, *u, _transpose );
599  }
600 
608  virtual real_type l1Norm () const = 0;
609 
620  virtual real_type linftyNorm () const = 0;
621 
626  virtual bool closed() const = 0;
627 
633  void print( std::ostream& os=std::cout ) const;
634 
639  template <typename U>
640  friend std::ostream& operator << ( std::ostream& os, const MatrixSparse<U>& m );
641 
646  virtual void printPersonal( std::ostream& /*os*/=std::cout ) const
647  {
648  std::cerr << "ERROR: Not Implemented in base class yet!" << std::endl;
649  FEELPP_ASSERT( 0 ).error( "invalid call" );
650  }
651 
658  virtual void printMatlab( const std::string name="NULL" ) const
659  {
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" );
663  }
664 
670  virtual void createSubmatrix( MatrixSparse<T>& submatrix,
671  const std::vector<size_type>& rows,
672  const std::vector<size_type>& cols ) const
673  {
674  this->_get_submatrix( submatrix,
675  rows,
676  cols,
677  false ); // false means DO NOT REUSE submatrix
678  }
679 
686  virtual void reinitSubmatrix( MatrixSparse<T>& submatrix,
687  const std::vector<size_type>& rows,
688  const std::vector<size_type>& cols ) const
689  {
690  this->_get_submatrix( submatrix,
691  rows,
692  cols,
693  true ); // true means REUSE submatrix
694  }
695 
703  virtual void zeroRows( std::vector<int> const& rows, Vector<value_type> const& values, Vector<value_type>& rhs, Context const& on_context ) = 0;
704 
708  virtual void updateBlockMat( boost::shared_ptr<MatrixSparse<T> > m, std::vector<size_type> start_i, std::vector<size_type> start_j ) = 0;
709 
710 
714  void setInitialized( bool _init )
715  {
716  M_is_initialized = _init;
717  }
718 #if 0
719  template<typename DomainSpace, typename ImageSpace>
720  void updateIndexSplit( DomainSpace const& dm, ImageSpace const& im )
721  {
722  auto nSpace = DomainSpace::element_type::nSpaces;
723 
724  if ( nSpace>1 )
725  {
726  //std::cout << "\n Debug : nSpace " << nSpace << "\n";
727  std::vector < std::vector<int> > is( nSpace );
728  uint cptSpaces=0;
729  uint start=0;
730  auto result = boost::make_tuple( cptSpaces,start );
731 
732  std::vector < std::vector<int> > indexSplit( nSpace );
733  //detail::computeNDofForEachSpace cndof(nSpace);
734  detail::computeNDofForEachSpace cndof( indexSplit );
735  boost::fusion::fold( dm->functionSpaces(), result, cndof );
736 
737  this->setIndexSplit( indexSplit );
738  }
739  }
740 #endif
741 protected:
748  const std::vector<size_type>& ,
749  const std::vector<size_type>& ,
750  const bool ) const
751  {
752  std::cerr << "Error! This function is not yet implemented in the base class!"
753  << std::endl;
754  FEELPP_ASSERT( 0 ).error( "invalid call" );
755  }
756 
758  //mpi::communicator M_comm;
759  WorldComm M_worldComm;
760 
766 
767  graph_ptrtype M_graph;
768 
769  Context M_mprop;
770 
771  std::vector < std::vector<size_type> > M_IndexSplit;
772 
776  datamap_ptrtype M_mapRow;
777  datamap_ptrtype M_mapCol;
778 
779 
780 };
781 
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;
785 
786 
787 //-----------------------------------------------------------------------
788 // MatrixSparse inline members
789 template <typename T>
790 inline
792  M_is_initialized( false ),
793  M_mprop( NON_HERMITIAN )
794 {}
795 
796 template <typename T>
797 inline
798 MatrixSparse<T>::MatrixSparse ( datamap_ptrtype const& dmRow, datamap_ptrtype const& dmCol, WorldComm const& worldComm ) :
799  M_worldComm( worldComm ),
800  M_is_initialized( false ),
801  M_mprop( NON_HERMITIAN ),
802  M_mapRow( dmRow ),
803  M_mapCol( dmCol )
804 {}
805 
806 template <typename T>
807 inline
809 {}
810 
811 
812 
813 template <typename T>
814 inline
815 void MatrixSparse<T>::print( std::ostream& os ) const
816 {
817  assert ( this->isInitialized() );
818 
819  for ( size_type i=0; i<this->size1(); i++ )
820  {
821  for ( size_type j=0; j<this->size2(); j++ )
822  os << std::setw( 8 ) << ( *this )( i,j ) << " ";
823 
824  os << std::endl;
825  }
826 }
827 
828 template <typename T>
830  Vector<T>& dest ) const
831 {
832  dest.zero();
833  this->multAddVector( arg,dest );
834 }
835 
836 
837 
838 template <typename T>
840  Vector<T>& dest ) const
841 {
842  /* This functionality is actually implemented in the \p
843  Vector class. */
844  dest.addVector( arg,*this );
845 }
846 
847 
848 #if 0
849 // Full specialization for Complex datatypes
850 template <>
851 inline
852 void MatrixSparse<std::complex<double> >::print( std::ostream& os ) const
853 {
854  // std::complex<>::operator<<() is defined, but use this form
855 
856  std::cout << "Real part:" << std::endl;
857 
858  for ( size_type i=0; i<this->size1(); i++ )
859  {
860  for ( size_type j=0; j<this->size2(); j++ )
861  os << std::setw( 8 ) << ( *this )( i,j ).real() << " ";
862 
863  os << std::endl;
864  }
865 
866  os << std::endl << "Imaginary part:" << std::endl;
867 
868  for ( size_type i=0; i<this->size1(); i++ )
869  {
870  for ( size_type j=0; j<this->size2(); j++ )
871  os << std::setw( 8 ) << ( *this )( i,j ).imag() << " ";
872 
873  os << std::endl;
874  }
875 }
876 #endif
877 
878 
879 // For SGI MIPSpro this implementation must occur after
880 // the partial specialization of the print() member.
881 template <typename T>
882 inline
883 std::ostream& operator << ( std::ostream& os, const MatrixSparse<T>& m )
884 {
885  m.print( os );
886  return os;
887 }
888 
889 } // Feel
890 
891 #endif // #ifndef __sparse_matrix_h__

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