1 #ifndef DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
2 #define DUNE_FEM_OPERATOR_MATRIX_DENSEMATRIX_HH
7 #include <dune/common/dynmatrixev.hh>
8 #include <dune/common/std/memory.hh>
40 : rows_( 0 ), cols_( 0 )
45 unsigned int rows ()
const {
return rows_; }
46 unsigned int cols ()
const {
return cols_; }
50 assert( (row <
rows()) && (col <
cols()) );
51 return fields_[ row*
cols() + col ];
56 assert( (row <
rows()) && (col <
cols()) );
57 return fields_[ row*
cols() + col ];
60 void add (
unsigned int row,
unsigned int col,
const Field &value )
62 assert( (row <
rows()) && (col <
cols()) );
63 fields_[ row*
cols() + col ] += value;
68 assert( row <
rows() );
69 return Row< const Field >(
cols(), fields_.get() + row*
cols() );
74 assert( row <
rows() );
75 return Row< Field >(
cols(), fields_.get() + row*
cols() );
82 for(
unsigned int row = 0; row <
rows(); ++row )
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 ];
103 DUNE_THROW( InvalidStateException,
"Requiring square matrix for eigenvalue computation" );
105 const long int N =
rows();
106 const char jobvl =
'n';
107 const char jobvr =
'n';
110 std::unique_ptr< double[] > work = Std::make_unique< double[] >( 5*N );
114 long int lwork = 3*N;
117 DynamicMatrixHelp::eigenValuesNonsymLapackCall( &jobvl, &jobvr, &N, fields_, &N, work.get(), work.get()+N, 0, &N, 0, &N, work.get()+2*N, &lwork, &
info );
120 DUNE_THROW( MathError,
"DenseRowMatrix: Eigenvalue computation failed" );
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 ); } );
129 if( (
rows != rows_) || (
cols != cols_) )
139 void print( std::ostream& s=std::cout )
const
142 for(
unsigned int row = 0; row <
rows(); ++row )
144 const Field *fields = fields_ + row*
cols();
145 for(
unsigned int col = 0; col <
cols(); ++col )
146 s << fields[ col ] <<
" ";
153 unsigned int rows_, cols_;
154 std::unique_ptr< Field[] > fields_;
168 template<
class >
friend class Row;
176 Row (
const Row< F > &row )
177 : cols_( row.cols_ ),
178 fields_( row.fields_ )
183 assert( col < size() );
184 return fields_[ col ];
189 assert( col < size() );
190 return fields_[ col ];
195 Field *
const end = fields_ + size();
196 for(
Field *it = fields_; it != end; ++it )
215 template<
class DomainSpace,
class RangeSpace >
224 typedef typename RangeSpaceType::RangeFieldType
Field;
233 typedef typename DomainSpace::GridType::template Codim< 0 >::Entity
ColEntityType;
234 typedef typename RangeSpace::GridType::template Codim< 0 >::Entity
RowEntityType;
239 class LocalMatrixTraits;
241 class LocalMatrixFactory;
254 domainSequence_( -1 ),
255 rangeSequence_( -1 ),
256 localMatrixFactory_( *this ),
257 localMatrixStack_( localMatrixFactory_ )
282 template<
class Stencil >
293 template<
class DomainFunction,
class RangeFunction >
294 void apply (
const DomainFunction &u, RangeFunction &w )
const
296 matrix_.
multiply( u.leakPointer(), w.leakPointer() );
303 RangeFunction vFunction(
"v (ddotOEM)",
rangeSpace(), v );
304 RangeFunction wFunction(
"w (ddotOEM)",
rangeSpace(), w );
305 return vFunction.scalarProductDofs( wFunction );
313 RangeFunction wFunction(
"w (multOEM)",
rangeSpace(), w );
317 template<
class DiscreteFunctionType >
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 )
325 assert( row < matrix_.
rows() );
326 (*dofIt) = matrix_( row, row );
331 const RangeSpace &
rangeSpace ()
const {
return rangeSpace_; }
344 LocalMatrixFactory localMatrixFactory_;
345 mutable LocalMatrixStack localMatrixStack_;
353 template<
class DomainSpace,
class RangeSpace >
373 template<
class DomainSpace,
class RangeSpace >
394 domainMapper_( domainMapper ),
395 rangeMapper_( rangeMapper )
399 ThisType &operator= (
const ThisType & ) =
delete;
403 BaseType::init( domainEntity, rangeEntity );
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_ ) );
411 int rows ()
const {
return rowIndices_.size(); }
412 int cols ()
const {
return colIndices_.size(); }
414 void add (
const int row,
const int col,
const DofType &value )
416 assert( (row >= 0) && (row < rows()) );
417 assert( (col >= 0) && (col < cols()) );
418 matrix_( rowIndices_[ row ], colIndices_[ col ] ) += value;
423 assert( (row >= 0) && (row < rows()) );
424 assert( (col >= 0) && (col < cols()) );
425 return matrix_( rowIndices_[ row ], colIndices_[ col ] );
428 void set (
const int row,
const int col,
const DofType &value )
430 assert( (row >= 0) && (row < rows()) );
431 assert( (col >= 0) && (col < cols()) );
432 matrix_( rowIndices_[ row ], colIndices_[ col ] ) = value;
437 assert( (row >= 0) && (row < rows()) );
438 const unsigned int rowIndex = rowIndices_[ row ];
439 matrix_[ rowIndex ].
clear();
450 typedef std::vector< unsigned int >::const_iterator Iterator;
451 const Iterator rowEnd = rowIndices_.end();
452 for( Iterator rowIt = rowIndices_.begin(); rowIt != rowEnd; ++rowIt )
454 const Iterator colEnd = colIndices_.end();
455 for( Iterator colIt = colIndices_.begin(); colIt != colEnd; ++colIt )
456 matrix_( *rowIt, *colIt ) =
DofType( 0 );
462 typedef std::vector< unsigned int >::const_iterator Iterator;
463 const Iterator rowEnd = rowIndices_.end();
464 for( Iterator rowIt = rowIndices_.begin(); rowIt != rowEnd; ++rowIt )
466 const Iterator colEnd = colIndices_.end();
467 for( Iterator colIt = colIndices_.begin(); colIt != colEnd; ++colIt )
468 matrix_( *rowIt, *colIt ) *= value;
473 using BaseType::domainSpace_;
474 using BaseType::rangeSpace_;
480 std::vector< unsigned int > rowIndices_;
481 std::vector< unsigned int > colIndices_;
489 template<
class DomainSpace,
class RangeSpace >
498 : matrixObject_( &matrixObject )
503 return new ObjectType( matrixObject_->matrix_, matrixObject_->domainSpace_, matrixObject_->rangeSpace_, matrixObject_->domainMapper_, matrixObject_->rangeMapper_ );
507 MatrixObject *matrixObject_;
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
Definition: densematrix.hh:355
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
LocalMatrix(const ThisType &)=delete
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
Definition: densematrix.hh:491
LocalMatrix ObjectType
Definition: densematrix.hh:495
ObjectType * newObject() const
Definition: densematrix.hh:501
LocalMatrixFactory(MatrixObject &matrixObject)
Definition: densematrix.hh:497