1#ifndef DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
2#define DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
7#include <dune/common/dynmatrixev.hh>
39 : rows_( 0 ), cols_( 0 )
44 unsigned int rows ()
const {
return rows_; }
45 unsigned int cols ()
const {
return cols_; }
49 assert( (row <
rows()) && (col <
cols()) );
50 return fields_[ row*
cols() + 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;
67 assert( row <
rows() );
73 assert( row <
rows() );
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 ];
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 ); } );
128 if( (
rows != rows_) || (
cols != 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_;
167 template<
class >
friend class Row;
176 : cols_( row.cols_ ),
177 fields_( row.fields_ )
182 assert( col < size() );
183 return fields_[ col ];
188 assert( col < size() );
189 return fields_[ col ];
194 Field *
const end = fields_ + size();
195 for(
Field *it = fields_; it != end; ++it )
214 template<
class DomainSpace,
class RangeSpace >
223 typedef typename RangeSpaceType::RangeFieldType
Field;
232 typedef typename DomainSpace::GridType::template Codim< 0 >::Entity
ColEntityType;
233 typedef typename RangeSpace::GridType::template Codim< 0 >::Entity
RowEntityType;
253 domainSequence_( -1 ),
254 rangeSequence_( -1 ),
255 localMatrixFactory_( *this ),
256 localMatrixStack_( localMatrixFactory_ )
281 template<
class Stencil >
292 template<
class DomainFunction,
class RangeFunction >
293 void apply (
const DomainFunction &u, RangeFunction &w )
const
295 matrix_.
multiply( u.leakPointer(), w.leakPointer() );
302 RangeFunction vFunction(
"v (ddotOEM)",
rangeSpace(), v );
303 RangeFunction wFunction(
"w (ddotOEM)",
rangeSpace(), w );
304 return vFunction.scalarProductDofs( wFunction );
312 RangeFunction wFunction(
"w (multOEM)",
rangeSpace(), w );
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 );
330 const RangeSpace &
rangeSpace ()
const {
return rangeSpace_; }
343 LocalMatrixFactory localMatrixFactory_;
344 mutable LocalMatrixStack localMatrixStack_;
352 template<
class DomainSpace,
class RangeSpace >
372 template<
class DomainSpace,
class RangeSpace >
393 domainMapper_( domainMapper ),
394 rangeMapper_( rangeMapper )
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;
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;
436 assert( (row >= 0) && (row < rows()) );
437 const unsigned int rowIndex = rowIndices_[ row ];
438 matrix_[ rowIndex ].
clear();
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 );
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_;
479 std::vector< unsigned int > rowIndices_;
480 std::vector< unsigned int > colIndices_;
488 template<
class DomainSpace,
class RangeSpace >
497 : matrixObject_( &matrixObject )
502 return new ObjectType( matrixObject_->matrix_, matrixObject_->domainSpace_, matrixObject_->rangeSpace_, matrixObject_->domainMapper_, matrixObject_->rangeMapper_ );
506 MatrixObject *matrixObject_;
Definition: bindguard.hh:11
Definition: adaptivefunction/adaptivefunction.hh:48
Definition: grcommon.hh:31
Definition: misc/functor.hh:31
Default implementation for local matrix classes.
Definition: localmatrix.hh:285
BaseType::RangeEntityType RangeEntityType
Definition: localmatrix.hh:299
BaseType::DomainEntityType DomainEntityType
Definition: localmatrix.hh:298
BaseType::DomainSpaceType DomainSpaceType
Definition: localmatrix.hh:292
BaseType::RangeSpaceType RangeSpaceType
Definition: localmatrix.hh:293
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
void print(std::ostream &s=std::cout) const
Definition: densematrix.hh:138
std::vector< std::complex< double > > eigenValues()
calculate eigenvalues
Definition: densematrix.hh:99
const Field & operator()(unsigned int row, unsigned int col) const
Definition: densematrix.hh:47
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
DomainSpace::GridType::template Codim< 0 >::Entity ColEntityType
Definition: densematrix.hh:232
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
const DomainSpace & domainSpace() const
Definition: densematrix.hh:329
LocalMatrixType localMatrix(const RowEntityType &rowEntity, const ColEntityType &colEntity)
Definition: densematrix.hh:264
MatrixType & matrix()
Definition: densematrix.hh:259
RangeSpace::GridType::template Codim< 0 >::Entity RowEntityType
Definition: densematrix.hh:233
DenseRowMatrixObject(const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace)
Definition: densematrix.hh:247
NonBlockMapper< DomainBlockMapperType, DomainSpaceType::localBlockSize > DomainMapperType
Definition: densematrix.hh:226
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
void reserve(const Stencil &stencil, bool verbose=false)
Definition: densematrix.hh:282
Definition: densematrix.hh:354
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
LocalMatrix(const ThisType &)=delete
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
void unitRow(const int row)
Definition: densematrix.hh:441
const DofType & get(const int row, const int col) const
Definition: densematrix.hh:420
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
Definition: densematrix.hh:490
ObjectType * newObject() const
Definition: densematrix.hh:500
LocalMatrix ObjectType
Definition: densematrix.hh:494
LocalMatrixFactory(MatrixObject &matrixObject)
Definition: densematrix.hh:496