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