dune-fem 2.8.0
Loading...
Searching...
No Matches
spmatrix.hh
Go to the documentation of this file.
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
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() :
59 values_(0), columns_(0), rows_(0), dim_({{0,0}}), maxNzPerRow_(0), compressed_( false )
60 {}
61
65 values_(0), columns_(0), rows_(0), dim_({{0,0}}), maxNzPerRow_(0), compressed_( false )
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 constexpr auto blockSize = ArgDFType::DiscreteFunctionSpaceType::localBlockSize;
151 auto ret_it = ret.dbegin();
152
153 for(size_type row = 0; row<dim_[0]; ++row)
154 {
155 const size_type endrow = endRow( row );
156 (*ret_it) = 0.0;
157 for(size_type col = startRow( row ); col<endrow; ++col)
158 {
159 const auto realCol = columns_[ col ];
160
161 if( ! compressed_ && ((realCol == defaultCol) || (realCol == zeroCol)) )
162 continue;
163
164 const auto blockNr = realCol / blockSize ;
165 const auto dofNr = realCol % blockSize ;
166 (*ret_it) += values_[ col ] * f.dofVector()[ blockNr ][ dofNr ];
167 }
168
169 ++ret_it;
170 }
171 }
172
174 field_type get(const size_type row, const size_type col) const
175 {
176 assert((col>=0) && (col < dim_[1]));
177 assert((row>=0) && (row < dim_[0]));
178
179 const size_type endrow = endRow( row );
180 for( size_type i = startRow( row ); i<endrow; ++i )
181 {
182 if(columns_[ i ] == col)
183 return values_[ i ];
184 }
185 return 0;
186 }
187
189 void clear()
190 {
191 std::fill( values_.begin(), values_.end(), 0 );
192 for (auto &c : columns_) c = defaultCol;
193 }
194
196 void clearRow(const size_type row)
197 {
198 assert((row>=0) && (row < dim_[0]));
199
200 const size_type endrow = endRow( row );
201 for(size_type idx = startRow( row ); idx<endrow; ++idx )
202 {
203 values_[ idx ] = 0;
204 // if ( !compressed_ )
205 // columns_[idx] = zeroCol;
206 }
207 }
208
210 void scale(const size_type row, const size_type col, const field_type val)
211 {
212 assert((row>=0) && (row < rows()) );
213 assert((col>=0) && (col < cols()) );
214
215 const size_type column = colIndex(row,col) ;
216 assert( column != defaultCol && column != zeroCol );
217
218 // scale value
219 values_ [ column ] *= val;
220 }
221
225 {
226 return maxNzPerRow_;
227 }
228
232 {
233 return dim_[0] > 0 ? rows_[ dim_[0] ] : 0;
234 }
235
239 {
240 assert( row >= 0 && row < dim_[0] );
241 return endRow( row ) - startRow( row );
242 }
243
246 std::pair<const field_type, size_type> realValue(size_type index) const
247 {
248 return std::pair<const field_type, size_type>(values_[index], columns_[index]);
249 }
250
252 void print(std::ostream& s=std::cout, unsigned int offset=0) const
253 {
254 for(size_type row=0; row<dim_[0]; ++row)
255 {
256 const size_type endrow = endRow( row );
257 for( size_type pos = startRow( row ); pos<endrow; ++pos )
258 {
259 const auto rv(realValue(pos));
260 const auto column(rv.second);
261 const auto value(rv.first);
262 if((std::abs(value) > 1.e-15) && (column != defaultCol))
263 s << row+offset << " " << column+offset << " " << value << std::endl;
264 }
265 }
266 }
267
268 template <class SizeT, class NumericT >
269 void fillCSRStorage( std::vector< std::map<SizeT, NumericT> >& matrix ) const
270 {
271 matrix.resize( dim_[0] );
272
273 size_type thisCol = 0;
274 for(size_type row = 0; row<dim_[ 0 ]; ++row )
275 {
276 auto& matRow = matrix[ row ];
277 const size_type endrow = endRow( row );
278 for(size_type col = startRow( row ); col<endrow; ++col)
279 {
280 const size_type realCol = columns_[ col ];
281
282 if( ! compressed_ && (realCol == defaultCol || realCol == zeroCol) )
283 continue;
284
285 matRow[ realCol ] = values_[ thisCol ];
286 }
287 }
288 }
289
290 void compress ()
291 {
292 if( ! compressed_ && (dim_[0] != 0) && (dim_[1] != 0))
293 {
294 // determine first row nonZeros
295 size_type newpos = 0 ;
296 for( newpos = startRow( 0 ); newpos < endRow( 0 ); ++newpos )
297 {
298 if( columns_[ newpos ] == defaultCol )
299 break;
300 }
301
302 for( size_type row = 1; row<dim_[0]; ++row )
303 {
304 const size_type endrow = endRow( row );
305 size_type col = startRow( row );
306 // update new row start position
307 rows_[ row ] = newpos;
308 for(; col<endrow; ++col, ++newpos )
309 {
310 if( columns_[ col ] == defaultCol )
311 break ;
312
313 assert( newpos <= col );
314 values_ [ newpos ] = values_ [ col ];
315 columns_[ newpos ] = columns_[ col ];
316 }
317 }
318 rows_[ dim_[0] ] = newpos ;
319
320 // values_.resize( newpos );
321 // columns_.resize( newpos );
322 compressed_ = true ;
323 }
324 }
325
326 size_type startRow ( const size_type row ) const
327 {
328 return rows_[ row ];
329 }
330
331 size_type endRow ( const size_type row ) const
332 {
333 return rows_[ row+1 ];
334 }
335
336 std::tuple< ValuesVector&, IndicesVector&, IndicesVector& > exportCRS()
337 {
338 // only return internal data in compressed status
339 compress();
340 return std::tie(values_,columns_,rows_);
341 }
342
343 protected:
346 {
347 constexpr auto colVal = defaultCol;
348 values_.resize( rows*nz , 0 );
349 columns_.resize( rows*nz , colVal );
350 rows_.resize( rows+1 , 0 );
351 rows_[ 0 ] = 0;
352 for( size_type i=1; i <= rows; ++i )
353 {
354 rows_[ i ] = rows_[ i-1 ] + nz ;
355 }
356 compressed_ = false;
357
358 dim_[0] = rows;
359 dim_[1] = cols;
361 }
362
365 {
366 assert((col>=0) && (col < dim_[1]));
367 assert((row>=0) && (row < dim_[0]));
368
369 const size_type endR = endRow( row );
370 size_type i = startRow( row );
371 // find local column or empty spot
372 for( ; i < endR; ++i )
373 {
374 if( columns_[ i ] == defaultCol || columns_[ i ] == zeroCol || columns_[ i ] == col )
375 {
376 return i;
377 }
378 }
379
380 assert(0);
381 DUNE_THROW( InvalidStateException, "Could not store entry in sparse matrix - no space available" );
382
383 // TODO: implement resize with 2*nz
384 std::abort();
385 return defaultCol;
386
387 /*
388 if(columns_[ i ] == col)
389 return i; // column already in matrix
390 else if( columns_[ i ] == defaultCol )
391 { // add this column at end of this row
392 ++nonZeros_[row];
393 return i;
394 }
395 else
396 {
397 std::abort();
398 // TODO re-implement
399 //
400 ++nonZeros_[row];
401 // must shift this row to add col at the position i
402 auto j = nz_-1; // last column
403 if (columns_[row*nz_+j] != defaultCol)
404 { // new space available - so resize
405 resize( rows(), cols(), (2 * nz_) );
406 j++;
407 }
408 for(;j>i;--j)
409 {
410 columns_[row*nz_+j] = columns_[row*nz_+j-1];
411 values_[row*nz_+j] = values_[row*nz_+j-1];
412 }
413 columns_[row*nz_+i] = col;
414 values_[row*nz_+i] = 0;
415 return i;
416 }
417 */
418 }
419
420 ValuesVector values_;
421 IndicesVector columns_;
422 IndicesVector rows_;
423
424 std::array<size_type,2> dim_;
427 };
428
429
430
432 template< class DomainSpace, class RangeSpace,
435 {
436 protected:
437 template< class MatrixObject >
438 struct LocalMatrixTraits;
439
440 template< class MatrixObject >
441 class LocalMatrix;
442
443 public:
444 typedef DomainSpace DomainSpaceType;
445 typedef RangeSpace RangeSpaceType;
446 typedef typename DomainSpaceType::EntityType DomainEntityType;
447 typedef typename RangeSpaceType::EntityType RangeEntityType;
448 typedef typename DomainSpaceType::EntityType ColumnEntityType;
449 typedef typename RangeSpaceType::EntityType RowEntityType;
450
451 typedef typename DomainSpaceType::BlockMapperType DomainBlockMapperType;
453 typedef typename RangeSpaceType::BlockMapperType RangeBlockMapperType;
455 typedef Matrix MatrixType;
456 typedef typename MatrixType::size_type size_type;
457 typedef typename MatrixType::field_type field_type;
458
460
461 static const size_type domainLocalBlockSize = DomainSpaceType::dimRange;
462 static const size_type rangeLocalBlockSize = RangeSpaceType::dimRange;
463
464 typedef Dune::FieldMatrix< field_type, rangeLocalBlockSize, domainLocalBlockSize > MatrixBlockType;
466
468
474
478 const SolverParameter& param = SolverParameter() )
481 domainMapper_( domainSpace_.blockMapper() ),
482 rangeMapper_( rangeSpace_.blockMapper() ),
483 sequence_( -1 ),
484 matrix_(),
486 param.preconditionMethod({SolverParameter::none,SolverParameter::jacobi})
488 localMatrixStack_( *this )
489 {}
490
493 {
494 return domainSpace_;
495 }
496
499 {
500 return rangeSpace_;
501 }
502
503 protected:
506 {
507 return matrix_;
508 }
509
510 void finalizeAssembly() const { const_cast< ThisType& > (*this).compress(); }
511
512 public:
515 {
517 return matrix_;
518 }
519
520
523 {
525 }
526
531 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
532 LocalMatrixType localMatrix( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const
533 {
534 return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
535 }
536
541 [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
543 {
545 }
546
549 {
550 return LocalColumnObjectType( *this, domainEntity );
551 }
552
553 void unitRow( const size_type row )
554 {
555 for( unsigned int i=0, r = row * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r )
556 {
557 matrix_.clearRow( r );
558 matrix_.set( r, r, 1.0 );
559 }
560 }
561
562 template <class LocalBlock>
563 void addBlock( const size_type row, const size_type col, const LocalBlock& block )
564 {
565 std::array< size_type, rangeLocalBlockSize > rows;
566 std::array< size_type, domainLocalBlockSize > cols;
567 for( unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
568 {
569 rows[ i ] = r;
570 cols[ i ] = c;
571 }
572
573 for( unsigned int i=0; i<domainLocalBlockSize; ++i )
574 {
575 for( unsigned int j=0; j<domainLocalBlockSize; ++j )
576 {
577 matrix_.add( rows[ i ], cols[ j ], block[ i ][ j ]);
578 }
579 }
580 }
581
582 template <class LocalBlock>
583 void setBlock( const size_type row, const size_type col, const LocalBlock& block )
584 {
585 std::array< size_type, rangeLocalBlockSize > rows;
586 std::array< size_type, domainLocalBlockSize > cols;
587 for( unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
588 {
589 rows[ i ] = r;
590 cols[ i ] = c;
591 }
592
593 for( unsigned int i=0; i<domainLocalBlockSize; ++i )
594 {
595 for( unsigned int j=0; j<domainLocalBlockSize; ++j )
596 {
597 matrix_.set( rows[ i ], cols[ j ], block[ i ][ j ]);
598 }
599 }
600 }
601
602 template< class LocalMatrix >
603 void addLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
604 {
605 auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
606 {
607 matrix_.add( global.first, global.second, localMat.get( local.first, local.second ) );
608 };
609
610 rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
611 }
612
613 template< class LocalMatrix, class Scalar >
614 void addScaledLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s )
615 {
616 auto functor = [ &localMat, &s, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
617 {
618 matrix_.add( global.first, global.second, s * localMat.get( local.first, local.second ) );
619 };
620
621 rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
622 }
623
624 template< class LocalMatrix >
625 void setLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
626 {
627 auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
628 {
629 matrix_.set( global.first, global.second, localMat.get( local.first, local.second ) );
630 };
631
632 rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
633 }
634
635 template< class LocalMatrix >
636 void getLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat ) const
637 {
638 auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
639 {
640 localMat.set( local.first, local.second, matrix_.get( global.first, global.second ) );
641 };
642
643 rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
644 }
645
647 void clear()
648 {
649 matrix_.clear();
650 }
651
653 void compress() { matrix_.compress(); }
654
655 template <class Set>
656 void reserve (const std::vector< Set >& sparsityPattern )
657 {
659 }
660
662 template <class Stencil>
663 void reserve(const Stencil &stencil, bool verbose = false )
664 {
665 if( sequence_ != domainSpace_.sequence() )
666 {
667 // if empty grid do nothing (can appear in parallel runs)
668 if( (domainSpace_.begin() != domainSpace_.end()) && (rangeSpace_.begin() != rangeSpace_.end()) )
669 {
670 // output some info
671 if( verbose )
672 {
673 std::cout << "Reserve Matrix with (" << rangeSpace_.size() << "," << domainSpace_.size()<< ")" << std::endl;
674 std::cout << "Max number of base functions = (" << rangeMapper_.maxNumDofs() << ","
675 << domainMapper_.maxNumDofs() << ")" << std::endl;
676 }
677
678 // reserve matrix
679 const auto nonZeros = std::max( static_cast<size_type>(stencil.maxNonZerosEstimate()*DomainSpaceType::localBlockSize),
680 matrix_.maxNzPerRow() );
681 matrix_.reserve( rangeSpace_.size(), domainSpace_.size(), nonZeros );
682 matrix_.fillPattern( stencil, RangeSpaceType::localBlockSize, DomainSpaceType::localBlockSize );
683 }
684 sequence_ = domainSpace_.sequence();
685 }
686 }
687
689 template< class DomainFunction, class RangeFunction >
690 void apply( const DomainFunction &arg, RangeFunction &dest ) const
691 {
692 // do matrix vector multiplication
693 matrix_.apply( arg, dest );
694 // communicate data
695 dest.communicate();
696 }
697
700 template < class DiscreteFunctionType >
701 void extractDiagonal( DiscreteFunctionType& diag ) const
702 {
703 assert( matrix_.rows() == matrix_.cols() );
704 const auto dofEnd = diag.dend();
705 size_type row = 0;
706 for( auto dofIt = diag.dbegin(); dofIt != dofEnd; ++dofIt, ++row )
707 {
708 assert( row < matrix_.rows() );
709 (*dofIt) = matrix_.get( row, row );
710 }
711 }
712
713 template <class Vector>
714 void setUnitRows( const Vector &rows )
715 {
716 const auto &auxiliaryDofs = domainSpace().auxiliaryDofs();
717 for (auto r : rows)
718 {
719 matrix_.clearRow(r);
720 matrix_.set(r,r,auxiliaryDofs.contains( r )? 0.0 : 1.0);
721 }
722 }
723
725 void resort()
726 {
727 DUNE_THROW(NotImplemented,"SpMatrixObject::resort is not implemented");
728 // this method does not even exist on SpMatrix!!!
729 // matrix_.resort();
730 }
731
732 protected:
741 };
742
743
744
746 template< class DomainSpace, class RangeSpace, class Matrix >
747 template< class MatrixObject >
748 struct SparseRowMatrixObject< DomainSpace, RangeSpace, Matrix >::LocalMatrixTraits
749 {
750 typedef DomainSpace DomainSpaceType;
751 typedef RangeSpace RangeSpaceType;
752
754
755 typedef typename SparseRowMatrixObjectType::template LocalMatrix< MatrixObject > LocalMatrixType;
756
757 typedef typename RangeSpaceType::RangeFieldType RangeFieldType;
759
762 };
763
764
765
767 template< class DomainSpace, class RangeSpace, class Matrix >
768 template< class MatrixObject >
769 class SparseRowMatrixObject< DomainSpace, RangeSpace, Matrix >::LocalMatrix
770 : public LocalMatrixDefault< LocalMatrixTraits< MatrixObject > >
771 {
772 public:
774 typedef MatrixObject MatrixObjectType;
775
778
779 private:
781
782 public:
784 typedef typename MatrixObjectType::MatrixType MatrixType;
785
788
791
794
799
800 typedef std::vector< typename RangeMapperType::SizeType > RowIndicesType;
801 typedef std::vector< typename DomainMapperType::SizeType > ColumnIndicesType;
802
804 LocalMatrix( const MatrixObjectType &matrixObject,
807 const DomainMapperType& domainMapper,
808 const RangeMapperType& rangeMapper )
810 matrix_( matrixObject.matrix() ),
811 domainMapper_( domainMapper ),
812 rangeMapper_( rangeMapper )
813 {}
814
815 LocalMatrix( const LocalMatrix & ) = delete;
816
817 void init( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
818 {
819 // initialize base functions sets
820 BaseType::init( domainEntity, rangeEntity );
821 // rows are determined by the range space
822 rowIndices_.resize( rangeMapper_.numDofs( rangeEntity ) );
823 rangeMapper_.mapEach( rangeEntity, AssignFunctor< RowIndicesType >( rowIndices_ ) );
824 // columns are determined by the domain space
825 columnIndices_.resize( domainMapper_.numDofs( domainEntity ) );
826 domainMapper_.mapEach( domainEntity, AssignFunctor< ColumnIndicesType >( columnIndices_ ) );
827 }
828
831 {
832 return rowIndices_.size();
833 }
834
837 {
838 return columnIndices_.size();
839 }
840
842 void add(size_type localRow, size_type localCol, DofType value)
843 {
844 assert( value == value );
845 assert( (localRow >= 0) && (localRow < rows()) );
846 assert( (localCol >= 0) && (localCol < columns()) );
847 matrix_.add( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
848 }
849
851 DofType get(size_type localRow, size_type localCol) const
852 {
853 assert( (localRow >= 0) && (localRow < rows()) );
854 assert( (localCol >= 0) && (localCol < columns()) );
855 return matrix_.get( rowIndices_[ localRow ], columnIndices_[ localCol ] );
856 }
857
859 void set(size_type localRow, size_type localCol, DofType value)
860 {
861 assert( (localRow >= 0) && (localRow < rows()) );
862 assert( (localCol >= 0) && (localCol < columns()) );
863 matrix_.set( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
864 }
865
867 void unitRow(size_type localRow)
868 {
869 assert( (localRow >= 0) && (localRow < rows()) );
870 matrix_.unitRow( rowIndices_[ localRow ] );
871 }
872
874 void clearRow( size_type localRow )
875 {
876 assert( (localRow >= 0) && (localRow < rows()) );
877 matrix_.clearRow( rowIndices_[localRow]);
878 }
879
881 void clearCol( size_type localCol )
882 {
883 assert( (localCol >= 0) && (localCol < columns()) );
884 matrix_.clearCol( columnIndices_[localCol] );
885 }
886
888 void clear()
889 {
890 const size_type nrows = rows();
891 for(size_type i=0; i < nrows; ++i )
892 matrix_.clearRow( rowIndices_[ i ] );
893 }
894
896 void resort()
897 {
898 DUNE_THROW(NotImplemented,"SpMatrixObject::LocalMatrix::resort is not implemented");
899 //const size_type nrows = rows();
900 //for(size_type i=0; i < nrows; ++i )
901 //matrix_.resortRow( rowIndices_[ i ] );
902 }
903
905 void scale( const DofType& value )
906 {
907 const size_type nrows = rows();
908 const size_type ncols = columns();
909 for(size_type i=0; i < nrows; ++i )
910 {
911 for( size_type j=0; j < ncols; ++j )
912 {
913 scale(i, j, value );
914 }
915 }
916 }
917
918 protected:
920 void scale(size_type localRow, size_type localCol, DofType value)
921 {
922 assert( (localRow >= 0) && (localRow < rows()) );
923 assert( (localCol >= 0) && (localCol < columns()) );
924 matrix_.scale( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
925 }
926
927 protected:
933 };
934
935 } // namespace Fem
936
937} // namespace Dune
938
939#endif // #ifndef DUNE_FEM_SPMATRIX_HH
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
Definition: bindguard.hh:11
PairFunctor< Mapper, Entity, Functor > makePairFunctor(const Mapper &mapper, const Entity &entity, Functor functor)
Definition: operator/matrix/functor.hh:65
Definition: misc/functor.hh:31
Default implementation for local matrix classes.
Definition: localmatrix.hh:285
BaseType::RangeEntityType RangeEntityType
Definition: localmatrix.hh:299
BaseType::DomainEntityType DomainEntityType
Definition: localmatrix.hh:298
Definition: localmatrixwrapper.hh:48
default implementation for a general operator stencil
Definition: stencil.hh:35
const GlobalStencilType & globalStencil() const
Return the full stencil.
Definition: stencil.hh:115
int maxNonZerosEstimate() const
Return an upper bound for the maximum number of non-zero entries in all rows.
Definition: stencil.hh:122
a simple wrapper class for sparsity patterns provide as vector< set< size_t > >
Definition: stencil.hh:233
Definition: columnobject.hh:12
SparseRowMatrix.
Definition: spmatrix.hh:40
void print(std::ostream &s=std::cout, unsigned int offset=0) const
print matrix
Definition: spmatrix.hh:252
size_type colIndex(size_type row, size_type col)
returns local col index for given global (row,col)
Definition: spmatrix.hh:364
std::array< size_type, 2 > dim_
Definition: spmatrix.hh:424
std::pair< const field_type, size_type > realValue(size_type index) const
Definition: spmatrix.hh:246
void fillCSRStorage(std::vector< std::map< SizeT, NumericT > > &matrix) const
Definition: spmatrix.hh:269
size_type endRow(const size_type row) const
Definition: spmatrix.hh:331
void clearRow(const size_type row)
set all entries in row to zero
Definition: spmatrix.hh:196
IndexT size_type
matrix index type
Definition: spmatrix.hh:45
size_type maxNzPerRow_
Definition: spmatrix.hh:425
void add(const size_type row, const size_type col, const field_type val)
add value to row,col entry
Definition: spmatrix.hh:134
std::tuple< ValuesVector &, IndicesVector &, IndicesVector & > exportCRS()
Definition: spmatrix.hh:336
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:238
IndicesVector columns_
Definition: spmatrix.hh:421
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:189
static const size_type zeroCol
Definition: spmatrix.hh:52
SparseRowMatrix< field_type, size_type, ValuesVector, IndicesVector > ThisType
Definition: spmatrix.hh:46
static const size_type defaultCol
Definition: spmatrix.hh:51
size_type numNonZeros() const
Definition: spmatrix.hh:231
field_type get(const size_type row, const size_type col) const
return value of entry (row,col)
Definition: spmatrix.hh:174
static const int firstCol
Definition: spmatrix.hh:53
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
size_type startRow(const size_type row) const
Definition: spmatrix.hh:326
void resize(size_type rows, size_type cols, size_type nz)
resize matrix
Definition: spmatrix.hh:345
size_type rows() const
return number of rows
Definition: spmatrix.hh:109
IndicesVector rows_
Definition: spmatrix.hh:422
size_type cols() const
return number of columns
Definition: spmatrix.hh:115
SparseRowMatrix(const size_type rows, const size_type cols, const size_type nz)
Definition: spmatrix.hh:64
ValuesVector values_
Definition: spmatrix.hh:420
bool compressed_
Definition: spmatrix.hh:426
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:210
void compress()
Definition: spmatrix.hh:290
SparseRowMatrix(const ThisType &)=delete
size_type maxNzPerRow() const
Definition: spmatrix.hh:224
SparseRowMatrix()
construct matrix of zero size
Definition: spmatrix.hh:58
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
SparseRowMatrixObject.
Definition: spmatrix.hh:435
NonBlockMapper< DomainBlockMapperType, DomainSpaceType::localBlockSize > DomainMapperType
Definition: spmatrix.hh:452
RangeSpaceType::EntityType RangeEntityType
Definition: spmatrix.hh:447
const DomainSpaceType & domainSpace() const
get domain space (i.e. space that builds the rows)
Definition: spmatrix.hh:492
const DomainSpaceType & domainSpace_
Definition: spmatrix.hh:733
DomainMapperType domainMapper_
Definition: spmatrix.hh:735
LocalMatrixType localMatrix() const
Definition: spmatrix.hh:542
DomainSpaceType::EntityType DomainEntityType
Definition: spmatrix.hh:446
RangeSpaceType::BlockMapperType RangeBlockMapperType
Definition: spmatrix.hh:453
Matrix MatrixType
Definition: spmatrix.hh:455
MatrixType::size_type size_type
Definition: spmatrix.hh:456
LocalMatrixStackType localMatrixStack_
Definition: spmatrix.hh:740
MatrixType & exportMatrix() const
get reference to storage object
Definition: spmatrix.hh:514
ThisType LocalMatrixFactoryType
Definition: spmatrix.hh:470
LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType
Definition: spmatrix.hh:472
void setUnitRows(const Vector &rows)
Definition: spmatrix.hh:714
void addScaledLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s)
Definition: spmatrix.hh:614
SparseRowMatrixObject< DomainSpaceType, RangeSpaceType, MatrixType > ThisType
Definition: spmatrix.hh:459
RangeSpace RangeSpaceType
Definition: spmatrix.hh:445
void setBlock(const size_type row, const size_type col, const LocalBlock &block)
Definition: spmatrix.hh:583
RangeSpaceType::EntityType RowEntityType
Definition: spmatrix.hh:449
static const size_type rangeLocalBlockSize
Definition: spmatrix.hh:462
void reserve(const std::vector< Set > &sparsityPattern)
Definition: spmatrix.hh:656
LocalMatrix< ThisType > ObjectType
Definition: spmatrix.hh:469
ColumnObject< ThisType > LocalColumnObjectType
Definition: spmatrix.hh:473
void apply(const DomainFunction &arg, RangeFunction &dest) const
apply matrix to discrete function
Definition: spmatrix.hh:690
Fem::ObjectStack< LocalMatrixFactoryType > LocalMatrixStackType
Definition: spmatrix.hh:471
int sequence_
Definition: spmatrix.hh:737
SparseRowMatrixObject(const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const SolverParameter &param=SolverParameter())
construct matrix object
Definition: spmatrix.hh:476
static const size_type domainLocalBlockSize
Definition: spmatrix.hh:461
LocalColumnObjectType localColumn(const DomainEntityType &domainEntity) const
get local column
Definition: spmatrix.hh:548
void extractDiagonal(DiscreteFunctionType &diag) const
Definition: spmatrix.hh:701
void setLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat)
Definition: spmatrix.hh:625
MatrixType::field_type field_type
Definition: spmatrix.hh:457
LocalMatrixType localMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity) const
Definition: spmatrix.hh:532
MatrixType PreconditionMatrixType
Definition: spmatrix.hh:467
DomainSpace DomainSpaceType
Definition: spmatrix.hh:444
DomainSpaceType::BlockMapperType DomainBlockMapperType
Definition: spmatrix.hh:451
DomainSpaceType::EntityType ColumnEntityType
Definition: spmatrix.hh:448
ObjectType * newObject() const
interface method from LocalMatrixFactory
Definition: spmatrix.hh:522
NonBlockMapper< RangeBlockMapperType, RangeSpaceType::localBlockSize > RangeMapperType
Definition: spmatrix.hh:454
Dune::FieldMatrix< field_type, rangeLocalBlockSize, domainLocalBlockSize > MatrixBlockType
Definition: spmatrix.hh:464
void reserve(const Stencil &stencil, bool verbose=false)
reserve memory
Definition: spmatrix.hh:663
void addLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat)
Definition: spmatrix.hh:603
MatrixType matrix_
Definition: spmatrix.hh:738
void compress()
compress matrix to a real CRS format
Definition: spmatrix.hh:653
void resort()
resort row numbering in matrix to have ascending numbering
Definition: spmatrix.hh:725
void getLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat) const
Definition: spmatrix.hh:636
void unitRow(const size_type row)
Definition: spmatrix.hh:553
void addBlock(const size_type row, const size_type col, const LocalBlock &block)
Definition: spmatrix.hh:563
RangeMapperType rangeMapper_
Definition: spmatrix.hh:736
MatrixBlockType block_type
Definition: spmatrix.hh:465
void clear()
clear matrix
Definition: spmatrix.hh:647
bool preconditioning_
Definition: spmatrix.hh:739
const RangeSpaceType & rangeSpace_
Definition: spmatrix.hh:734
const RangeSpaceType & rangeSpace() const
get range space (i.e. space that builds the columns)
Definition: spmatrix.hh:498
MatrixType & matrix() const
get reference to storage object, for internal use
Definition: spmatrix.hh:505
void finalizeAssembly() const
Definition: spmatrix.hh:510
LocalMatrixTraits.
Definition: spmatrix.hh:749
SparseRowMatrixObject< DomainSpaceType, RangeSpaceType, Matrix > SparseRowMatrixObjectType
Definition: spmatrix.hh:753
RangeFieldType LittleBlockType
Definition: spmatrix.hh:758
RangeSpaceType::RangeFieldType RangeFieldType
Definition: spmatrix.hh:757
DomainSpace DomainSpaceType
Definition: spmatrix.hh:750
RangeSpace RangeSpaceType
Definition: spmatrix.hh:751
SparseRowMatrixObjectType::template LocalMatrix< MatrixObject > LocalMatrixType
Definition: spmatrix.hh:755
SparseRowMatrixObjectType::RangeMapperType RangeMapperType
Definition: spmatrix.hh:761
SparseRowMatrixObjectType::DomainMapperType DomainMapperType
Definition: spmatrix.hh:760
LocalMatrix.
Definition: spmatrix.hh:771
MatrixType & matrix_
Definition: spmatrix.hh:928
LocalMatrix(const LocalMatrix &)=delete
void clearRow(size_type localRow)
set matrix row to zero
Definition: spmatrix.hh:874
const DomainMapperType & domainMapper_
Definition: spmatrix.hh:929
void add(size_type localRow, size_type localCol, DofType value)
add value to matrix entry
Definition: spmatrix.hh:842
void init(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity)
Definition: spmatrix.hh:817
MatrixObjectType::MatrixType MatrixType
type of matrix
Definition: spmatrix.hh:784
void clearCol(size_type localCol)
set matrix column to zero
Definition: spmatrix.hh:881
void resort()
resort all global rows of matrix to have ascending numbering
Definition: spmatrix.hh:896
size_type columns() const
return number of columns
Definition: spmatrix.hh:836
ColumnIndicesType columnIndices_
Definition: spmatrix.hh:932
size_type rows() const
return number of rows
Definition: spmatrix.hh:830
RowIndicesType rowIndices_
Definition: spmatrix.hh:931
std::vector< typename DomainMapperType::SizeType > ColumnIndicesType
Definition: spmatrix.hh:801
Traits::RangeFieldType RangeFieldType
type of entries of little blocks
Definition: spmatrix.hh:787
DofType get(size_type localRow, size_type localCol) const
get matrix entry
Definition: spmatrix.hh:851
MatrixObject MatrixObjectType
type of matrix object
Definition: spmatrix.hh:774
std::vector< typename RangeMapperType::SizeType > RowIndicesType
Definition: spmatrix.hh:800
Traits::DomainMapperType DomainMapperType
type of nonblocked domain mapper
Definition: spmatrix.hh:796
Traits::RangeMapperType RangeMapperType
type of nonblocked domain mapper
Definition: spmatrix.hh:798
LocalMatrixTraits< MatrixObjectType > Traits
type of the traits
Definition: spmatrix.hh:777
LocalMatrix(const MatrixObjectType &matrixObject, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const DomainMapperType &domainMapper, const RangeMapperType &rangeMapper)
constructor
Definition: spmatrix.hh:804
const RangeMapperType & rangeMapper_
Definition: spmatrix.hh:930
void scale(const DofType &value)
scale local matrix with a certain value
Definition: spmatrix.hh:905
RangeFieldType DofType
type of the DoFs
Definition: spmatrix.hh:790
void scale(size_type localRow, size_type localCol, DofType value)
scale matrix entry with value
Definition: spmatrix.hh:920
void clear()
clear all entries belonging to local matrix
Definition: spmatrix.hh:888
Traits::LittleBlockType LittleBlockType
type of little blocks
Definition: spmatrix.hh:793
void set(size_type localRow, size_type localCol, DofType value)
set matrix entry to value
Definition: spmatrix.hh:859
void unitRow(size_type localRow)
set matrix row to zero except diagonla entry
Definition: spmatrix.hh:867
Definition: solver/parameter.hh:15
static const int none
Definition: solver/parameter.hh:40
static const int jacobi
Definition: solver/parameter.hh:45