DUNE-FEM (unstable)

spmatrix.hh
1#ifndef DUNE_FEM_SPMATRIX_HH
2#define DUNE_FEM_SPMATRIX_HH
3
4// system includes
5#include <algorithm>
6#include <array>
7#include <cstdlib>
8#include <iostream>
9#include <limits>
10#include <string>
11#include <utility>
12#include <vector>
13
14// local includes
15#include <dune/fem/function/adaptivefunction/adaptivefunction.hh>
16#include <dune/fem/misc/functor.hh>
17#include <dune/fem/operator/common/localmatrix.hh>
18#include <dune/fem/operator/common/localmatrixwrapper.hh>
19#include <dune/fem/io/file/asciiparser.hh>
20#include <dune/fem/io/parameter.hh>
21#include <dune/fem/operator/common/operator.hh>
22#include <dune/fem/operator/matrix/columnobject.hh>
23#include <dune/fem/space/mapper/nonblockmapper.hh>
24#include <dune/fem/storage/objectstack.hh>
25#include <dune/fem/operator/common/stencil.hh>
26#include <dune/fem/operator/matrix/functor.hh>
27#include <dune/fem/solver/parameter.hh>
28
29namespace Dune
30{
31
32 namespace Fem
33 {
34
36 template <class T, class IndexT = std::size_t,
37 class ValuesVector = std::vector< T >,
38 class IndicesVector = std::vector< IndexT > >
40 {
41 public:
43 typedef T field_type;
45 typedef IndexT size_type;
50
51 static const size_type defaultCol = std::numeric_limits<size_type>::max();
52 static const size_type zeroCol = std::numeric_limits<size_type>::max()-1;
53 static const int firstCol = 0;
54
55 SparseRowMatrix(const ThisType& ) = delete;
56
58 explicit SparseRowMatrix( const bool threading = true ) :
59 values_(0), columns_(0), rows_(0), dim_({{0,0}}), maxNzPerRow_(0), compressed_( false ), threading_( threading )
60 {}
61
64 SparseRowMatrix(const size_type rows, const size_type cols, const size_type nz, const bool threading = true ) :
65 values_(0), columns_(0), rows_(0), dim_({{0,0}}), maxNzPerRow_(0), compressed_( false ), threading_( threading )
66 {
67 reserve(rows,cols,nz);
68 }
69
71 void reserve(const size_type rows, const size_type cols, const size_type nz)
72 {
73 // if( (rows != dim_[0]) || (cols != dim_[1]) || (nz != maxNzPerRow_))
74 resize(rows,cols,nz);
75 clear();
76 }
77
79 template <class Stencil>
80 void fillPattern(const Stencil& stencil,
81 const size_type rowBlockSize,
82 const size_type colBlockSize )
83 {
84 const auto& sparsityPattern = stencil.globalStencil();
85 for( const auto& entry : sparsityPattern )
86 {
87 const auto& blockRow = entry.first;
88 const auto& blockColumnSet = entry.second;
89
90 // blocking of rows
91 const size_type nextRow = (blockRow + 1) * rowBlockSize;
92 for( size_type row = blockRow * rowBlockSize; row < nextRow; ++row )
93 {
94 size_type column = startRow( row );
95 for( const auto& blockColEntry : blockColumnSet )
96 {
97 size_type col = blockColEntry * colBlockSize;
98 for( size_type c = 0; c<colBlockSize; ++c, ++col, ++column )
99 {
100 assert( column < endRow( row ) );
101 columns_[ column ] = col ;
102 }
103 }
104 }
105 }
106 }
107
110 {
111 return dim_[0];
112 }
113
116 {
117 return dim_[1];
118 }
119
121 void set(const size_type row, const size_type col, const field_type val)
122 {
123 assert((col>=0) && (col < dim_[1]));
124 assert((row>=0) && (row < dim_[0]));
125
126 const size_type column = colIndex(row,col) ;
127 assert( column != defaultCol && column != zeroCol );
128
129 values_ [ column ] = val;
130 columns_[ column ] = col;
131 }
132
134 void add(const size_type row, const size_type col, const field_type val)
135 {
136 assert((col>=0) && (col < dim_[1]));
137 assert((row>=0) && (row < dim_[0]));
138
139 const size_type column = colIndex(row,col) ;
140 assert( column != defaultCol && column != zeroCol );
141
142 values_ [ column ] += val;
143 columns_[ column ] = col;
144 }
145
147 template<class ArgDFType, class DestDFType>
148 void apply(const ArgDFType& f, DestDFType& ret ) const
149 {
150 bool runParallel = threading_;
151
152 auto doApply = [this, &f, &ret]()
153 {
154 constexpr auto blockSize = ArgDFType::DiscreteFunctionSpaceType::localBlockSize;
155
156 // compute slice of rows to be worked on
157 const auto slice = sliceBeginEnd( MPIManager::thread(), MPIManager::numThreads(), std::true_type() );
158
159 // same as begin just with a row not necessarily zero
160 auto ret_it = ret.dofVector().find( slice.first );
161
162 const auto& fDofVec = f.dofVector();
163 for(size_type row = slice.first; row<slice.second; ++row)
164 {
165 const size_type endrow = endRow( row );
166 (*ret_it) = 0.0;
167 for(size_type col = startRow( row ); col<endrow; ++col)
168 {
169 const auto realCol = columns_[ col ];
170
171 if( ! compressed_ && ((realCol == defaultCol) || (realCol == zeroCol)) )
172 continue;
173
174 const auto blockNr = realCol / blockSize ;
175 const auto dofNr = realCol % blockSize ;
176 (*ret_it) += values_[ col ] * fDofVec[ blockNr ][ dofNr ];
177 }
178
179 ++ret_it;
180 }
181 };
182
183 if( runParallel )
184 {
185 try {
186 // execute in parallel
187 MPIManager :: run ( doApply );
188 }
189 catch ( const SingleThreadModeError& e )
190 {
191 runParallel = false;
192 }
193 }
194
195 // run serial if threading disabled or something else went wrong
196 if( ! runParallel )
197 {
198 doApply();
199 }
200 }
201
203 field_type get(const size_type row, const size_type col) const
204 {
205 assert((col>=0) && (col < dim_[1]));
206 assert((row>=0) && (row < dim_[0]));
207
208 const size_type endrow = endRow( row );
209 for( size_type i = startRow( row ); i<endrow; ++i )
210 {
211 if(columns_[ i ] == col)
212 return values_[ i ];
213 }
214 return 0;
215 }
216
218 void clear()
219 {
220 std::fill( values_.begin(), values_.end(), 0 );
221 for (auto &c : columns_) c = defaultCol;
222 }
223
225 void clearRow(const size_type row)
226 {
227 assert((row>=0) && (row < dim_[0]));
228
229 const size_type endrow = endRow( row );
230 for(size_type idx = startRow( row ); idx<endrow; ++idx )
231 {
232 values_[ idx ] = 0;
233 // if ( !compressed_ )
234 // columns_[idx] = zeroCol;
235 }
236 }
237
239 void scale(const size_type row, const size_type col, const field_type val)
240 {
241 assert((row>=0) && (row < rows()) );
242 assert((col>=0) && (col < cols()) );
243
244 const size_type column = colIndex(row,col) ;
245 assert( column != defaultCol && column != zeroCol );
246
247 // scale value
248 values_ [ column ] *= val;
249 }
250
254 {
255 return maxNzPerRow_;
256 }
257
261 {
262 return dim_[0] > 0 ? rows_[ dim_[0] ] : 0;
263 }
264
268 {
269 assert( row >= 0 && row < dim_[0] );
270 return endRow( row ) - startRow( row );
271 }
272
275 std::pair<const field_type, size_type> realValue(size_type index) const
276 {
277 return std::pair<const field_type, size_type>(values_[index], columns_[index]);
278 }
279
281 void print(std::ostream& s=std::cout, unsigned int offset=0) const
282 {
283 for(size_type row=0; row<dim_[0]; ++row)
284 {
285 const size_type endrow = endRow( row );
286 for( size_type pos = startRow( row ); pos<endrow; ++pos )
287 {
288 const auto rv(realValue(pos));
289 const auto column(rv.second);
290 const auto value(rv.first);
291 if((std::abs(value) > 1.e-15) && (column != defaultCol))
292 s << row+offset << " " << column+offset << " " << value << std::endl;
293 }
294 }
295 }
296
297 template <class SizeT, class NumericT >
298 void fillCSRStorage( std::vector< std::map<SizeT, NumericT> >& matrix ) const
299 {
300 matrix.resize( dim_[0] );
301
302 size_type thisCol = 0;
303 for(size_type row = 0; row<dim_[ 0 ]; ++row )
304 {
305 auto& matRow = matrix[ row ];
306 const size_type endrow = endRow( row );
307 for(size_type col = startRow( row ); col<endrow; ++col)
308 {
309 const size_type realCol = columns_[ col ];
310
311 if( ! compressed_ && (realCol == defaultCol || realCol == zeroCol) )
312 continue;
313
314 matRow[ realCol ] = values_[ thisCol ];
315 }
316 }
317 }
318
319 void compress ()
320 {
321 if( ! compressed_ && (dim_[0] != 0) && (dim_[1] != 0))
322 {
323 // determine first row nonZeros
324 size_type newpos = 0 ;
325 for( newpos = startRow( 0 ); newpos < endRow( 0 ); ++newpos )
326 {
327 if( columns_[ newpos ] == defaultCol )
328 break;
329 }
330
331 for( size_type row = 1; row<dim_[0]; ++row )
332 {
333 const size_type endrow = endRow( row );
334 size_type col = startRow( row );
335 // update new row start position
336 rows_[ row ] = newpos;
337 for(; col<endrow; ++col, ++newpos )
338 {
339 if( columns_[ col ] == defaultCol )
340 break ;
341
342 assert( newpos <= col );
343 values_ [ newpos ] = values_ [ col ];
344 columns_[ newpos ] = columns_[ col ];
345 }
346 }
347 rows_[ dim_[0] ] = newpos ;
348
349 // values_.resize( newpos );
350 // columns_.resize( newpos );
351 compressed_ = true ;
352 }
353 }
354
355 size_type startRow ( const size_type row ) const
356 {
357 return rows_[ row ];
358 }
359
360 size_type endRow ( const size_type row ) const
361 {
362 return rows_[ row+1 ];
363 }
364
365 std::tuple< ValuesVector&, IndicesVector&, IndicesVector& > exportCRS()
366 {
367 // only return internal data in compressed status
368 compress();
369 return std::tie(values_,columns_,rows_);
370 }
371
373 template<class DiagType, class ArgDFType, class DestDFType, class WType>
374 void forwardIterative(const DiagType& diagInv, const ArgDFType& b, const DestDFType& xold, DestDFType& xnew, const WType& w ) const
375 {
376 parallelIterative(diagInv, b, xold, xnew, w, std::true_type() );
377 }
378
380 template<class DiagType, class ArgDFType, class DestDFType, class WType>
381 void backwardIterative(const DiagType& diagInv, const ArgDFType& b, const DestDFType& xold, DestDFType& xnew, const WType& w ) const
382 {
383 parallelIterative(diagInv, b, xold, xnew, w, std::false_type() );
384 }
385
386 protected:
387 // return slice of rows for given thread in forward direction
388 std::pair< size_type, size_type > sliceBeginEnd(const size_type thread, const size_type numThreads, std::true_type ) const
389 {
390 const size_type sliceSize = this->rows() / numThreads;
391 const size_type sliceBegin = (thread * sliceSize) ;
392 const size_type sliceEnd = (thread == numThreads-1 ) ? this->rows(): (sliceBegin + sliceSize);
393 return std::make_pair( sliceBegin, sliceEnd );
394 }
395
396 // return slice of rows for given thread in backward direction
397 std::pair< size_type, size_type > sliceBeginEnd(const size_type thread, const size_type numThreads, std::false_type ) const
398 {
399 std::pair< size_type, size_type > beginEnd = sliceBeginEnd( thread, numThreads, std::true_type() );
400 return std::make_pair( beginEnd.second-1, beginEnd.first-1 );
401 }
402
403 template<class DiagType, class ArgDFType, class DestDFType, class WType, bool forward >
404 void parallelIterative(const DiagType& diagInv, const ArgDFType& b, const DestDFType& xold, DestDFType& xnew,
405 const WType& w, std::integral_constant<bool, forward> direction ) const
406 {
407 bool runParallel = threading_ && (&xold != &xnew) ; // only for Jacobi
408
409 auto doIterate = [this, &diagInv, &b, &xold, &xnew, &w, &direction] ()
410 {
411 // compute slice to be worked on depending on direction
412 const auto slice = sliceBeginEnd( MPIManager::thread(), MPIManager::numThreads(), direction );
413
414 doParallelIterative( diagInv.dofVector().find( slice.first ), // still O(1) operation just like for begin()
415 b.dofVector().find( slice.first ),
416 xold,
417 xnew.dofVector().find( slice.first ),
418 w,
419 slice.first, // row begin
420 slice.second, // row end
421 direction );
422 };
423
424 if( runParallel )
425 {
426 try {
427 // execute in parallel
428 MPIManager :: run ( doIterate );
429 }
430 catch ( const SingleThreadModeError& e )
431 {
432 runParallel = false;
433 }
434 }
435
436 // run serial if threading disabled or something else went wrong
437 if( ! runParallel )
438 {
439 doIterate();
440 }
441 }
442
444 template<class DiagIt, class ArgDFIt, class DestDFType, class DestDFIt,
445 class WType, bool forward>
446 void doParallelIterative(DiagIt diag, ArgDFIt bit, const DestDFType& xold, DestDFIt xit,
447 const WType& w,
448 size_type row, const size_type end,
449 std::integral_constant<bool, forward> ) const
450 {
451 constexpr auto blockSize = DestDFType::DiscreteFunctionSpaceType::localBlockSize;
452
453 const auto nextRow = [&diag, &xit, &bit](size_type &row)
454 {
455 if constexpr ( forward )
456 {
457 ++diag; ++xit; ++bit; ++row;
458 }
459 else
460 {
461 --diag; --xit; --bit; --row;
462 }
463 };
464
465 const auto& xOldVec = xold.dofVector();
466 for(; row != end; nextRow(row))
467 {
468 auto rhs = (*bit);
469
470 const size_type endrow = endRow( row );
471 for(size_type col = startRow( row ); col<endrow; ++col)
472 {
473 const auto realCol = columns_[ col ];
474
475 if( (realCol == row ) || (! compressed_ && ((realCol == defaultCol) || (realCol == zeroCol))) )
476 continue;
477
478 const auto blockNr = realCol / blockSize ;
479 const auto dofNr = realCol % blockSize ;
480
481 rhs -= values_[ col ] * xOldVec[ blockNr ][ dofNr ] ;
482 }
483
484 (*xit) = w * (rhs * (*diag));
485 }
486 }
489 {
490 constexpr auto colVal = defaultCol;
491 values_.resize( rows*nz , 0 );
492 columns_.resize( rows*nz , colVal );
493 rows_.resize( rows+1 , 0 );
494 rows_[ 0 ] = 0;
495 for( size_type i=1; i <= rows; ++i )
496 {
497 rows_[ i ] = rows_[ i-1 ] + nz ;
498 }
499 compressed_ = false;
500
501 dim_[0] = rows;
502 dim_[1] = cols;
503 maxNzPerRow_ = nz+firstCol;
504 }
505
508 {
509 assert((col>=0) && (col < dim_[1]));
510 assert((row>=0) && (row < dim_[0]));
511
512 const size_type endR = endRow( row );
513 size_type i = startRow( row );
514 // find local column or empty spot
515 for( ; i < endR; ++i )
516 {
517 if( columns_[ i ] == defaultCol || columns_[ i ] == zeroCol || columns_[ i ] == col )
518 {
519 return i;
520 }
521 }
522
523 assert(0);
524 DUNE_THROW( InvalidStateException, "Could not store entry in sparse matrix - no space available" );
525
526 // TODO: implement resize with 2*nz
527 std::abort();
528 return defaultCol;
529
530 /*
531 if(columns_[ i ] == col)
532 return i; // column already in matrix
533 else if( columns_[ i ] == defaultCol )
534 { // add this column at end of this row
535 ++nonZeros_[row];
536 return i;
537 }
538 else
539 {
540 std::abort();
541 // TODO re-implement
542 //
543 ++nonZeros_[row];
544 // must shift this row to add col at the position i
545 auto j = nz_-1; // last column
546 if (columns_[row*nz_+j] != defaultCol)
547 { // new space available - so resize
548 resize( rows(), cols(), (2 * nz_) );
549 j++;
550 }
551 for(;j>i;--j)
552 {
553 columns_[row*nz_+j] = columns_[row*nz_+j-1];
554 values_[row*nz_+j] = values_[row*nz_+j-1];
555 }
556 columns_[row*nz_+i] = col;
557 values_[row*nz_+i] = 0;
558 return i;
559 }
560 */
561 }
562
563 ValuesVector values_;
564 IndicesVector columns_;
565 IndicesVector rows_;
566
567 std::array<size_type,2> dim_;
568 size_type maxNzPerRow_;
569 bool compressed_;
570 const bool threading_;
571 };
572
573
574
576 template< class DomainSpace, class RangeSpace,
579 {
580 protected:
581 template< class MatrixObject >
582 struct LocalMatrixTraits;
583
584 template< class MatrixObject >
585 class LocalMatrix;
586
587 public:
588 typedef DomainSpace DomainSpaceType;
589 typedef RangeSpace RangeSpaceType;
590 typedef typename DomainSpaceType::EntityType DomainEntityType;
591 typedef typename RangeSpaceType::EntityType RangeEntityType;
592 typedef typename DomainSpaceType::EntityType ColumnEntityType;
593 typedef typename RangeSpaceType::EntityType RowEntityType;
594
595 typedef typename DomainSpaceType::BlockMapperType DomainBlockMapperType;
597 typedef typename RangeSpaceType::BlockMapperType RangeBlockMapperType;
599 typedef Matrix MatrixType;
600 typedef typename MatrixType::size_type size_type;
601 typedef typename MatrixType::field_type field_type;
602
604
605 static const size_type domainLocalBlockSize = DomainSpaceType::dimRange;
606 static const size_type rangeLocalBlockSize = RangeSpaceType::dimRange;
607
610
612
616 typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
617 typedef ColumnObject< ThisType > LocalColumnObjectType;
618
620 SparseRowMatrixObject( const DomainSpaceType &domainSpace,
621 const RangeSpaceType &rangeSpace,
622 const SolverParameter& param = SolverParameter() )
623 : domainSpace_( domainSpace ),
624 rangeSpace_( rangeSpace ),
625 domainMapper_( domainSpace_.blockMapper() ),
626 rangeMapper_( rangeSpace_.blockMapper() ),
627 sequence_( -1 ),
628 matrix_( param.threading() ),
629 localMatrixStack_( *this )
630 {}
631
633 const DomainSpaceType& domainSpace() const
634 {
635 return domainSpace_;
636 }
637
639 const RangeSpaceType& rangeSpace() const
640 {
641 return rangeSpace_;
642 }
643
644 protected:
647 {
648 return matrix_;
649 }
650
651 void finalizeAssembly() const { const_cast< ThisType& > (*this).compress(); }
652
653 public:
656 {
657 finalizeAssembly();
658 return matrix_;
659 }
660
661
664 {
665 return new ObjectType( *this, domainSpace_, rangeSpace_, domainMapper_, rangeMapper_ );
666 }
667
672 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
673 LocalMatrixType localMatrix( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const
674 {
675 return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
676 }
677
682 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
683 LocalMatrixType localMatrix() const
684 {
685 return LocalMatrixType( localMatrixStack_ );
686 }
687
689 LocalColumnObjectType localColumn( const DomainEntityType &domainEntity ) const
690 {
691 return LocalColumnObjectType( *this, domainEntity );
692 }
693
694 void unitRow( const size_type row )
695 {
696 for( unsigned int i=0, r = row * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r )
697 {
698 matrix_.clearRow( r );
699 matrix_.set( r, r, 1.0 );
700 }
701 }
702
703 template <class LocalBlock>
704 void addBlock( const size_type row, const size_type col, const LocalBlock& block )
705 {
706 std::array< size_type, rangeLocalBlockSize > rows;
707 std::array< size_type, domainLocalBlockSize > cols;
708 for( unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
709 {
710 rows[ i ] = r;
711 cols[ i ] = c;
712 }
713
714 for( unsigned int i=0; i<domainLocalBlockSize; ++i )
715 {
716 for( unsigned int j=0; j<domainLocalBlockSize; ++j )
717 {
718 matrix_.add( rows[ i ], cols[ j ], block[ i ][ j ]);
719 }
720 }
721 }
722
723 template <class LocalBlock>
724 void setBlock( const size_type row, const size_type col, const LocalBlock& block )
725 {
726 std::array< size_type, rangeLocalBlockSize > rows;
727 std::array< size_type, domainLocalBlockSize > cols;
728 for( unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
729 {
730 rows[ i ] = r;
731 cols[ i ] = c;
732 }
733
734 for( unsigned int i=0; i<domainLocalBlockSize; ++i )
735 {
736 for( unsigned int j=0; j<domainLocalBlockSize; ++j )
737 {
738 matrix_.set( rows[ i ], cols[ j ], block[ i ][ j ]);
739 }
740 }
741 }
742
743 template< class LocalMatrix >
744 void addLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
745 {
746 auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
747 {
748 matrix_.add( global.first, global.second, localMat.get( local.first, local.second ) );
749 };
750
751 rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
752 }
753
754 template< class LocalMatrix, class Scalar >
755 void addScaledLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s )
756 {
757 auto functor = [ &localMat, &s, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
758 {
759 matrix_.add( global.first, global.second, s * localMat.get( local.first, local.second ) );
760 };
761
762 rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
763 }
764
765 template< class LocalMatrix >
766 void setLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
767 {
768 auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
769 {
770 matrix_.set( global.first, global.second, localMat.get( local.first, local.second ) );
771 };
772
773 rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
774 }
775
776 template< class LocalMatrix >
777 void getLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat ) const
778 {
779 auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
780 {
781 localMat.set( local.first, local.second, matrix_.get( global.first, global.second ) );
782 };
783
784 rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
785 }
786
788 void clear()
789 {
790 matrix_.clear();
791 }
792
794 void compress() { matrix_.compress(); }
795
796 template <class Set>
797 void reserve (const std::vector< Set >& sparsityPattern )
798 {
799 reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ) );
800 }
801
803 template <class Stencil>
804 void reserve(const Stencil &stencil, bool verbose = false )
805 {
806 if( sequence_ != domainSpace_.sequence() )
807 {
808 // if empty grid do nothing (can appear in parallel runs)
809 if( (domainSpace_.begin() != domainSpace_.end()) && (rangeSpace_.begin() != rangeSpace_.end()) )
810 {
811 // output some info
812 if( verbose )
813 {
814 std::cout << "Reserve Matrix with (" << rangeSpace_.size() << "," << domainSpace_.size()<< ")" << std::endl;
815 std::cout << "Max number of base functions = (" << rangeMapper_.maxNumDofs() << ","
816 << domainMapper_.maxNumDofs() << ")" << std::endl;
817 }
818
819 // reserve matrix
820 const auto nonZeros = std::max( static_cast<size_type>(stencil.maxNonZerosEstimate()*DomainSpaceType::localBlockSize),
821 matrix_.maxNzPerRow() );
822 matrix_.reserve( rangeSpace_.size(), domainSpace_.size(), nonZeros );
823 matrix_.fillPattern( stencil, RangeSpaceType::localBlockSize, DomainSpaceType::localBlockSize );
824 }
825 sequence_ = domainSpace_.sequence();
826 }
827 }
828
830 template< class DomainFunction, class RangeFunction >
831 void apply( const DomainFunction &arg, RangeFunction &dest ) const
832 {
833 // do matrix vector multiplication
834 matrix_.apply( arg, dest );
835 // communicate data
836 dest.communicate();
837 }
838
841 template < class DiscreteFunctionType >
843 {
844 assert( matrix_.rows() == matrix_.cols() );
845 const auto dofEnd = diag.dend();
846 size_type row = 0;
847 for( auto dofIt = diag.dbegin(); dofIt != dofEnd; ++dofIt, ++row )
848 {
849 assert( row < matrix_.rows() );
850 (*dofIt) = matrix_.get( row, row );
851 }
852 }
853
854 template <class Container>
855 void setUnitRows( const Container& unitRows, const Container& auxRows )
856 {
857 for (const auto& r : unitRows )
858 {
859 matrix_.clearRow(r);
860 matrix_.set(r, r, 1.0);
861 }
862
863 for (const auto& r : auxRows )
864 {
865 matrix_.clearRow(r);
866 // not sure if this is really needed,
867 // but for consistency with previous code
868 // matrix_.set(r, r, 0.0);
869 }
870 }
871
873 void resort()
874 {
875 DUNE_THROW(NotImplemented,"SpMatrixObject::resort is not implemented");
876 // this method does not even exist on SpMatrix!!!
877 // matrix_.resort();
878 }
879
880 protected:
881 const DomainSpaceType &domainSpace_;
882 const RangeSpaceType &rangeSpace_;
883 DomainMapperType domainMapper_ ;
884 RangeMapperType rangeMapper_ ;
885 int sequence_;
886 mutable MatrixType matrix_;
887 mutable LocalMatrixStackType localMatrixStack_;
888 };
889
890
891
893 template< class DomainSpace, class RangeSpace, class Matrix >
894 template< class MatrixObject >
895 struct SparseRowMatrixObject< DomainSpace, RangeSpace, Matrix >::LocalMatrixTraits
896 {
897 typedef DomainSpace DomainSpaceType;
898 typedef RangeSpace RangeSpaceType;
899
901
902 typedef typename SparseRowMatrixObjectType::template LocalMatrix< MatrixObject > LocalMatrixType;
903
904 typedef typename RangeSpaceType::RangeFieldType RangeFieldType;
905 typedef RangeFieldType LittleBlockType;
906
909 };
910
911
912
914 template< class DomainSpace, class RangeSpace, class Matrix >
915 template< class MatrixObject >
916 class SparseRowMatrixObject< DomainSpace, RangeSpace, Matrix >::LocalMatrix
917 : public LocalMatrixDefault< LocalMatrixTraits< MatrixObject > >
918 {
919 public:
921 typedef MatrixObject MatrixObjectType;
922
925
926 private:
928
929 public:
931 typedef typename MatrixObjectType::MatrixType MatrixType;
932
934 typedef typename Traits::RangeFieldType RangeFieldType;
935
938
940 typedef typename Traits::LittleBlockType LittleBlockType;
941
946
947 typedef std::vector< typename RangeMapperType::SizeType > RowIndicesType;
948 typedef std::vector< typename DomainMapperType::SizeType > ColumnIndicesType;
949
951 LocalMatrix( const MatrixObjectType &matrixObject,
952 const DomainSpaceType &domainSpace,
953 const RangeSpaceType &rangeSpace,
954 const DomainMapperType& domainMapper,
955 const RangeMapperType& rangeMapper )
957 matrix_( matrixObject.matrix() ),
958 domainMapper_( domainMapper ),
959 rangeMapper_( rangeMapper )
960 {}
961
962 LocalMatrix( const LocalMatrix & ) = delete;
963
964 void init( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
965 {
966 // initialize base functions sets
967 BaseType::init( domainEntity, rangeEntity );
968 // rows are determined by the range space
969 rowIndices_.resize( rangeMapper_.numDofs( rangeEntity ) );
970 rangeMapper_.mapEach( rangeEntity, AssignFunctor< RowIndicesType >( rowIndices_ ) );
971 // columns are determined by the domain space
972 columnIndices_.resize( domainMapper_.numDofs( domainEntity ) );
973 domainMapper_.mapEach( domainEntity, AssignFunctor< ColumnIndicesType >( columnIndices_ ) );
974 }
975
977 size_type rows() const
978 {
979 return rowIndices_.size();
980 }
981
983 size_type columns() const
984 {
985 return columnIndices_.size();
986 }
987
989 void add(size_type localRow, size_type localCol, DofType value)
990 {
991 assert( value == value );
992 assert( (localRow >= 0) && (localRow < rows()) );
993 assert( (localCol >= 0) && (localCol < columns()) );
994 matrix_.add( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
995 }
996
998 DofType get(size_type localRow, size_type localCol) const
999 {
1000 assert( (localRow >= 0) && (localRow < rows()) );
1001 assert( (localCol >= 0) && (localCol < columns()) );
1002 return matrix_.get( rowIndices_[ localRow ], columnIndices_[ localCol ] );
1003 }
1004
1006 void set(size_type localRow, size_type localCol, DofType value)
1007 {
1008 assert( (localRow >= 0) && (localRow < rows()) );
1009 assert( (localCol >= 0) && (localCol < columns()) );
1010 matrix_.set( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
1011 }
1012
1014 void unitRow(size_type localRow)
1015 {
1016 assert( (localRow >= 0) && (localRow < rows()) );
1017 matrix_.unitRow( rowIndices_[ localRow ] );
1018 }
1019
1021 void clearRow( size_type localRow )
1022 {
1023 assert( (localRow >= 0) && (localRow < rows()) );
1024 matrix_.clearRow( rowIndices_[localRow]);
1025 }
1026
1028 void clearCol( size_type localCol )
1029 {
1030 assert( (localCol >= 0) && (localCol < columns()) );
1031 matrix_.clearCol( columnIndices_[localCol] );
1032 }
1033
1035 void clear()
1036 {
1037 const size_type nrows = rows();
1038 for(size_type i=0; i < nrows; ++i )
1039 matrix_.clearRow( rowIndices_[ i ] );
1040 }
1041
1043 void resort()
1044 {
1045 DUNE_THROW(NotImplemented,"SpMatrixObject::LocalMatrix::resort is not implemented");
1046 //const size_type nrows = rows();
1047 //for(size_type i=0; i < nrows; ++i )
1048 //matrix_.resortRow( rowIndices_[ i ] );
1049 }
1050
1052 void scale( const DofType& value )
1053 {
1054 const size_type nrows = rows();
1055 const size_type ncols = columns();
1056 for(size_type i=0; i < nrows; ++i )
1057 {
1058 for( size_type j=0; j < ncols; ++j )
1059 {
1060 scale(i, j, value );
1061 }
1062 }
1063 }
1064
1065 protected:
1067 void scale(size_type localRow, size_type localCol, DofType value)
1068 {
1069 assert( (localRow >= 0) && (localRow < rows()) );
1070 assert( (localCol >= 0) && (localCol < columns()) );
1071 matrix_.scale( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
1072 }
1073
1074 protected:
1075 MatrixType &matrix_;
1076 const DomainMapperType& domainMapper_;
1077 const RangeMapperType& rangeMapper_;
1078 RowIndicesType rowIndices_;
1079 ColumnIndicesType columnIndices_;
1080 };
1081
1082 } // namespace Fem
1083
1084} // namespace Dune
1085
1086#endif // #ifndef DUNE_FEM_SPMATRIX_HH
ConstDofIteratorType dbegin() const
Obtain the constant iterator pointing to the first dof.
Definition: discretefunction.hh:761
ConstDofIteratorType dend() const
Obtain the constant iterator pointing to the last dof.
Definition: discretefunction.hh:773
Default implementation for local matrix classes.
Definition: localmatrix.hh:287
Exception thrown when a code segment that is supposed to be only accessed in single thread mode is ac...
Definition: mpimanager.hh:43
LocalMatrix.
Definition: spmatrix.hh:918
void clearRow(size_type localRow)
set matrix row to zero
Definition: spmatrix.hh:1021
void add(size_type localRow, size_type localCol, DofType value)
add value to matrix entry
Definition: spmatrix.hh:989
MatrixObjectType::MatrixType MatrixType
type of matrix
Definition: spmatrix.hh:931
void clearCol(size_type localCol)
set matrix column to zero
Definition: spmatrix.hh:1028
void resort()
resort all global rows of matrix to have ascending numbering
Definition: spmatrix.hh:1043
size_type columns() const
return number of columns
Definition: spmatrix.hh:983
size_type rows() const
return number of rows
Definition: spmatrix.hh:977
Traits::RangeFieldType RangeFieldType
type of entries of little blocks
Definition: spmatrix.hh:934
DofType get(size_type localRow, size_type localCol) const
get matrix entry
Definition: spmatrix.hh:998
MatrixObject MatrixObjectType
type of matrix object
Definition: spmatrix.hh:921
Traits::DomainMapperType DomainMapperType
type of nonblocked domain mapper
Definition: spmatrix.hh:943
Traits::RangeMapperType RangeMapperType
type of nonblocked domain mapper
Definition: spmatrix.hh:945
LocalMatrixTraits< MatrixObjectType > Traits
type of the traits
Definition: spmatrix.hh:924
LocalMatrix(const MatrixObjectType &matrixObject, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const DomainMapperType &domainMapper, const RangeMapperType &rangeMapper)
constructor
Definition: spmatrix.hh:951
void scale(const DofType &value)
scale local matrix with a certain value
Definition: spmatrix.hh:1052
RangeFieldType DofType
type of the DoFs
Definition: spmatrix.hh:937
void scale(size_type localRow, size_type localCol, DofType value)
scale matrix entry with value
Definition: spmatrix.hh:1067
void clear()
clear all entries belonging to local matrix
Definition: spmatrix.hh:1035
Traits::LittleBlockType LittleBlockType
type of little blocks
Definition: spmatrix.hh:940
void set(size_type localRow, size_type localCol, DofType value)
set matrix entry to value
Definition: spmatrix.hh:1006
void unitRow(size_type localRow)
set matrix row to zero except diagonla entry
Definition: spmatrix.hh:1014
SparseRowMatrixObject.
Definition: spmatrix.hh:579
const DomainSpaceType & domainSpace() const
get domain space (i.e. space that builds the rows)
Definition: spmatrix.hh:633
LocalMatrixType localMatrix() const
Definition: spmatrix.hh:683
MatrixType & exportMatrix() const
get reference to storage object
Definition: spmatrix.hh:655
void apply(const DomainFunction &arg, RangeFunction &dest) const
apply matrix to discrete function
Definition: spmatrix.hh:831
SparseRowMatrixObject(const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const SolverParameter &param=SolverParameter())
construct matrix object
Definition: spmatrix.hh:620
LocalColumnObjectType localColumn(const DomainEntityType &domainEntity) const
get local column
Definition: spmatrix.hh:689
void extractDiagonal(DiscreteFunctionType &diag) const
Definition: spmatrix.hh:842
LocalMatrixType localMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity) const
Definition: spmatrix.hh:673
ObjectType * newObject() const
interface method from LocalMatrixFactory
Definition: spmatrix.hh:663
void reserve(const Stencil &stencil, bool verbose=false)
reserve memory
Definition: spmatrix.hh:804
void compress()
compress matrix to a real CRS format
Definition: spmatrix.hh:794
void resort()
resort row numbering in matrix to have ascending numbering
Definition: spmatrix.hh:873
void clear()
clear matrix
Definition: spmatrix.hh:788
const RangeSpaceType & rangeSpace() const
get range space (i.e. space that builds the columns)
Definition: spmatrix.hh:639
MatrixType & matrix() const
get reference to storage object, for internal use
Definition: spmatrix.hh:646
SparseRowMatrix.
Definition: spmatrix.hh:40
void print(std::ostream &s=std::cout, unsigned int offset=0) const
print matrix
Definition: spmatrix.hh:281
size_type colIndex(size_type row, size_type col)
returns local col index for given global (row,col)
Definition: spmatrix.hh:507
std::pair< const field_type, size_type > realValue(size_type index) const
Definition: spmatrix.hh:275
SparseRowMatrix(const bool threading=true)
construct matrix of zero size
Definition: spmatrix.hh:58
void clearRow(const size_type row)
set all entries in row to zero
Definition: spmatrix.hh:225
IndexT size_type
matrix index type
Definition: spmatrix.hh:45
void add(const size_type row, const size_type col, const field_type val)
add value to row,col entry
Definition: spmatrix.hh:134
void apply(const ArgDFType &f, DestDFType &ret) const
ret = A*f
Definition: spmatrix.hh:148
size_type numNonZeros(size_type row) const
Definition: spmatrix.hh:267
void reserve(const size_type rows, const size_type cols, const size_type nz)
reserve memory for given rows, columns and number of non zeros
Definition: spmatrix.hh:71
void clear()
set all matrix entries to zero
Definition: spmatrix.hh:218
size_type numNonZeros() const
Definition: spmatrix.hh:260
void forwardIterative(const DiagType &diagInv, const ArgDFType &b, const DestDFType &xold, DestDFType &xnew, const WType &w) const
Apply Jacobi/SOR method.
Definition: spmatrix.hh:374
field_type get(const size_type row, const size_type col) const
return value of entry (row,col)
Definition: spmatrix.hh:203
void doParallelIterative(DiagIt diag, ArgDFIt bit, const DestDFType &xold, DestDFIt xit, const WType &w, size_type row, const size_type end, std::integral_constant< bool, forward >) const
Apply Jacobi/SOR method.
Definition: spmatrix.hh:446
void set(const size_type row, const size_type col, const field_type val)
set entry to value (also setting 0 will result in an entry)
Definition: spmatrix.hh:121
void resize(size_type rows, size_type cols, size_type nz)
resize matrix
Definition: spmatrix.hh:488
size_type rows() const
return number of rows
Definition: spmatrix.hh:109
size_type cols() const
return number of columns
Definition: spmatrix.hh:115
void scale(const size_type row, const size_type col, const field_type val)
scale all entries in row with a given value
Definition: spmatrix.hh:239
size_type maxNzPerRow() const
Definition: spmatrix.hh:253
void backwardIterative(const DiagType &diagInv, const ArgDFType &b, const DestDFType &xold, DestDFType &xnew, const WType &w) const
Apply Jacobi/SOR method.
Definition: spmatrix.hh:381
SparseRowMatrix(const size_type rows, const size_type cols, const size_type nz, const bool threading=true)
Definition: spmatrix.hh:64
ThisType MatrixBaseType
Definition: spmatrix.hh:49
void fillPattern(const Stencil &stencil, const size_type rowBlockSize, const size_type colBlockSize)
reserve memory for given rows, columns and number of non zeros
Definition: spmatrix.hh:80
T field_type
matrix field type
Definition: spmatrix.hh:43
a simple wrapper class for sparsity patterns provide as vector< set< size_t > >
Definition: stencil.hh:234
default implementation for a general operator stencil
Definition: stencil.hh:35
const GlobalStencilType & globalStencil() const
Return the full stencil.
Definition: stencil.hh:116
int maxNonZerosEstimate() const
Return an upper bound for the maximum number of non-zero entries in all rows.
Definition: stencil.hh:123
forward declaration
Definition: discretefunction.hh:51
A dense n x m matrix.
Definition: fmatrix.hh:117
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
A generic dynamic dense matrix.
Definition: matrix.hh:561
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:577
typename Imp::BlockTraits< T >::field_type field_type
Export the type representing the underlying field.
Definition: matrix.hh:565
Default exception for dummy implementations.
Definition: exceptions.hh:263
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
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
LocalMatrixTraits.
Definition: spmatrix.hh:896
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)