dune-fem  2.8-git
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
16 #include <dune/fem/misc/functor.hh>
20 #include <dune/fem/io/parameter.hh>
28 
29 namespace 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 
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  size_type rows () const
80  {
81  return dim_[0];
82  }
83 
85  size_type cols () const
86  {
87  return dim_[1];
88  }
89 
91  void set(const size_type row, const size_type col, const field_type val)
92  {
93  assert((col>=0) && (col < dim_[1]));
94  assert((row>=0) && (row < dim_[0]));
95 
96  const size_type column = colIndex(row,col) ;
97  assert( column != defaultCol && column != zeroCol );
98 
99  values_ [ column ] = val;
100  columns_[ column ] = col;
101  }
102 
104  void add(const size_type row, const size_type col, const field_type val)
105  {
106  assert((col>=0) && (col < dim_[1]));
107  assert((row>=0) && (row < dim_[0]));
108 
109  const size_type column = colIndex(row,col) ;
110  assert( column != defaultCol && column != zeroCol );
111 
112  values_ [ column ] += val;
113  columns_[ column ] = col;
114  }
115 
117  template<class ArgDFType, class DestDFType>
118  void apply(const ArgDFType& f, DestDFType& ret ) const
119  {
120  constexpr auto blockSize = ArgDFType::DiscreteFunctionSpaceType::localBlockSize;
121  auto ret_it = ret.dbegin();
122 
123  for(size_type row = 0; row<dim_[0]; ++row)
124  {
125  const size_type endrow = endRow( row );
126  (*ret_it) = 0.0;
127  for(size_type col = startRow( row ); col<endrow; ++col)
128  {
129  const auto realCol = columns_[ col ];
130 
131  if( ! compressed_ && ((realCol == defaultCol) || (realCol == zeroCol)) )
132  continue;
133 
134  const auto blockNr = realCol / blockSize ;
135  const auto dofNr = realCol % blockSize ;
136  (*ret_it) += values_[ col ] * f.dofVector()[ blockNr ][ dofNr ];
137  }
138 
139  ++ret_it;
140  }
141  }
142 
144  field_type get(const size_type row, const size_type col) const
145  {
146  assert((col>=0) && (col < dim_[1]));
147  assert((row>=0) && (row < dim_[0]));
148 
149  const size_type endrow = endRow( row );
150  for( size_type i = startRow( row ); i<endrow; ++i )
151  {
152  if(columns_[ i ] == col)
153  return values_[ i ];
154  }
155  return 0;
156  }
157 
159  void clear()
160  {
161  std::fill( values_.begin(), values_.end(), 0 );
162  for (auto &c : columns_) c = defaultCol;
163  }
164 
166  void clearRow(const size_type row)
167  {
168  assert((row>=0) && (row < dim_[0]));
169 
170  const size_type endrow = endRow( row );
171  for(size_type idx = startRow( row ); idx<endrow; ++idx )
172  {
173  values_[ idx ] = 0;
174  // if ( !compressed_ )
175  // columns_[idx] = zeroCol;
176  }
177  }
178 
180  void scale(const size_type row, const size_type col, const field_type val)
181  {
182  assert((row>=0) && (row < rows()) );
183  assert((col>=0) && (col < cols()) );
184 
185  const size_type column = colIndex(row,col) ;
186  assert( column != defaultCol && column != zeroCol );
187 
188  // scale value
189  values_ [ column ] *= val;
190  }
191 
195  {
196  return maxNzPerRow_;
197  }
198 
202  {
203  return dim_[0] > 0 ? rows_[ dim_[0] ] : 0;
204  }
205 
209  {
210  assert( row >= 0 && row < dim_[0] );
211  return endRow( row ) - startRow( row );
212  }
213 
216  std::pair<const field_type, size_type> realValue(size_type index) const
217  {
218  return std::pair<const field_type, size_type>(values_[index], columns_[index]);
219  }
220 
222  void print(std::ostream& s=std::cout, unsigned int offset=0) const
223  {
224  for(size_type row=0; row<dim_[0]; ++row)
225  {
226  const size_type endrow = endRow( row );
227  for( size_type pos = startRow( row ); pos<endrow; ++pos )
228  {
229  const auto rv(realValue(pos));
230  const auto column(rv.second);
231  const auto value(rv.first);
232  if((std::abs(value) > 1.e-15) && (column != defaultCol))
233  s << row+offset << " " << column+offset << " " << value << std::endl;
234  }
235  }
236  }
237 
238  template <class SizeT, class NumericT >
239  void fillCSRStorage( std::vector< std::map<SizeT, NumericT> >& matrix ) const
240  {
241  matrix.resize( dim_[0] );
242 
243  size_type thisCol = 0;
244  for(size_type row = 0; row<dim_[ 0 ]; ++row )
245  {
246  auto& matRow = matrix[ row ];
247  const size_type endrow = endRow( row );
248  for(size_type col = startRow( row ); col<endrow; ++col)
249  {
250  const size_type realCol = columns_[ col ];
251 
252  if( ! compressed_ && (realCol == defaultCol || realCol == zeroCol) )
253  continue;
254 
255  matRow[ realCol ] = values_[ thisCol ];
256  }
257  }
258  }
259 
260  void compress ()
261  {
262  if( ! compressed_ && (dim_[0] != 0) && (dim_[1] != 0))
263  {
264  // determine first row nonZeros
265  size_type newpos = 0 ;
266  for( newpos = startRow( 0 ); newpos < endRow( 0 ); ++newpos )
267  {
268  if( columns_[ newpos ] == defaultCol )
269  break;
270  }
271 
272  for( size_type row = 1; row<dim_[0]; ++row )
273  {
274  const size_type endrow = endRow( row );
275  size_type col = startRow( row );
276  // update new row start position
277  rows_[ row ] = newpos;
278  for(; col<endrow; ++col, ++newpos )
279  {
280  if( columns_[ col ] == defaultCol )
281  break ;
282 
283  assert( newpos <= col );
284  values_ [ newpos ] = values_ [ col ];
285  columns_[ newpos ] = columns_[ col ];
286  }
287  }
288  rows_[ dim_[0] ] = newpos ;
289 
290  // values_.resize( newpos );
291  // columns_.resize( newpos );
292  compressed_ = true ;
293  }
294  }
295 
296  size_type startRow ( const size_type row ) const
297  {
298  return rows_[ row ];
299  }
300 
301  size_type endRow ( const size_type row ) const
302  {
303  return rows_[ row+1 ];
304  }
305 
306  std::tuple< ValuesVector&, IndicesVector&, IndicesVector& > exportCRS()
307  {
308  // only return internal data in compressed status
309  compress();
310  return std::tie(values_,columns_,rows_);
311  }
312 
313  protected:
316  {
317  constexpr auto colVal = defaultCol;
318  values_.resize( rows*nz , 0 );
319  columns_.resize( rows*nz , colVal );
320  rows_.resize( rows+1 , 0 );
321  rows_[ 0 ] = 0;
322  for( size_type i=1; i <= rows; ++i )
323  {
324  rows_[ i ] = rows_[ i-1 ] + nz ;
325  }
326  compressed_ = false;
327 
328  dim_[0] = rows;
329  dim_[1] = cols;
330  maxNzPerRow_ = nz+firstCol;
331  }
332 
335  {
336  assert((col>=0) && (col < dim_[1]));
337  assert((row>=0) && (row < dim_[0]));
338 
339  const size_type endR = endRow( row );
340  size_type i = startRow( row );
341  // find local column or empty spot
342  for( ; i < endR; ++i )
343  {
344  if( columns_[ i ] == defaultCol || columns_[ i ] == zeroCol || columns_[ i ] == col )
345  {
346  return i;
347  }
348  }
349 
350  assert(0);
351  DUNE_THROW( InvalidStateException, "Could not store entry in sparse matrix - no space available" );
352 
353  // TODO: implement resize with 2*nz
354  std::abort();
355  return defaultCol;
356 
357  /*
358  if(columns_[ i ] == col)
359  return i; // column already in matrix
360  else if( columns_[ i ] == defaultCol )
361  { // add this column at end of this row
362  ++nonZeros_[row];
363  return i;
364  }
365  else
366  {
367  std::abort();
368  // TODO re-implement
369  //
370  ++nonZeros_[row];
371  // must shift this row to add col at the position i
372  auto j = nz_-1; // last column
373  if (columns_[row*nz_+j] != defaultCol)
374  { // new space available - so resize
375  resize( rows(), cols(), (2 * nz_) );
376  j++;
377  }
378  for(;j>i;--j)
379  {
380  columns_[row*nz_+j] = columns_[row*nz_+j-1];
381  values_[row*nz_+j] = values_[row*nz_+j-1];
382  }
383  columns_[row*nz_+i] = col;
384  values_[row*nz_+i] = 0;
385  return i;
386  }
387  */
388  }
389 
390  ValuesVector values_;
391  IndicesVector columns_;
392  IndicesVector rows_;
393 
394  std::array<size_type,2> dim_;
397  };
398 
399 
400 
402  template< class DomainSpace, class RangeSpace,
405  {
406  protected:
407  template< class MatrixObject >
408  struct LocalMatrixTraits;
409 
410  template< class MatrixObject >
411  class LocalMatrix;
412 
413  public:
414  typedef DomainSpace DomainSpaceType;
415  typedef RangeSpace RangeSpaceType;
416  typedef typename DomainSpaceType::EntityType DomainEntityType;
417  typedef typename RangeSpaceType::EntityType RangeEntityType;
418  typedef typename DomainSpaceType::EntityType ColumnEntityType;
419  typedef typename RangeSpaceType::EntityType RowEntityType;
420 
421  typedef typename DomainSpaceType::BlockMapperType DomainBlockMapperType;
423  typedef typename RangeSpaceType::BlockMapperType RangeBlockMapperType;
425  typedef Matrix MatrixType;
426  typedef typename MatrixType::size_type size_type;
427  typedef typename MatrixType::field_type field_type;
428 
430 
431  static const size_type domainLocalBlockSize = DomainSpaceType::dimRange;
432  static const size_type rangeLocalBlockSize = RangeSpaceType::dimRange;
433 
434  typedef Dune::FieldMatrix< field_type, rangeLocalBlockSize, domainLocalBlockSize > MatrixBlockType;
436 
438 
444 
447  const RangeSpaceType &rangeSpace,
448  const SolverParameter& param = SolverParameter() )
451  domainMapper_( domainSpace_.blockMapper() ),
452  rangeMapper_( rangeSpace_.blockMapper() ),
453  sequence_( -1 ),
454  matrix_(),
456  param.preconditionMethod({SolverParameter::none,SolverParameter::jacobi})
458  localMatrixStack_( *this )
459  {}
460 
463  {
464  return domainSpace_;
465  }
466 
468  const RangeSpaceType& rangeSpace() const
469  {
470  return rangeSpace_;
471  }
472 
473  protected:
476  {
477  return matrix_;
478  }
479 
480  void finalizeAssembly() const { const_cast< ThisType& > (*this).compress(); }
481 
482  public:
485  {
487  return matrix_;
488  }
489 
490 
493  {
495  }
496 
501  [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
502  LocalMatrixType localMatrix( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const
503  {
504  return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
505  }
506 
511  [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
513  {
515  }
516 
519  {
520  return LocalColumnObjectType( *this, domainEntity );
521  }
522 
523  void unitRow( const size_type row )
524  {
525  for( unsigned int i=0, r = row * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r )
526  {
527  matrix_.clearRow( r );
528  matrix_.set( r, r, 1.0 );
529  }
530  }
531 
532  template <class LocalBlock>
533  void addBlock( const size_type row, const size_type col, const LocalBlock& block )
534  {
535  std::array< size_type, rangeLocalBlockSize > rows;
536  std::array< size_type, domainLocalBlockSize > cols;
537  for( unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
538  {
539  rows[ i ] = r;
540  cols[ i ] = c;
541  }
542 
543  for( unsigned int i=0; i<domainLocalBlockSize; ++i )
544  {
545  for( unsigned int j=0; j<domainLocalBlockSize; ++j )
546  {
547  matrix_.add( rows[ i ], cols[ j ], block[ i ][ j ]);
548  }
549  }
550  }
551 
552  template <class LocalBlock>
553  void setBlock( const size_type row, const size_type col, const LocalBlock& block )
554  {
555  std::array< size_type, rangeLocalBlockSize > rows;
556  std::array< size_type, domainLocalBlockSize > cols;
557  for( unsigned int i=0, r = row * domainLocalBlockSize, c = col * domainLocalBlockSize; i<domainLocalBlockSize; ++i, ++r, ++c )
558  {
559  rows[ i ] = r;
560  cols[ i ] = c;
561  }
562 
563  for( unsigned int i=0; i<domainLocalBlockSize; ++i )
564  {
565  for( unsigned int j=0; j<domainLocalBlockSize; ++j )
566  {
567  matrix_.set( rows[ i ], cols[ j ], block[ i ][ j ]);
568  }
569  }
570  }
571 
572  template< class LocalMatrix >
573  void addLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
574  {
575  auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
576  {
577  matrix_.add( global.first, global.second, localMat.get( local.first, local.second ) );
578  };
579 
580  rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
581  }
582 
583  template< class LocalMatrix, class Scalar >
584  void addScaledLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s )
585  {
586  auto functor = [ &localMat, &s, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
587  {
588  matrix_.add( global.first, global.second, s * localMat.get( local.first, local.second ) );
589  };
590 
591  rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
592  }
593 
594  template< class LocalMatrix >
595  void setLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
596  {
597  auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
598  {
599  matrix_.set( global.first, global.second, localMat.get( local.first, local.second ) );
600  };
601 
602  rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
603  }
604 
605  template< class LocalMatrix >
606  void getLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat ) const
607  {
608  auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< size_type, size_type >& global )
609  {
610  localMat.set( local.first, local.second, matrix_.get( global.first, global.second ) );
611  };
612 
613  rangeMapper_.mapEach( rangeEntity, makePairFunctor( domainMapper_, domainEntity, functor ) );
614  }
615 
617  void clear()
618  {
619  matrix_.clear();
620  }
621 
623  void compress() { matrix_.compress(); }
624 
625  template <class Set>
626  void reserve (const std::vector< Set >& sparsityPattern )
627  {
629  }
630 
632  template <class Stencil>
633  void reserve(const Stencil &stencil, bool verbose = false )
634  {
635  if( sequence_ != domainSpace_.sequence() )
636  {
637  // if empty grid do nothing (can appear in parallel runs)
638  if( (domainSpace_.begin() != domainSpace_.end()) && (rangeSpace_.begin() != rangeSpace_.end()) )
639  {
640  // output some info
641  if( verbose )
642  {
643  std::cout << "Reserve Matrix with (" << rangeSpace_.size() << "," << domainSpace_.size()<< ")" << std::endl;
644  std::cout << "Max number of base functions = (" << rangeMapper_.maxNumDofs() << ","
645  << domainMapper_.maxNumDofs() << ")" << std::endl;
646  }
647  // reserve matrix
648  const auto nonZeros = std::max( static_cast<size_type>(stencil.maxNonZerosEstimate()*DomainSpaceType::localBlockSize),
649  matrix_.maxNzPerRow() );
650  matrix_.reserve( rangeSpace_.size(), domainSpace_.size(), nonZeros );
651  }
652  sequence_ = domainSpace_.sequence();
653  }
654  }
655 
657  template< class DomainFunction, class RangeFunction >
658  void apply( const DomainFunction &arg, RangeFunction &dest ) const
659  {
660  // do matrix vector multiplication
661  matrix_.apply( arg, dest );
662  // communicate data
663  dest.communicate();
664  }
665 
668  template < class DiscreteFunctionType >
669  void extractDiagonal( DiscreteFunctionType& diag ) const
670  {
671  assert( matrix_.rows() == matrix_.cols() );
672  const auto dofEnd = diag.dend();
673  size_type row = 0;
674  for( auto dofIt = diag.dbegin(); dofIt != dofEnd; ++dofIt, ++row )
675  {
676  assert( row < matrix_.rows() );
677  (*dofIt) = matrix_.get( row, row );
678  }
679  }
680 
681  template <class Vector>
682  void setUnitRows( const Vector &rows )
683  {
684  const auto &auxiliaryDofs = domainSpace().auxiliaryDofs();
685  for (auto r : rows)
686  {
687  matrix_.clearRow(r);
688  matrix_.set(r,r,auxiliaryDofs.contains( r )? 0.0 : 1.0);
689  }
690  }
691 
693  void resort()
694  {
695  DUNE_THROW(NotImplemented,"SpMatrixObject::resort is not implemented");
696  // this method does not even exist on SpMatrix!!!
697  // matrix_.resort();
698  }
699 
700  protected:
709  };
710 
711 
712 
714  template< class DomainSpace, class RangeSpace, class Matrix >
715  template< class MatrixObject >
716  struct SparseRowMatrixObject< DomainSpace, RangeSpace, Matrix >::LocalMatrixTraits
717  {
718  typedef DomainSpace DomainSpaceType;
719  typedef RangeSpace RangeSpaceType;
720 
722 
723  typedef typename SparseRowMatrixObjectType::template LocalMatrix< MatrixObject > LocalMatrixType;
724 
725  typedef typename RangeSpaceType::RangeFieldType RangeFieldType;
727 
730  };
731 
732 
733 
735  template< class DomainSpace, class RangeSpace, class Matrix >
736  template< class MatrixObject >
737  class SparseRowMatrixObject< DomainSpace, RangeSpace, Matrix >::LocalMatrix
738  : public LocalMatrixDefault< LocalMatrixTraits< MatrixObject > >
739  {
740  public:
742  typedef MatrixObject MatrixObjectType;
743 
746 
747  private:
749 
750  public:
752  typedef typename MatrixObjectType::MatrixType MatrixType;
753 
755  typedef typename Traits::RangeFieldType RangeFieldType;
756 
759 
761  typedef typename Traits::LittleBlockType LittleBlockType;
762 
764  typedef typename Traits::DomainMapperType DomainMapperType;
766  typedef typename Traits::RangeMapperType RangeMapperType;
767 
768  typedef std::vector< typename RangeMapperType::SizeType > RowIndicesType;
769  typedef std::vector< typename DomainMapperType::SizeType > ColumnIndicesType;
770 
772  LocalMatrix( const MatrixObjectType &matrixObject,
774  const RangeSpaceType &rangeSpace,
775  const DomainMapperType& domainMapper,
776  const RangeMapperType& rangeMapper )
778  matrix_( matrixObject.matrix() ),
779  domainMapper_( domainMapper ),
780  rangeMapper_( rangeMapper )
781  {}
782 
783  LocalMatrix( const LocalMatrix & ) = delete;
784 
785  void init( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
786  {
787  // initialize base functions sets
788  BaseType::init( domainEntity, rangeEntity );
789  // rows are determined by the range space
790  rowIndices_.resize( rangeMapper_.numDofs( rangeEntity ) );
791  rangeMapper_.mapEach( rangeEntity, AssignFunctor< RowIndicesType >( rowIndices_ ) );
792  // columns are determined by the domain space
793  columnIndices_.resize( domainMapper_.numDofs( domainEntity ) );
794  domainMapper_.mapEach( domainEntity, AssignFunctor< ColumnIndicesType >( columnIndices_ ) );
795  }
796 
798  size_type rows() const
799  {
800  return rowIndices_.size();
801  }
802 
805  {
806  return columnIndices_.size();
807  }
808 
810  void add(size_type localRow, size_type localCol, DofType value)
811  {
812  assert( value == value );
813  assert( (localRow >= 0) && (localRow < rows()) );
814  assert( (localCol >= 0) && (localCol < columns()) );
815  matrix_.add( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
816  }
817 
819  DofType get(size_type localRow, size_type localCol) const
820  {
821  assert( (localRow >= 0) && (localRow < rows()) );
822  assert( (localCol >= 0) && (localCol < columns()) );
823  return matrix_.get( rowIndices_[ localRow ], columnIndices_[ localCol ] );
824  }
825 
827  void set(size_type localRow, size_type localCol, DofType value)
828  {
829  assert( (localRow >= 0) && (localRow < rows()) );
830  assert( (localCol >= 0) && (localCol < columns()) );
831  matrix_.set( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
832  }
833 
835  void unitRow(size_type localRow)
836  {
837  assert( (localRow >= 0) && (localRow < rows()) );
838  matrix_.unitRow( rowIndices_[ localRow ] );
839  }
840 
842  void clearRow( size_type localRow )
843  {
844  assert( (localRow >= 0) && (localRow < rows()) );
845  matrix_.clearRow( rowIndices_[localRow]);
846  }
847 
849  void clearCol( size_type localCol )
850  {
851  assert( (localCol >= 0) && (localCol < columns()) );
852  matrix_.clearCol( columnIndices_[localCol] );
853  }
854 
856  void clear()
857  {
858  const size_type nrows = rows();
859  for(size_type i=0; i < nrows; ++i )
860  matrix_.clearRow( rowIndices_[ i ] );
861  }
862 
864  void resort()
865  {
866  DUNE_THROW(NotImplemented,"SpMatrixObject::LocalMatrix::resort is not implemented");
867  //const size_type nrows = rows();
868  //for(size_type i=0; i < nrows; ++i )
869  //matrix_.resortRow( rowIndices_[ i ] );
870  }
871 
873  void scale( const DofType& value )
874  {
875  const size_type nrows = rows();
876  const size_type ncols = columns();
877  for(size_type i=0; i < nrows; ++i )
878  {
879  for( size_type j=0; j < ncols; ++j )
880  {
881  scale(i, j, value );
882  }
883  }
884  }
885 
886  protected:
888  void scale(size_type localRow, size_type localCol, DofType value)
889  {
890  assert( (localRow >= 0) && (localRow < rows()) );
891  assert( (localCol >= 0) && (localCol < columns()) );
892  matrix_.scale( rowIndices_[ localRow ], columnIndices_[ localCol ], value );
893  }
894 
895  protected:
901  };
902 
903  } // namespace Fem
904 
905 } // namespace Dune
906 
907 #endif // #ifndef DUNE_FEM_SPMATRIX_HH
Definition: bindguard.hh:11
Double abs(const Double &a)
Definition: double.hh:871
PairFunctor< Mapper, Entity, Functor > makePairFunctor(const Mapper &mapper, const Entity &entity, Functor functor)
Definition: operator/matrix/functor.hh:65
static constexpr T max(T a)
Definition: utility.hh:77
Definition: bartonnackmaninterface.hh:17
Definition: misc/functor.hh:31
Traits ::RangeFieldType RangeFieldType
type of range field
Definition: localmatrix.hh:49
Traits ::RangeSpaceType RangeSpaceType
type of range discrete function space
Definition: localmatrix.hh:55
Traits ::DomainSpaceType DomainSpaceType
type of domain discrete function space
Definition: localmatrix.hh:52
RangeSpaceType::EntityType RangeEntityType
Definition: localmatrix.hh:66
DomainSpaceType::EntityType DomainEntityType
Definition: localmatrix.hh:65
Default implementation for local matrix classes.
Definition: localmatrix.hh:285
Definition: localmatrixwrapper.hh:48
default implementation for a general operator stencil
Definition: stencil.hh:35
int maxNonZerosEstimate() const
Return an upper bound for the maximum number of non-zero entries in all row.
Definition: stencil.hh:104
a simple wrapper class for sparsity patterns provide as vector< set< size_t > >
Definition: stencil.hh:191
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:222
size_type colIndex(size_type row, size_type col)
returns local col index for given global (row,col)
Definition: spmatrix.hh:334
std::array< size_type, 2 > dim_
Definition: spmatrix.hh:394
void fillCSRStorage(std::vector< std::map< SizeT, NumericT > > &matrix) const
Definition: spmatrix.hh:239
size_type endRow(const size_type row) const
Definition: spmatrix.hh:301
void clearRow(const size_type row)
set all entries in row to zero
Definition: spmatrix.hh:166
IndexT size_type
matrix index type
Definition: spmatrix.hh:45
size_type maxNzPerRow_
Definition: spmatrix.hh:395
void add(const size_type row, const size_type col, const field_type val)
add value to row,col entry
Definition: spmatrix.hh:104
void apply(const ArgDFType &f, DestDFType &ret) const
ret = A*f
Definition: spmatrix.hh:118
size_type numNonZeros(size_type row) const
Definition: spmatrix.hh:208
IndicesVector columns_
Definition: spmatrix.hh:391
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:159
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:201
field_type get(const size_type row, const size_type col) const
return value of entry (row,col)
Definition: spmatrix.hh:144
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:91
size_type startRow(const size_type row) const
Definition: spmatrix.hh:296
void resize(size_type rows, size_type cols, size_type nz)
resize matrix
Definition: spmatrix.hh:315
size_type rows() const
return number of rows
Definition: spmatrix.hh:79
IndicesVector rows_
Definition: spmatrix.hh:392
size_type cols() const
return number of columns
Definition: spmatrix.hh:85
SparseRowMatrix(const size_type rows, const size_type cols, const size_type nz)
Definition: spmatrix.hh:64
ValuesVector values_
Definition: spmatrix.hh:390
bool compressed_
Definition: spmatrix.hh:396
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:180
void compress()
Definition: spmatrix.hh:260
SparseRowMatrix(const ThisType &)=delete
std::tuple< ValuesVector &, IndicesVector &, IndicesVector & > exportCRS()
Definition: spmatrix.hh:306
size_type maxNzPerRow() const
Definition: spmatrix.hh:194
SparseRowMatrix()
construct matrix of zero size
Definition: spmatrix.hh:58
ThisType MatrixBaseType
Definition: spmatrix.hh:49
T field_type
matrix field type
Definition: spmatrix.hh:43
std::pair< const field_type, size_type > realValue(size_type index) const
Definition: spmatrix.hh:216
SparseRowMatrixObject.
Definition: spmatrix.hh:405
NonBlockMapper< DomainBlockMapperType, DomainSpaceType::localBlockSize > DomainMapperType
Definition: spmatrix.hh:422
RangeSpaceType::EntityType RangeEntityType
Definition: spmatrix.hh:417
const DomainSpaceType & domainSpace_
Definition: spmatrix.hh:701
DomainMapperType domainMapper_
Definition: spmatrix.hh:703
LocalMatrixType localMatrix() const
Definition: spmatrix.hh:512
DomainSpaceType::EntityType DomainEntityType
Definition: spmatrix.hh:416
ObjectType * newObject() const
interface method from LocalMatrixFactory
Definition: spmatrix.hh:492
RangeSpaceType::BlockMapperType RangeBlockMapperType
Definition: spmatrix.hh:423
Matrix MatrixType
Definition: spmatrix.hh:425
MatrixType::size_type size_type
Definition: spmatrix.hh:426
LocalMatrixStackType localMatrixStack_
Definition: spmatrix.hh:708
ThisType LocalMatrixFactoryType
Definition: spmatrix.hh:440
LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType
Definition: spmatrix.hh:442
void setUnitRows(const Vector &rows)
Definition: spmatrix.hh:682
void addScaledLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s)
Definition: spmatrix.hh:584
SparseRowMatrixObject< DomainSpaceType, RangeSpaceType, MatrixType > ThisType
Definition: spmatrix.hh:429
RangeSpace RangeSpaceType
Definition: spmatrix.hh:415
MatrixType & exportMatrix() const
get reference to storage object
Definition: spmatrix.hh:484
void setBlock(const size_type row, const size_type col, const LocalBlock &block)
Definition: spmatrix.hh:553
RangeSpaceType::EntityType RowEntityType
Definition: spmatrix.hh:419
static const size_type rangeLocalBlockSize
Definition: spmatrix.hh:432
void reserve(const std::vector< Set > &sparsityPattern)
Definition: spmatrix.hh:626
LocalMatrix< ThisType > ObjectType
Definition: spmatrix.hh:439
ColumnObject< ThisType > LocalColumnObjectType
Definition: spmatrix.hh:443
void apply(const DomainFunction &arg, RangeFunction &dest) const
apply matrix to discrete function
Definition: spmatrix.hh:658
Fem::ObjectStack< LocalMatrixFactoryType > LocalMatrixStackType
Definition: spmatrix.hh:441
int sequence_
Definition: spmatrix.hh:705
SparseRowMatrixObject(const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const SolverParameter &param=SolverParameter())
construct matrix object
Definition: spmatrix.hh:446
static const size_type domainLocalBlockSize
Definition: spmatrix.hh:431
LocalColumnObjectType localColumn(const DomainEntityType &domainEntity) const
get local column
Definition: spmatrix.hh:518
void extractDiagonal(DiscreteFunctionType &diag) const
Definition: spmatrix.hh:669
void setLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat)
Definition: spmatrix.hh:595
MatrixType::field_type field_type
Definition: spmatrix.hh:427
LocalMatrixType localMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity) const
Definition: spmatrix.hh:502
MatrixType PreconditionMatrixType
Definition: spmatrix.hh:437
DomainSpace DomainSpaceType
Definition: spmatrix.hh:411
DomainSpaceType::BlockMapperType DomainBlockMapperType
Definition: spmatrix.hh:421
DomainSpaceType::EntityType ColumnEntityType
Definition: spmatrix.hh:418
NonBlockMapper< RangeBlockMapperType, RangeSpaceType::localBlockSize > RangeMapperType
Definition: spmatrix.hh:424
Dune::FieldMatrix< field_type, rangeLocalBlockSize, domainLocalBlockSize > MatrixBlockType
Definition: spmatrix.hh:434
void reserve(const Stencil &stencil, bool verbose=false)
reserve memory
Definition: spmatrix.hh:633
void addLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat)
Definition: spmatrix.hh:573
MatrixType matrix_
Definition: spmatrix.hh:706
void compress()
compress matrix to a real CRS format
Definition: spmatrix.hh:623
void resort()
resort row numbering in matrix to have ascending numbering
Definition: spmatrix.hh:693
void getLocalMatrix(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat) const
Definition: spmatrix.hh:606
void unitRow(const size_type row)
Definition: spmatrix.hh:523
const DomainSpaceType & domainSpace() const
get domain space (i.e. space that builds the rows)
Definition: spmatrix.hh:462
void addBlock(const size_type row, const size_type col, const LocalBlock &block)
Definition: spmatrix.hh:533
RangeMapperType rangeMapper_
Definition: spmatrix.hh:704
MatrixBlockType block_type
Definition: spmatrix.hh:435
MatrixType & matrix() const
get reference to storage object, for internal use
Definition: spmatrix.hh:475
void clear()
clear matrix
Definition: spmatrix.hh:617
bool preconditioning_
Definition: spmatrix.hh:707
const RangeSpaceType & rangeSpace_
Definition: spmatrix.hh:702
const RangeSpaceType & rangeSpace() const
get range space (i.e. space that builds the columns)
Definition: spmatrix.hh:468
void finalizeAssembly() const
Definition: spmatrix.hh:480
LocalMatrixTraits.
Definition: spmatrix.hh:717
SparseRowMatrixObject< DomainSpaceType, RangeSpaceType, Matrix > SparseRowMatrixObjectType
Definition: spmatrix.hh:721
RangeFieldType LittleBlockType
Definition: spmatrix.hh:726
RangeSpaceType::RangeFieldType RangeFieldType
Definition: spmatrix.hh:725
DomainSpace DomainSpaceType
Definition: spmatrix.hh:718
RangeSpace RangeSpaceType
Definition: spmatrix.hh:719
SparseRowMatrixObjectType::template LocalMatrix< MatrixObject > LocalMatrixType
Definition: spmatrix.hh:723
SparseRowMatrixObjectType::RangeMapperType RangeMapperType
Definition: spmatrix.hh:729
SparseRowMatrixObjectType::DomainMapperType DomainMapperType
Definition: spmatrix.hh:728
LocalMatrix.
Definition: spmatrix.hh:739
MatrixType & matrix_
Definition: spmatrix.hh:896
LocalMatrix(const LocalMatrix &)=delete
void clearRow(size_type localRow)
set matrix row to zero
Definition: spmatrix.hh:842
const DomainMapperType & domainMapper_
Definition: spmatrix.hh:897
void add(size_type localRow, size_type localCol, DofType value)
add value to matrix entry
Definition: spmatrix.hh:810
void init(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity)
Definition: spmatrix.hh:785
MatrixObjectType::MatrixType MatrixType
type of matrix
Definition: spmatrix.hh:752
void clearCol(size_type localCol)
set matrix column to zero
Definition: spmatrix.hh:849
void resort()
resort all global rows of matrix to have ascending numbering
Definition: spmatrix.hh:864
size_type columns() const
return number of columns
Definition: spmatrix.hh:804
ColumnIndicesType columnIndices_
Definition: spmatrix.hh:900
size_type rows() const
return number of rows
Definition: spmatrix.hh:798
RowIndicesType rowIndices_
Definition: spmatrix.hh:899
std::vector< typename DomainMapperType::SizeType > ColumnIndicesType
Definition: spmatrix.hh:769
Traits::RangeFieldType RangeFieldType
type of entries of little blocks
Definition: spmatrix.hh:755
DofType get(size_type localRow, size_type localCol) const
get matrix entry
Definition: spmatrix.hh:819
MatrixObject MatrixObjectType
type of matrix object
Definition: spmatrix.hh:742
std::vector< typename RangeMapperType::SizeType > RowIndicesType
Definition: spmatrix.hh:768
Traits::DomainMapperType DomainMapperType
type of nonblocked domain mapper
Definition: spmatrix.hh:764
Traits::RangeMapperType RangeMapperType
type of nonblocked domain mapper
Definition: spmatrix.hh:766
LocalMatrixTraits< MatrixObjectType > Traits
type of the traits
Definition: spmatrix.hh:745
LocalMatrix(const MatrixObjectType &matrixObject, const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const DomainMapperType &domainMapper, const RangeMapperType &rangeMapper)
constructor
Definition: spmatrix.hh:772
const RangeMapperType & rangeMapper_
Definition: spmatrix.hh:898
void scale(const DofType &value)
scale local matrix with a certain value
Definition: spmatrix.hh:873
RangeFieldType DofType
type of the DoFs
Definition: spmatrix.hh:758
void scale(size_type localRow, size_type localCol, DofType value)
scale matrix entry with value
Definition: spmatrix.hh:888
void clear()
clear all entries belonging to local matrix
Definition: spmatrix.hh:856
Traits::LittleBlockType LittleBlockType
type of little blocks
Definition: spmatrix.hh:761
void set(size_type localRow, size_type localCol, DofType value)
set matrix entry to value
Definition: spmatrix.hh:827
void unitRow(size_type localRow)
set matrix row to zero except diagonla entry
Definition: spmatrix.hh:835
Definition: solver/parameter.hh:15
static const int none
Definition: solver/parameter.hh:40
static const int jacobi
Definition: solver/parameter.hh:45