1#ifndef DUNE_FEM_ISTLMATRIXWRAPPER_HH
2#define DUNE_FEM_ISTLMATRIXWRAPPER_HH
25#include <dune/fem/function/blockvectorfunction.hh>
26#include <dune/fem/operator/common/localmatrix.hh>
27#include <dune/fem/operator/common/localmatrixwrapper.hh>
28#include <dune/fem/operator/common/operator.hh>
29#include <dune/fem/operator/common/stencil.hh>
30#include <dune/fem/function/common/scalarproducts.hh>
31#include <dune/fem/io/parameter.hh>
32#include <dune/fem/operator/matrix/columnobject.hh>
33#include <dune/fem/storage/objectstack.hh>
35#include <dune/fem/operator/matrix/istlmatrixadapter.hh>
36#include <dune/fem/operator/matrix/istlpreconditioner.hh>
37#include <dune/fem/operator/matrix/functor.hh>
47 template<
class MatrixObject >
48 class ISTLLocalMatrix;
50 template <
class DomainSpaceImp,
class RangeSpaceImp,
53 class ISTLMatrixObject;
58 template <
class LittleBlockType,
class RowDiscreteFunctionImp,
class ColDiscreteFunctionImp = RowDiscreteFunctionImp>
59 class ImprovedBCRSMatrix :
public BCRSMatrix<LittleBlockType>
61 friend struct MatrixDimension<ImprovedBCRSMatrix>;
63 typedef RowDiscreteFunctionImp RowDiscreteFunctionType;
64 typedef ColDiscreteFunctionImp ColDiscreteFunctionType;
66 typedef BCRSMatrix<LittleBlockType> BaseType;
68 typedef BaseType MatrixBaseType;
69 typedef typename BaseType :: RowIterator RowIteratorType ;
70 typedef typename BaseType :: ColIterator ColIteratorType ;
72 typedef ImprovedBCRSMatrix< LittleBlockType, RowDiscreteFunctionImp, ColDiscreteFunctionImp > ThisType;
74 typedef typename BaseType :: size_type size_type;
79 typedef typename BaseType::field_type field_type;
82 typedef typename BaseType::block_type block_type;
85 typedef typename BaseType:: allocator_type allocator_type;
88 typedef typename BaseType :: row_type row_type;
91 typedef typename BaseType :: ColIterator ColIterator;
94 typedef typename BaseType :: ConstColIterator ConstColIterator;
97 typedef typename BaseType :: RowIterator RowIterator;
100 typedef typename BaseType :: ConstRowIterator ConstRowIterator;
103 typedef typename ColDiscreteFunctionType :: DiscreteFunctionSpaceType RangeSpaceType;
106 typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
109 typedef typename ColDiscreteFunctionType :: DofStorageType ColBlockVectorType;
112 typedef typename RangeSpaceType :: GridType :: Traits :: Communication CommunicationType ;
114 typedef typename BaseType :: BuildMode BuildMode ;
118 ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nnz,
double overflowFraction) :
119 BaseType (rows, cols, nnz, overflowFraction, BaseType::implicit)
123 ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nz = 0 ) :
124 BaseType (rows, cols, BaseType :: row_wise)
128 ImprovedBCRSMatrix( ) :
133 ImprovedBCRSMatrix(
const ImprovedBCRSMatrix& org) :
137 ConstRowIterator slicedBegin(
const size_type row )
const
139 return ConstRowIterator( this->r, row );
142 ConstRowIterator slicedEnd(
const size_type row )
const
144 return ConstRowIterator( this->r, row );
147 std::pair< size_type, size_type > sliceBeginEnd(
const size_type thread,
const size_type numThreads )
const
149 const size_type sliceSize = this->N() / numThreads ;
150 const size_type sliceBegin = thread * sliceSize ;
151 const size_type sliceEnd = ( thread == numThreads-1 ) ? this->N() : (sliceBegin + sliceSize) ;
152 return std::make_pair( sliceBegin, sliceEnd );
156 template <
class X,
class Y,
bool isMV>
157 void mvThreadedImpl(
const field_type& alpha,
158 const X& x, Y& y, std::integral_constant<bool, isMV> )
const
160 auto doMV = [
this, &alpha, &x, &y] ()
162 const auto slice = sliceBeginEnd( MPIManager::thread(), MPIManager::numThreads() );
163 const ConstRowIterator endi = slicedEnd( slice.second );
164 for (ConstRowIterator i = slicedBegin(slice.first); i!=endi; ++i)
166 if constexpr ( isMV )
169 ConstColIterator endj = (*i).end();
170 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
172 auto&& xj = Dune::Impl::asVector(x[j.index()]);
173 auto&& yi = Dune::Impl::asVector(y[i.index()]);
174 if constexpr ( isMV )
175 Dune::Impl::asMatrix(*j).umv(xj, yi);
177 Dune::Impl::asMatrix(*j).usmv(alpha, xj, yi);
184 MPIManager :: run ( doMV );
186 catch (
const SingleThreadModeError& e )
189 if constexpr ( isMV )
193 DUNE_THROW(SingleThreadModeError,
"ImprovedBCRSMatrix::usmvThreaded: Cannot recover from threading error. Disable threading!");
200 template <
class X,
class Y>
201 void usmvThreaded(
const field_type& alpha,
const X& x, Y& y )
const
203 mvThreadedImpl(alpha, x, y, std::false_type() );
206 template <
class X,
class Y>
207 void mvThreaded(
const X& x, Y& y )
const
209 mvThreadedImpl(1.0, x, y, std::true_type() );
212 template <
class SparsityPattern >
213 void createEntries(
const SparsityPattern& sparsityPattern)
215 const auto spEnd = sparsityPattern.end();
218 auto endcreate = this->createend();
219 for(
auto create = this->createbegin(); create != endcreate; ++create)
221 const auto row = sparsityPattern.find( create.index() );
223 if( row == spEnd )
continue;
226 for (
const auto& col : row->second )
227 create.insert( col );
234 for (
auto& row : *
this)
235 for (
auto& entry : row)
240 void unitRow(
const size_t row )
242 block_type idBlock( 0 );
243 for (
int i = 0; i < idBlock.rows; ++i)
246 auto& matRow = (*this)[ row ];
247 auto colIt = matRow.begin();
248 const auto& colEndIt = matRow.end();
249 for (; colIt != colEndIt; ++colIt)
251 if( colIt.index() == row )
259 template <
class HangingNodesType>
260 void setup(ThisType& oldMatrix,
const HangingNodesType& hangingNodes)
263 typedef std::set< std::pair<int, block_type> > LocalEntryType;
264 typedef std::map< int , LocalEntryType > EntriesType;
268 std::map< int , std::set<int> > indices;
270 auto rowend = oldMatrix.end();
271 for(
auto it = oldMatrix.begin(); it != rowend; ++it)
273 const auto row = it.index();
274 auto& localIndices = indices[ row ];
276 if( hangingNodes.isHangingNode( row ) )
279 const auto& cols = hangingNodes.associatedDofs( row );
280 const auto colSize = cols.size();
281 for(
auto i=0; i<colSize; ++i)
283 assert( ! hangingNodes.isHangingNode( cols[i].first ) );
286 auto& localColIndices = indices[ cols[i].first ];
287 auto& localEntry = entries[ cols[i].first ];
290 auto endj = (*it).end();
291 for (
auto j= (*it).begin(); j!=endj; ++j)
293 localColIndices.insert( j.index () );
294 localEntry.insert( std::make_pair( j.index(), (cols[i].second * (*j)) ));
299 localIndices.insert( row );
300 for(
auto i=0; i<colSize; ++i)
301 localIndices.insert( cols[i].first );
306 auto endj = (*it).end();
307 for (
auto j= (*it).begin(); j!=endj; ++j)
308 localIndices.insert( j.index () );
313 createEntries( indices );
316 auto rowit = oldMatrix.begin();
318 auto endcreate = this->end();
319 for(
auto create = this->begin(); create != endcreate; ++create, ++rowit )
321 assert( rowit != oldMatrix.end() );
323 const auto row = create.index();
324 if( hangingNodes.isHangingNode( row ) )
326 const auto& cols = hangingNodes.associatedDofs( row );
328 std::map< const int , block_type > colMap;
330 assert( block_type :: rows == 1 );
332 const auto colSize = cols.size();
333 for(
auto i=0; i<colSize; ++i)
334 colMap[ cols[i].first ] = -cols[i].second;
338 auto endj = (*create).end();
339 for (
auto j= (*create).begin(); j!=endj; ++j)
341 assert( colMap.find( j.index() ) != colMap.end() );
342 (*j) = colMap[ j.index() ];
346 else if ( entries.find( row ) == entries.end() )
348 auto colit = (*rowit).begin();
349 auto endj = (*create).end();
350 for (
auto j= (*create).begin(); j!=endj; ++j, ++colit )
352 assert( colit != (*rowit).end() );
358 std::map< int , block_type > oldCols;
361 auto colend = (*rowit).end();
362 for(
auto colit = (*rowit).begin(); colit != colend; ++colit)
363 oldCols[ colit.index() ] = 0;
366 auto entry = entries.find( row );
367 assert( entry != entries.end ());
370 auto endcol = (*entry).second.end();
371 for(
auto co = (*entry).second.begin(); co != endcol; ++co)
372 oldCols[ (*co).first ] = 0;
376 auto colend = (*rowit).end();
377 for(
auto colit = (*rowit).begin(); colit != colend; ++colit)
378 oldCols[ colit.index() ] += (*colit);
382 auto endcol = (*entry).second.end();
383 for(
auto co = (*entry).second.begin(); co != endcol; ++co)
384 oldCols[ (*co).first ] += (*co).second;
387 auto endj = (*create).end();
388 for (
auto j= (*create).begin(); j!=endj; ++j )
390 auto colEntry = oldCols.find( j.index() );
391 if( colEntry != oldCols.end() )
392 (*j) = (*colEntry).second;
401 void extractDiagonal( ColBlockVectorType& diag )
const
403 const auto endi = this->end();
404 for (
auto i = this->begin(); i!=endi; ++i)
407 const auto row = i.index();
408 auto entry = (*i).find( row );
409 const LittleBlockType& block = (*entry);
410 enum { blockSize = LittleBlockType :: rows };
411 for(
auto l=0; l<blockSize; ++l )
412 diag[ row ][ l ] = block[ l ][ l ];
418 field_type operator()(
const std::size_t row,
const std::size_t col)
const
420 const std::size_t blockRow(row/(LittleBlockType :: rows));
421 const std::size_t localRowIdx(row%(LittleBlockType :: rows));
422 const std::size_t blockCol(col/(LittleBlockType :: cols));
423 const std::size_t localColIdx(col%(LittleBlockType :: cols));
425 const auto& matrixRow(this->
operator[](blockRow));
426 auto entry = matrixRow.find( blockCol );
427 const LittleBlockType& block = (*entry);
428 return block[localRowIdx][localColIdx];
433 void set(
const std::size_t row,
const std::size_t col, field_type value)
435 const std::size_t blockRow(row/(LittleBlockType :: rows));
436 const std::size_t localRowIdx(row%(LittleBlockType :: rows));
437 const std::size_t blockCol(col/(LittleBlockType :: cols));
438 const std::size_t localColIdx(col%(LittleBlockType :: cols));
440 auto& matrixRow(this->
operator[](blockRow));
441 auto entry = matrixRow.find( blockCol );
442 LittleBlockType& block = (*entry);
443 block[localRowIdx][localColIdx] = value;
447 void print(std::ostream& s=std::cout,
unsigned int offset=0)
const
450 const auto endi=this->end();
451 for (
auto i=this->begin(); i!=endi; ++i)
453 const auto endj = (*i).end();
454 for (
auto j=(*i).begin(); j!=endj; ++j)
455 if( (*j).infinity_norm() > 1.e-15)
456 s << i.index()+offset <<
" " << j.index()+offset <<
" " << *j << std::endl;
466 template<
class MatrixObject >
467 struct ISTLLocalMatrixTraits
469 typedef typename MatrixObject::DomainSpaceType DomainSpaceType;
470 typedef typename MatrixObject::RangeSpaceType RangeSpaceType;
471 typedef typename DomainSpaceType::RangeFieldType RangeFieldType;
473 typedef ISTLLocalMatrix< MatrixObject > LocalMatrixType;
474 typedef typename MatrixObject::MatrixType::block_type LittleBlockType;
481 template<
class MatrixObject >
482 class ISTLLocalMatrix
483 :
public LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > >
487 typedef LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > > BaseType;
490 typedef MatrixObject MatrixObjectType;
492 typedef typename MatrixObjectType::MatrixType MatrixType;
494 typedef typename MatrixType::block_type LittleBlockType;
496 typedef typename MatrixType :: size_type size_type;
498 typedef typename MatrixObjectType::DomainSpaceType DomainSpaceType;
499 typedef typename MatrixObjectType::RangeSpaceType RangeSpaceType;
501 typedef typename MatrixObjectType::DomainEntityType DomainEntityType;
502 typedef typename MatrixObjectType::RangeEntityType RangeEntityType;
505 typedef typename DomainSpaceType::RangeFieldType DofType;
506 typedef typename MatrixType::row_type RowType;
509 typedef typename DomainSpaceType::BlockMapperType ColMapperType;
511 typedef typename RangeSpaceType::BlockMapperType RowMapperType;
513 static const int littleCols = MatrixObjectType::littleCols;
514 static const int littleRows = MatrixObjectType::littleRows;
516 typedef typename MatrixType::size_type Index;
521 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
522 ISTLLocalMatrix (
const MatrixObjectType& mObj,
const DomainSpaceType& domainSpace,
const RangeSpaceType& rangeSpace )
523 : BaseType( domainSpace, rangeSpace ),
524 rowMapper_( rangeSpace.blockMapper() ),
525 colMapper_( domainSpace.blockMapper() ),
526 numRows_( rowMapper_.maxNumDofs() ),
527 numCols_( colMapper_.maxNumDofs() ),
534 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
535 ISTLLocalMatrix (
const ISTLLocalMatrix& org )
537 rowMapper_(org.rowMapper_),
538 colMapper_(org.colMapper_),
539 numRows_( org.numRows_ ),
540 numCols_( org.numCols_ ),
541 matrices_(org.matrices_),
542 matrixObj_(org.matrixObj_)
546 void init (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
549 BaseType :: init ( domainEntity, rangeEntity );
551 numRows_ = rowMapper_.numDofs( rangeEntity );
552 numCols_ = colMapper_.numDofs( domainEntity );
555 matrices_.resize( numRows_ );
556 for(
auto& row : matrices_ )
558 row.resize( numCols_,
nullptr );
561 if( matrixObj_.implicitModeActive() )
563 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
565 return matrixObj_.matrix().entry( index.first, index.second );
567 initBlocks( blockAccess, domainEntity, rangeEntity );
571 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
573 return matrixObj_.matrix()[ index.first ][ index.second ];
575 initBlocks( blockAccess, domainEntity, rangeEntity );
579 template <
class BlockAccess>
580 void initBlocks( BlockAccess& blockAccess,
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
582 auto functor = [
this, &blockAccess ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
584 matrices_[ local.first ][ local.second ] = &blockAccess( index );
587 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
592 void check(
int localRow,
int localCol)
const
595 const std::size_t row = (int) localRow / littleRows;
596 const std::size_t col = (int) localCol / littleCols;
597 const int lRow = localRow%littleRows;
598 const int lCol = localCol%littleCols;
599 assert( row < matrices_.size() ) ;
600 assert( col < matrices_[row].
size() );
601 assert( lRow < littleRows );
602 assert( lCol < littleCols );
606 DofType& getValue(
const int localRow,
const int localCol)
608 const int row = (int) localRow / littleRows;
609 const int col = (int) localCol / littleCols;
610 const int lRow = localRow%littleRows;
611 const int lCol = localCol%littleCols;
612 return (*matrices_[row][col])[lRow][lCol];
616 const DofType
get(
const int localRow,
const int localCol)
const
618 const int row = (int) localRow / littleRows;
619 const int col = (int) localCol / littleCols;
620 const int lRow = localRow%littleRows;
621 const int lCol = localCol%littleCols;
622 return (*matrices_[row][col])[lRow][lCol];
625 void scale (
const DofType& scalar)
627 const size_type nrows = matrices_.size();
628 for(size_type i=0; i<nrows; ++i)
630 const size_type ncols = matrices_[i].size();
631 for(size_type j=0; j<ncols; ++j)
633 (*matrices_[i][j]) *= scalar;
638 void add(
const int localRow,
const int localCol ,
const DofType value)
641 check(localRow,localCol);
643 getValue(localRow,localCol) += value;
646 void set(
const int localRow,
const int localCol ,
const DofType value)
649 check(localRow,localCol);
651 getValue(localRow,localCol) = value;
655 void unitRow(
const int localRow)
657 const int row = (int) localRow / littleRows;
658 const int lRow = localRow%littleRows;
661 doClearRow( row, lRow );
664 (*matrices_[row][row])[lRow][lRow] = 1;
670 const size_type nrows = matrices_.size();
671 for(size_type i=0; i<nrows; ++i)
673 const size_type ncols = matrices_[i].size();
674 for(size_type j=0; j<ncols; ++j)
676 (*matrices_[i][j]) = (DofType) 0;
682 void clearRow (
const int localRow )
684 const int row = (int) localRow / littleRows;
685 const int lRow = localRow%littleRows;
688 doClearRow( row, lRow );
697 void doClearRow (
const int row,
const int lRow )
700 const auto col = this->columns();
701 for(
auto localCol=0; localCol<col; ++localCol)
703 const int col = (int) localCol / littleCols;
704 const int lCol = localCol%littleCols;
705 (*matrices_[row][col])[lRow][lCol] = 0;
711 const RowMapperType& rowMapper_;
712 const ColMapperType& colMapper_;
719 std::vector< std::vector< LittleBlockType* > > matrices_;
722 const MatrixObjectType& matrixObj_;
727 template <
class DomainSpaceImp,
class RangeSpaceImp,
class DomainBlock,
class RangeBlock >
728 class ISTLMatrixObject
732 typedef DomainSpaceImp DomainSpaceType;
734 typedef RangeSpaceImp RangeSpaceType;
737 typedef ISTLMatrixObject< DomainSpaceImp, RangeSpaceImp, DomainBlock, RangeBlock > ThisType;
738 typedef ThisType PreconditionMatrixType;
740 typedef typename DomainSpaceType::GridType GridType;
742 typedef typename RangeSpaceType :: EntityType RangeEntityType ;
743 typedef typename DomainSpaceType :: EntityType DomainEntityType ;
745 enum { littleCols = DomainSpaceType :: localBlockSize };
746 enum { littleRows = RangeSpaceType :: localBlockSize };
748 typedef FieldMatrix<typename DomainSpaceType :: RangeFieldType, littleRows, littleCols> LittleBlockType;
749 typedef LittleBlockType block_type;
750 typedef LittleBlockType MatrixBlockType;
752 typedef ISTLBlockVectorDiscreteFunction< RangeSpaceType, RangeBlock > RowDiscreteFunctionType;
753 typedef ISTLBlockVectorDiscreteFunction< DomainSpaceType, DomainBlock > ColumnDiscreteFunctionType;
755 typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
756 typedef typename ColumnDiscreteFunctionType :: DofStorageType ColumnBlockVectorType;
758 typedef typename RangeSpaceType :: BlockMapperType RowMapperType;
759 typedef typename DomainSpaceType :: BlockMapperType ColMapperType;
762 typedef ImprovedBCRSMatrix< LittleBlockType , ColumnDiscreteFunctionType , RowDiscreteFunctionType > MatrixType;
763 typedef ISTLParallelMatrixAdapterInterface< MatrixType > MatrixAdapterType;
766 typedef ISTLLocalMatrix<ThisType> ObjectType;
767 friend class ISTLLocalMatrix<ThisType>;
769 typedef ThisType LocalMatrixFactoryType;
772 typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
773 typedef ColumnObject< ThisType > LocalColumnObjectType;
776 const DomainSpaceType& domainSpace_;
777 const RangeSpaceType& rangeSpace_;
780 RowMapperType& rowMapper_;
782 ColMapperType& colMapper_;
787 mutable std::unique_ptr< MatrixType > matrix_;
789 mutable LocalMatrixStackType localMatrixStack_;
791 mutable std::unique_ptr< MatrixAdapterType > matrixAdap_;
793 const double overflowFraction_;
794 const bool threading_ ;
796 ISTLSolverParameter param_;
798 ISTLMatrixObject(
const ISTLMatrixObject&) =
delete;
803 ISTLMatrixObject (
const DomainSpaceType &domainSpace,
const RangeSpaceType &rangeSpace,
const ISTLSolverParameter& param = ISTLSolverParameter() ) :
804 domainSpace_(domainSpace)
805 , rangeSpace_(rangeSpace)
806 , rowMapper_( rangeSpace.blockMapper() )
807 , colMapper_( domainSpace.blockMapper() )
810 , localMatrixStack_( *this )
811 , overflowFraction_( param.overflowFraction() )
812 , threading_( param.threading() )
819 MatrixType& matrix()
const
826 bool threading()
const {
return threading_; }
828 void printTexInfo(std::ostream& out)
const
830 out <<
"ISTL MatrixObj: ";
835 std::string preconditionName()
const
840 void createMatrixAdapter (
const ISTLSolverParameter& param )
const
844 matrixAdap_ = ISTLMatrixAdapterFactory< ThisType > :: matrixAdapter( *
this, param );
846 assert( matrixAdap_ );
850 MatrixAdapterType& matrixAdapter(
const ISTLSolverParameter& parameter )
const
852 const ISTLSolverParameter* param =
dynamic_cast< const ISTLSolverParameter*
> (¶meter);
854 createMatrixAdapter( *param );
857 ISTLSolverParameter newparam( parameter );
858 createMatrixAdapter( newparam );
866 MatrixAdapterType& matrixAdapter()
const
868 return matrixAdapter( param_ );
872 MatrixType& exportMatrix()
const
878 bool implicitModeActive()
const
884 if( matrix().buildMode() == MatrixType::implicit && matrix().buildStage() == MatrixType::building )
890 void finalizeAssembly()
const
893 const_cast< ThisType&
> (*this).compress();
899 if( implicitModeActive() )
913 void unitRow(
const size_t row )
915 matrix().unitRow( row );
919 template <
class Container>
920 void setUnitRowImpl(
const Container &rows,
const double diagonal )
922 for (
const auto& r : rows)
924 const std::size_t blockRow( r/(LittleBlockType :: rows) );
925 const std::size_t localRowIdx( r%(LittleBlockType :: rows) );
926 auto& row = matrix()[blockRow];
927 const auto endcol = row.end();
931 for (
auto col=row.begin(); col!=endcol; ++col)
933 for (
auto& entry : (*col)[localRowIdx])
935 if (col.index() == blockRow)
937 (*col)[localRowIdx][localRowIdx] = diagonal;
948 template <
class Container>
949 void setUnitRows(
const Container& unitRows,
const Container& auxRows )
951 setUnitRowImpl( unitRows, 1.0 );
952 setUnitRowImpl( auxRows, 0.0 );
958 reserve( Stencil<DomainSpaceType,RangeSpaceType>(domainSpace_, rangeSpace_),
true );
963 void reserve (
const std::vector< Set >& sparsityPattern )
965 reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ),
false );
969 template <
class Stencil>
970 void reserve(
const Stencil &stencil,
const bool proposeImplicit =
true )
973 if(sequence_ != domainSpace().sequence())
978 const bool implicit = proposeImplicit && MPIManager::numThreads() == 1;
981 auto nnz = stencil.maxNonZerosEstimate();
984 if (domainSpace_.begin() != domainSpace_.end() && rangeSpace_.begin() != rangeSpace_.end())
986 Stencil tmpStencil( stencil );
987 tmpStencil.fill( *(domainSpace_.begin()), *(rangeSpace_.begin()) );
988 nnz = tmpStencil.maxNonZerosEstimate();
991 matrix_.reset(
new MatrixType( rowMapper_.size(), colMapper_.size(), nnz, overflowFraction_ ) );
995 matrix_.reset(
new MatrixType( rowMapper_.size(), colMapper_.size() ) );
996 matrix().createEntries( stencil.globalStencil() );
999 sequence_ = domainSpace().sequence();
1004 template <
class HangingNodesType>
1005 void changeHangingNodes(
const HangingNodesType& hangingNodes)
1008 MatrixType* newMatrix =
new MatrixType(rowMapper_.size(), colMapper_.size());
1011 newMatrix->setup( *matrix_ , hangingNodes );
1017 matrix_.reset( newMatrix );
1022 void extractDiagonal( ColumnDiscreteFunctionType& diag )
const
1025 matrix().extractDiagonal( diag.blockVector() );
1029 bool rightPrecondition()
const
1035 void apply(
const ColumnDiscreteFunctionType& arg, RowDiscreteFunctionType& dest)
const
1037 matrixAdapter().apply( arg.blockVector(), dest.blockVector() );
1041 template <
class CDF,
class RDF>
1042 void apply(
const CDF& arg, RDF& dest)
const
1046 DUNE_THROW(NotImplemented,
"ISTLMatrixObj::apply called for DiscreteFunctions not specified in the template list");
1055 void createPreconditionMatrix()
1059 void print(std::ostream & s)
const
1064 const DomainSpaceType& domainSpace()
const
1066 return domainSpace_;
1068 const RangeSpaceType& rangeSpace()
const
1073 const RowMapperType& rowMapper()
const
1077 const ColMapperType& colMapper()
const
1083 ObjectType* newObject()
const
1085 return new ObjectType(*
this, domainSpace(), rangeSpace());
1092 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1093 LocalMatrixType localMatrix(
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
const
1095 return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
1102 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1103 LocalMatrixType localMatrix()
const
1105 return LocalMatrixType( localMatrixStack_ );
1108 LocalColumnObjectType localColumn(
const DomainEntityType &domainEntity )
const
1110 return LocalColumnObjectType ( *
this, domainEntity );
1114 template<
class LocalBlock,
class Operation >
1115 void applyToBlock (
const size_t row,
const size_t col,
1116 const LocalBlock &localBlock,
1117 Operation& operation )
1119 LittleBlockType& block = ( implicitModeActive() ) ? matrix().entry( row, col ) : matrix()[ row ][ col ];
1120 for(
int i = 0; i < littleRows; ++i )
1121 for(
int j = 0; j < littleCols; ++j )
1122 operation( block[ i ][ j ], localBlock[ i ][ j ] );
1125 template<
class LocalBlock >
1126 void setBlock (
const size_t row,
const size_t col,
1127 const LocalBlock &localBlock )
1129 typedef typename DomainSpaceType :: RangeFieldType Field;
1130 auto copy = [] ( Field& a,
const typename LocalBlock::field_type& b ) { a = b; };
1131 applyToBlock( row, col, localBlock, copy );
1134 template<
class LocalBlock >
1135 void addBlock (
const size_t row,
const size_t col,
1136 const LocalBlock &localBlock )
1138 typedef typename DomainSpaceType :: RangeFieldType Field;
1139 auto add = [] ( Field& a,
const typename LocalBlock::field_type& b ) { a += b; };
1140 applyToBlock( row, col, localBlock, add );
1144 template<
class LocalMatrix,
class Operation >
1145 void applyToLocalMatrix (
const DomainEntityType &domainEntity,
1146 const RangeEntityType &rangeEntity,
1147 const LocalMatrix &localMat,
1148 Operation& operation )
1150 typedef typename MatrixType::size_type Index;
1151 if( implicitModeActive() )
1153 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
1155 return matrix().entry( index.first, index.second );
1157 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
1159 LittleBlockType& block = blockAccess( index );
1160 for(
int i = 0; i < littleRows; ++i )
1161 for(
int j = 0; j < littleCols; ++j )
1162 operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
1164 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1168 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
1170 return matrix()[ index.first][ index.second ];
1172 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
1174 LittleBlockType& block = blockAccess( index );
1175 for(
int i = 0; i < littleRows; ++i )
1176 for(
int j = 0; j < littleCols; ++j )
1177 operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
1179 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1183 template<
class LocalMatrix >
1184 void setLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
1186 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1187 auto operation = [] ( RangeFieldType& a,
const RangeFieldType& b )
1191 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1194 template<
class LocalMatrix >
1195 void addLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
1197 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1198 auto operation = [] ( RangeFieldType& a,
const RangeFieldType& b )
1202 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1205 template<
class LocalMatrix,
class Scalar >
1206 void addScaledLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat,
const Scalar &s )
1208 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1209 auto operation = [ &s ] ( RangeFieldType& a,
const RangeFieldType& b )
1213 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1216 template<
class LocalMatrix >
1217 void getLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity, LocalMatrix &localMat )
const
1222 typedef typename MatrixType::size_type Index;
1223 auto functor = [ &localMat, this ] ( std::pair< int, int > local,
const std::pair< Index, Index > &global )
1225 for( std::size_t i = 0; i < littleRows; ++i )
1226 for( std::size_t j = 0; j < littleCols; ++j )
1227 localMat.set( local.first * littleRows + i, local.second *littleCols + j,
1228 matrix()[ global.first ][ global.second ][i][j] );
1231 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1235 void preConErrorMsg(
int preCon)
const
1242 matrixAdap_.reset(
nullptr );
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
vector space out of a tensor product of fields.
Definition: fvector.hh:92
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_THROW(E,...)
Definition: exceptions.hh:312
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
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
Define general preconditioner interface.