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/operator/matrix/spmatrix.hh>
32#include <dune/fem/io/parameter.hh>
33#include <dune/fem/operator/matrix/columnobject.hh>
34#include <dune/fem/storage/objectstack.hh>
36#include <dune/fem/operator/matrix/istlmatrixadapter.hh>
37#include <dune/fem/operator/matrix/istlpreconditioner.hh>
38#include <dune/fem/operator/matrix/functor.hh>
48 template<
class MatrixObject >
49 class ISTLLocalMatrix;
51 template <
class DomainSpaceImp,
class RangeSpaceImp,
54 class ISTLMatrixObject;
59 template <
class LittleBlockType,
class RowDiscreteFunctionImp,
class ColDiscreteFunctionImp = RowDiscreteFunctionImp>
60 class ImprovedBCRSMatrix :
public BCRSMatrix<LittleBlockType>
62 friend struct MatrixDimension<ImprovedBCRSMatrix>;
64 typedef RowDiscreteFunctionImp RowDiscreteFunctionType;
65 typedef ColDiscreteFunctionImp ColDiscreteFunctionType;
67 typedef BCRSMatrix<LittleBlockType> BaseType;
69 typedef BaseType MatrixBaseType;
70 typedef typename BaseType :: RowIterator RowIteratorType ;
71 typedef typename BaseType :: ColIterator ColIteratorType ;
73 typedef ImprovedBCRSMatrix< LittleBlockType, RowDiscreteFunctionImp, ColDiscreteFunctionImp > ThisType;
75 typedef typename BaseType :: size_type size_type;
80 typedef typename BaseType::field_type field_type;
83 typedef typename BaseType::block_type block_type;
86 typedef typename BaseType:: allocator_type allocator_type;
89 typedef typename BaseType :: row_type row_type;
92 typedef typename BaseType :: ColIterator ColIterator;
95 typedef typename BaseType :: ConstColIterator ConstColIterator;
98 typedef typename BaseType :: RowIterator RowIterator;
101 typedef typename BaseType :: ConstRowIterator ConstRowIterator;
104 typedef typename ColDiscreteFunctionType :: DiscreteFunctionSpaceType RangeSpaceType;
107 typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
110 typedef typename ColDiscreteFunctionType :: DofStorageType ColBlockVectorType;
113 typedef typename RangeSpaceType :: GridType :: Traits :: Communication CommunicationType ;
115 typedef typename BaseType :: BuildMode BuildMode ;
119 ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nnz,
double overflowFraction) :
120 BaseType (rows, cols, nnz, overflowFraction, BaseType::implicit)
124 ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nz = 0 ) :
125 BaseType (rows, cols, BaseType :: row_wise)
129 ImprovedBCRSMatrix( ) :
134 ImprovedBCRSMatrix(
const ImprovedBCRSMatrix& org) :
138 ConstRowIterator slicedBegin(
const size_type row )
const
140 return ConstRowIterator( this->r, row );
143 ConstRowIterator slicedEnd(
const size_type row )
const
145 return ConstRowIterator( this->r, row );
148 std::pair< size_type, size_type > sliceBeginEnd(
const size_type thread,
const size_type numThreads )
const
150 const size_type sliceSize = this->N() / numThreads ;
151 const size_type sliceBegin = thread * sliceSize ;
152 const size_type sliceEnd = ( thread == numThreads-1 ) ? this->N() : (sliceBegin + sliceSize) ;
153 return std::make_pair( sliceBegin, sliceEnd );
157 template <
class X,
class Y,
bool isMV>
158 void mvThreadedImpl(
const field_type& alpha,
159 const X& x, Y& y, std::integral_constant<bool, isMV> )
const
161 auto doMV = [
this, &alpha, &x, &y] ()
163 const auto slice = sliceBeginEnd( MPIManager::thread(), MPIManager::numThreads() );
164 const ConstRowIterator endi = slicedEnd( slice.second );
165 for (ConstRowIterator i = slicedBegin(slice.first); i!=endi; ++i)
167 if constexpr ( isMV )
170 ConstColIterator endj = (*i).end();
171 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
173 auto&& xj = Dune::Impl::asVector(x[j.index()]);
174 auto&& yi = Dune::Impl::asVector(y[i.index()]);
175 if constexpr ( isMV )
176 Dune::Impl::asMatrix(*j).umv(xj, yi);
178 Dune::Impl::asMatrix(*j).usmv(alpha, xj, yi);
185 MPIManager :: run ( doMV );
187 catch (
const SingleThreadModeError& e )
190 if constexpr ( isMV )
194 DUNE_THROW(SingleThreadModeError,
"ImprovedBCRSMatrix::usmvThreaded: Cannot recover from threading error. Disable threading!");
201 template <
class X,
class Y>
202 void usmvThreaded(
const field_type& alpha,
const X& x, Y& y )
const
204 mvThreadedImpl(alpha, x, y, std::false_type() );
207 template <
class X,
class Y>
208 void mvThreaded(
const X& x, Y& y )
const
210 mvThreadedImpl(1.0, x, y, std::true_type() );
213 template <
class SparsityPattern >
214 void createEntries(
const SparsityPattern& sparsityPattern)
216 const auto spEnd = sparsityPattern.end();
219 auto endcreate = this->createend();
220 for(
auto create = this->createbegin(); create != endcreate; ++create)
222 const auto row = sparsityPattern.find( create.index() );
224 if( row == spEnd )
continue;
227 for (
const auto& col : row->second )
228 create.insert( col );
235 for (
auto& row : *
this)
236 for (
auto& entry : row)
241 void unitRow(
const size_t row )
243 block_type idBlock( 0 );
244 for (
int i = 0; i < idBlock.rows; ++i)
247 auto& matRow = (*this)[ row ];
248 auto colIt = matRow.begin();
249 const auto& colEndIt = matRow.end();
250 for (; colIt != colEndIt; ++colIt)
252 if( colIt.index() == row )
260 template <
class HangingNodesType>
261 void setup(ThisType& oldMatrix,
const HangingNodesType& hangingNodes)
264 typedef std::set< std::pair<int, block_type> > LocalEntryType;
265 typedef std::map< int , LocalEntryType > EntriesType;
269 std::map< int , std::set<int> > indices;
271 auto rowend = oldMatrix.end();
272 for(
auto it = oldMatrix.begin(); it != rowend; ++it)
274 const auto row = it.index();
275 auto& localIndices = indices[ row ];
277 if( hangingNodes.isHangingNode( row ) )
280 const auto& cols = hangingNodes.associatedDofs( row );
281 const auto colSize = cols.size();
282 for(
auto i=0; i<colSize; ++i)
284 assert( ! hangingNodes.isHangingNode( cols[i].first ) );
287 auto& localColIndices = indices[ cols[i].first ];
288 auto& localEntry = entries[ cols[i].first ];
291 auto endj = (*it).end();
292 for (
auto j= (*it).begin(); j!=endj; ++j)
294 localColIndices.insert( j.index () );
295 localEntry.insert( std::make_pair( j.index(), (cols[i].second * (*j)) ));
300 localIndices.insert( row );
301 for(
auto i=0; i<colSize; ++i)
302 localIndices.insert( cols[i].first );
307 auto endj = (*it).end();
308 for (
auto j= (*it).begin(); j!=endj; ++j)
309 localIndices.insert( j.index () );
314 createEntries( indices );
317 auto rowit = oldMatrix.begin();
319 auto endcreate = this->end();
320 for(
auto create = this->begin(); create != endcreate; ++create, ++rowit )
322 assert( rowit != oldMatrix.end() );
324 const auto row = create.index();
325 if( hangingNodes.isHangingNode( row ) )
327 const auto& cols = hangingNodes.associatedDofs( row );
329 std::map< const int , block_type > colMap;
331 assert( block_type :: rows == 1 );
333 const auto colSize = cols.size();
334 for(
auto i=0; i<colSize; ++i)
335 colMap[ cols[i].first ] = -cols[i].second;
339 auto endj = (*create).end();
340 for (
auto j= (*create).begin(); j!=endj; ++j)
342 assert( colMap.find( j.index() ) != colMap.end() );
343 (*j) = colMap[ j.index() ];
347 else if ( entries.find( row ) == entries.end() )
349 auto colit = (*rowit).begin();
350 auto endj = (*create).end();
351 for (
auto j= (*create).begin(); j!=endj; ++j, ++colit )
353 assert( colit != (*rowit).end() );
359 std::map< int , block_type > oldCols;
362 auto colend = (*rowit).end();
363 for(
auto colit = (*rowit).begin(); colit != colend; ++colit)
364 oldCols[ colit.index() ] = 0;
367 auto entry = entries.find( row );
368 assert( entry != entries.end ());
371 auto endcol = (*entry).second.end();
372 for(
auto co = (*entry).second.begin(); co != endcol; ++co)
373 oldCols[ (*co).first ] = 0;
377 auto colend = (*rowit).end();
378 for(
auto colit = (*rowit).begin(); colit != colend; ++colit)
379 oldCols[ colit.index() ] += (*colit);
383 auto endcol = (*entry).second.end();
384 for(
auto co = (*entry).second.begin(); co != endcol; ++co)
385 oldCols[ (*co).first ] += (*co).second;
388 auto endj = (*create).end();
389 for (
auto j= (*create).begin(); j!=endj; ++j )
391 auto colEntry = oldCols.find( j.index() );
392 if( colEntry != oldCols.end() )
393 (*j) = (*colEntry).second;
402 void extractDiagonal( ColBlockVectorType& diag )
const
404 const auto endi = this->end();
405 for (
auto i = this->begin(); i!=endi; ++i)
408 const auto row = i.index();
409 auto entry = (*i).find( row );
410 const LittleBlockType& block = (*entry);
411 enum { blockSize = LittleBlockType :: rows };
412 for(
auto l=0; l<blockSize; ++l )
413 diag[ row ][ l ] = block[ l ][ l ];
419 field_type operator()(
const std::size_t row,
const std::size_t col)
const
421 const std::size_t blockRow(row/(LittleBlockType :: rows));
422 const std::size_t localRowIdx(row%(LittleBlockType :: rows));
423 const std::size_t blockCol(col/(LittleBlockType :: cols));
424 const std::size_t localColIdx(col%(LittleBlockType :: cols));
426 const auto& matrixRow(this->
operator[](blockRow));
427 auto entry = matrixRow.find( blockCol );
428 const LittleBlockType& block = (*entry);
429 return block[localRowIdx][localColIdx];
434 void set(
const std::size_t row,
const std::size_t col, field_type value)
436 const std::size_t blockRow(row/(LittleBlockType :: rows));
437 const std::size_t localRowIdx(row%(LittleBlockType :: rows));
438 const std::size_t blockCol(col/(LittleBlockType :: cols));
439 const std::size_t localColIdx(col%(LittleBlockType :: cols));
441 auto& matrixRow(this->
operator[](blockRow));
442 auto entry = matrixRow.find( blockCol );
443 LittleBlockType& block = (*entry);
444 block[localRowIdx][localColIdx] = value;
448 void print(std::ostream& s=std::cout,
unsigned int offset=0)
const
451 const auto endi=this->end();
452 for (
auto i=this->begin(); i!=endi; ++i)
454 const auto endj = (*i).end();
455 for (
auto j=(*i).begin(); j!=endj; ++j)
456 if( (*j).infinity_norm() > 1.e-15)
457 s << i.index()+offset <<
" " << j.index()+offset <<
" " << *j << std::endl;
467 template<
class MatrixObject >
468 struct ISTLLocalMatrixTraits
470 typedef typename MatrixObject::DomainSpaceType DomainSpaceType;
471 typedef typename MatrixObject::RangeSpaceType RangeSpaceType;
472 typedef typename DomainSpaceType::RangeFieldType RangeFieldType;
474 typedef ISTLLocalMatrix< MatrixObject > LocalMatrixType;
475 typedef typename MatrixObject::MatrixType::block_type LittleBlockType;
482 template<
class MatrixObject >
483 class ISTLLocalMatrix
484 :
public LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > >
488 typedef LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > > BaseType;
491 typedef MatrixObject MatrixObjectType;
493 typedef typename MatrixObjectType::MatrixType MatrixType;
495 typedef typename MatrixType::block_type LittleBlockType;
497 typedef typename MatrixType :: size_type size_type;
499 typedef typename MatrixObjectType::DomainSpaceType DomainSpaceType;
500 typedef typename MatrixObjectType::RangeSpaceType RangeSpaceType;
502 typedef typename MatrixObjectType::DomainEntityType DomainEntityType;
503 typedef typename MatrixObjectType::RangeEntityType RangeEntityType;
506 typedef typename DomainSpaceType::RangeFieldType DofType;
507 typedef typename MatrixType::row_type RowType;
510 typedef typename DomainSpaceType::BlockMapperType ColMapperType;
512 typedef typename RangeSpaceType::BlockMapperType RowMapperType;
514 static const int littleCols = MatrixObjectType::littleCols;
515 static const int littleRows = MatrixObjectType::littleRows;
517 typedef typename MatrixType::size_type Index;
522 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
523 ISTLLocalMatrix (
const MatrixObjectType& mObj,
const DomainSpaceType& domainSpace,
const RangeSpaceType& rangeSpace )
524 : BaseType( domainSpace, rangeSpace ),
525 rowMapper_( rangeSpace.blockMapper() ),
526 colMapper_( domainSpace.blockMapper() ),
527 numRows_( rowMapper_.maxNumDofs() ),
528 numCols_( colMapper_.maxNumDofs() ),
535 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
536 ISTLLocalMatrix (
const ISTLLocalMatrix& org )
538 rowMapper_(org.rowMapper_),
539 colMapper_(org.colMapper_),
540 numRows_( org.numRows_ ),
541 numCols_( org.numCols_ ),
542 matrices_(org.matrices_),
543 matrixObj_(org.matrixObj_)
547 void init (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
550 BaseType :: init ( domainEntity, rangeEntity );
552 numRows_ = rowMapper_.numDofs( rangeEntity );
553 numCols_ = colMapper_.numDofs( domainEntity );
556 matrices_.resize( numRows_ );
557 for(
auto& row : matrices_ )
559 row.resize( numCols_,
nullptr );
562 if( matrixObj_.implicitModeActive() )
564 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
566 return matrixObj_.matrix().entry( index.first, index.second );
568 initBlocks( blockAccess, domainEntity, rangeEntity );
572 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
574 return matrixObj_.matrix()[ index.first ][ index.second ];
576 initBlocks( blockAccess, domainEntity, rangeEntity );
580 template <
class BlockAccess>
581 void initBlocks( BlockAccess& blockAccess,
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
583 auto functor = [
this, &blockAccess ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
585 matrices_[ local.first ][ local.second ] = &blockAccess( index );
588 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
593 void check(
int localRow,
int localCol)
const
596 const std::size_t row = (int) localRow / littleRows;
597 const std::size_t col = (int) localCol / littleCols;
598 const int lRow = localRow%littleRows;
599 const int lCol = localCol%littleCols;
600 assert( row < matrices_.size() ) ;
601 assert( col < matrices_[row].
size() );
602 assert( lRow < littleRows );
603 assert( lCol < littleCols );
607 DofType& getValue(
const int localRow,
const int localCol)
609 const int row = (int) localRow / littleRows;
610 const int col = (int) localCol / littleCols;
611 const int lRow = localRow%littleRows;
612 const int lCol = localCol%littleCols;
613 return (*matrices_[row][col])[lRow][lCol];
617 const DofType
get(
const int localRow,
const int localCol)
const
619 const int row = (int) localRow / littleRows;
620 const int col = (int) localCol / littleCols;
621 const int lRow = localRow%littleRows;
622 const int lCol = localCol%littleCols;
623 return (*matrices_[row][col])[lRow][lCol];
626 void scale (
const DofType& scalar)
628 const size_type nrows = matrices_.size();
629 for(size_type i=0; i<nrows; ++i)
631 const size_type ncols = matrices_[i].size();
632 for(size_type j=0; j<ncols; ++j)
634 (*matrices_[i][j]) *= scalar;
639 void add(
const int localRow,
const int localCol ,
const DofType value)
642 check(localRow,localCol);
644 getValue(localRow,localCol) += value;
647 void set(
const int localRow,
const int localCol ,
const DofType value)
650 check(localRow,localCol);
652 getValue(localRow,localCol) = value;
656 void unitRow(
const int localRow)
658 const int row = (int) localRow / littleRows;
659 const int lRow = localRow%littleRows;
662 doClearRow( row, lRow );
665 (*matrices_[row][row])[lRow][lRow] = 1;
671 const size_type nrows = matrices_.size();
672 for(size_type i=0; i<nrows; ++i)
674 const size_type ncols = matrices_[i].size();
675 for(size_type j=0; j<ncols; ++j)
677 (*matrices_[i][j]) = (DofType) 0;
683 void clearRow (
const int localRow )
685 const int row = (int) localRow / littleRows;
686 const int lRow = localRow%littleRows;
689 doClearRow( row, lRow );
698 void doClearRow (
const int row,
const int lRow )
701 const auto col = this->columns();
702 for(
auto localCol=0; localCol<col; ++localCol)
704 const int col = (int) localCol / littleCols;
705 const int lCol = localCol%littleCols;
706 (*matrices_[row][col])[lRow][lCol] = 0;
712 const RowMapperType& rowMapper_;
713 const ColMapperType& colMapper_;
720 std::vector< std::vector< LittleBlockType* > > matrices_;
723 const MatrixObjectType& matrixObj_;
728 template <
class DomainSpaceImp,
class RangeSpaceImp,
class DomainBlock,
class RangeBlock >
729 class ISTLMatrixObject
733 typedef DomainSpaceImp DomainSpaceType;
735 typedef RangeSpaceImp RangeSpaceType;
738 typedef ISTLMatrixObject< DomainSpaceImp, RangeSpaceImp, DomainBlock, RangeBlock > ThisType;
739 typedef ThisType PreconditionMatrixType;
741 typedef typename DomainSpaceType::GridType GridType;
743 typedef typename RangeSpaceType :: EntityType RangeEntityType ;
744 typedef typename DomainSpaceType :: EntityType DomainEntityType ;
746 enum { littleCols = DomainSpaceType :: localBlockSize };
747 enum { littleRows = RangeSpaceType :: localBlockSize };
749 typedef FieldMatrix<typename DomainSpaceType :: RangeFieldType, littleRows, littleCols> LittleBlockType;
750 typedef LittleBlockType block_type;
751 typedef LittleBlockType MatrixBlockType;
753 typedef ISTLBlockVectorDiscreteFunction< RangeSpaceType, RangeBlock > RowDiscreteFunctionType;
754 typedef ISTLBlockVectorDiscreteFunction< DomainSpaceType, DomainBlock > ColumnDiscreteFunctionType;
756 typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
757 typedef typename ColumnDiscreteFunctionType :: DofStorageType ColumnBlockVectorType;
759 typedef typename RangeSpaceType :: BlockMapperType RowMapperType;
760 typedef typename DomainSpaceType :: BlockMapperType ColMapperType;
763 typedef ImprovedBCRSMatrix< LittleBlockType , ColumnDiscreteFunctionType , RowDiscreteFunctionType > MatrixType;
764 typedef ISTLParallelMatrixAdapterInterface< MatrixType > MatrixAdapterType;
767 typedef ISTLLocalMatrix<ThisType> ObjectType;
768 friend class ISTLLocalMatrix<ThisType>;
770 typedef ThisType LocalMatrixFactoryType;
773 typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
774 typedef ColumnObject< ThisType > LocalColumnObjectType;
777 const DomainSpaceType& domainSpace_;
778 const RangeSpaceType& rangeSpace_;
781 RowMapperType& rowMapper_;
783 ColMapperType& colMapper_;
788 mutable std::unique_ptr< MatrixType > matrix_;
790 mutable LocalMatrixStackType localMatrixStack_;
792 mutable std::unique_ptr< MatrixAdapterType > matrixAdap_;
794 const double overflowFraction_;
795 const bool threading_ ;
797 ISTLSolverParameter param_;
799 ISTLMatrixObject(
const ISTLMatrixObject&) =
delete;
804 ISTLMatrixObject (
const DomainSpaceType &domainSpace,
const RangeSpaceType &rangeSpace,
const ISTLSolverParameter& param = ISTLSolverParameter() ) :
805 domainSpace_(domainSpace)
806 , rangeSpace_(rangeSpace)
807 , rowMapper_( rangeSpace.blockMapper() )
808 , colMapper_( domainSpace.blockMapper() )
811 , localMatrixStack_( *this )
812 , overflowFraction_( param.overflowFraction() )
813 , threading_( param.threading() )
820 MatrixType& matrix()
const
827 bool threading()
const {
return threading_; }
829 void printTexInfo(std::ostream& out)
const
831 out <<
"ISTL MatrixObj: ";
836 std::string preconditionName()
const
841 void createMatrixAdapter (
const ISTLSolverParameter& param )
const
845 matrixAdap_ = ISTLMatrixAdapterFactory< ThisType > :: matrixAdapter( *
this, param );
847 assert( matrixAdap_ );
851 MatrixAdapterType& matrixAdapter(
const ISTLSolverParameter& parameter )
const
853 const ISTLSolverParameter* param =
dynamic_cast< const ISTLSolverParameter*
> (¶meter);
855 createMatrixAdapter( *param );
858 ISTLSolverParameter newparam( parameter );
859 createMatrixAdapter( newparam );
867 MatrixAdapterType& matrixAdapter()
const
869 return matrixAdapter( param_ );
873 MatrixType& exportMatrix()
const
879 bool implicitModeActive()
const
885 if( matrix().buildMode() == MatrixType::implicit && matrix().buildStage() == MatrixType::building )
891 void finalizeAssembly()
const
894 const_cast< ThisType&
> (*this).compress();
900 if( implicitModeActive() )
914 void unitRow(
const size_t row )
916 matrix().unitRow( row );
920 template <
class Container>
921 void setUnitRowImpl(
const Container &rows,
const double diagonal )
923 for (
const auto& r : rows)
925 const std::size_t blockRow( r/(LittleBlockType :: rows) );
926 const std::size_t localRowIdx( r%(LittleBlockType :: rows) );
927 auto& row = matrix()[blockRow];
928 const auto endcol = row.end();
932 for (
auto col=row.begin(); col!=endcol; ++col)
934 for (
auto& entry : (*col)[localRowIdx])
936 if (col.index() == blockRow)
938 (*col)[localRowIdx][localRowIdx] = diagonal;
949 template <
class Container>
950 void setUnitRows(
const Container& unitRows,
const Container& auxRows )
952 setUnitRowImpl( unitRows, 1.0 );
953 setUnitRowImpl( auxRows, 0.0 );
959 reserve( Stencil<DomainSpaceType,RangeSpaceType>(domainSpace_, rangeSpace_),
true );
964 void reserve (
const std::vector< Set >& sparsityPattern )
966 reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ),
false );
970 template <
class Stencil>
971 void reserve(
const Stencil &stencil,
const bool proposeImplicit =
true )
974 if(sequence_ != domainSpace().sequence())
979 const bool implicit = proposeImplicit && MPIManager::numThreads() == 1;
982 auto nnz = stencil.maxNonZerosEstimate();
985 if (domainSpace_.begin() != domainSpace_.end() && rangeSpace_.begin() != rangeSpace_.end())
987 Stencil tmpStencil( stencil );
988 tmpStencil.fill( *(domainSpace_.begin()), *(rangeSpace_.begin()) );
989 nnz = tmpStencil.maxNonZerosEstimate();
992 matrix_.reset(
new MatrixType( rowMapper_.size(), colMapper_.size(), nnz, overflowFraction_ ) );
996 matrix_.reset(
new MatrixType( rowMapper_.size(), colMapper_.size() ) );
997 matrix().createEntries( stencil.globalStencil() );
1000 sequence_ = domainSpace().sequence();
1005 template <
class HangingNodesType>
1006 void changeHangingNodes(
const HangingNodesType& hangingNodes)
1009 MatrixType* newMatrix =
new MatrixType(rowMapper_.size(), colMapper_.size());
1012 newMatrix->setup( *matrix_ , hangingNodes );
1018 matrix_.reset( newMatrix );
1023 void extractDiagonal( ColumnDiscreteFunctionType& diag )
const
1026 matrix().extractDiagonal( diag.blockVector() );
1030 bool rightPrecondition()
const
1036 void apply(
const ColumnDiscreteFunctionType& arg, RowDiscreteFunctionType& dest)
const
1038 matrixAdapter().apply( arg.blockVector(), dest.blockVector() );
1042 template <
class CDF,
class RDF>
1043 void apply(
const CDF& arg, RDF& dest)
const
1047 DUNE_THROW(NotImplemented,
"ISTLMatrixObj::apply called for DiscreteFunctions not specified in the template list");
1056 void createPreconditionMatrix()
1060 void print(std::ostream & s)
const
1065 const DomainSpaceType& domainSpace()
const
1067 return domainSpace_;
1069 const RangeSpaceType& rangeSpace()
const
1074 const RowMapperType& rowMapper()
const
1078 const ColMapperType& colMapper()
const
1084 ObjectType* newObject()
const
1086 return new ObjectType(*
this, domainSpace(), rangeSpace());
1093 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1094 LocalMatrixType localMatrix(
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity )
const
1096 return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
1103 [[deprecated(
"Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1104 LocalMatrixType localMatrix()
const
1106 return LocalMatrixType( localMatrixStack_ );
1109 LocalColumnObjectType localColumn(
const DomainEntityType &domainEntity )
const
1111 return LocalColumnObjectType ( *
this, domainEntity );
1115 template<
class LocalBlock,
class Operation >
1116 void applyToBlock (
const size_t row,
const size_t col,
1117 const LocalBlock &localBlock,
1118 Operation& operation )
1120 LittleBlockType& block = ( implicitModeActive() ) ? matrix().entry( row, col ) : matrix()[ row ][ col ];
1121 for(
int i = 0; i < littleRows; ++i )
1122 for(
int j = 0; j < littleCols; ++j )
1123 operation( block[ i ][ j ], localBlock[ i ][ j ] );
1126 template<
class LocalBlock >
1127 void setBlock (
const size_t row,
const size_t col,
1128 const LocalBlock &localBlock )
1130 typedef typename DomainSpaceType :: RangeFieldType Field;
1131 auto copy = [] ( Field& a,
const typename LocalBlock::field_type& b ) { a = b; };
1132 applyToBlock( row, col, localBlock, copy );
1135 template<
class LocalBlock >
1136 void addBlock (
const size_t row,
const size_t col,
1137 const LocalBlock &localBlock )
1139 typedef typename DomainSpaceType :: RangeFieldType Field;
1140 auto add = [] ( Field& a,
const typename LocalBlock::field_type& b ) { a += b; };
1141 applyToBlock( row, col, localBlock, add );
1145 template<
class LocalMatrix,
class Operation >
1146 void applyToLocalMatrix (
const DomainEntityType &domainEntity,
1147 const RangeEntityType &rangeEntity,
1148 const LocalMatrix &localMat,
1149 Operation& operation )
1151 typedef typename MatrixType::size_type Index;
1152 if( implicitModeActive() )
1154 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
1156 return matrix().entry( index.first, index.second );
1158 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
1160 LittleBlockType& block = blockAccess( index );
1161 for(
int i = 0; i < littleRows; ++i )
1162 for(
int j = 0; j < littleCols; ++j )
1163 operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
1165 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1169 auto blockAccess = [ this ] (
const std::pair< Index, Index > &index ) -> LittleBlockType&
1171 return matrix()[ index.first][ index.second ];
1173 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local,
const std::pair< Index, Index > &index )
1175 LittleBlockType& block = blockAccess( index );
1176 for(
int i = 0; i < littleRows; ++i )
1177 for(
int j = 0; j < littleCols; ++j )
1178 operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
1180 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1184 template<
class LocalMatrix >
1185 void setLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
1187 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1188 auto operation = [] ( RangeFieldType& a,
const RangeFieldType& b )
1192 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1195 template<
class LocalMatrix >
1196 void addLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat )
1198 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1199 auto operation = [] ( RangeFieldType& a,
const RangeFieldType& b )
1203 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1206 template<
class LocalMatrix,
class Scalar >
1207 void addScaledLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity,
const LocalMatrix &localMat,
const Scalar &s )
1209 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1210 auto operation = [ &s ] ( RangeFieldType& a,
const RangeFieldType& b )
1214 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1217 template<
class LocalMatrix >
1218 void getLocalMatrix (
const DomainEntityType &domainEntity,
const RangeEntityType &rangeEntity, LocalMatrix &localMat )
const
1223 typedef typename MatrixType::size_type Index;
1224 auto functor = [ &localMat, this ] ( std::pair< int, int > local,
const std::pair< Index, Index > &global )
1226 for( std::size_t i = 0; i < littleRows; ++i )
1227 for( std::size_t j = 0; j < littleCols; ++j )
1228 localMat.set( local.first * littleRows + i, local.second *littleCols + j,
1229 matrix()[ global.first ][ global.second ][i][j] );
1232 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1236 void preConErrorMsg(
int preCon)
const
1243 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:91
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, m)
Definition: exceptions.hh:218
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.