dune-fem  2.6-git
densematrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
2 #define DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
3 
4 #include <memory>
5 #include <iostream>
6 
7 #include <dune/common/dynmatrixev.hh>
8 #include <dune/common/std/memory.hh>
9 
11 #include <dune/fem/misc/functor.hh>
16 
17 namespace Dune
18 {
19 
20  namespace Fem
21  {
22 
23  // DenseRowMatrix
24  // --------------
25 
26  template< class F >
28  {
30 
31  public:
32  typedef F Field;
33 
34  template< class RF >
35  class Row;
36 
37  DenseRowMatrix () : rows_( 0 ), cols_( 0 ) {}
38 
39  DenseRowMatrix ( unsigned int rows, unsigned int cols )
40  : rows_( 0 ), cols_( 0 )
41  {
42  reserve( rows, cols );
43  }
44 
45  unsigned int rows () const { return rows_; }
46  unsigned int cols () const { return cols_; }
47 
48  const Field &operator() ( unsigned int row, unsigned int col ) const
49  {
50  assert( (row < rows()) && (col < cols()) );
51  return fields_[ row*cols() + col ];
52  }
53 
54  Field &operator() ( unsigned int row, unsigned int col )
55  {
56  assert( (row < rows()) && (col < cols()) );
57  return fields_[ row*cols() + col ];
58  }
59 
60  void add ( unsigned int row, unsigned int col, const Field &value )
61  {
62  assert( (row < rows()) && (col < cols()) );
63  fields_[ row*cols() + col ] += value;
64  }
65 
66  Row< const Field > operator[] ( unsigned int row ) const
67  {
68  assert( row < rows() );
69  return Row< const Field >( cols(), fields_.get() + row*cols() );
70  }
71 
72  Row< Field > operator[] ( unsigned int row )
73  {
74  assert( row < rows() );
75  return Row< Field >( cols(), fields_.get() + row*cols() );
76  }
77 
78  void clear () { std::fill( fields_.get(), fields_.get() + (rows()*cols()), Field( 0 ) ); }
79 
80  void multiply ( const Field *x, Field *y ) const
81  {
82  for( unsigned int row = 0; row < rows(); ++row )
83  {
84  const Field *fields = fields_.get() + row*cols();
85  y[ row ] = Field( 0 );
86  for( unsigned int col = 0; col < cols(); ++col )
87  y[ row ] += fields[ col ] * x[ col ];
88  }
89  }
90 
100  std::vector< std::complex< double > > eigenValues ()
101  {
102  if( rows() != cols() )
103  DUNE_THROW( InvalidStateException, "Requiring square matrix for eigenvalue computation" );
104 
105  const long int N = rows();
106  const char jobvl = 'n';
107  const char jobvr = 'n';
108 
109  // working memory
110  std::unique_ptr< double[] > work = Std::make_unique< double[] >( 5*N );
111 
112  // return value information
113  long int info = 0;
114  long int lwork = 3*N;
115 
116  // call LAPACK routine (see fmatrixev_ext.cc)
117  DynamicMatrixHelp::eigenValuesNonsymLapackCall( &jobvl, &jobvr, &N, fields_, &N, work.get(), work.get()+N, 0, &N, 0, &N, work.get()+2*N, &lwork, &info );
118 
119  if( info != 0 )
120  DUNE_THROW( MathError, "DenseRowMatrix: Eigenvalue computation failed" );
121 
122  std::vector< std::complex< double > > eigenValues( N );
123  std::transform( work.get(), work.get()+N, work.get()+N, eigenValues.begin(), [] ( double r, double i ) { return std::complex< double >( r, i ); } );
124  return eigenValues;
125  }
126 
127  void reserve ( unsigned int rows, unsigned int cols )
128  {
129  if( (rows != rows_) || (cols != cols_) )
130  {
131  fields_.reset( new Field[ rows*cols ] );
132  rows_ = rows;
133  cols_ = cols;
134  }
135  // Martin: Is clearing required, here?
136  clear();
137  }
138 
139  void print( std::ostream& s=std::cout ) const
140  {
141  s.precision( 6 );
142  for( unsigned int row = 0; row < rows(); ++row )
143  {
144  const Field *fields = fields_ + row*cols();
145  for( unsigned int col = 0; col < cols(); ++col )
146  s << fields[ col ] << " ";
147 
148  s << std::endl;
149  }
150  }
151 
152  private:
153  unsigned int rows_, cols_;
154  std::unique_ptr< Field[] > fields_;
155  };
156 
157 
158 
159  // DenseRowMatrix::Row
160  // -------------------
161 
162  template< class F >
163  template< class RF >
164  class DenseRowMatrix< F >::Row
165  {
166  typedef Row< RF > ThisType;
167 
168  template< class > friend class Row;
169 
170  public:
171  Row ( unsigned int cols, RF *fields )
172  : cols_( cols ),
173  fields_( fields )
174  {}
175 
176  Row ( const Row< F > &row )
177  : cols_( row.cols_ ),
178  fields_( row.fields_ )
179  {}
180 
181  const RF &operator[] ( const unsigned int col ) const
182  {
183  assert( col < size() );
184  return fields_[ col ];
185  }
186 
187  RF &operator[] ( const unsigned int col )
188  {
189  assert( col < size() );
190  return fields_[ col ];
191  }
192 
193  void clear ()
194  {
195  Field *const end = fields_ + size();
196  for( Field *it = fields_; it != end; ++it )
197  *it = Field( 0 );
198  }
199 
200  unsigned int size () const
201  {
202  return cols_;
203  }
204 
205  private:
206  unsigned int cols_;
207  RF *fields_;
208  };
209 
210 
211 
212  // DenseRowMatrixObject
213  // --------------------
214 
215  template< class DomainSpace, class RangeSpace >
217  {
219 
220  public:
221  typedef DomainSpace DomainSpaceType;
222  typedef RangeSpace RangeSpaceType;
223 
224  typedef typename RangeSpaceType::RangeFieldType Field;
225 
226  typedef typename DomainSpaceType::BlockMapperType DomainBlockMapperType;
228  typedef typename RangeSpaceType::BlockMapperType RangeBlockMapperType;
230 
231  typedef typename DomainSpaceType::EntityType DomainEntityType;
232  typedef typename RangeSpaceType::EntityType RangeEntityType;
233  typedef typename DomainSpace::GridType::template Codim< 0 >::Entity ColEntityType;
234  typedef typename RangeSpace::GridType::template Codim< 0 >::Entity RowEntityType;
235 
237 
238  private:
239  class LocalMatrixTraits;
240  class LocalMatrix;
241  class LocalMatrixFactory;
242 
244 
245  public:
247 
249  const RangeSpaceType &rangeSpace )
250  : domainSpace_( domainSpace ),
251  rangeSpace_( rangeSpace ),
252  domainMapper_( domainSpace.blockMapper() ),
253  rangeMapper_( rangeSpace.blockMapper() ),
254  domainSequence_( -1 ),
255  rangeSequence_( -1 ),
256  localMatrixFactory_( *this ),
257  localMatrixStack_( localMatrixFactory_ )
258  {}
259 
261  {
262  return matrix_;
263  }
264 
266  const ColEntityType &colEntity )
267  {
268  return LocalMatrixType( localMatrixStack_, rowEntity, colEntity );
269  }
270 
271 
273  {
274  return LocalMatrixType( localMatrixStack_ );
275  }
276 
277  void clear ()
278  {
279  matrix_.clear();
280  }
281 
282  template< class Stencil >
283  void reserve ( const Stencil &stencil, bool verbose = false )
284  {
285  if( (domainSequence_ != domainSpace().sequence()) || (rangeSequence_ != rangeSpace().sequence()) )
286  {
287  matrix_.reserve( rangeSpace().size(), domainSpace().size() );
288  domainSequence_ = domainSpace().sequence();
289  rangeSequence_ = rangeSpace().sequence();
290  }
291  }
292 
293  template< class DomainFunction, class RangeFunction >
294  void apply ( const DomainFunction &u, RangeFunction &w ) const
295  {
296  matrix_.multiply( u.leakPointer(), w.leakPointer() );
297  rangeSpace().communicate( w );
298  }
299 
300  Field ddotOEM ( const Field *v, const Field *w ) const
301  {
302  typedef AdaptiveDiscreteFunction< RangeSpaceType > RangeFunction;
303  RangeFunction vFunction( "v (ddotOEM)", rangeSpace(), v );
304  RangeFunction wFunction( "w (ddotOEM)", rangeSpace(), w );
305  return vFunction.scalarProductDofs( wFunction );
306  }
307 
308  void multOEM ( const Field *u, Field *w ) const
309  {
310  matrix_.multiply( u, w );
311 
312  typedef AdaptiveDiscreteFunction< RangeSpaceType > RangeFunction;
313  RangeFunction wFunction( "w (multOEM)", rangeSpace(), w );
314  rangeSpace().communicate( wFunction );
315  }
316 
317  template< class DiscreteFunctionType >
318  void extractDiagonal( DiscreteFunctionType &diag ) const
319  {
320  assert( matrix_.rows() == matrix_.cols() );
321  const auto dofEnd = diag.dend();
322  unsigned int row = 0;
323  for( auto dofIt = diag.dbegin(); dofIt != dofEnd; ++dofIt, ++row )
324  {
325  assert( row < matrix_.rows() );
326  (*dofIt) = matrix_( row, row );
327  }
328  }
329 
330  const DomainSpace &domainSpace () const { return domainSpace_; }
331  const RangeSpace &rangeSpace () const { return rangeSpace_; }
332 
333  private:
334  const DomainSpaceType &domainSpace_;
335  const RangeSpaceType &rangeSpace_;
336  DomainMapperType domainMapper_;
337  RangeMapperType rangeMapper_;
338 
339  int domainSequence_;
340  int rangeSequence_;
341 
342  MatrixType matrix_;
343 
344  LocalMatrixFactory localMatrixFactory_;
345  mutable LocalMatrixStack localMatrixStack_;
346  };
347 
348 
349 
350  // DenseRowMatrixObject::LocalMatrixTraits
351  // ---------------------------------------
352 
353  template< class DomainSpace, class RangeSpace >
354  class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrixTraits
355  {
357 
358  public:
361 
364 
365  typedef typename MatrixObject::LocalMatrix LocalMatrixType;
366  };
367 
368 
369 
370  // DenseRowMatrixObject::LocalMatrix
371  // ---------------------------------
372 
373  template< class DomainSpace, class RangeSpace >
374  class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrix
375  : public LocalMatrixDefault< LocalMatrixTraits >
376  {
378 
379  typedef LocalMatrix ThisType;
381 
382  public:
383  typedef LocalMatrixTraits Traits;
384 
386 
387  typedef typename Traits::RangeFieldType RangeFieldType;
388  typedef typename Traits::LittleBlockType LittleBlockType;
390 
393  matrix_( matrix ),
394  domainMapper_( domainMapper ),
395  rangeMapper_( rangeMapper )
396  {}
397 
398  LocalMatrix ( const ThisType & ) = delete;
399  ThisType &operator= ( const ThisType & ) = delete;
400 
401  void init ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
402  {
403  BaseType::init( domainEntity, rangeEntity );
404 
405  rowIndices_.resize( rangeMapper_.numDofs( rangeEntity ) );
406  rangeMapper_.mapEach( rangeEntity, Fem::AssignFunctor< std::vector< unsigned int > >( rowIndices_ ) );
407  colIndices_.resize( domainMapper_.numDofs( domainEntity ) );
408  domainMapper_.mapEach( domainEntity, Fem::AssignFunctor< std::vector< unsigned int > >( colIndices_ ) );
409  }
410 
411  int rows () const { return rowIndices_.size(); }
412  int cols () const { return colIndices_.size(); }
413 
414  void add ( const int row, const int col, const DofType &value )
415  {
416  assert( (row >= 0) && (row < rows()) );
417  assert( (col >= 0) && (col < cols()) );
418  matrix_( rowIndices_[ row ], colIndices_[ col ] ) += value;
419  }
420 
421  const DofType &get ( const int row, const int col ) const
422  {
423  assert( (row >= 0) && (row < rows()) );
424  assert( (col >= 0) && (col < cols()) );
425  return matrix_( rowIndices_[ row ], colIndices_[ col ] );
426  }
427 
428  void set ( const int row, const int col, const DofType &value )
429  {
430  assert( (row >= 0) && (row < rows()) );
431  assert( (col >= 0) && (col < cols()) );
432  matrix_( rowIndices_[ row ], colIndices_[ col ] ) = value;
433  }
434 
435  void clearRow ( const int row )
436  {
437  assert( (row >= 0) && (row < rows()) );
438  const unsigned int rowIndex = rowIndices_[ row ];
439  matrix_[ rowIndex ].clear();
440  }
441 
442  void unitRow ( const int row )
443  {
444  clearRow( row );
445  set( row, row, DofType( 1 ) );
446  }
447 
448  void clear ()
449  {
450  typedef std::vector< unsigned int >::const_iterator Iterator;
451  const Iterator rowEnd = rowIndices_.end();
452  for( Iterator rowIt = rowIndices_.begin(); rowIt != rowEnd; ++rowIt )
453  {
454  const Iterator colEnd = colIndices_.end();
455  for( Iterator colIt = colIndices_.begin(); colIt != colEnd; ++colIt )
456  matrix_( *rowIt, *colIt ) = DofType( 0 );
457  }
458  }
459 
460  void scale ( const DofType &value )
461  {
462  typedef std::vector< unsigned int >::const_iterator Iterator;
463  const Iterator rowEnd = rowIndices_.end();
464  for( Iterator rowIt = rowIndices_.begin(); rowIt != rowEnd; ++rowIt )
465  {
466  const Iterator colEnd = colIndices_.end();
467  for( Iterator colIt = colIndices_.begin(); colIt != colEnd; ++colIt )
468  matrix_( *rowIt, *colIt ) *= value;
469  }
470  }
471 
472  protected:
473  using BaseType::domainSpace_;
474  using BaseType::rangeSpace_;
475 
476  private:
477  MatrixType &matrix_;
478  const DomainMapperType &domainMapper_;
479  const RangeMapperType &rangeMapper_;
480  std::vector< unsigned int > rowIndices_;
481  std::vector< unsigned int > colIndices_;
482  };
483 
484 
485 
486  // DenseRowMatrixObject::LocalMatrixFactory
487  // ----------------------------------------
488 
489  template< class DomainSpace, class RangeSpace >
490  class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrixFactory
491  {
493 
494  public:
495  typedef LocalMatrix ObjectType;
496 
498  : matrixObject_( &matrixObject )
499  {}
500 
502  {
503  return new ObjectType( matrixObject_->matrix_, matrixObject_->domainSpace_, matrixObject_->rangeSpace_, matrixObject_->domainMapper_, matrixObject_->rangeMapper_ );
504  }
505 
506  private:
507  MatrixObject *matrixObject_;
508  };
509 
510  } // namespace Fem
511 
512 } // namespace Dune
513 
514 #endif // #ifndef DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
Definition: bindguard.hh:11
Definition: function/adaptivefunction/adaptivefunction.hh:45
Definition: grcommon.hh:31
Definition: bartonnackmaninterface.hh:17
Definition: misc/functor.hh:31
Interface for local matrix classes.
Definition: localmatrix.hh:32
Default implementation for local matrix classes.
Definition: localmatrix.hh:269
Definition: localmatrixwrapper.hh:48
default implementation for a general operator stencil
Definition: stencil.hh:32
Definition: densematrix.hh:28
DenseRowMatrix()
Definition: densematrix.hh:37
void reserve(unsigned int rows, unsigned int cols)
Definition: densematrix.hh:127
unsigned int cols() const
Definition: densematrix.hh:46
unsigned int rows() const
Definition: densematrix.hh:45
F Field
Definition: densematrix.hh:32
void add(unsigned int row, unsigned int col, const Field &value)
Definition: densematrix.hh:60
Row< const Field > operator[](unsigned int row) const
Definition: densematrix.hh:66
DenseRowMatrix(unsigned int rows, unsigned int cols)
Definition: densematrix.hh:39
void clear()
Definition: densematrix.hh:78
void multiply(const Field *x, Field *y) const
Definition: densematrix.hh:80
const Field & operator()(unsigned int row, unsigned int col) const
Definition: densematrix.hh:48
void print(std::ostream &s=std::cout) const
Definition: densematrix.hh:139
std::vector< std::complex< double > > eigenValues()
calculate eigenvalues
Definition: densematrix.hh:100
Definition: densematrix.hh:165
Row(const Row< F > &row)
Definition: densematrix.hh:176
unsigned int size() const
Definition: densematrix.hh:200
void clear()
Definition: densematrix.hh:193
Row(unsigned int cols, RF *fields)
Definition: densematrix.hh:171
Definition: densematrix.hh:217
RangeSpaceType::EntityType RangeEntityType
Definition: densematrix.hh:232
DenseRowMatrix< Field > MatrixType
Definition: densematrix.hh:236
NonBlockMapper< RangeBlockMapperType, RangeSpaceType::localBlockSize > RangeMapperType
Definition: densematrix.hh:229
const RangeSpace & rangeSpace() const
Definition: densematrix.hh:331
DomainSpace DomainSpaceType
Definition: densematrix.hh:221
void apply(const DomainFunction &u, RangeFunction &w) const
Definition: densematrix.hh:294
void extractDiagonal(DiscreteFunctionType &diag) const
Definition: densematrix.hh:318
LocalMatrixType localMatrix() const
Definition: densematrix.hh:272
LocalMatrixWrapper< LocalMatrixStack > LocalMatrixType
Definition: densematrix.hh:246
RangeSpace RangeSpaceType
Definition: densematrix.hh:222
RangeSpaceType::RangeFieldType Field
Definition: densematrix.hh:224
DomainSpaceType::EntityType DomainEntityType
Definition: densematrix.hh:231
void clear()
Definition: densematrix.hh:277
LocalMatrixType localMatrix(const RowEntityType &rowEntity, const ColEntityType &colEntity)
Definition: densematrix.hh:265
DenseRowMatrixObject(const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace)
Definition: densematrix.hh:248
NonBlockMapper< DomainBlockMapperType, DomainSpaceType::localBlockSize > DomainMapperType
Definition: densematrix.hh:227
RangeSpace::GridType::template Codim< 0 >::Entity RowEntityType
Definition: densematrix.hh:234
const DomainSpace & domainSpace() const
Definition: densematrix.hh:330
RangeSpaceType::BlockMapperType RangeBlockMapperType
Definition: densematrix.hh:228
DomainSpaceType::BlockMapperType DomainBlockMapperType
Definition: densematrix.hh:226
Field ddotOEM(const Field *v, const Field *w) const
Definition: densematrix.hh:300
void multOEM(const Field *u, Field *w) const
Definition: densematrix.hh:308
DomainSpace::GridType::template Codim< 0 >::Entity ColEntityType
Definition: densematrix.hh:233
MatrixType & matrix()
Definition: densematrix.hh:260
void reserve(const Stencil &stencil, bool verbose=false)
Definition: densematrix.hh:283
MatrixObject::LocalMatrix LocalMatrixType
Definition: densematrix.hh:365
RangeFieldType LittleBlockType
Definition: densematrix.hh:363
MatrixObject::DomainSpaceType DomainSpaceType
Definition: densematrix.hh:359
MatrixObject::RangeSpaceType RangeSpaceType
Definition: densematrix.hh:360
MatrixObject::Field RangeFieldType
Definition: densematrix.hh:362
Definition: densematrix.hh:376
Traits::LittleBlockType LittleBlockType
Definition: densematrix.hh:388
LocalMatrix(MatrixType &matrix, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const DomainMapperType &domainMapper, const RangeMapperType &rangeMapper)
Definition: densematrix.hh:391
MatrixObject::MatrixType MatrixType
Definition: densematrix.hh:385
int rows() const
Definition: densematrix.hh:411
void scale(const DofType &value)
Definition: densematrix.hh:460
Traits::RangeFieldType RangeFieldType
Definition: densematrix.hh:387
void clear()
Definition: densematrix.hh:448
void init(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity)
Definition: densematrix.hh:401
int cols() const
Definition: densematrix.hh:412
RangeFieldType DofType
Definition: densematrix.hh:389
const DofType & get(const int row, const int col) const
Definition: densematrix.hh:421
void unitRow(const int row)
Definition: densematrix.hh:442
LocalMatrixTraits Traits
Definition: densematrix.hh:383
void add(const int row, const int col, const DofType &value)
Definition: densematrix.hh:414
void set(const int row, const int col, const DofType &value)
Definition: densematrix.hh:428
void clearRow(const int row)
Definition: densematrix.hh:435
LocalMatrix ObjectType
Definition: densematrix.hh:495
ObjectType * newObject() const
Definition: densematrix.hh:501
LocalMatrixFactory(MatrixObject &matrixObject)
Definition: densematrix.hh:497