1#ifndef DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
2#define DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
9#include <dune/fem/function/adaptivefunction.hh>
10#include <dune/fem/misc/functor.hh>
11#include <dune/fem/operator/common/localmatrix.hh>
12#include <dune/fem/operator/common/localmatrixwrapper.hh>
13#include <dune/fem/solver/krylovinverseoperators.hh>
14#include <dune/fem/storage/objectstack.hh>
28 typedef DenseRowMatrix< F > ThisType;
36 DenseRowMatrix () : rows_( 0 ), cols_( 0 ) {}
38 DenseRowMatrix (
unsigned int rows,
unsigned int cols )
39 : rows_( 0 ), cols_( 0 )
41 reserve( rows, cols );
44 unsigned int rows ()
const {
return rows_; }
45 unsigned int cols ()
const {
return cols_; }
47 const Field &operator() (
unsigned int row,
unsigned int col )
const
49 assert( (row < rows()) && (col < cols()) );
50 return fields_[ row*cols() + col ];
53 Field &operator() (
unsigned int row,
unsigned int col )
55 assert( (row < rows()) && (col < cols()) );
56 return fields_[ row*cols() + col ];
59 void add (
unsigned int row,
unsigned int col,
const Field &value )
61 assert( (row < rows()) && (col < cols()) );
62 fields_[ row*cols() + col ] += value;
65 Row< const Field > operator[] (
unsigned int row )
const
67 assert( row < rows() );
68 return Row< const Field >( cols(), fields_.get() + row*cols() );
71 Row< Field > operator[] (
unsigned int row )
73 assert( row < rows() );
74 return Row< Field >( cols(), fields_.get() + row*cols() );
77 void clear () { std::fill( fields_.get(), fields_.get() + (rows()*cols()), Field( 0 ) ); }
79 void multiply (
const Field *x, Field *y )
const
81 for(
unsigned int row = 0; row < rows(); ++row )
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 ];
99 std::vector< std::complex< double > > eigenValues ()
101 if( rows() != cols() )
102 DUNE_THROW( InvalidStateException,
"Requiring square matrix for eigenvalue computation" );
104 const long int N = rows();
105 const char jobvl =
'n';
106 const char jobvr =
'n';
109 std::unique_ptr< double[] > work = std::make_unique< double[] >( 5*N );
113 long int lwork = 3*N;
116 DynamicMatrixHelp::eigenValuesNonsymLapackCall( &jobvl, &jobvr, &N, fields_, &N, work.get(), work.get()+N, 0, &N, 0, &N, work.get()+2*N, &lwork, &info );
119 DUNE_THROW( MathError,
"DenseRowMatrix: Eigenvalue computation failed" );
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 ); } );
126 void reserve (
unsigned int rows,
unsigned int cols )
128 if( (rows != rows_) || (cols != cols_) )
130 fields_.reset(
new Field[ rows*cols ] );
138 void print( std::ostream& s=std::cout )
const
141 for(
unsigned int row = 0; row < rows(); ++row )
143 const Field *fields = fields_ + row*cols();
144 for(
unsigned int col = 0; col < cols(); ++col )
145 s << fields[ col ] <<
" ";
152 unsigned int rows_, cols_;
153 std::unique_ptr< Field[] > fields_;
163 class DenseRowMatrix< F >::Row
165 typedef Row< RF > ThisType;
167 template<
class >
friend class Row;
170 Row (
unsigned int cols, RF *fields )
175 Row (
const Row< F > &row )
176 : cols_( row.cols_ ),
177 fields_( row.fields_ )
180 const RF &operator[] (
const unsigned int col )
const
182 assert( col <
size() );
183 return fields_[ col ];
186 RF &operator[] (
const unsigned int col )
188 assert( col <
size() );
189 return fields_[ col ];
194 Field *
const end = fields_ +
size();
195 for( Field *it = fields_; it != end; ++it )
199 unsigned int size ()
const
214 template<
class DomainSpace,
class RangeSpace >
215 class DenseRowMatrixObject
217 typedef DenseRowMatrixObject< DomainSpace, RangeSpace > ThisType;
220 typedef DomainSpace DomainSpaceType;
221 typedef RangeSpace RangeSpaceType;
223 typedef typename RangeSpaceType::RangeFieldType Field;
225 typedef typename DomainSpaceType::BlockMapperType DomainBlockMapperType;
227 typedef typename RangeSpaceType::BlockMapperType RangeBlockMapperType;
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;
235 typedef DenseRowMatrix< Field > MatrixType;
238 class LocalMatrixTraits;
240 class LocalMatrixFactory;
242 typedef Fem::ObjectStack< LocalMatrixFactory > LocalMatrixStack;
245 typedef LocalMatrixWrapper< LocalMatrixStack > LocalMatrixType;
247 DenseRowMatrixObject (
const DomainSpaceType &domainSpace,
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_ )
259 MatrixType &matrix ()
264 LocalMatrixType localMatrix (
const RowEntityType &rowEntity,
265 const ColEntityType &colEntity )
267 return LocalMatrixType( localMatrixStack_, rowEntity, colEntity );
271 LocalMatrixType localMatrix ()
const
273 return LocalMatrixType( localMatrixStack_ );
281 template<
class Stencil >
282 void reserve (
const Stencil &stencil,
bool verbose =
false )
284 if( (domainSequence_ != domainSpace().sequence()) || (rangeSequence_ != rangeSpace().sequence()) )
286 matrix_.reserve( rangeSpace().
size(), domainSpace().
size() );
287 domainSequence_ = domainSpace().sequence();
288 rangeSequence_ = rangeSpace().sequence();
292 template<
class DomainFunction,
class RangeFunction >
293 void apply (
const DomainFunction &u, RangeFunction &w )
const
295 matrix_.multiply( u.leakPointer(), w.leakPointer() );
296 rangeSpace().communicate( w );
299 Field ddotOEM (
const Field *v,
const Field *w )
const
301 typedef AdaptiveDiscreteFunction< RangeSpaceType > RangeFunction;
302 RangeFunction vFunction(
"v (ddotOEM)", rangeSpace(), v );
303 RangeFunction wFunction(
"w (ddotOEM)", rangeSpace(), w );
304 return vFunction.scalarProductDofs( wFunction );
307 void multOEM (
const Field *u, Field *w )
const
309 matrix_.multiply( u, w );
311 typedef AdaptiveDiscreteFunction< RangeSpaceType > RangeFunction;
312 RangeFunction wFunction(
"w (multOEM)", rangeSpace(), w );
313 rangeSpace().communicate( wFunction );
316 template<
class DiscreteFunctionType >
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 )
324 assert( row < matrix_.rows() );
325 (*dofIt) = matrix_( row, row );
329 const DomainSpace &domainSpace ()
const {
return domainSpace_; }
330 const RangeSpace &rangeSpace ()
const {
return rangeSpace_; }
333 const DomainSpaceType &domainSpace_;
334 const RangeSpaceType &rangeSpace_;
335 DomainMapperType domainMapper_;
336 RangeMapperType rangeMapper_;
343 LocalMatrixFactory localMatrixFactory_;
344 mutable LocalMatrixStack localMatrixStack_;
352 template<
class DomainSpace,
class RangeSpace >
353 class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrixTraits
355 typedef DenseRowMatrixObject< DomainSpace, RangeSpace > MatrixObject;
358 typedef typename MatrixObject::DomainSpaceType DomainSpaceType;
359 typedef typename MatrixObject::RangeSpaceType RangeSpaceType;
361 typedef typename MatrixObject::Field RangeFieldType;
362 typedef RangeFieldType LittleBlockType;
364 typedef typename MatrixObject::LocalMatrix LocalMatrixType;
372 template<
class DomainSpace,
class RangeSpace >
373 class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrix
374 :
public LocalMatrixDefault< LocalMatrixTraits >
376 typedef DenseRowMatrixObject< DomainSpace, RangeSpace > MatrixObject;
378 typedef LocalMatrix ThisType;
379 typedef LocalMatrixDefault< LocalMatrixTraits > BaseType;
382 typedef LocalMatrixTraits Traits;
384 typedef typename MatrixObject::MatrixType MatrixType;
386 typedef typename Traits::RangeFieldType RangeFieldType;
387 typedef typename Traits::LittleBlockType LittleBlockType;
388 typedef RangeFieldType DofType;
390 LocalMatrix ( MatrixType &matrix,
const DomainSpaceType &domainSpace,
const RangeSpaceType &rangeSpace,
const DomainMapperType &domainMapper,
const RangeMapperType &rangeMapper )
391 : BaseType( domainSpace, rangeSpace ),
393 domainMapper_( domainMapper ),
394 rangeMapper_( rangeMapper )
397 LocalMatrix (
const ThisType & ) =
delete;
398 ThisType &operator= (
const ThisType & ) =
delete;
400 void init (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
402 BaseType::init( domainEntity, rangeEntity );
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_ ) );
410 int rows ()
const {
return rowIndices_.size(); }
411 int cols ()
const {
return colIndices_.size(); }
413 void add (
const int row,
const int col,
const DofType &value )
415 assert( (row >= 0) && (row < rows()) );
416 assert( (col >= 0) && (col < cols()) );
417 matrix_( rowIndices_[ row ], colIndices_[ col ] ) += value;
420 const DofType &
get (
const int row,
const int col )
const
422 assert( (row >= 0) && (row < rows()) );
423 assert( (col >= 0) && (col < cols()) );
424 return matrix_( rowIndices_[ row ], colIndices_[ col ] );
427 void set (
const int row,
const int col,
const DofType &value )
429 assert( (row >= 0) && (row < rows()) );
430 assert( (col >= 0) && (col < cols()) );
431 matrix_( rowIndices_[ row ], colIndices_[ col ] ) = value;
434 void clearRow (
const int row )
436 assert( (row >= 0) && (row < rows()) );
437 const unsigned int rowIndex = rowIndices_[ row ];
438 matrix_[ rowIndex ].clear();
441 void unitRow (
const int row )
444 set( row, row, DofType( 1 ) );
449 typedef std::vector< unsigned int >::const_iterator Iterator;
450 const Iterator rowEnd = rowIndices_.end();
451 for( Iterator rowIt = rowIndices_.begin(); rowIt != rowEnd; ++rowIt )
453 const Iterator colEnd = colIndices_.end();
454 for( Iterator colIt = colIndices_.begin(); colIt != colEnd; ++colIt )
455 matrix_( *rowIt, *colIt ) = DofType( 0 );
459 void scale (
const DofType &value )
461 typedef std::vector< unsigned int >::const_iterator Iterator;
462 const Iterator rowEnd = rowIndices_.end();
463 for( Iterator rowIt = rowIndices_.begin(); rowIt != rowEnd; ++rowIt )
465 const Iterator colEnd = colIndices_.end();
466 for( Iterator colIt = colIndices_.begin(); colIt != colEnd; ++colIt )
467 matrix_( *rowIt, *colIt ) *= value;
472 using BaseType::domainSpace_;
473 using BaseType::rangeSpace_;
477 const DomainMapperType &domainMapper_;
478 const RangeMapperType &rangeMapper_;
479 std::vector< unsigned int > rowIndices_;
480 std::vector< unsigned int > colIndices_;
488 template<
class DomainSpace,
class RangeSpace >
489 class DenseRowMatrixObject< DomainSpace, RangeSpace >::LocalMatrixFactory
491 typedef DenseRowMatrixObject< DomainSpace, RangeSpace > MatrixObject;
494 typedef LocalMatrix ObjectType;
496 LocalMatrixFactory ( MatrixObject &matrixObject )
497 : matrixObject_( &matrixObject )
500 ObjectType *newObject ()
const
502 return new ObjectType( matrixObject_->matrix_, matrixObject_->domainSpace_, matrixObject_->rangeSpace_, matrixObject_->domainMapper_, matrixObject_->rangeMapper_ );
506 MatrixObject *matrixObject_;
ConstDofIteratorType dbegin() const
Obtain the constant iterator pointing to the first dof.
Definition: discretefunction.hh:761
ConstDofIteratorType dend() const
Obtain the constant iterator pointing to the last dof.
Definition: discretefunction.hh:773
forward declaration
Definition: discretefunction.hh:51
utility functions to compute eigenvalues for dense matrices.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22