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
matrixblock.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): Vincent Chabannes <vincent.chabannes@imag.fr>
6  Date: 2008-01-03
7 
8  Copyright (C) 2011 Université Joseph Fourier (Grenoble I)
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 */
30 #ifndef __MatrixBlock_H
31 #define __MatrixBlock_H 1
32 
33 
34 #include <boost/spirit/home/phoenix.hpp>
35 #include <boost/spirit/home/phoenix/core/argument.hpp>
36 
38 #include <feel/feelalg/backend.hpp>
39 #include <feel/feelvf/block.hpp>
40 
41 
42 namespace Feel
43 {
44 
45 
46 template<typename T> class Backend;
47 
48 template <typename T=double>
49 class BlocksBaseSparseMatrix : public vf::BlocksBase<boost::shared_ptr<MatrixSparse<T> > >
50 {
51 public :
52  typedef vf::BlocksBase<boost::shared_ptr<MatrixSparse<T> > > super_type;
53  typedef BlocksBaseSparseMatrix<T> self_type;
54  typedef boost::shared_ptr<MatrixSparse<T> > matrix_sparse_ptrtype;
55 
56  BlocksBaseSparseMatrix(uint16_type nr,uint16_type nc)
57  :
58  super_type(nr,nc)
59  {}
60 
61  BlocksBaseSparseMatrix(super_type const & b)
62  :
63  super_type(b)
64  {}
65 
66  self_type
67  operator<<( matrix_sparse_ptrtype const& m ) const
68  {
69  return super_type::operator<<( m );
70  }
71 
72 };
73 
74 template <int NR, int NC, typename T=double>
75 class BlocksSparseMatrix : public BlocksBaseSparseMatrix<T>
76 {
77 public :
78  static const uint16_type NBLOCKROWS = NR;
79  static const uint16_type NBLOCKCOLS = NC;
80 
81  typedef BlocksBaseSparseMatrix<T> super_type;
82 
83  BlocksSparseMatrix()
84  :
85  super_type(NBLOCKROWS,NBLOCKCOLS)
86  {}
87 
88 };
89 
90 
91 
92 
107 template< typename T>
108 class MatrixBlockBase : public MatrixSparse<T>
109 {
110  typedef MatrixSparse<T> super;
111 public:
112 
116 
117  typedef MatrixBlockBase<T> self_type;
118 
119  typedef typename super::value_type value_type;
120  typedef typename super::real_type real_type;
121 
122  typedef Backend<value_type> backend_type;
123  typedef boost::shared_ptr<backend_type> backend_ptrtype;
124 
125  typedef super matrix_type;
126  typedef boost::shared_ptr<matrix_type> matrix_ptrtype;
127 
128  typedef std::vector<matrix_ptrtype> vector_matrix_ptrtype;
129 
130  typedef typename super::graph_type graph_type;
131  typedef typename super::graph_ptrtype graph_ptrtype;
132 
134 
138 
139  MatrixBlockBase( vf::BlocksBase<matrix_ptrtype > const & blockSet,
140  backend_type &backend,
141  bool copy_values=true,
142  bool diag_is_nonzero=true );
143 
144  MatrixBlockBase( vf::BlocksBase<graph_ptrtype> const & graph,
145  backend_type &backend,
146  bool diag_is_nonzero=true );
147 
148  MatrixBlockBase( MatrixBlockBase const & mb )
149  :
150  super( mb ),
151  M_mat( mb.M_mat )
152  {}
153 
154  ~MatrixBlockBase()
155  {}
156 
158 
162 
163  MatrixBlockBase operator=( MatrixBlockBase const& mb )
164  {
165  if ( this != &mb )
166  {
167  M_mat = mb.M_mat;
168  }
169 
170  return *this;
171  }
172 
174 
178  matrix_ptrtype getSparseMatrix()
179  {
180  return M_mat;
181  }
183 
187 
188 
190 
194 
203  void init ( const size_type m,
204  const size_type n,
205  const size_type m_l,
206  const size_type n_l,
207  const size_type nnz=30,
208  const size_type noz=10 );
209 
213  void init ( const size_type m,
214  const size_type n,
215  const size_type m_l,
216  const size_type n_l,
217  graph_ptrtype const& graph );
218 
219 
224  void clear ();
225 
229  void zero ();
230 
234  void zero ( size_type start1, size_type size1,
235  size_type start2, size_type size2 );
236 
237 
242  void close () const;
243 
248  size_type size1 () const;
249 
254  size_type size2 () const;
255 
260  size_type rowStart () const;
261 
266  size_type rowStop () const;
267 
274  void set ( const size_type i,
275  const size_type j,
276  const value_type& value );
277 
286  void add ( const size_type i,
287  const size_type j,
288  const value_type& value );
289 
296  void addMatrix ( const ublas::matrix<value_type> &dm,
297  const std::vector<size_type> &rows,
298  const std::vector<size_type> &cols );
299 
306  void addMatrix ( int* rows, int nrows,
307  int* cols, int ncols,
308  value_type* data );
309 
314  void addMatrix ( const ublas::matrix<value_type> &dm,
315  const std::vector<size_type> &dof_indices )
316  {
317  this->addMatrix ( dm, dof_indices, dof_indices );
318  }
319 
320 
326  void addMatrix ( const value_type, MatrixSparse<value_type> & );
327 
328  void scale ( const value_type );
329 
341  value_type operator () ( const size_type i,
342  const size_type j ) const;
343 
347  self_type & operator = ( MatrixSparse<value_type> const& M );
348 
354  void diagonal( Vector<value_type>& out ) const;
355 
361  void transpose( MatrixSparse<value_type>& Mt ) const;
362 
366  value_type
367  energy( Vector<value_type> const& __v,
368  Vector<value_type> const& __u,
369  bool transpose = false ) const;
370 
378  real_type l1Norm () const;
379 
390  real_type linftyNorm () const;
391 
396  bool closed() const;
397 
403  void print( std::ostream& os=std::cout ) const;
404 
409  template <typename U>
410  friend std::ostream& operator << ( std::ostream& os, const MatrixSparse<U>& m );
411 
416  void printPersonal( std::ostream& /*os*/=std::cout ) const
417  {
418  std::cerr << "ERROR: Not Implemented in base class yet!" << std::endl;
419  FEELPP_ASSERT( 0 ).error( "invalid call" );
420  }
421 
428  void printMatlab( const std::string name="NULL" ) const;
429 
436  const std::vector<size_type>& rows,
437  const std::vector<size_type>& cols ) const
438  {
439  this->_get_submatrix( submatrix,
440  rows,
441  cols,
442  false ); // false means DO NOT REUSE submatrix
443  }
444 
452  const std::vector<size_type>& rows,
453  const std::vector<size_type>& cols ) const
454  {
455  this->_get_submatrix( submatrix,
456  rows,
457  cols,
458  true ); // true means REUSE submatrix
459  }
460 
468  void zeroRows( std::vector<int> const& rows, Vector<value_type> const& values, Vector<value_type>& rhs, Context const& on_context );
469 
470  void updateBlockMat( boost::shared_ptr<MatrixSparse<value_type> > m, std::vector<size_type> start_i, std::vector<size_type> start_j );
471 
473 
474 
475 
476 protected:
477 
478 private:
479 
480  //vector_matrix_ptrtype M_v;
481 
482  boost::shared_ptr<MatrixSparse<value_type> > M_mat;
483 };
484 
485 
486 template<int NR, int NC, typename T>
487 class MatrixBlock : public MatrixBlockBase<T>
488 {
489  typedef MatrixBlockBase<T> super_type;
490 
491 public:
492 
493  static const uint16_type NBLOCKROWS = NR;
494  static const uint16_type NBLOCKCOLS = NC;
495  static const uint16_type NBLOCKSIZE = NR * NC;
496 
497  typedef typename super_type::value_type value_type;
498  typedef typename super_type::matrix_ptrtype matrix_ptrtype;
499  typedef typename super_type::backend_type backend_type;
500  typedef vf::Blocks<NBLOCKROWS,NBLOCKCOLS,matrix_ptrtype > blocks_type;
501  typedef vf::BlocksBase<matrix_ptrtype> blocksbase_type;
502  MatrixBlock( blocksbase_type const & blockSet,
503  backend_type &backend,
504  bool copy_values=true,
505  bool diag_is_nonzero=true )
506  :
507  super_type( blockSet,backend,copy_values,diag_is_nonzero )
508  {}
509 
510  MatrixBlock( MatrixBlock const & mb )
511  :
512  super_type( mb )
513  {}
514 
515  MatrixBlock operator=( MatrixBlock const& mb )
516  {
517  super_type::operator=( mb );
518  return *this;
519  }
520 
521  MatrixBlock & operator = ( matrix_ptrtype const& M )
522  {
523  super_type::operator=( M );
524  return *this;
525  }
526 
527 }; // MatrixBlock
528 
529 
530 
531 } // Feel
532 
533 
534 //#include <feel/feelalg/matrixblock.cpp>
535 
536 #endif /* __MatrixBlock_H */

Generated on Sun Oct 20 2013 08:25:01 for Feel++ by doxygen 1.8.4