dune-istl  3.0-git
bcrsmatrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_ISTL_BCRSMATRIX_HH
5 #define DUNE_ISTL_BCRSMATRIX_HH
6 
7 #include <cmath>
8 #include <complex>
9 #include <set>
10 #include <iostream>
11 #include <algorithm>
12 #include <numeric>
13 #include <vector>
14 #include <map>
15 
16 #include "istlexception.hh"
17 #include "bvector.hh"
18 #include "matrixutils.hh"
19 #include <dune/common/stdstreams.hh>
20 #include <dune/common/iteratorfacades.hh>
21 #include <dune/common/typetraits.hh>
22 #include <dune/common/ftraits.hh>
23 
28 namespace Dune {
29 
69  template<typename M>
71 
73 
79  template<typename size_type>
81  {
83  double avg;
85  size_type maximum;
87  size_type overflow_total;
89 
92  double mem_ratio;
93  };
94 
96 
108  template<class M_>
110  {
111 
112  public:
113 
115  typedef M_ Matrix;
116 
118  typedef typename Matrix::block_type block_type;
119 
121  typedef typename Matrix::size_type size_type;
122 
124 
130  {
131 
132  public:
133 
135  block_type& operator[](size_type j) const
136  {
137  return _m.entry(_i,j);
138  }
139 
140 #ifndef DOXYGEN
141 
142  row_object(Matrix& m, size_type i)
143  : _m(m)
144  , _i(i)
145  {}
146 
147 #endif
148 
149  private:
150 
151  Matrix& _m;
152  size_type _i;
153 
154  };
155 
157 
164  : _m(m)
165  {
166  if (m.buildMode() != Matrix::implicit)
167  DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix in implicit build mode");
168  if (m.buildStage() != Matrix::building)
169  DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix with set size that has not been compressed() yet");
170  }
171 
173 
187  ImplicitMatrixBuilder(Matrix& m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
188  : _m(m)
189  {
190  if (m.buildStage() != Matrix::notAllocated)
191  DUNE_THROW(BCRSMatrixError,"You can only set up a matrix for this ImplicitBuilder if it has no memory allocated yet");
192  m.setBuildMode(Matrix::implicit);
193  m.setImplicitBuildModeParameters(avg_cols_per_row,overflow_fraction);
194  m.setSize(rows,cols);
195  }
196 
198  row_object operator[](size_type i) const
199  {
200  return row_object(_m,i);
201  }
202 
204  size_type N() const
205  {
206  return _m.N();
207  }
208 
210  size_type M() const
211  {
212  return _m.M();
213  }
214 
215  private:
216 
217  Matrix& _m;
218 
219  };
220 
410  template<class B, class A=std::allocator<B> >
412  {
413  friend struct MatrixDimension<BCRSMatrix>;
414  public:
415  enum BuildStage {
417  notbuilt=0,
419  notAllocated=0,
421  building=1,
426  rowSizesBuilt=2,
428  built=3
429  };
430 
431  //===== type definitions and constants
432 
434  typedef typename B::field_type field_type;
435 
437  typedef B block_type;
438 
440  typedef A allocator_type;
441 
444 
446  typedef typename A::size_type size_type;
447 
449  typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
450 
452  enum {
454  blocklevel = B::blocklevel+1
455  };
456 
458  enum BuildMode {
491  unknown
492  };
493 
494  //===== random access interface to rows of the matrix
495 
497  row_type& operator[] (size_type i)
498  {
499 #ifdef DUNE_ISTL_WITH_CHECKING
500  if (build_mode == implicit && ready != built)
501  DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
502  if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
503  if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
504 #endif
505  return r[i];
506  }
507 
509  const row_type& operator[] (size_type i) const
510  {
511 #ifdef DUNE_ISTL_WITH_CHECKING
512  if (build_mode == implicit && ready != built)
513  DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
514  if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
515  if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
516 #endif
517  return r[i];
518  }
519 
520 
521  //===== iterator interface to rows of the matrix
522 
524  template<class T>
526  : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
527  {
528 
529  public:
531  typedef typename std::remove_const<T>::type ValueType;
532 
533  friend class RandomAccessIteratorFacade<RealRowIterator<const ValueType>, const ValueType>;
534  friend class RandomAccessIteratorFacade<RealRowIterator<ValueType>, ValueType>;
535  friend class RealRowIterator<const ValueType>;
536  friend class RealRowIterator<ValueType>;
537 
539  RealRowIterator (row_type* _p, size_type _i)
540  : p(_p), i(_i)
541  {}
542 
545  : p(0), i(0)
546  {}
547 
549  : p(it.p), i(it.i)
550  {}
551 
552 
554  size_type index () const
555  {
556  return i;
557  }
558 
559  std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
560  {
561  assert(other.p==p);
562  return (other.i-i);
563  }
564 
565  std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
566  {
567  assert(other.p==p);
568  return (other.i-i);
569  }
570 
572  bool equals (const RealRowIterator<ValueType>& other) const
573  {
574  assert(other.p==p);
575  return i==other.i;
576  }
577 
579  bool equals (const RealRowIterator<const ValueType>& other) const
580  {
581  assert(other.p==p);
582  return i==other.i;
583  }
584 
585  private:
587  void increment()
588  {
589  ++i;
590  }
591 
593  void decrement()
594  {
595  --i;
596  }
597 
598  void advance(std::ptrdiff_t diff)
599  {
600  i+=diff;
601  }
602 
603  T& elementAt(std::ptrdiff_t diff) const
604  {
605  return p[i+diff];
606  }
607 
609  row_type& dereference () const
610  {
611  return p[i];
612  }
613 
614  row_type* p;
615  size_type i;
616  };
617 
621 
623  Iterator begin ()
624  {
625  return Iterator(r,0);
626  }
627 
629  Iterator end ()
630  {
631  return Iterator(r,n);
632  }
633 
636  Iterator beforeEnd ()
637  {
638  return Iterator(r,n-1);
639  }
640 
643  Iterator beforeBegin ()
644  {
645  return Iterator(r,-1);
646  }
647 
649  typedef Iterator RowIterator;
650 
652  typedef typename row_type::Iterator ColIterator;
653 
657 
658 
660  ConstIterator begin () const
661  {
662  return ConstIterator(r,0);
663  }
664 
666  ConstIterator end () const
667  {
668  return ConstIterator(r,n);
669  }
670 
673  ConstIterator beforeEnd() const
674  {
675  return ConstIterator(r,n-1);
676  }
677 
680  ConstIterator beforeBegin () const
681  {
682  return ConstIterator(r,-1);
683  }
684 
686  typedef ConstIterator ConstRowIterator;
687 
690 
691  //===== constructors & resizers
692 
693  // we use a negative overflowsize to indicate that the implicit
694  // mode parameters have not been set yet
695 
698  : build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
699  allocationSize_(0), r(0), a(0),
700  avg(0), overflowsize(-1.0)
701  {}
702 
704  BCRSMatrix (size_type _n, size_type _m, size_type _nnz, BuildMode bm)
705  : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
706  allocationSize_(0), r(0), a(0),
707  avg(0), overflowsize(-1.0)
708  {
709  allocate(_n, _m, _nnz,true,false);
710  }
711 
713  BCRSMatrix (size_type _n, size_type _m, BuildMode bm)
714  : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
715  allocationSize_(0), r(0), a(0),
716  avg(0), overflowsize(-1.0)
717  {
718  allocate(_n, _m,0,true,false);
719  }
720 
722 
732  BCRSMatrix (size_type _n, size_type _m, size_type _avg, double _overflowsize, BuildMode bm)
733  : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
734  allocationSize_(0), r(0), a(0),
735  avg(_avg), overflowsize(_overflowsize)
736  {
737  if (bm != implicit)
738  DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
739  // Prevent user from setting a negative overflowsize:
740  // 1) It doesn't make sense
741  // 2) We use a negative overflow value to indicate that the parameters
742  // have not been set yet
743  if (_overflowsize < 0.0)
744  DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
745  implicit_allocate(_n,_m);
746  }
747 
753  BCRSMatrix (const BCRSMatrix& Mat)
754  : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
755  allocationSize_(0), r(0), a(0),
756  avg(Mat.avg), overflowsize(Mat.overflowsize)
757  {
758  if (!(Mat.ready == notAllocated || Mat.ready == built))
759  DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
760 
761  // deep copy in global array
762  size_type _nnz = Mat.nnz_;
763 
764  // in case of row-wise allocation
765  if (_nnz<=0)
766  {
767  _nnz = 0;
768  for (size_type i=0; i<Mat.n; i++)
769  _nnz += Mat.r[i].getsize();
770  }
771 
772  j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
773  allocate(Mat.n, Mat.m, _nnz, true, true);
774 
775  // build window structure
776  copyWindowStructure(Mat);
777  }
778 
781  {
782  deallocate();
783  }
784 
790  {
791  if (ready == notAllocated)
792  {
793  build_mode = bm;
794  return;
795  }
796  if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
797  build_mode = bm;
798  else
799  DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
800  }
801 
817  void setSize(size_type rows, size_type columns, size_type nnz=0)
818  {
819  // deallocate already setup memory
820  deallocate();
821 
822  if (build_mode == implicit)
823  {
824  if (nnz>0)
825  DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");
826 
827  // implicit allocates differently
828  implicit_allocate(rows,columns);
829  }
830  else
831  {
832  // allocate matrix memory
833  allocate(rows, columns, nnz, true, false);
834  }
835  }
836 
845  void setImplicitBuildModeParameters(size_type _avg, double _overflow)
846  {
847  // Prevent user from setting a negative overflowsize:
848  // 1) It doesn't make sense
849  // 2) We use a negative overflow value to indicate that the parameters
850  // have not been set yet
851  if (_overflow < 0.0)
852  DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
853 
854  // make sure the parameters aren't changed after memory has been allocated
855  if (ready != notAllocated)
856  DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
857  avg = _avg;
858  overflowsize = _overflow;
859  }
860 
867  BCRSMatrix& operator= (const BCRSMatrix& Mat)
868  {
869  // return immediately when self-assignment
870  if (&Mat==this) return *this;
871 
872  if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
873  DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");
874 
875  // make it simple: ALWAYS throw away memory for a and j_
876  // and deallocate rows only if n != Mat.n
877  deallocate(n!=Mat.n);
878 
879  // reallocate the rows if required
880  if (n>0 && n!=Mat.n) {
881  // free rows
882  for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
883  rowAllocator_.destroy(riter);
884  rowAllocator_.deallocate(r,n);
885  }
886 
887  nnz_ = Mat.nnz_;
888  if (nnz_ <= 0)
889  {
890  for (size_type i=0; i<Mat.n; i++)
891  nnz_ += Mat.r[i].getsize();
892  }
893 
894  // allocate a, share j_
895  j_ = Mat.j_;
896  allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);
897 
898  // build window structure
899  copyWindowStructure(Mat);
900  return *this;
901  }
902 
904  BCRSMatrix& operator= (const field_type& k)
905  {
906 
907  if (!(ready == notAllocated || ready == built))
908  DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");
909 
910  for (size_type i=0; i<n; i++) r[i] = k;
911  return *this;
912  }
913 
914  //===== row-wise creation interface
915 
918  {
919  public:
921  CreateIterator (BCRSMatrix& _Mat, size_type _i)
922  : Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
923  {
924  if (Mat.build_mode == unknown && Mat.ready == building)
925  {
926  Mat.build_mode = row_wise;
927  }
928  if (i==0 && Mat.ready != building)
929  DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
930  if(Mat.build_mode!=row_wise)
931  DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
932  }
933 
936  {
937  // this should only be called if matrix is in creation
938  if (Mat.ready != building)
939  DUNE_THROW(BCRSMatrixError,"matrix already built up");
940 
941  // row i is defined through the pattern
942  // get memory for the row and initialize the j_ array
943  // this depends on the allocation mode
944 
945  // compute size of the row
946  size_type s = pattern.size();
947 
948  if(s>0) {
949  // update number of nonzeroes including this row
950  nnz += s;
951 
952  // alloc memory / set window
953  if (Mat.nnz_ > 0)
954  {
955  // memory is allocated in one long array
956 
957  // check if that memory is sufficient
958  if (nnz > Mat.nnz_)
959  DUNE_THROW(BCRSMatrixError,"allocated nnz too small");
960 
961  // set row i
962  Mat.r[i].set(s,nullptr,current_row.getindexptr());
963  current_row.setindexptr(current_row.getindexptr()+s);
964  }else{
965  // memory is allocated individually per row
966  // allocate and set row i
967  B* b = Mat.allocator_.allocate(s);
968  // use placement new to call constructor that allocates
969  // additional memory.
970  new (b) B[s];
971  size_type* j = Mat.sizeAllocator_.allocate(s);
972  Mat.r[i].set(s,b,j);
973  }
974  }else
975  // setup empty row
976  Mat.r[i].set(0,nullptr,nullptr);
977 
978  // initialize the j array for row i from pattern
979  std::copy(pattern.cbegin(), pattern.cend(), Mat.r[i].getindexptr());
980 
981  // now go to next row
982  i++;
983  pattern.clear();
984 
985  // check if this was last row
986  if (i==Mat.n)
987  {
988  Mat.ready = built;
989  if(Mat.nnz_ > 0)
990  {
991  // Set nnz to the exact number of nonzero blocks inserted
992  // as some methods rely on it
993  Mat.nnz_ = nnz;
994  // allocate data array
995  Mat.allocateData();
996  Mat.setDataPointers();
997  }
998  }
999  // done
1000  return *this;
1001  }
1002 
1004  bool operator!= (const CreateIterator& it) const
1005  {
1006  return (i!=it.i) || (&Mat!=&it.Mat);
1007  }
1008 
1010  bool operator== (const CreateIterator& it) const
1011  {
1012  return (i==it.i) && (&Mat==&it.Mat);
1013  }
1014 
1016  size_type index () const
1017  {
1018  return i;
1019  }
1020 
1022  void insert (size_type j)
1023  {
1024  pattern.insert(j);
1025  }
1026 
1028  bool contains (size_type j)
1029  {
1030  return pattern.find(j) != pattern.end();
1031  }
1037  size_type size() const
1038  {
1039  return pattern.size();
1040  }
1041 
1042  private:
1043  BCRSMatrix& Mat; // the matrix we are defining
1044  size_type i; // current row to be defined
1045  size_type nnz; // count total number of nonzeros
1046  typedef std::set<size_type,std::less<size_type> > PatternType;
1047  PatternType pattern; // used to compile entries in a row
1048  row_type current_row; // row pointing to the current row to setup
1049  };
1050 
1052  friend class CreateIterator;
1053 
1056  {
1057  return CreateIterator(*this,0);
1058  }
1059 
1062  {
1063  return CreateIterator(*this,n);
1064  }
1065 
1066 
1067  //===== random creation interface
1068 
1075  void setrowsize (size_type i, size_type s)
1076  {
1077  if (build_mode!=random)
1078  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1079  if (ready != building)
1080  DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1081 
1082  r[i].setsize(s);
1083  }
1084 
1086  size_type getrowsize (size_type i) const
1087  {
1088 #ifdef DUNE_ISTL_WITH_CHECKING
1089  if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
1090  if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
1091 #endif
1092  return r[i].getsize();
1093  }
1094 
1096  void incrementrowsize (size_type i, size_type s = 1)
1097  {
1098  if (build_mode!=random)
1099  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1100  if (ready != building)
1101  DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1102 
1103  r[i].setsize(r[i].getsize()+s);
1104  }
1105 
1107  void endrowsizes ()
1108  {
1109  if (build_mode!=random)
1110  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1111  if (ready != building)
1112  DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1113 
1114  // compute total size, check positivity
1115  size_type total=0;
1116  for (size_type i=0; i<n; i++)
1117  {
1118  total += r[i].getsize();
1119  }
1120 
1121  if(nnz_ == 0)
1122  // allocate/check memory
1123  allocate(n,m,total,false,false);
1124  else if(nnz_ < total)
1125  DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
1126  <<"sufficient for calculated nonzeros ("<<total<<"! ");
1127 
1128  // set the window pointers correctly
1129  setColumnPointers(begin());
1130 
1131  // initialize j_ array with m (an invalid column index)
1132  // this indicates an unused entry
1133  for (size_type k=0; k<nnz_; k++)
1134  j_.get()[k] = m;
1135  ready = rowSizesBuilt;
1136  }
1137 
1139 
1149  void addindex (size_type row, size_type col)
1150  {
1151  if (build_mode!=random)
1152  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1153  if (ready==built)
1154  DUNE_THROW(BCRSMatrixError,"matrix already built up");
1155  if (ready==building)
1156  DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
1157  if (ready==notAllocated)
1158  DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1159 
1160  if (col >= m)
1161  DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1162 
1163  // get row range
1164  size_type* const first = r[row].getindexptr();
1165  size_type* const last = first + r[row].getsize();
1166 
1167  // find correct insertion position for new column index
1168  size_type* pos = std::lower_bound(first,last,col);
1169 
1170  // check if index is already in row
1171  if (pos!=last && *pos == col) return;
1172 
1173  // find end of already inserted column indices
1174  size_type* end = std::lower_bound(pos,last,m);
1175  if (end==last)
1176  DUNE_THROW(BCRSMatrixError,"row is too small");
1177 
1178  // insert new column index at correct position
1179  std::copy_backward(pos,end,end+1);
1180  *pos = col;
1181  }
1182 
1184 
1191  template<typename It>
1192  void setIndices(size_type row, It begin, It end)
1193  {
1194  size_type row_size = r[row].size();
1195  size_type* col_begin = r[row].getindexptr();
1196  size_type* col_end;
1197  // consistency check between allocated row size and number of passed column indices
1198  if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
1199  DUNE_THROW(BCRSMatrixError,"Given size of row " << row
1200  << " (" << row_size
1201  << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
1202  std::sort(col_begin,col_end);
1203  }
1204 
1206  void endindices ()
1207  {
1208  if (build_mode!=random)
1209  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1210  if (ready==built)
1211  DUNE_THROW(BCRSMatrixError,"matrix already built up");
1212  if (ready==building)
1213  DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
1214  if (ready==notAllocated)
1215  DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1216 
1217  // check if there are undefined indices
1218  RowIterator endi=end();
1219  for (RowIterator i=begin(); i!=endi; ++i)
1220  {
1221  ColIterator endj = (*i).end();
1222  for (ColIterator j=(*i).begin(); j!=endj; ++j) {
1223  if (j.index() >= m) {
1224  dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
1225  <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
1226  r[i.index()].setsize(j.offset());
1227  break;
1228  }
1229  }
1230  }
1231 
1232  allocateData();
1233  setDataPointers();
1234 
1235  // if not, set matrix to built
1236  ready = built;
1237  }
1238 
1239  //===== implicit creation interface
1240 
1242 
1253  B& entry(size_type row, size_type col)
1254  {
1255 #ifdef DUNE_ISTL_WITH_CHECKING
1256  if (build_mode!=implicit)
1257  DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1258  if (ready==built)
1259  DUNE_THROW(BCRSMatrixError,"matrix already built up, use operator[] for entry access now");
1260  if (ready==notAllocated)
1261  DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1262  if (ready!=building)
1263  DUNE_THROW(InvalidStateException,"You may only use entry() during the 'building' stage");
1264 
1265  if (row >= n)
1266  DUNE_THROW(BCRSMatrixError,"row index exceeds matrix size");
1267  if (col >= m)
1268  DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1269 #endif
1270 
1271  size_type* begin = r[row].getindexptr();
1272  size_type* end = begin + r[row].getsize();
1273 
1274  size_type* pos = std::find(begin, end, col);
1275 
1276  //treat the case that there was a match in the array
1277  if (pos != end)
1278  if (*pos == col)
1279  {
1280  std::ptrdiff_t offset = pos - r[row].getindexptr();
1281  B* aptr = r[row].getptr() + offset;
1282 
1283  return *aptr;
1284  }
1285 
1286  //determine whether overflow has to be taken into account or not
1287  if (r[row].getsize() == avg)
1288  return overflow[std::make_pair(row,col)];
1289  else
1290  {
1291  //modify index array
1292  *end = col;
1293 
1294  //do simulatenous operations on data array a
1295  std::ptrdiff_t offset = end - r[row].getindexptr();
1296  B* apos = r[row].getptr() + offset;
1297 
1298  //increase rowsize
1299  r[row].setsize(r[row].getsize()+1);
1300 
1301  //return reference to the newly created entry
1302  return *apos;
1303  }
1304  }
1305 
1307 
1317  CompressionStatistics compress()
1318  {
1319  if (build_mode!=implicit)
1320  DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1321  if (ready==built)
1322  DUNE_THROW(BCRSMatrixError,"matrix already built up, no more need for compression");
1323  if (ready==notAllocated)
1324  DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1325  if (ready!=building)
1326  DUNE_THROW(InvalidStateException,"You may only call compress() at the end of the 'building' stage");
1327 
1328  //calculate statistics
1329  CompressionStatistics stats;
1330  stats.overflow_total = overflow.size();
1331  stats.maximum = 0;
1332 
1333  //get insertion iterators pointing to one before start (for later use of ++it)
1334  size_type* jiit = j_.get();
1335  B* aiit = a;
1336 
1337  //get iterator to the smallest overflow element
1338  typename OverflowType::iterator oit = overflow.begin();
1339 
1340  //store a copy of index pointers on which to perform sortation
1341  std::vector<size_type*> perm;
1342 
1343  //iterate over all rows and copy elements into their position in the compressed array
1344  for (size_type i=0; i<n; i++)
1345  {
1346  //get old pointers into a and j and size without overflow changes
1347  size_type* begin = r[i].getindexptr();
1348  //B* apos = r[i].getptr();
1349  size_type size = r[i].getsize();
1350 
1351  perm.resize(size);
1352 
1353  typename std::vector<size_type*>::iterator it = perm.begin();
1354  for (size_type* iit = begin; iit < begin + size; ++iit, ++it)
1355  *it = iit;
1356 
1357  //sort permutation array
1358  std::sort(perm.begin(),perm.end(),PointerCompare<size_type>());
1359 
1360  //change row window pointer to their new positions
1361  r[i].setindexptr(jiit);
1362  r[i].setptr(aiit);
1363 
1364  for (it = perm.begin(); it != perm.end(); ++it)
1365  {
1366  //check whether there are elements in the overflow area which take precedence
1367  while ((oit!=overflow.end()) && (oit->first < std::make_pair(i,**it)))
1368  {
1369  //check whether there is enough memory to write to
1370  if (jiit > begin)
1372  "Allocated memory for BCRSMatrix exhausted during compress()!"
1373  "Please increase either the average number of entries per row or the overflow fraction."
1374  );
1375  //copy an element from the overflow area to the insertion position in a and j
1376  *jiit = oit->first.second;
1377  ++jiit;
1378  *aiit = oit->second;
1379  ++aiit;
1380  ++oit;
1381  r[i].setsize(r[i].getsize()+1);
1382  }
1383 
1384  //check whether there is enough memory to write to
1385  if (jiit > begin)
1387  "Allocated memory for BCRSMatrix exhausted during compress()!"
1388  "Please increase either the average number of entries per row or the overflow fraction."
1389  );
1390 
1391  //copy element from array
1392  *jiit = **it;
1393  ++jiit;
1394  B* apos = *it - j_.get() + a;
1395  *aiit = *apos;
1396  ++aiit;
1397  }
1398 
1399  //copy remaining elements from the overflow area
1400  while ((oit!=overflow.end()) && (oit->first.first == i))
1401  {
1402  //check whether there is enough memory to write to
1403  if (jiit > begin)
1405  "Allocated memory for BCRSMatrix exhausted during compress()!"
1406  "Please increase either the average number of entries per row or the overflow fraction."
1407  );
1408 
1409  //copy and element from the overflow area to the insertion position in a and j
1410  *jiit = oit->first.second;
1411  ++jiit;
1412  *aiit = oit->second;
1413  ++aiit;
1414  ++oit;
1415  r[i].setsize(r[i].getsize()+1);
1416  }
1417 
1418  // update maximum row size
1419  if (r[i].getsize()>stats.maximum)
1420  stats.maximum = r[i].getsize();
1421  }
1422 
1423  // overflow area may be cleared
1424  overflow.clear();
1425 
1426  //determine average number of entries and memory usage
1427  std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
1428  nnz_ = diff;
1429  stats.avg = (double) (nnz_) / (double) n;
1430  stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
1431 
1432  //matrix is now built
1433  ready = built;
1434 
1435  return stats;
1436  }
1437 
1438  //===== vector space arithmetic
1439 
1441  BCRSMatrix& operator*= (const field_type& k)
1442  {
1443 #ifdef DUNE_ISTL_WITH_CHECKING
1444  if (ready != built)
1445  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1446 #endif
1447 
1448  if (nnz_ > 0)
1449  {
1450  // process 1D array
1451  for (size_type i=0; i<nnz_; i++)
1452  a[i] *= k;
1453  }
1454  else
1455  {
1456  RowIterator endi=end();
1457  for (RowIterator i=begin(); i!=endi; ++i)
1458  {
1459  ColIterator endj = (*i).end();
1460  for (ColIterator j=(*i).begin(); j!=endj; ++j)
1461  (*j) *= k;
1462  }
1463  }
1464 
1465  return *this;
1466  }
1467 
1469  BCRSMatrix& operator/= (const field_type& k)
1470  {
1471 #ifdef DUNE_ISTL_WITH_CHECKING
1472  if (ready != built)
1473  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1474 #endif
1475 
1476  if (nnz_ > 0)
1477  {
1478  // process 1D array
1479  for (size_type i=0; i<nnz_; i++)
1480  a[i] /= k;
1481  }
1482  else
1483  {
1484  RowIterator endi=end();
1485  for (RowIterator i=begin(); i!=endi; ++i)
1486  {
1487  ColIterator endj = (*i).end();
1488  for (ColIterator j=(*i).begin(); j!=endj; ++j)
1489  (*j) /= k;
1490  }
1491  }
1492 
1493  return *this;
1494  }
1495 
1496 
1502  BCRSMatrix& operator+= (const BCRSMatrix& b)
1503  {
1504 #ifdef DUNE_ISTL_WITH_CHECKING
1505  if (ready != built || b.ready != built)
1506  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1507  if(N()!=b.N() || M() != b.M())
1508  DUNE_THROW(RangeError, "Matrix sizes do not match!");
1509 #endif
1510  RowIterator endi=end();
1511  ConstRowIterator j=b.begin();
1512  for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1513  i->operator+=(*j);
1514  }
1515 
1516  return *this;
1517  }
1518 
1524  BCRSMatrix& operator-= (const BCRSMatrix& b)
1525  {
1526 #ifdef DUNE_ISTL_WITH_CHECKING
1527  if (ready != built || b.ready != built)
1528  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1529  if(N()!=b.N() || M() != b.M())
1530  DUNE_THROW(RangeError, "Matrix sizes do not match!");
1531 #endif
1532  RowIterator endi=end();
1533  ConstRowIterator j=b.begin();
1534  for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1535  i->operator-=(*j);
1536  }
1537 
1538  return *this;
1539  }
1540 
1549  BCRSMatrix& axpy(field_type alpha, const BCRSMatrix& b)
1550  {
1551 #ifdef DUNE_ISTL_WITH_CHECKING
1552  if (ready != built || b.ready != built)
1553  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1554  if(N()!=b.N() || M() != b.M())
1555  DUNE_THROW(RangeError, "Matrix sizes do not match!");
1556 #endif
1557  RowIterator endi=end();
1558  ConstRowIterator j=b.begin();
1559  for(RowIterator i=begin(); i!=endi; ++i, ++j)
1560  i->axpy(alpha, *j);
1561 
1562  return *this;
1563  }
1564 
1565  //===== linear maps
1566 
1568  template<class X, class Y>
1569  void mv (const X& x, Y& y) const
1570  {
1571 #ifdef DUNE_ISTL_WITH_CHECKING
1572  if (ready != built)
1573  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1574  if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,
1575  "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
1576  if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,
1577  "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
1578 #endif
1579  ConstRowIterator endi=end();
1580  for (ConstRowIterator i=begin(); i!=endi; ++i)
1581  {
1582  y[i.index()]=0;
1583  ConstColIterator endj = (*i).end();
1584  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1585  (*j).umv(x[j.index()],y[i.index()]);
1586  }
1587  }
1588 
1590  template<class X, class Y>
1591  void umv (const X& x, Y& y) const
1592  {
1593 #ifdef DUNE_ISTL_WITH_CHECKING
1594  if (ready != built)
1595  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1596  if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1597  if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1598 #endif
1599  ConstRowIterator endi=end();
1600  for (ConstRowIterator i=begin(); i!=endi; ++i)
1601  {
1602  ConstColIterator endj = (*i).end();
1603  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1604  (*j).umv(x[j.index()],y[i.index()]);
1605  }
1606  }
1607 
1609  template<class X, class Y>
1610  void mmv (const X& x, Y& y) const
1611  {
1612 #ifdef DUNE_ISTL_WITH_CHECKING
1613  if (ready != built)
1614  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1615  if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1616  if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1617 #endif
1618  ConstRowIterator endi=end();
1619  for (ConstRowIterator i=begin(); i!=endi; ++i)
1620  {
1621  ConstColIterator endj = (*i).end();
1622  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1623  (*j).mmv(x[j.index()],y[i.index()]);
1624  }
1625  }
1626 
1628  template<class X, class Y, class F>
1629  void usmv (F&& alpha, const X& x, Y& y) const
1630  {
1631 #ifdef DUNE_ISTL_WITH_CHECKING
1632  if (ready != built)
1633  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1634  if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1635  if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1636 #endif
1637  ConstRowIterator endi=end();
1638  for (ConstRowIterator i=begin(); i!=endi; ++i)
1639  {
1640  ConstColIterator endj = (*i).end();
1641  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1642  (*j).usmv(alpha,x[j.index()],y[i.index()]);
1643  }
1644  }
1645 
1647  template<class X, class Y>
1648  void mtv (const X& x, Y& y) const
1649  {
1650 #ifdef DUNE_ISTL_WITH_CHECKING
1651  if (ready != built)
1652  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1653  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1654  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1655 #endif
1656  for(size_type i=0; i<y.N(); ++i)
1657  y[i]=0;
1658  umtv(x,y);
1659  }
1660 
1662  template<class X, class Y>
1663  void umtv (const X& x, Y& y) const
1664  {
1665 #ifdef DUNE_ISTL_WITH_CHECKING
1666  if (ready != built)
1667  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1668  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1669  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1670 #endif
1671  ConstRowIterator endi=end();
1672  for (ConstRowIterator i=begin(); i!=endi; ++i)
1673  {
1674  ConstColIterator endj = (*i).end();
1675  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1676  (*j).umtv(x[i.index()],y[j.index()]);
1677  }
1678  }
1679 
1681  template<class X, class Y>
1682  void mmtv (const X& x, Y& y) const
1683  {
1684 #ifdef DUNE_ISTL_WITH_CHECKING
1685  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1686  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1687 #endif
1688  ConstRowIterator endi=end();
1689  for (ConstRowIterator i=begin(); i!=endi; ++i)
1690  {
1691  ConstColIterator endj = (*i).end();
1692  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1693  (*j).mmtv(x[i.index()],y[j.index()]);
1694  }
1695  }
1696 
1698  template<class X, class Y>
1699  void usmtv (const field_type& alpha, const X& x, Y& y) const
1700  {
1701 #ifdef DUNE_ISTL_WITH_CHECKING
1702  if (ready != built)
1703  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1704  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1705  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1706 #endif
1707  ConstRowIterator endi=end();
1708  for (ConstRowIterator i=begin(); i!=endi; ++i)
1709  {
1710  ConstColIterator endj = (*i).end();
1711  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1712  (*j).usmtv(alpha,x[i.index()],y[j.index()]);
1713  }
1714  }
1715 
1717  template<class X, class Y>
1718  void umhv (const X& x, Y& y) const
1719  {
1720 #ifdef DUNE_ISTL_WITH_CHECKING
1721  if (ready != built)
1722  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1723  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1724  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1725 #endif
1726  ConstRowIterator endi=end();
1727  for (ConstRowIterator i=begin(); i!=endi; ++i)
1728  {
1729  ConstColIterator endj = (*i).end();
1730  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1731  (*j).umhv(x[i.index()],y[j.index()]);
1732  }
1733  }
1734 
1736  template<class X, class Y>
1737  void mmhv (const X& x, Y& y) const
1738  {
1739 #ifdef DUNE_ISTL_WITH_CHECKING
1740  if (ready != built)
1741  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1742  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1743  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1744 #endif
1745  ConstRowIterator endi=end();
1746  for (ConstRowIterator i=begin(); i!=endi; ++i)
1747  {
1748  ConstColIterator endj = (*i).end();
1749  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1750  (*j).mmhv(x[i.index()],y[j.index()]);
1751  }
1752  }
1753 
1755  template<class X, class Y>
1756  void usmhv (const field_type& alpha, const X& x, Y& y) const
1757  {
1758 #ifdef DUNE_ISTL_WITH_CHECKING
1759  if (ready != built)
1760  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1761  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1762  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1763 #endif
1764  ConstRowIterator endi=end();
1765  for (ConstRowIterator i=begin(); i!=endi; ++i)
1766  {
1767  ConstColIterator endj = (*i).end();
1768  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1769  (*j).usmhv(alpha,x[i.index()],y[j.index()]);
1770  }
1771  }
1772 
1773 
1774  //===== norms
1775 
1777  typename FieldTraits<field_type>::real_type frobenius_norm2 () const
1778  {
1779 #ifdef DUNE_ISTL_WITH_CHECKING
1780  if (ready != built)
1781  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1782 #endif
1783 
1784  double sum=0;
1785 
1786  ConstRowIterator endi=end();
1787  for (ConstRowIterator i=begin(); i!=endi; ++i)
1788  {
1789  ConstColIterator endj = (*i).end();
1790  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1791  sum += (*j).frobenius_norm2();
1792  }
1793 
1794  return sum;
1795  }
1796 
1798  typename FieldTraits<field_type>::real_type frobenius_norm () const
1799  {
1800  return sqrt(frobenius_norm2());
1801  }
1802 
1804  template <typename ft = field_type,
1805  typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
1806  typename FieldTraits<ft>::real_type infinity_norm() const {
1807  if (ready != built)
1808  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1809 
1810  using real_type = typename FieldTraits<ft>::real_type;
1811  using std::max;
1812 
1813  real_type norm = 0;
1814  for (auto const &x : *this) {
1815  real_type sum = 0;
1816  for (auto const &y : x)
1817  sum += y.infinity_norm();
1818  norm = max(sum, norm);
1819  }
1820  return norm;
1821  }
1822 
1824  template <typename ft = field_type,
1825  typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
1826  typename FieldTraits<ft>::real_type infinity_norm_real() const {
1827  if (ready != built)
1828  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1829 
1830  using real_type = typename FieldTraits<ft>::real_type;
1831  using std::max;
1832 
1833  real_type norm = 0;
1834  for (auto const &x : *this) {
1835  real_type sum = 0;
1836  for (auto const &y : x)
1837  sum += y.infinity_norm_real();
1838  norm = max(sum, norm);
1839  }
1840  return norm;
1841  }
1842 
1844  template <typename ft = field_type,
1845  typename std::enable_if<has_nan<ft>::value, int>::type = 0>
1846  typename FieldTraits<ft>::real_type infinity_norm() const {
1847  if (ready != built)
1848  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1849 
1850  using real_type = typename FieldTraits<ft>::real_type;
1851  using std::max;
1852 
1853  real_type norm = 0;
1854  real_type isNaN = 1;
1855  for (auto const &x : *this) {
1856  real_type sum = 0;
1857  for (auto const &y : x)
1858  sum += y.infinity_norm();
1859  norm = max(sum, norm);
1860  isNaN += sum;
1861  }
1862  isNaN /= isNaN;
1863  return norm * isNaN;
1864  }
1865 
1867  template <typename ft = field_type,
1868  typename std::enable_if<has_nan<ft>::value, int>::type = 0>
1869  typename FieldTraits<ft>::real_type infinity_norm_real() const {
1870  if (ready != built)
1871  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1872 
1873  using real_type = typename FieldTraits<ft>::real_type;
1874  using std::max;
1875 
1876  real_type norm = 0;
1877  real_type isNaN = 1;
1878  for (auto const &x : *this) {
1879  real_type sum = 0;
1880  for (auto const &y : x)
1881  sum += y.infinity_norm_real();
1882  norm = max(sum, norm);
1883  isNaN += sum;
1884  }
1885  isNaN /= isNaN;
1886  return norm * isNaN;
1887  }
1888 
1889  //===== sizes
1890 
1892  size_type N () const
1893  {
1894  return n;
1895  }
1896 
1898  size_type M () const
1899  {
1900  return m;
1901  }
1902 
1904  size_type nonzeroes () const
1905  {
1906  return nnz_;
1907  }
1908 
1911  {
1912  return ready;
1913  }
1914 
1917  {
1918  return build_mode;
1919  }
1920 
1921  //===== query
1922 
1924  bool exists (size_type i, size_type j) const
1925  {
1926 #ifdef DUNE_ISTL_WITH_CHECKING
1927  if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
1928  if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
1929 #endif
1930  return (r[i].size() && r[i].find(j) != r[i].end());
1931  }
1932 
1933 
1934  protected:
1935  // state information
1936  BuildMode build_mode; // row wise or whole matrix
1937  BuildStage ready; // indicate the stage the matrix building is in
1938 
1939  // The allocator used for memory management
1940  typename A::template rebind<B>::other allocator_;
1941 
1942  typename A::template rebind<row_type>::other rowAllocator_;
1943 
1944  typename A::template rebind<size_type>::other sizeAllocator_;
1945 
1946  // size of the matrix
1947  size_type n; // number of rows
1948  size_type m; // number of columns
1949  size_type nnz_; // number of nonzeroes contained in the matrix
1950  size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
1951  // zero means that memory is allocated separately for each row.
1952 
1953  // the rows are dynamically allocated
1954  row_type* r; // [n] the individual rows having pointers into a,j arrays
1955 
1956  // dynamically allocated memory
1957  B* a; // [allocationSize] non-zero entries of the matrix in row-wise ordering
1958  // If a single array of column indices is used, it can be shared
1959  // between different matrices with the same sparsity pattern
1960  std::shared_ptr<size_type> j_; // [allocationSize] column indices of entries
1961 
1962  // additional data is needed in implicit buildmode
1963  size_type avg;
1965 
1966  typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
1967  OverflowType overflow;
1968 
1969  void setWindowPointers(ConstRowIterator row)
1970  {
1971  row_type current_row(a,j_.get(),0); // Pointers to current row data
1972  for (size_type i=0; i<n; i++, ++row) {
1973  // set row i
1974  size_type s = row->getsize();
1975 
1976  if (s>0) {
1977  // setup pointers and size
1978  r[i].set(s,current_row.getptr(), current_row.getindexptr());
1979  // update pointer for next row
1980  current_row.setptr(current_row.getptr()+s);
1981  current_row.setindexptr(current_row.getindexptr()+s);
1982  } else{
1983  // empty row
1984  r[i].set(0,nullptr,nullptr);
1985  }
1986  }
1987  }
1988 
1990 
1994  void setColumnPointers(ConstRowIterator row)
1995  {
1996  size_type* jptr = j_.get();
1997  for (size_type i=0; i<n; ++i, ++row) {
1998  // set row i
1999  size_type s = row->getsize();
2000 
2001  if (s>0) {
2002  // setup pointers and size
2003  r[i].setsize(s);
2004  r[i].setindexptr(jptr);
2005  } else{
2006  // empty row
2007  r[i].set(0,nullptr,nullptr);
2008  }
2009 
2010  // advance position in global array
2011  jptr += s;
2012  }
2013  }
2014 
2016 
2021  {
2022  B* aptr = a;
2023  for (size_type i=0; i<n; ++i) {
2024  // set row i
2025  if (r[i].getsize() > 0) {
2026  // setup pointers and size
2027  r[i].setptr(aptr);
2028  } else{
2029  // empty row
2030  r[i].set(0,nullptr,nullptr);
2031  }
2032 
2033  // advance position in global array
2034  aptr += r[i].getsize();
2035  }
2036  }
2037 
2040  {
2041  setWindowPointers(Mat.begin());
2042 
2043  // copy data
2044  for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
2045 
2046  // finish off
2047  build_mode = row_wise; // dummy
2048  ready = built;
2049  }
2050 
2056  void deallocate(bool deallocateRows=true)
2057  {
2058 
2059  if (notAllocated)
2060  return;
2061 
2062  if (allocationSize_>0)
2063  {
2064  // a,j_ have been allocated as one long vector
2065  j_.reset();
2066  if (a)
2067  {
2068  for(B *aiter=a+(allocationSize_-1), *aend=a-1; aiter!=aend; --aiter)
2069  allocator_.destroy(aiter);
2070  allocator_.deallocate(a,allocationSize_);
2071  a = nullptr;
2072  }
2073  }
2074  else if (r)
2075  {
2076  // check if memory for rows have been allocated individually
2077  for (size_type i=0; i<n; i++)
2078  if (r[i].getsize()>0)
2079  {
2080  for (B *col=r[i].getptr()+(r[i].getsize()-1),
2081  *colend = r[i].getptr()-1; col!=colend; --col) {
2082  allocator_.destroy(col);
2083  }
2084  sizeAllocator_.deallocate(r[i].getindexptr(),1);
2085  allocator_.deallocate(r[i].getptr(),1);
2086  // clear out row data in case we don't want to deallocate the rows
2087  // otherwise we might run into a double free problem here later
2088  r[i].set(0,nullptr,nullptr);
2089  }
2090  }
2091 
2092  // deallocate the rows
2093  if (n>0 && deallocateRows && r) {
2094  for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
2095  rowAllocator_.destroy(riter);
2096  rowAllocator_.deallocate(r,n);
2097  r = nullptr;
2098  }
2099 
2100  // Mark matrix as not built at all.
2101  ready=notAllocated;
2102 
2103  }
2104 
2107  {
2108  typename A::template rebind<size_type>::other& sizeAllocator_;
2109 
2110  public:
2111  Deallocator(typename A::template rebind<size_type>::other& sizeAllocator)
2112  : sizeAllocator_(sizeAllocator)
2113  {}
2114 
2115  void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); }
2116  };
2117 
2118 
2136  void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
2137  {
2138  // Store size
2139  n = rows;
2140  m = columns;
2141  nnz_ = allocationSize;
2142  allocationSize_ = allocationSize;
2143 
2144  // allocate rows
2145  if(allocateRows) {
2146  if (n>0) {
2147  if (r)
2148  DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
2149  r = rowAllocator_.allocate(rows);
2150  }else{
2151  r = 0;
2152  }
2153  }
2154 
2155  // allocate a and j_ array
2156  if (allocate_data)
2157  allocateData();
2158  if (allocationSize_>0) {
2159  // allocate column indices only if not yet present (enable sharing)
2160  if (!j_.get())
2161  j_.reset(sizeAllocator_.allocate(allocationSize_),Deallocator(sizeAllocator_));
2162  }else{
2163  j_.reset();
2164  for(row_type* ri=r; ri!=r+rows; ++ri)
2165  rowAllocator_.construct(ri, row_type());
2166  }
2167 
2168  // Mark the matrix as not built.
2169  ready = building;
2170  }
2171 
2173  {
2174  if (a)
2175  DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
2176  if (allocationSize_>0) {
2177  a = allocator_.allocate(allocationSize_);
2178  // use placement new to call constructor that allocates
2179  // additional memory.
2180  new (a) B[allocationSize_];
2181  } else {
2182  a = nullptr;
2183  }
2184  }
2185 
2191  void implicit_allocate(size_type _n, size_type _m)
2192  {
2193  if (build_mode != implicit)
2194  DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
2195  if (ready != notAllocated)
2196  DUNE_THROW(InvalidStateException,"memory has already been allocated");
2197 
2198  // check to make sure the user has actually set the parameters
2199  if (overflowsize < 0)
2200  DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
2201  //calculate size of overflow region, add buffer for row sort!
2202  size_type osize = (size_type) (_n*avg)*overflowsize + 4*avg;
2203  allocationSize_ = _n*avg + osize;
2204 
2205  allocate(_n, _m, allocationSize_,true,true);
2206 
2207  //set row pointers correctly
2208  size_type* jptr = j_.get() + osize;
2209  B* aptr = a + osize;
2210  for (size_type i=0; i<n; i++)
2211  {
2212  r[i].set(0,aptr,jptr);
2213  jptr = jptr + avg;
2214  aptr = aptr + avg;
2215  }
2216 
2217  ready = building;
2218  }
2219  };
2220 
2221 
2224 } // end namespace
2225 
2226 #endif
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1022
Class used by shared_ptr to deallocate memory using the proper allocator.
Definition: bcrsmatrix.hh:2106
BuildMode
we support two modes
Definition: bcrsmatrix.hh:458
Col col
Definition: matrixmatrix.hh:347
ConstIterator begin() const
Get const iterator to first row.
Definition: bcrsmatrix.hh:660
OverflowType overflow
Definition: bcrsmatrix.hh:1967
BCRSMatrix(size_type _n, size_type _m, BuildMode bm)
matrix with unknown number of nonzeroes
Definition: bcrsmatrix.hh:713
::Dune::CompressionStatistics< size_type > CompressionStatistics
The type for the statistics object returned by compress()
Definition: bcrsmatrix.hh:449
BCRSMatrix(size_type _n, size_type _m, size_type _nnz, BuildMode bm)
matrix with known number of nonzeroes
Definition: bcrsmatrix.hh:704
void mtv(const X &x, Y &y) const
y = A^T x
Definition: bcrsmatrix.hh:1648
void setImplicitBuildModeParameters(size_type _avg, double _overflow)
Set parameters needed for creation in implicit build mode.
Definition: bcrsmatrix.hh:845
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:571
Definition: bvector.hh:1049
void umhv(const X &x, Y &y) const
y += A^H x
Definition: bcrsmatrix.hh:1718
void setsize(size_type _n)
set size only
Definition: bvector.hh:1157
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: bcrsmatrix.hh:1798
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1569
std::ptrdiff_t distanceTo(const RealRowIterator< ValueType > &other) const
Definition: bcrsmatrix.hh:559
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:446
size_type allocationSize_
Definition: bcrsmatrix.hh:1950
size_type M() const
The number of columns in the matrix.
Definition: bcrsmatrix.hh:210
BuildStage
Definition: bcrsmatrix.hh:415
ConstIterator beforeBegin() const
Definition: bcrsmatrix.hh:680
void setrowsize(size_type i, size_type s)
Set number of indices in row i to s.
Definition: bcrsmatrix.hh:1075
~BCRSMatrix()
destructor
Definition: bcrsmatrix.hh:780
size_type m
Definition: bcrsmatrix.hh:1948
ImplicitMatrixBuilder(Matrix &m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
Sets up matrix m for implicit construction using the given parameters and creates an ImplicitBmatrixu...
Definition: bcrsmatrix.hh:187
size_type nonzeroes() const
number of blocks that are stored (the number of blocks that possibly are nonzero) ...
Definition: bcrsmatrix.hh:1904
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1206
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:629
void set(size_type _n, B *_p, size_type *_j)
set size and pointer
Definition: bvector.hh:1149
Definition: matrixutils.hh:511
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:623
BCRSMatrix()
an empty matrix
Definition: bcrsmatrix.hh:697
B::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:434
size_type n
Definition: bcrsmatrix.hh:1947
std::remove_const< T >::type ValueType
The unqualified value type.
Definition: bcrsmatrix.hh:531
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:789
double avg
average number of non-zeroes per row.
Definition: bcrsmatrix.hh:83
void setWindowPointers(ConstRowIterator row)
Definition: bcrsmatrix.hh:1969
ConstIterator end() const
Get const iterator to one beyond last row.
Definition: bcrsmatrix.hh:666
bool equals(const RealRowIterator< ValueType > &other) const
equality
Definition: bcrsmatrix.hh:572
Definition: bcrsmatrix.hh:70
Error specific to BCRSMatrix.
Definition: istlexception.hh:19
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: bcrsmatrix.hh:1826
void operator()(size_type *p)
Definition: bcrsmatrix.hh:2115
RealRowIterator< row_type > iterator
The iterator over the (mutable matrix rows.
Definition: bcrsmatrix.hh:619
double mem_ratio
fraction of wasted memory resulting from non-used overflow area.
Definition: bcrsmatrix.hh:92
T block_type
Export the type representing the components.
Definition: matrix.hh:562
B & entry(size_type row, size_type col)
Returns reference to entry (row,col) of the matrix.
Definition: bcrsmatrix.hh:1253
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:292
Definition: matrixmarket.hh:258
Iterator beforeBegin()
Definition: bcrsmatrix.hh:643
block_type & operator[](size_type j) const
Returns entry in column j.
Definition: bcrsmatrix.hh:135
void copyWindowStructure(const BCRSMatrix &Mat)
Copy the window structure from another matrix.
Definition: bcrsmatrix.hh:2039
size_type avg
Definition: bcrsmatrix.hh:1963
bool contains(size_type j)
return true if column index is in row
Definition: bcrsmatrix.hh:1028
CompressedBlockVectorWindow< B, A > row_type
implement row_type with compressed vector
Definition: bcrsmatrix.hh:443
CreateIterator(BCRSMatrix &_Mat, size_type _i)
constructor
Definition: bcrsmatrix.hh:921
BuildStage ready
Definition: bcrsmatrix.hh:1937
void setindexptr(size_type *_j)
set pointer only
Definition: bvector.hh:1169
Build entries randomly with an educated guess for the number of entries per row.
Definition: bcrsmatrix.hh:487
B * a
Definition: bcrsmatrix.hh:1957
BuildStage buildStage() const
The current build stage of the matrix.
Definition: bcrsmatrix.hh:1910
Iterator RowIterator
rename the iterators for easier access
Definition: bcrsmatrix.hh:649
Matrix::block_type block_type
The block_type of the underlying matrix.
Definition: bcrsmatrix.hh:118
Matrix::size_type size_type
The size_type of the underlying matrix.
Definition: bcrsmatrix.hh:121
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: bcrsmatrix.hh:1699
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1061
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: bcrsmatrix.hh:1806
BuildMode build_mode
Definition: bcrsmatrix.hh:1936
size_type getsize() const
get size
Definition: bvector.hh:1198
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:411
void usmv(F &&alpha, const X &x, Y &y) const
y += alpha A x
Definition: bcrsmatrix.hh:1629
std::map< std::pair< size_type, size_type >, B > OverflowType
Definition: bcrsmatrix.hh:1966
void mmv(const X &x, Y &y) const
y -= A x
Definition: bcrsmatrix.hh:1610
BCRSMatrix & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition: bcrsmatrix.hh:1549
void setptr(B *_p)
set pointer only
Definition: bvector.hh:1163
Statistics about compression achieved in implicit mode.
Definition: bcrsmatrix.hh:80
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:437
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:817
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: bcrsmatrix.hh:1756
void setIndices(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition: bcrsmatrix.hh:1192
Definition: basearray.hh:19
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1107
RealRowIterator()
empty constructor, use with care!
Definition: bcrsmatrix.hh:544
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition: bcrsmatrix.hh:753
A allocator_type
export the allocator type
Definition: bcrsmatrix.hh:440
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Iterator class for sequential creation of blocks
Definition: bcrsmatrix.hh:917
size_type overflow_total
total number of elements written to the overflow area during construction.
Definition: bcrsmatrix.hh:87
The overflow error used during implicit BCRSMatrix construction was exhausted.
Definition: istlexception.hh:32
void setDataPointers()
Set data pointers for all rows.
Definition: bcrsmatrix.hh:2020
Build in a row-wise manner.
Definition: bcrsmatrix.hh:469
M_ Matrix
The underlying matrix.
Definition: bcrsmatrix.hh:115
ImplicitMatrixBuilder(Matrix &m)
Creates an ImplicitMatrixBuilder for matrix m.
Definition: bcrsmatrix.hh:163
RealRowIterator(const RealRowIterator< ValueType > &it)
Definition: bcrsmatrix.hh:548
size_type N() const
The number of rows in the matrix.
Definition: bcrsmatrix.hh:204
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: bcrsmatrix.hh:1924
A::template rebind< size_type >::other sizeAllocator_
Definition: bcrsmatrix.hh:1944
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: bcrsmatrix.hh:1682
Iterator beforeEnd()
Definition: bcrsmatrix.hh:636
void implicit_allocate(size_type _n, size_type _m)
organizes allocation implicit mode calculates correct array size to be allocated and sets the the win...
Definition: bcrsmatrix.hh:2191
void umv(const X &x, Y &y) const
y += A x
Definition: bcrsmatrix.hh:1591
BCRSMatrix(size_type _n, size_type _m, size_type _avg, double _overflowsize, BuildMode bm)
construct matrix with a known average number of entries per row
Definition: bcrsmatrix.hh:732
size_type index() const
return index
Definition: bcrsmatrix.hh:554
RealRowIterator(row_type *_p, size_type _i)
constructor
Definition: bcrsmatrix.hh:539
Build entries randomly.
Definition: bcrsmatrix.hh:478
std::shared_ptr< size_type > j_
Definition: bcrsmatrix.hh:1960
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1892
row_type * r
Definition: bcrsmatrix.hh:1954
double overflowsize
Definition: bcrsmatrix.hh:1964
size_type index() const
The number of the row that the iterator currently points to.
Definition: bcrsmatrix.hh:1016
size_type maximum
maximum number of non-zeroes per row.
Definition: bcrsmatrix.hh:85
RealRowIterator< const row_type > const_iterator
The const iterator over the matrix rows.
Definition: bcrsmatrix.hh:655
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition: bcrsmatrix.hh:686
ConstIterator beforeEnd() const
Definition: bcrsmatrix.hh:673
row_object operator[](size_type i) const
Returns a proxy for entries in row i.
Definition: bcrsmatrix.hh:198
size_type size() const
Get the current row size.
Definition: bcrsmatrix.hh:1037
void setColumnPointers(ConstRowIterator row)
Copy row sizes from iterator range starting at row and set column index pointers for all rows...
Definition: bcrsmatrix.hh:1994
void deallocate(bool deallocateRows=true)
deallocate memory of the matrix.
Definition: bcrsmatrix.hh:2056
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:652
RealRowIterator< row_type > Iterator
Definition: bcrsmatrix.hh:620
Proxy row object for entry access.
Definition: bcrsmatrix.hh:129
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1055
A::template rebind< row_type >::other rowAllocator_
Definition: bcrsmatrix.hh:1942
Deallocator(typename A::template rebind< size_type >::other &sizeAllocator)
Definition: bcrsmatrix.hh:2111
RealRowIterator< const row_type > ConstIterator
Definition: bcrsmatrix.hh:656
std::ptrdiff_t distanceTo(const RealRowIterator< const ValueType > &other) const
Definition: bcrsmatrix.hh:565
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1898
Some handy generic functions for ISTL matrices.
void umtv(const X &x, Y &y) const
y += A^T x
Definition: bcrsmatrix.hh:1663
A::template rebind< B >::other allocator_
Definition: bcrsmatrix.hh:1940
iterator class for sequential access
Definition: basearray.hh:584
size_type getrowsize(size_type i) const
get current number of indices in row i
Definition: bcrsmatrix.hh:1086
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: bcrsmatrix.hh:1777
void incrementrowsize(size_type i, size_type s=1)
increment size of row i by s (1 by default)
Definition: bcrsmatrix.hh:1096
Iterator access to matrix rows
Definition: bcrsmatrix.hh:525
A wrapper for uniform access to the BCRSMatrix during and after the build stage in implicit build mod...
Definition: bcrsmatrix.hh:109
void allocateData()
Definition: bcrsmatrix.hh:2172
CreateIterator & operator++()
prefix increment
Definition: bcrsmatrix.hh:935
bool equals(const RealRowIterator< const ValueType > &other) const
equality
Definition: bcrsmatrix.hh:579
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: bcrsmatrix.hh:1737
size_type nnz_
Definition: bcrsmatrix.hh:1949
BuildMode buildMode() const
The currently selected build mode of the matrix.
Definition: bcrsmatrix.hh:1916
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:689
void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
Allocate memory for the matrix structure.
Definition: bcrsmatrix.hh:2136
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1149
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition: bcrsmatrix.hh:1317