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