1#ifndef DUNE_FEM_SPMATRIX_HH
2#define DUNE_FEM_SPMATRIX_HH
36 template <
class T,
class IndexT = std::size_t,
37 class ValuesVector = std::vector< T >,
38 class IndicesVector = std::vector< IndexT > >
79 template <
class Stencil>
85 for(
const auto& entry : sparsityPattern )
87 const auto& blockRow = entry.first;
88 const auto& blockColumnSet = entry.second;
91 const size_type nextRow = (blockRow + 1) * rowBlockSize;
92 for(
size_type row = blockRow * rowBlockSize; row < nextRow; ++row )
95 for(
const auto& blockColEntry : blockColumnSet )
97 size_type col = blockColEntry * colBlockSize;
98 for(
size_type c = 0; c<colBlockSize; ++c, ++col, ++column )
100 assert( column <
endRow( row ) );
123 assert((col>=0) && (col <
dim_[1]));
124 assert((row>=0) && (row <
dim_[0]));
136 assert((col>=0) && (col <
dim_[1]));
137 assert((row>=0) && (row <
dim_[0]));
147 template<
class ArgDFType,
class DestDFType>
148 void apply(
const ArgDFType& f, DestDFType& ret )
const
150 constexpr auto blockSize = ArgDFType::DiscreteFunctionSpaceType::localBlockSize;
151 auto ret_it = ret.dbegin();
159 const auto realCol =
columns_[ col ];
164 const auto blockNr = realCol / blockSize ;
165 const auto dofNr = realCol % blockSize ;
166 (*ret_it) +=
values_[ col ] * f.dofVector()[ blockNr ][ dofNr ];
176 assert((col>=0) && (col <
dim_[1]));
177 assert((row>=0) && (row <
dim_[0]));
198 assert((row>=0) && (row <
dim_[0]));
212 assert((row>=0) && (row <
rows()) );
213 assert((col>=0) && (col <
cols()) );
240 assert( row >= 0 && row <
dim_[0] );
248 return std::pair<const field_type, size_type>(
values_[index],
columns_[index]);
252 void print(std::ostream& s=std::cout,
unsigned int offset=0)
const
260 const auto column(rv.second);
261 const auto value(rv.first);
263 s << row+offset <<
" " << column+offset <<
" " << value << std::endl;
268 template <
class SizeT,
class NumericT >
271 matrix.resize(
dim_[0] );
276 auto& matRow = matrix[ row ];
285 matRow[ realCol ] =
values_[ thisCol ];
307 rows_[ row ] = newpos;
308 for(; col<endrow; ++col, ++newpos )
313 assert( newpos <= col );
333 return rows_[ row+1 ];
336 std::tuple< ValuesVector&, IndicesVector&, IndicesVector& >
exportCRS()
366 assert((col>=0) && (col <
dim_[1]));
367 assert((row>=0) && (row <
dim_[0]));
372 for( ; i < endR; ++i )
381 DUNE_THROW( InvalidStateException,
"Could not store entry in sparse matrix - no space available" );
432 template<
class DomainSpace,
class RangeSpace,
437 template<
class MatrixObject >
440 template<
class MatrixObject >
464 typedef Dune::FieldMatrix< field_type, rangeLocalBlockSize, domainLocalBlockSize >
MatrixBlockType;
531 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
541 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
562 template <
class LocalBlock>
565 std::array< size_type, rangeLocalBlockSize > rows;
566 std::array< size_type, domainLocalBlockSize > cols;
577 matrix_.add( rows[ i ], cols[ j ], block[ i ][ j ]);
582 template <
class LocalBlock>
585 std::array< size_type, rangeLocalBlockSize > rows;
586 std::array< size_type, domainLocalBlockSize > cols;
597 matrix_.set( rows[ i ], cols[ j ], block[ i ][ j ]);
602 template<
class LocalMatrix >
605 auto functor = [ &localMat, this ] ( std::pair< int, int > local,
const std::pair< size_type, size_type >& global )
607 matrix_.add( global.first, global.second, localMat.
get( local.first, local.second ) );
613 template<
class LocalMatrix,
class Scalar >
616 auto functor = [ &localMat, &s, this ] ( std::pair< int, int > local,
const std::pair< size_type, size_type >& global )
618 matrix_.add( global.first, global.second, s * localMat.
get( local.first, local.second ) );
624 template<
class LocalMatrix >
627 auto functor = [ &localMat, this ] ( std::pair< int, int > local,
const std::pair< size_type, size_type >& global )
629 matrix_.set( global.first, global.second, localMat.
get( local.first, local.second ) );
635 template<
class LocalMatrix >
638 auto functor = [ &localMat, this ] ( std::pair< int, int > local,
const std::pair< size_type, size_type >& global )
640 localMat.
set( local.first, local.second,
matrix_.get( global.first, global.second ) );
656 void reserve (
const std::vector< Set >& sparsityPattern )
662 template <
class Stencil>
674 std::cout <<
"Max number of base functions = (" <<
rangeMapper_.maxNumDofs() <<
","
682 matrix_.fillPattern( stencil, RangeSpaceType::localBlockSize, DomainSpaceType::localBlockSize );
689 template<
class DomainFunction,
class RangeFunction >
690 void apply(
const DomainFunction &arg, RangeFunction &dest )
const
700 template <
class DiscreteFunctionType >
704 const auto dofEnd = diag.dend();
706 for(
auto dofIt = diag.dbegin(); dofIt != dofEnd; ++dofIt, ++row )
708 assert( row <
matrix_.rows() );
709 (*dofIt) =
matrix_.get( row, row );
713 template <
class Vector>
716 const auto &auxiliaryDofs =
domainSpace().auxiliaryDofs();
720 matrix_.set(r,r,auxiliaryDofs.contains( r )? 0.0 : 1.0);
727 DUNE_THROW(NotImplemented,
"SpMatrixObject::resort is not implemented");
746 template<
class DomainSpace,
class RangeSpace,
class Matrix >
747 template<
class MatrixObject >
767 template<
class DomainSpace,
class RangeSpace,
class Matrix >
768 template<
class MatrixObject >
820 BaseType::init( domainEntity, rangeEntity );
822 rowIndices_.resize(
rangeMapper_.numDofs( rangeEntity ) );
825 columnIndices_.resize(
domainMapper_.numDofs( domainEntity ) );
832 return rowIndices_.size();
838 return columnIndices_.size();
844 assert( value == value );
845 assert( (localRow >= 0) && (localRow < rows()) );
846 assert( (localCol >= 0) && (localCol < columns()) );
847 matrix_.add( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
853 assert( (localRow >= 0) && (localRow < rows()) );
854 assert( (localCol >= 0) && (localCol < columns()) );
855 return matrix_.get( rowIndices_[ localRow ], columnIndices_[ localCol ] );
861 assert( (localRow >= 0) && (localRow < rows()) );
862 assert( (localCol >= 0) && (localCol < columns()) );
863 matrix_.set( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
869 assert( (localRow >= 0) && (localRow < rows()) );
870 matrix_.unitRow( rowIndices_[ localRow ] );
876 assert( (localRow >= 0) && (localRow < rows()) );
877 matrix_.clearRow( rowIndices_[localRow]);
883 assert( (localCol >= 0) && (localCol < columns()) );
884 matrix_.clearCol( columnIndices_[localCol] );
892 matrix_.clearRow( rowIndices_[ i ] );
898 DUNE_THROW(NotImplemented,
"SpMatrixObject::LocalMatrix::resort is not implemented");
922 assert( (localRow >= 0) && (localRow < rows()) );
923 assert( (localCol >= 0) && (localCol < columns()) );
924 matrix_.scale( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
Definition: bindguard.hh:11
PairFunctor< Mapper, Entity, Functor > makePairFunctor(const Mapper &mapper, const Entity &entity, Functor functor)
Definition: operator/matrix/functor.hh:65
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
Definition: localmatrixwrapper.hh:48
default implementation for a general operator stencil
Definition: stencil.hh:35
const GlobalStencilType & globalStencil() const
Return the full stencil.
Definition: stencil.hh:115
int maxNonZerosEstimate() const
Return an upper bound for the maximum number of non-zero entries in all rows.
Definition: stencil.hh:122
a simple wrapper class for sparsity patterns provide as vector< set< size_t > >
Definition: stencil.hh:233
Definition: columnobject.hh:12
SparseRowMatrix.
Definition: spmatrix.hh:40
void print(std::ostream &s=std::cout, unsigned int offset=0) const
print matrix
Definition: spmatrix.hh:252
size_type colIndex(size_type row, size_type col)
returns local col index for given global (row,col)
Definition: spmatrix.hh:364
std::array< size_type, 2 > dim_
Definition: spmatrix.hh:424
std::pair< const field_type, size_type > realValue(size_type index) const
Definition: spmatrix.hh:246
void fillCSRStorage(std::vector< std::map< SizeT, NumericT > > &matrix) const
Definition: spmatrix.hh:269
size_type endRow(const size_type row) const
Definition: spmatrix.hh:331
void clearRow(const size_type row)
set all entries in row to zero
Definition: spmatrix.hh:196
IndexT size_type
matrix index type
Definition: spmatrix.hh:45
size_type maxNzPerRow_
Definition: spmatrix.hh:425
void add(const size_type row, const size_type col, const field_type val)
add value to row,col entry
Definition: spmatrix.hh:134
std::tuple< ValuesVector &, IndicesVector &, IndicesVector & > exportCRS()
Definition: spmatrix.hh:336
void apply(const ArgDFType &f, DestDFType &ret) const
ret = A*f
Definition: spmatrix.hh:148
size_type numNonZeros(size_type row) const
Definition: spmatrix.hh:238
IndicesVector columns_
Definition: spmatrix.hh:421
void reserve(const size_type rows, const size_type cols, const size_type nz)
reserve memory for given rows, columns and number of non zeros
Definition: spmatrix.hh:71
void clear()
set all matrix entries to zero
Definition: spmatrix.hh:189
static const size_type zeroCol
Definition: spmatrix.hh:52
SparseRowMatrix< field_type, size_type, ValuesVector, IndicesVector > ThisType
Definition: spmatrix.hh:46
static const size_type defaultCol
Definition: spmatrix.hh:51
size_type numNonZeros() const
Definition: spmatrix.hh:231
field_type get(const size_type row, const size_type col) const
return value of entry (row,col)
Definition: spmatrix.hh:174
static const int firstCol
Definition: spmatrix.hh:53
void set(const size_type row, const size_type col, const field_type val)
set entry to value (also setting 0 will result in an entry)
Definition: spmatrix.hh:121
size_type startRow(const size_type row) const
Definition: spmatrix.hh:326
void resize(size_type rows, size_type cols, size_type nz)
resize matrix
Definition: spmatrix.hh:345
size_type rows() const
return number of rows
Definition: spmatrix.hh:109
IndicesVector rows_
Definition: spmatrix.hh:422
size_type cols() const
return number of columns
Definition: spmatrix.hh:115
SparseRowMatrix(const size_type rows, const size_type cols, const size_type nz)
Definition: spmatrix.hh:64
ValuesVector values_
Definition: spmatrix.hh:420
bool compressed_
Definition: spmatrix.hh:426
void scale(const size_type row, const size_type col, const field_type val)
scale all entries in row with a given value
Definition: spmatrix.hh:210
void compress()
Definition: spmatrix.hh:290
SparseRowMatrix(const ThisType &)=delete
size_type maxNzPerRow() const
Definition: spmatrix.hh:224
SparseRowMatrix()
construct matrix of zero size
Definition: spmatrix.hh:58
ThisType MatrixBaseType
Definition: spmatrix.hh:49
void fillPattern(const Stencil &stencil, const size_type rowBlockSize, const size_type colBlockSize)
reserve memory for given rows, columns and number of non zeros
Definition: spmatrix.hh:80
T field_type
matrix field type
Definition: spmatrix.hh:43
SparseRowMatrixObject.
Definition: spmatrix.hh:435
NonBlockMapper< DomainBlockMapperType, DomainSpaceType::localBlockSize > DomainMapperType
Definition: spmatrix.hh:452
RangeSpaceType::EntityType RangeEntityType
Definition: spmatrix.hh:447
const DomainSpaceType & domainSpace() const
get domain space (i.e. space that builds the rows)
Definition: spmatrix.hh:492
const DomainSpaceType & domainSpace_
Definition: spmatrix.hh:733
DomainMapperType domainMapper_
Definition: spmatrix.hh:735
LocalMatrixType localMatrix() const
Definition: spmatrix.hh:542
DomainSpaceType::EntityType DomainEntityType
Definition: spmatrix.hh:446
RangeSpaceType::BlockMapperType RangeBlockMapperType
Definition: spmatrix.hh:453
Matrix MatrixType
Definition: spmatrix.hh:455
MatrixType::size_type size_type
Definition: spmatrix.hh:456
LocalMatrixStackType localMatrixStack_
Definition: spmatrix.hh:740
MatrixType & exportMatrix() const
get reference to storage object
Definition: spmatrix.hh:514
ThisType LocalMatrixFactoryType
Definition: spmatrix.hh:470
LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType
Definition: spmatrix.hh:472
void setUnitRows(const Vector &rows)
Definition: spmatrix.hh:714
void addScaledLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s)
Definition: spmatrix.hh:614
SparseRowMatrixObject< DomainSpaceType, RangeSpaceType, MatrixType > ThisType
Definition: spmatrix.hh:459
RangeSpace RangeSpaceType
Definition: spmatrix.hh:445
void setBlock(const size_type row, const size_type col, const LocalBlock &block)
Definition: spmatrix.hh:583
RangeSpaceType::EntityType RowEntityType
Definition: spmatrix.hh:449
static const size_type rangeLocalBlockSize
Definition: spmatrix.hh:462
void reserve(const std::vector< Set > &sparsityPattern)
Definition: spmatrix.hh:656
LocalMatrix< ThisType > ObjectType
Definition: spmatrix.hh:469
ColumnObject< ThisType > LocalColumnObjectType
Definition: spmatrix.hh:473
void apply(const DomainFunction &arg, RangeFunction &dest) const
apply matrix to discrete function
Definition: spmatrix.hh:690
Fem::ObjectStack< LocalMatrixFactoryType > LocalMatrixStackType
Definition: spmatrix.hh:471
int sequence_
Definition: spmatrix.hh:737
SparseRowMatrixObject(const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const SolverParameter ¶m=SolverParameter())
construct matrix object
Definition: spmatrix.hh:476
static const size_type domainLocalBlockSize
Definition: spmatrix.hh:461
LocalColumnObjectType localColumn(const DomainEntityType &domainEntity) const
get local column
Definition: spmatrix.hh:548
void extractDiagonal(DiscreteFunctionType &diag) const
Definition: spmatrix.hh:701
void setLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat)
Definition: spmatrix.hh:625
MatrixType::field_type field_type
Definition: spmatrix.hh:457
LocalMatrixType localMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity) const
Definition: spmatrix.hh:532
MatrixType PreconditionMatrixType
Definition: spmatrix.hh:467
DomainSpace DomainSpaceType
Definition: spmatrix.hh:444
DomainSpaceType::BlockMapperType DomainBlockMapperType
Definition: spmatrix.hh:451
DomainSpaceType::EntityType ColumnEntityType
Definition: spmatrix.hh:448
ObjectType * newObject() const
interface method from LocalMatrixFactory
Definition: spmatrix.hh:522
NonBlockMapper< RangeBlockMapperType, RangeSpaceType::localBlockSize > RangeMapperType
Definition: spmatrix.hh:454
Dune::FieldMatrix< field_type, rangeLocalBlockSize, domainLocalBlockSize > MatrixBlockType
Definition: spmatrix.hh:464
void reserve(const Stencil &stencil, bool verbose=false)
reserve memory
Definition: spmatrix.hh:663
void addLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat)
Definition: spmatrix.hh:603
MatrixType matrix_
Definition: spmatrix.hh:738
void compress()
compress matrix to a real CRS format
Definition: spmatrix.hh:653
void resort()
resort row numbering in matrix to have ascending numbering
Definition: spmatrix.hh:725
void getLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat) const
Definition: spmatrix.hh:636
void unitRow(const size_type row)
Definition: spmatrix.hh:553
void addBlock(const size_type row, const size_type col, const LocalBlock &block)
Definition: spmatrix.hh:563
RangeMapperType rangeMapper_
Definition: spmatrix.hh:736
MatrixBlockType block_type
Definition: spmatrix.hh:465
void clear()
clear matrix
Definition: spmatrix.hh:647
bool preconditioning_
Definition: spmatrix.hh:739
const RangeSpaceType & rangeSpace_
Definition: spmatrix.hh:734
const RangeSpaceType & rangeSpace() const
get range space (i.e. space that builds the columns)
Definition: spmatrix.hh:498
MatrixType & matrix() const
get reference to storage object, for internal use
Definition: spmatrix.hh:505
void finalizeAssembly() const
Definition: spmatrix.hh:510
LocalMatrixTraits.
Definition: spmatrix.hh:749
SparseRowMatrixObject< DomainSpaceType, RangeSpaceType, Matrix > SparseRowMatrixObjectType
Definition: spmatrix.hh:753
RangeFieldType LittleBlockType
Definition: spmatrix.hh:758
RangeSpaceType::RangeFieldType RangeFieldType
Definition: spmatrix.hh:757
DomainSpace DomainSpaceType
Definition: spmatrix.hh:750
RangeSpace RangeSpaceType
Definition: spmatrix.hh:751
SparseRowMatrixObjectType::template LocalMatrix< MatrixObject > LocalMatrixType
Definition: spmatrix.hh:755
SparseRowMatrixObjectType::RangeMapperType RangeMapperType
Definition: spmatrix.hh:761
SparseRowMatrixObjectType::DomainMapperType DomainMapperType
Definition: spmatrix.hh:760
LocalMatrix.
Definition: spmatrix.hh:771
MatrixType & matrix_
Definition: spmatrix.hh:928
LocalMatrix(const LocalMatrix &)=delete
void clearRow(size_type localRow)
set matrix row to zero
Definition: spmatrix.hh:874
const DomainMapperType & domainMapper_
Definition: spmatrix.hh:929
void add(size_type localRow, size_type localCol, DofType value)
add value to matrix entry
Definition: spmatrix.hh:842
void init(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity)
Definition: spmatrix.hh:817
MatrixObjectType::MatrixType MatrixType
type of matrix
Definition: spmatrix.hh:784
void clearCol(size_type localCol)
set matrix column to zero
Definition: spmatrix.hh:881
void resort()
resort all global rows of matrix to have ascending numbering
Definition: spmatrix.hh:896
size_type columns() const
return number of columns
Definition: spmatrix.hh:836
ColumnIndicesType columnIndices_
Definition: spmatrix.hh:932
size_type rows() const
return number of rows
Definition: spmatrix.hh:830
RowIndicesType rowIndices_
Definition: spmatrix.hh:931
std::vector< typename DomainMapperType::SizeType > ColumnIndicesType
Definition: spmatrix.hh:801
Traits::RangeFieldType RangeFieldType
type of entries of little blocks
Definition: spmatrix.hh:787
DofType get(size_type localRow, size_type localCol) const
get matrix entry
Definition: spmatrix.hh:851
MatrixObject MatrixObjectType
type of matrix object
Definition: spmatrix.hh:774
std::vector< typename RangeMapperType::SizeType > RowIndicesType
Definition: spmatrix.hh:800
Traits::DomainMapperType DomainMapperType
type of nonblocked domain mapper
Definition: spmatrix.hh:796
Traits::RangeMapperType RangeMapperType
type of nonblocked domain mapper
Definition: spmatrix.hh:798
LocalMatrixTraits< MatrixObjectType > Traits
type of the traits
Definition: spmatrix.hh:777
LocalMatrix(const MatrixObjectType &matrixObject, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const DomainMapperType &domainMapper, const RangeMapperType &rangeMapper)
constructor
Definition: spmatrix.hh:804
const RangeMapperType & rangeMapper_
Definition: spmatrix.hh:930
void scale(const DofType &value)
scale local matrix with a certain value
Definition: spmatrix.hh:905
RangeFieldType DofType
type of the DoFs
Definition: spmatrix.hh:790
void scale(size_type localRow, size_type localCol, DofType value)
scale matrix entry with value
Definition: spmatrix.hh:920
void clear()
clear all entries belonging to local matrix
Definition: spmatrix.hh:888
Traits::LittleBlockType LittleBlockType
type of little blocks
Definition: spmatrix.hh:793
void set(size_type localRow, size_type localCol, DofType value)
set matrix entry to value
Definition: spmatrix.hh:859
void unitRow(size_type localRow)
set matrix row to zero except diagonla entry
Definition: spmatrix.hh:867
Definition: solver/parameter.hh:15
static const int none
Definition: solver/parameter.hh:40
static const int jacobi
Definition: solver/parameter.hh:45