DUNE-FEM (unstable)

istlmatrix.hh
1#ifndef DUNE_FEM_ISTLMATRIXWRAPPER_HH
2#define DUNE_FEM_ISTLMATRIXWRAPPER_HH
3
4#if HAVE_DUNE_ISTL
5
6//- system includes
7#include <vector>
8#include <iostream>
9#include <fstream>
10#include <algorithm>
11#include <set>
12#include <map>
13#include <string>
14
15//- Dune common includes
18
19//- Dune istl includes
20#include <dune/istl/bvector.hh>
23
24//- Dune fem includes
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>
34
35#include <dune/fem/operator/matrix/istlmatrixadapter.hh>
36#include <dune/fem/operator/matrix/istlpreconditioner.hh>
37#include <dune/fem/operator/matrix/functor.hh>
38
39namespace Dune
40{
41
42 namespace Fem
43 {
44
45
47 template< class MatrixObject >
48 class ISTLLocalMatrix;
49
50 template <class DomainSpaceImp, class RangeSpaceImp,
53 class ISTLMatrixObject;
54
56 // --ISTLMatrixHandle
58 template <class LittleBlockType, class RowDiscreteFunctionImp, class ColDiscreteFunctionImp = RowDiscreteFunctionImp>
59 class ImprovedBCRSMatrix : public BCRSMatrix<LittleBlockType>
60 {
61 friend struct MatrixDimension<ImprovedBCRSMatrix>;
62 public:
63 typedef RowDiscreteFunctionImp RowDiscreteFunctionType;
64 typedef ColDiscreteFunctionImp ColDiscreteFunctionType;
65
66 typedef BCRSMatrix<LittleBlockType> BaseType;
68 typedef BaseType MatrixBaseType;
69 typedef typename BaseType :: RowIterator RowIteratorType ;
70 typedef typename BaseType :: ColIterator ColIteratorType ;
71
72 typedef ImprovedBCRSMatrix< LittleBlockType, RowDiscreteFunctionImp, ColDiscreteFunctionImp > ThisType;
73
74 typedef typename BaseType :: size_type size_type;
75
76 //===== type definitions and constants
77
79 typedef typename BaseType::field_type field_type;
80
82 typedef typename BaseType::block_type block_type;
83
85 typedef typename BaseType:: allocator_type allocator_type;
86
88 typedef typename BaseType :: row_type row_type;
89
91 typedef typename BaseType :: ColIterator ColIterator;
92
94 typedef typename BaseType :: ConstColIterator ConstColIterator;
95
97 typedef typename BaseType :: RowIterator RowIterator;
98
100 typedef typename BaseType :: ConstRowIterator ConstRowIterator;
101
103 typedef typename ColDiscreteFunctionType :: DiscreteFunctionSpaceType RangeSpaceType;
104
106 typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
107
109 typedef typename ColDiscreteFunctionType :: DofStorageType ColBlockVectorType;
110
112 typedef typename RangeSpaceType :: GridType :: Traits :: Communication CommunicationType ;
113
114 typedef typename BaseType :: BuildMode BuildMode ;
115
116 public:
118 ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nnz, double overflowFraction) :
119 BaseType (rows, cols, nnz, overflowFraction, BaseType::implicit)
120 {}
121
123 ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nz = 0 ) :
124 BaseType (rows, cols, BaseType :: row_wise)
125 {}
126
128 ImprovedBCRSMatrix( ) :
129 BaseType ()
130 {}
131
133 ImprovedBCRSMatrix(const ImprovedBCRSMatrix& org) :
134 BaseType(org)
135 {}
136
137 ConstRowIterator slicedBegin( const size_type row ) const
138 {
139 return ConstRowIterator( this->r, row );
140 }
141
142 ConstRowIterator slicedEnd( const size_type row ) const
143 {
144 return ConstRowIterator( this->r, row );
145 }
146
147 std::pair< size_type, size_type > sliceBeginEnd( const size_type thread, const size_type numThreads ) const
148 {
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 );
153 }
154
155 protected:
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
159 {
160 auto doMV = [this, &alpha, &x, &y] ()
161 {
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)
165 {
166 if constexpr ( isMV )
167 y[i.index()] = 0;
168
169 ConstColIterator endj = (*i).end();
170 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
171 {
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);
176 else
177 Dune::Impl::asMatrix(*j).usmv(alpha, xj, yi);
178 }
179 }
180 };
181
182 try {
183 // execute in parallel
184 MPIManager :: run ( doMV );
185 }
186 catch ( const SingleThreadModeError& e )
187 {
188 // serial version
189 if constexpr ( isMV )
190 this->mv( x, y );
191 else
192 {
193 DUNE_THROW(SingleThreadModeError,"ImprovedBCRSMatrix::usmvThreaded: Cannot recover from threading error. Disable threading!");
194 // this->usmv( alpha, x, y );
195 }
196 }
197 }
198
199 public:
200 template <class X, class Y>
201 void usmvThreaded( const field_type& alpha, const X& x, Y& y ) const
202 {
203 mvThreadedImpl(alpha, x, y, std::false_type() );
204 }
205
206 template <class X, class Y>
207 void mvThreaded( const X& x, Y& y ) const
208 {
209 mvThreadedImpl(1.0, x, y, std::true_type() );
210 }
211
212 template < class SparsityPattern >
213 void createEntries(const SparsityPattern& sparsityPattern)
214 {
215 const auto spEnd = sparsityPattern.end();
216
217 // insert map of indices into matrix
218 auto endcreate = this->createend();
219 for(auto create = this->createbegin(); create != endcreate; ++create)
220 {
221 const auto row = sparsityPattern.find( create.index() );
222 // if a row is empty then do nothing
223 if( row == spEnd ) continue;
224
225 // insert all indices for this row
226 for ( const auto& col : row->second )
227 create.insert( col );
228 }
229 }
230
232 void clear()
233 {
234 for (auto& row : *this)
235 for (auto& entry : row)
236 entry = 0;
237 }
238
240 void unitRow( const size_t row )
241 {
242 block_type idBlock( 0 );
243 for (int i = 0; i < idBlock.rows; ++i)
244 idBlock[i][i] = 1.0;
245
246 auto& matRow = (*this)[ row ];
247 auto colIt = matRow.begin();
248 const auto& colEndIt = matRow.end();
249 for (; colIt != colEndIt; ++colIt)
250 {
251 if( colIt.index() == row )
252 *colIt = idBlock;
253 else
254 *colIt = 0.0;
255 }
256 }
257
259 template <class HangingNodesType>
260 void setup(ThisType& oldMatrix, const HangingNodesType& hangingNodes)
261 {
262 // necessary because element traversal not necessarily is in ascending order
263 typedef std::set< std::pair<int, block_type> > LocalEntryType;
264 typedef std::map< int , LocalEntryType > EntriesType;
265 EntriesType entries;
266
267 // map of indices
268 std::map< int , std::set<int> > indices;
269 // not insert map of indices into matrix
270 auto rowend = oldMatrix.end();
271 for(auto it = oldMatrix.begin(); it != rowend; ++it)
272 {
273 const auto row = it.index();
274 auto& localIndices = indices[ row ];
275
276 if( hangingNodes.isHangingNode( row ) )
277 {
278 // insert columns into other columns
279 const auto& cols = hangingNodes.associatedDofs( row );
280 const auto colSize = cols.size();
281 for(auto i=0; i<colSize; ++i)
282 {
283 assert( ! hangingNodes.isHangingNode( cols[i].first ) );
284
285 // get local indices of col
286 auto& localColIndices = indices[ cols[i].first ];
287 auto& localEntry = entries[ cols[i].first ];
288
289 // copy from old matrix
290 auto endj = (*it).end();
291 for (auto j= (*it).begin(); j!=endj; ++j)
292 {
293 localColIndices.insert( j.index () );
294 localEntry.insert( std::make_pair( j.index(), (cols[i].second * (*j)) ));
295 }
296 }
297
298 // insert diagonal and hanging columns
299 localIndices.insert( row );
300 for(auto i=0; i<colSize; ++i)
301 localIndices.insert( cols[i].first );
302 }
303 else
304 {
305 // copy from old matrix
306 auto endj = (*it).end();
307 for (auto j= (*it).begin(); j!=endj; ++j)
308 localIndices.insert( j.index () );
309 }
310 }
311
312 // create matrix from entry map
313 createEntries( indices );
314
315 // not insert map of indices into matrix
316 auto rowit = oldMatrix.begin();
317
318 auto endcreate = this->end();
319 for(auto create = this->begin(); create != endcreate; ++create, ++rowit )
320 {
321 assert( rowit != oldMatrix.end() );
322
323 const auto row = create.index();
324 if( hangingNodes.isHangingNode( row ) )
325 {
326 const auto& cols = hangingNodes.associatedDofs( row );
327
328 std::map< const int , block_type > colMap;
329 // only working for block size 1 ath the moment
330 assert( block_type :: rows == 1 );
331 // insert columns into map
332 const auto colSize = cols.size();
333 for( auto i=0; i<colSize; ++i)
334 colMap[ cols[i].first ] = -cols[i].second;
335 // insert diagonal into map
336 colMap[ row ] = 1;
337
338 auto endj = (*create).end();
339 for (auto j= (*create).begin(); j!=endj; ++j)
340 {
341 assert( colMap.find( j.index() ) != colMap.end() );
342 (*j) = colMap[ j.index() ];
343 }
344 }
345 // if entries are equal, just copy
346 else if ( entries.find( row ) == entries.end() )
347 {
348 auto colit = (*rowit).begin();
349 auto endj = (*create).end();
350 for (auto j= (*create).begin(); j!=endj; ++j, ++colit )
351 {
352 assert( colit != (*rowit).end() );
353 (*j) = (*colit);
354 }
355 }
356 else
357 {
358 std::map< int , block_type > oldCols;
359
360 {
361 auto colend = (*rowit).end();
362 for(auto colit = (*rowit).begin(); colit != colend; ++colit)
363 oldCols[ colit.index() ] = 0;
364 }
365
366 auto entry = entries.find( row );
367 assert( entry != entries.end ());
368
369 {
370 auto endcol = (*entry).second.end();
371 for( auto co = (*entry).second.begin(); co != endcol; ++co)
372 oldCols[ (*co).first ] = 0;
373 }
374
375 {
376 auto colend = (*rowit).end();
377 for(auto colit = (*rowit).begin(); colit != colend; ++colit)
378 oldCols[ colit.index() ] += (*colit);
379 }
380
381 {
382 auto endcol = (*entry).second.end();
383 for( auto co = (*entry).second.begin(); co != endcol; ++co)
384 oldCols[ (*co).first ] += (*co).second;
385 }
386
387 auto endj = (*create).end();
388 for (auto j= (*create).begin(); j!=endj; ++j )
389 {
390 auto colEntry = oldCols.find( j.index() );
391 if( colEntry != oldCols.end() )
392 (*j) = (*colEntry).second;
393 else
394 abort();
395 }
396 }
397 }
398 }
399
401 void extractDiagonal( ColBlockVectorType& diag ) const
402 {
403 const auto endi = this->end();
404 for (auto i = this->begin(); i!=endi; ++i)
405 {
406 // get diagonal entry of matrix
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 ];
413 }
414 }
415
418 field_type operator()(const std::size_t row, const std::size_t col) const
419 {
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));
424
425 const auto& matrixRow(this->operator[](blockRow));
426 auto entry = matrixRow.find( blockCol );
427 const LittleBlockType& block = (*entry);
428 return block[localRowIdx][localColIdx];
429 }
430
433 void set(const std::size_t row, const std::size_t col, field_type value)
434 {
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));
439
440 auto& matrixRow(this->operator[](blockRow));
441 auto entry = matrixRow.find( blockCol );
442 LittleBlockType& block = (*entry);
443 block[localRowIdx][localColIdx] = value;
444 }
445
447 void print(std::ostream& s=std::cout, unsigned int offset=0) const
448 {
449 s.precision( 6 );
450 const auto endi=this->end();
451 for (auto i=this->begin(); i!=endi; ++i)
452 {
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;
457 }
458 }
459 };
460
461
462
463 // ISTLLocalMatrixTraits
464 // ---------------------
465
466 template< class MatrixObject >
467 struct ISTLLocalMatrixTraits
468 {
469 typedef typename MatrixObject::DomainSpaceType DomainSpaceType;
470 typedef typename MatrixObject::RangeSpaceType RangeSpaceType;
471 typedef typename DomainSpaceType::RangeFieldType RangeFieldType;
472
473 typedef ISTLLocalMatrix< MatrixObject > LocalMatrixType;
474 typedef typename MatrixObject::MatrixType::block_type LittleBlockType;
475 };
476
477
478 // ISTLLocalMatrix
479 // ---------------
480
481 template< class MatrixObject >
482 class ISTLLocalMatrix
483 : public LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > >
484 {
485 public:
487 typedef LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > > BaseType;
488
490 typedef MatrixObject MatrixObjectType;
492 typedef typename MatrixObjectType::MatrixType MatrixType;
494 typedef typename MatrixType::block_type LittleBlockType;
495 // type of row and column indices
496 typedef typename MatrixType :: size_type size_type;
497
498 typedef typename MatrixObjectType::DomainSpaceType DomainSpaceType;
499 typedef typename MatrixObjectType::RangeSpaceType RangeSpaceType;
500
501 typedef typename MatrixObjectType::DomainEntityType DomainEntityType;
502 typedef typename MatrixObjectType::RangeEntityType RangeEntityType;
503
505 typedef typename DomainSpaceType::RangeFieldType DofType;
506 typedef typename MatrixType::row_type RowType;
507
509 typedef typename DomainSpaceType::BlockMapperType ColMapperType;
511 typedef typename RangeSpaceType::BlockMapperType RowMapperType;
512
513 static const int littleCols = MatrixObjectType::littleCols;
514 static const int littleRows = MatrixObjectType::littleRows;
515
516 typedef typename MatrixType::size_type Index;
517
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() ),
528 matrixObj_( mObj )
529 {}
530
534 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
535 ISTLLocalMatrix ( const ISTLLocalMatrix& org )
536 : BaseType( 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_)
543 {}
544
546 void init ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
547 {
548 // initialize base functions sets
549 BaseType :: init ( domainEntity, rangeEntity );
550
551 numRows_ = rowMapper_.numDofs( rangeEntity );
552 numCols_ = colMapper_.numDofs( domainEntity );
553
554 // resize matrix pointer storage for new row/col numbers
555 matrices_.resize( numRows_ );
556 for( auto& row : matrices_ )
557 {
558 row.resize( numCols_, nullptr );
559 }
560
561 if( matrixObj_.implicitModeActive() )
562 {
563 auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
564 {
565 return matrixObj_.matrix().entry( index.first, index.second );
566 };
567 initBlocks( blockAccess, domainEntity, rangeEntity );
568 }
569 else
570 {
571 auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
572 {
573 return matrixObj_.matrix()[ index.first ][ index.second ];
574 };
575 initBlocks( blockAccess, domainEntity, rangeEntity );
576 }
577 }
578
579 template <class BlockAccess>
580 void initBlocks( BlockAccess& blockAccess, const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
581 {
582 auto functor = [ this, &blockAccess ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
583 {
584 matrices_[ local.first ][ local.second ] = &blockAccess( index );
585 };
586
587 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
588 }
589
590 private:
591 // check whether given (row,col) pair is valid
592 void check(int localRow, int localCol) const
593 {
594#ifndef NDEBUG
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 );
603#endif
604 }
605
606 DofType& getValue(const int localRow, const int localCol)
607 {
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];
613 }
614
615 public:
616 const DofType get(const int localRow, const int localCol) const
617 {
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];
623 }
624
625 void scale (const DofType& scalar)
626 {
627 const size_type nrows = matrices_.size();
628 for(size_type i=0; i<nrows; ++i)
629 {
630 const size_type ncols = matrices_[i].size();
631 for(size_type j=0; j<ncols; ++j)
632 {
633 (*matrices_[i][j]) *= scalar;
634 }
635 }
636 }
637
638 void add(const int localRow, const int localCol , const DofType value)
639 {
640#ifndef NDEBUG
641 check(localRow,localCol);
642#endif
643 getValue(localRow,localCol) += value;
644 }
645
646 void set(const int localRow, const int localCol , const DofType value)
647 {
648#ifndef NDEBUG
649 check(localRow,localCol);
650#endif
651 getValue(localRow,localCol) = value;
652 }
653
655 void unitRow(const int localRow)
656 {
657 const int row = (int) localRow / littleRows;
658 const int lRow = localRow%littleRows;
659
660 // clear row
661 doClearRow( row, lRow );
662
663 // set diagonal entry to 1
664 (*matrices_[row][row])[lRow][lRow] = 1;
665 }
666
668 void clear ()
669 {
670 const size_type nrows = matrices_.size();
671 for(size_type i=0; i<nrows; ++i)
672 {
673 const size_type ncols = matrices_[i].size();
674 for(size_type j=0; j<ncols; ++j)
675 {
676 (*matrices_[i][j]) = (DofType) 0;
677 }
678 }
679 }
680
682 void clearRow ( const int localRow )
683 {
684 const int row = (int) localRow / littleRows;
685 const int lRow = localRow%littleRows;
686
687 // clear the row
688 doClearRow( row, lRow );
689 }
690
692 void resort ()
693 {}
694
695 protected:
697 void doClearRow ( const int row, const int lRow )
698 {
699 // get number of columns
700 const auto col = this->columns();
701 for(auto localCol=0; localCol<col; ++localCol)
702 {
703 const int col = (int) localCol / littleCols;
704 const int lCol = localCol%littleCols;
705 (*matrices_[row][col])[lRow][lCol] = 0;
706 }
707 }
708
709 private:
710 // special mapper omitting block size
711 const RowMapperType& rowMapper_;
712 const ColMapperType& colMapper_;
713
714 // number of local matrices
715 int numRows_;
716 int numCols_;
717
718 // dynamic matrix with pointers to block matrices
719 std::vector< std::vector< LittleBlockType* > > matrices_;
720
721 // matrix to build
722 const MatrixObjectType& matrixObj_;
723 };
724
725
727 template <class DomainSpaceImp, class RangeSpaceImp, class DomainBlock, class RangeBlock >
728 class ISTLMatrixObject
729 {
730 public:
732 typedef DomainSpaceImp DomainSpaceType;
734 typedef RangeSpaceImp RangeSpaceType;
735
737 typedef ISTLMatrixObject< DomainSpaceImp, RangeSpaceImp, DomainBlock, RangeBlock > ThisType;
738 typedef ThisType PreconditionMatrixType;
739
740 typedef typename DomainSpaceType::GridType GridType;
741
742 typedef typename RangeSpaceType :: EntityType RangeEntityType ;
743 typedef typename DomainSpaceType :: EntityType DomainEntityType ;
744
745 enum { littleCols = DomainSpaceType :: localBlockSize };
746 enum { littleRows = RangeSpaceType :: localBlockSize };
747
748 typedef FieldMatrix<typename DomainSpaceType :: RangeFieldType, littleRows, littleCols> LittleBlockType;
749 typedef LittleBlockType block_type;
750 typedef LittleBlockType MatrixBlockType;
751
752 typedef ISTLBlockVectorDiscreteFunction< RangeSpaceType, RangeBlock > RowDiscreteFunctionType;
753 typedef ISTLBlockVectorDiscreteFunction< DomainSpaceType, DomainBlock > ColumnDiscreteFunctionType;
754
755 typedef typename RowDiscreteFunctionType :: DofStorageType RowBlockVectorType;
756 typedef typename ColumnDiscreteFunctionType :: DofStorageType ColumnBlockVectorType;
757
758 typedef typename RangeSpaceType :: BlockMapperType RowMapperType;
759 typedef typename DomainSpaceType :: BlockMapperType ColMapperType;
760
762 typedef ImprovedBCRSMatrix< LittleBlockType , ColumnDiscreteFunctionType , RowDiscreteFunctionType > MatrixType;
763 typedef ISTLParallelMatrixAdapterInterface< MatrixType > MatrixAdapterType;
764
766 typedef ISTLLocalMatrix<ThisType> ObjectType;
767 friend class ISTLLocalMatrix<ThisType>;
768
769 typedef ThisType LocalMatrixFactoryType;
770 typedef ObjectStack< LocalMatrixFactoryType > LocalMatrixStackType;
772 typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
773 typedef ColumnObject< ThisType > LocalColumnObjectType;
774
775 protected:
776 const DomainSpaceType& domainSpace_;
777 const RangeSpaceType& rangeSpace_;
778
779 // special row mapper
780 RowMapperType& rowMapper_;
781 // special col mapper
782 ColMapperType& colMapper_;
783
784 int size_;
785 int sequence_;
786
787 mutable std::unique_ptr< MatrixType > matrix_;
788
789 mutable LocalMatrixStackType localMatrixStack_;
790
791 mutable std::unique_ptr< MatrixAdapterType > matrixAdap_;
792 // overflow fraction for implicit build mode
793 const double overflowFraction_;
794 const bool threading_ ;
795
796 ISTLSolverParameter param_;
797 public:
798 ISTLMatrixObject(const ISTLMatrixObject&) = delete;
799
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() )
808 , size_(-1)
809 , sequence_(-1)
810 , localMatrixStack_( *this )
811 , overflowFraction_( param.overflowFraction() )
812 , threading_( param.threading() )
813 , param_( param )
814 {}
815
816 protected:
819 MatrixType& matrix() const
820 {
821 assert( matrix_ );
822 return *matrix_;
823 }
824
825 public:
826 bool threading() const { return threading_; }
827
828 void printTexInfo(std::ostream& out) const
829 {
830 out << "ISTL MatrixObj: ";
831 out << "\\\\ \n";
832 }
833
835 std::string preconditionName() const
836 {
837 return "";
838 }
839
840 void createMatrixAdapter ( const ISTLSolverParameter& param ) const
841 {
842 if( !matrixAdap_ )
843 {
844 matrixAdap_ = ISTLMatrixAdapterFactory< ThisType > :: matrixAdapter( *this, param );
845 }
846 assert( matrixAdap_ );
847 }
848
850 MatrixAdapterType& matrixAdapter( const ISTLSolverParameter& parameter ) const
851 {
852 const ISTLSolverParameter* param = dynamic_cast< const ISTLSolverParameter* > (&parameter);
853 if( param )
854 createMatrixAdapter( *param );
855 else
856 {
857 ISTLSolverParameter newparam( parameter );
858 createMatrixAdapter( newparam );
859 }
860
861 finalizeAssembly();
862 return *matrixAdap_;
863 }
864
866 MatrixAdapterType& matrixAdapter() const
867 {
868 return matrixAdapter( param_ );
869 }
870
871 public:
872 MatrixType& exportMatrix() const
873 {
874 finalizeAssembly();
875 return matrix();
876 }
877
878 bool implicitModeActive() const
879 {
880 // implicit build mode is only active when the
881 // build mode of the matrix is implicit and the
882 // matrix is currently being build
883
884 if( matrix().buildMode() == MatrixType::implicit && matrix().buildStage() == MatrixType::building )
885 return true;
886 else
887 return false;
888 }
889
890 void finalizeAssembly() const
891 {
892 // finalize build mode
893 const_cast< ThisType& > (*this).compress();
894 }
895
896 // compress matrix if not already done before and only in implicit build mode
897 void compress( )
898 {
899 if( implicitModeActive() )
900 {
901 matrix().compress();
902 }
903 }
904
906 void clear()
907 {
908 matrix().clear();
909 // clean matrix adapter and other helper classes
910 removeObj();
911 }
912
913 void unitRow( const size_t row )
914 {
915 matrix().unitRow( row );
916 }
917
918 protected:
919 template <class Container> // could be std::set or std::vector
920 void setUnitRowImpl( const Container &rows, const double diagonal )
921 {
922 for (const auto& r : rows)
923 {
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();
928#ifndef NDEBUG
929 bool set = false;
930#endif
931 for (auto col=row.begin(); col!=endcol; ++col)
932 {
933 for (auto& entry : (*col)[localRowIdx])
934 entry = 0;
935 if (col.index() == blockRow)
936 {
937 (*col)[localRowIdx][localRowIdx] = diagonal;
938#ifndef NDEBUG
939 set = true;
940#endif
941 }
942 }
943 assert(set);
944 }
945 }
946
947 public:
948 template <class Container>
949 void setUnitRows( const Container& unitRows, const Container& auxRows )
950 {
951 setUnitRowImpl( unitRows, 1.0 );
952 setUnitRowImpl( auxRows, 0.0 );
953 }
954
956 void reserve()
957 {
958 reserve( Stencil<DomainSpaceType,RangeSpaceType>(domainSpace_, rangeSpace_), true );
959 }
960
962 template <class Set>
963 void reserve (const std::vector< Set >& sparsityPattern )
964 {
965 reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ), false );
966 }
967
969 template <class Stencil>
970 void reserve(const Stencil &stencil, const bool proposeImplicit = true )
971 {
972 // if grid sequence number changed, rebuild matrix
973 if(sequence_ != domainSpace().sequence())
974 {
975 removeObj();
976
977 // do not use implicit build mode when multi threading is enabled
978 const bool implicit = proposeImplicit && MPIManager::numThreads() == 1;
979 if( implicit )
980 {
981 auto nnz = stencil.maxNonZerosEstimate();
982 if( nnz == 0 )
983 {
984 if (domainSpace_.begin() != domainSpace_.end() && rangeSpace_.begin() != rangeSpace_.end())
985 {
986 Stencil tmpStencil( stencil );
987 tmpStencil.fill( *(domainSpace_.begin()), *(rangeSpace_.begin()) );
988 nnz = tmpStencil.maxNonZerosEstimate();
989 }
990 }
991 matrix_.reset( new MatrixType( rowMapper_.size(), colMapper_.size(), nnz, overflowFraction_ ) );
992 }
993 else
994 {
995 matrix_.reset( new MatrixType( rowMapper_.size(), colMapper_.size() ) );
996 matrix().createEntries( stencil.globalStencil() );
997 }
998
999 sequence_ = domainSpace().sequence();
1000 }
1001 }
1002
1004 template <class HangingNodesType>
1005 void changeHangingNodes(const HangingNodesType& hangingNodes)
1006 {
1007 // create new matrix
1008 MatrixType* newMatrix = new MatrixType(rowMapper_.size(), colMapper_.size());
1009
1010 // setup with hanging rows
1011 newMatrix->setup( *matrix_ , hangingNodes );
1012
1013 // remove old matrix
1014 removeObj();
1015
1016 // store new matrix
1017 matrix_.reset( newMatrix );
1018 }
1019
1022 void extractDiagonal( ColumnDiscreteFunctionType& diag ) const
1023 {
1024 // extract diagonal entries
1025 matrix().extractDiagonal( diag.blockVector() );
1026 }
1027
1029 bool rightPrecondition() const
1030 {
1031 return false;
1032 }
1033
1035 void apply(const ColumnDiscreteFunctionType& arg, RowDiscreteFunctionType& dest) const
1036 {
1037 matrixAdapter().apply( arg.blockVector(), dest.blockVector() );
1038 }
1039
1041 template <class CDF, class RDF>
1042 void apply(const CDF& arg, RDF& dest) const
1043 {
1044 // this can happen when different storages are used for discrete
1045 // function and matrix/solver objects
1046 DUNE_THROW(NotImplemented,"ISTLMatrixObj::apply called for DiscreteFunctions not specified in the template list");
1047 }
1048
1050 void resort()
1051 {}
1052
1055 void createPreconditionMatrix()
1056 {}
1057
1059 void print(std::ostream & s) const
1060 {
1061 matrix().print(s);
1062 }
1063
1064 const DomainSpaceType& domainSpace() const
1065 {
1066 return domainSpace_;
1067 }
1068 const RangeSpaceType& rangeSpace() const
1069 {
1070 return rangeSpace_;
1071 }
1072
1073 const RowMapperType& rowMapper() const
1074 {
1075 return rowMapper_;
1076 }
1077 const ColMapperType& colMapper() const
1078 {
1079 return colMapper_;
1080 }
1081
1083 ObjectType* newObject() const
1084 {
1085 return new ObjectType(*this, domainSpace(), rangeSpace());
1086 }
1087
1092 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1093 LocalMatrixType localMatrix( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const
1094 {
1095 return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
1096 }
1097
1102 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
1103 LocalMatrixType localMatrix() const
1104 {
1105 return LocalMatrixType( localMatrixStack_ );
1106 }
1107
1108 LocalColumnObjectType localColumn( const DomainEntityType &domainEntity ) const
1109 {
1110 return LocalColumnObjectType ( *this, domainEntity );
1111 }
1112
1113 protected:
1114 template< class LocalBlock, class Operation >
1115 void applyToBlock ( const size_t row, const size_t col,
1116 const LocalBlock &localBlock,
1117 Operation& operation )
1118 {
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 ] );
1123 }
1124
1125 template< class LocalBlock >
1126 void setBlock ( const size_t row, const size_t col,
1127 const LocalBlock &localBlock )
1128 {
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 );
1132 }
1133
1134 template< class LocalBlock >
1135 void addBlock ( const size_t row, const size_t col,
1136 const LocalBlock &localBlock )
1137 {
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 );
1141 }
1142
1143 public:
1144 template< class LocalMatrix, class Operation >
1145 void applyToLocalMatrix ( const DomainEntityType &domainEntity,
1146 const RangeEntityType &rangeEntity,
1147 const LocalMatrix &localMat,
1148 Operation& operation )
1149 {
1150 typedef typename MatrixType::size_type Index;
1151 if( implicitModeActive() )
1152 {
1153 auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
1154 {
1155 return matrix().entry( index.first, index.second );
1156 };
1157 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
1158 {
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 ) );
1163 };
1164 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1165 }
1166 else
1167 {
1168 auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
1169 {
1170 return matrix()[ index.first][ index.second ];
1171 };
1172 auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
1173 {
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 ) );
1178 };
1179 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1180 }
1181 }
1182
1183 template< class LocalMatrix >
1184 void setLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
1185 {
1186 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1187 auto operation = [] ( RangeFieldType& a, const RangeFieldType& b )
1188 {
1189 a = b;
1190 };
1191 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1192 }
1193
1194 template< class LocalMatrix >
1195 void addLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
1196 {
1197 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1198 auto operation = [] ( RangeFieldType& a, const RangeFieldType& b )
1199 {
1200 a += b;
1201 };
1202 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1203 }
1204
1205 template< class LocalMatrix, class Scalar >
1206 void addScaledLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s )
1207 {
1208 typedef typename DomainSpaceType :: RangeFieldType RangeFieldType;
1209 auto operation = [ &s ] ( RangeFieldType& a, const RangeFieldType& b )
1210 {
1211 a += s * b;
1212 };
1213 applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
1214 }
1215
1216 template< class LocalMatrix >
1217 void getLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat ) const
1218 {
1219 // make sure that matrix is in compressed state if build mode is implicit
1220 finalizeAssembly();
1221
1222 typedef typename MatrixType::size_type Index;
1223 auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< Index, Index > &global )
1224 {
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] );
1229 };
1230
1231 rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
1232 }
1233
1234 protected:
1235 void preConErrorMsg(int preCon) const
1236 {
1237 exit(1);
1238 }
1239
1240 void removeObj ()
1241 {
1242 matrixAdap_.reset( nullptr );
1243 }
1244 };
1245
1246 } // namespace Fem
1247
1248} // namespace Dune
1249
1250#endif // #if HAVE_DUNE_ISTL
1251
1252#endif // #ifndef DUNE_FEM_ISTLMATRIXWRAPPER_HH
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.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Feb 17, 23:30, 2025)