dune-istl  2.8.0
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 #include <memory>
16 
17 #include "istlexception.hh"
18 #include "bvector.hh"
19 #include "matrixutils.hh"
20 #include <dune/common/stdstreams.hh>
21 #include <dune/common/iteratorfacades.hh>
22 #include <dune/common/typetraits.hh>
23 #include <dune/common/ftraits.hh>
24 #include <dune/common/scalarvectorview.hh>
25 #include <dune/common/scalarmatrixview.hh>
26 
27 #include <dune/istl/blocklevel.hh>
28 
33 namespace Dune {
34 
74  template<typename M>
75  struct MatrixDimension;
76 
78 
84  template<typename size_type>
86  {
88  double avg;
90  size_type maximum;
92  size_type overflow_total;
94 
97  double mem_ratio;
98  };
99 
101 
113  template<class M_>
115  {
116 
117  public:
118 
120  typedef M_ Matrix;
121 
123  typedef typename Matrix::block_type block_type;
124 
126  typedef typename Matrix::size_type size_type;
127 
129 
135  {
136 
137  public:
138 
141  {
142  return _m.entry(_i,j);
143  }
144 
145 #ifndef DOXYGEN
146 
148  : _m(m)
149  , _i(i)
150  {}
151 
152 #endif
153 
154  private:
155 
156  Matrix& _m;
158 
159  };
160 
162 
169  : _m(m)
170  {
171  if (m.buildMode() != Matrix::implicit)
172  DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix in implicit build mode");
173  if (m.buildStage() != Matrix::building)
174  DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix with set size that has not been compressed() yet");
175  }
176 
178 
192  ImplicitMatrixBuilder(Matrix& m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
193  : _m(m)
194  {
195  if (m.buildStage() != Matrix::notAllocated)
196  DUNE_THROW(BCRSMatrixError,"You can only set up a matrix for this ImplicitBuilder if it has no memory allocated yet");
197  m.setBuildMode(Matrix::implicit);
198  m.setImplicitBuildModeParameters(avg_cols_per_row,overflow_fraction);
199  m.setSize(rows,cols);
200  }
201 
204  {
205  return row_object(_m,i);
206  }
207 
209  size_type N() const
210  {
211  return _m.N();
212  }
213 
215  size_type M() const
216  {
217  return _m.M();
218  }
219 
220  private:
221 
222  Matrix& _m;
223 
224  };
225 
462  template<class B, class A=std::allocator<B> >
464  {
465  friend struct MatrixDimension<BCRSMatrix>;
466  public:
467  enum BuildStage {
480  built=3
481  };
482 
483  //===== type definitions and constants
484 
486  using field_type = typename Imp::BlockTraits<B>::field_type;
487 
489  typedef B block_type;
490 
492  typedef A allocator_type;
493 
495  typedef Imp::CompressedBlockVectorWindow<B,A> row_type;
496 
498  typedef typename A::size_type size_type;
499 
501  typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
502 
504  [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
505  static constexpr unsigned int blocklevel = blockLevel<B>()+1;
506 
508  enum BuildMode {
541  unknown
542  };
543 
544  //===== random access interface to rows of the matrix
545 
548  {
549 #ifdef DUNE_ISTL_WITH_CHECKING
550  if (build_mode == implicit && ready != built)
551  DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
552  if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
553  if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
554 #endif
555  return r[i];
556  }
557 
560  {
561 #ifdef DUNE_ISTL_WITH_CHECKING
562  if (build_mode == implicit && ready != built)
563  DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
564  if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
565  if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
566 #endif
567  return r[i];
568  }
569 
570 
571  //===== iterator interface to rows of the matrix
572 
574  template<class T>
576  : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
577  {
578 
579  public:
581  typedef typename std::remove_const<T>::type ValueType;
582 
583  friend class RandomAccessIteratorFacade<RealRowIterator<const ValueType>, const ValueType>;
584  friend class RandomAccessIteratorFacade<RealRowIterator<ValueType>, ValueType>;
585  friend class RealRowIterator<const ValueType>;
586  friend class RealRowIterator<ValueType>;
587 
590  : p(_p), i(_i)
591  {}
592 
595  : p(0), i(0)
596  {}
597 
599  : p(it.p), i(it.i)
600  {}
601 
602 
604  size_type index () const
605  {
606  return i;
607  }
608 
609  std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
610  {
611  assert(other.p==p);
612  return (other.i-i);
613  }
614 
615  std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
616  {
617  assert(other.p==p);
618  return (other.i-i);
619  }
620 
622  bool equals (const RealRowIterator<ValueType>& other) const
623  {
624  assert(other.p==p);
625  return i==other.i;
626  }
627 
629  bool equals (const RealRowIterator<const ValueType>& other) const
630  {
631  assert(other.p==p);
632  return i==other.i;
633  }
634 
635  private:
637  void increment()
638  {
639  ++i;
640  }
641 
643  void decrement()
644  {
645  --i;
646  }
647 
648  void advance(std::ptrdiff_t diff)
649  {
650  i+=diff;
651  }
652 
653  T& elementAt(std::ptrdiff_t diff) const
654  {
655  return p[i+diff];
656  }
657 
659  row_type& dereference () const
660  {
661  return p[i];
662  }
663 
664  row_type* p;
665  size_type i;
666  };
667 
671 
674  {
675  return Iterator(r,0);
676  }
677 
680  {
681  return Iterator(r,n);
682  }
683 
687  {
688  return Iterator(r,n-1);
689  }
690 
694  {
695  return Iterator(r,-1);
696  }
697 
700 
702  typedef typename row_type::Iterator ColIterator;
703 
707 
708 
711  {
712  return ConstIterator(r,0);
713  }
714 
717  {
718  return ConstIterator(r,n);
719  }
720 
724  {
725  return ConstIterator(r,n-1);
726  }
727 
731  {
732  return ConstIterator(r,-1);
733  }
734 
737 
739  typedef typename row_type::ConstIterator ConstColIterator;
740 
741  //===== constructors & resizers
742 
743  // we use a negative compressionBufferSize to indicate that the implicit
744  // mode parameters have not been set yet
745 
748  : build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
749  allocationSize_(0), r(0), a(0),
750  avg(0), compressionBufferSize_(-1.0)
751  {}
752 
755  : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
756  allocationSize_(0), r(0), a(0),
757  avg(0), compressionBufferSize_(-1.0)
758  {
759  allocate(_n, _m, _nnz,true,false);
760  }
761 
764  : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
765  allocationSize_(0), r(0), a(0),
766  avg(0), compressionBufferSize_(-1.0)
767  {
768  allocate(_n, _m,0,true,false);
769  }
770 
772 
782  BCRSMatrix (size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
783  : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
784  allocationSize_(0), r(0), a(0),
785  avg(_avg), compressionBufferSize_(compressionBufferSize)
786  {
787  if (bm != implicit)
788  DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
789  // Prevent user from setting a negative compression buffer size:
790  // 1) It doesn't make sense
791  // 2) We use a negative value to indicate that the parameters
792  // have not been set yet
793  if (compressionBufferSize_ < 0.0)
794  DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
795  implicit_allocate(_n,_m);
796  }
797 
803  BCRSMatrix (const BCRSMatrix& Mat)
804  : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
805  allocationSize_(0), r(0), a(0),
807  {
808  if (!(Mat.ready == notAllocated || Mat.ready == built))
809  DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
810 
811  // deep copy in global array
812  size_type _nnz = Mat.nonzeroes();
813 
814  j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
815  allocate(Mat.n, Mat.m, _nnz, true, true);
816 
817  // build window structure
818  copyWindowStructure(Mat);
819  }
820 
823  {
824  deallocate();
825  }
826 
832  {
833  if (ready == notAllocated)
834  {
835  build_mode = bm;
836  return;
837  }
838  if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
839  build_mode = bm;
840  else
841  DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
842  }
843 
859  void setSize(size_type rows, size_type columns, size_type nnz=0)
860  {
861  // deallocate already setup memory
862  deallocate();
863 
864  if (build_mode == implicit)
865  {
866  if (nnz>0)
867  DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");
868 
869  // implicit allocates differently
870  implicit_allocate(rows,columns);
871  }
872  else
873  {
874  // allocate matrix memory
875  allocate(rows, columns, nnz, true, false);
876  }
877  }
878 
887  void setImplicitBuildModeParameters(size_type _avg, double compressionBufferSize)
888  {
889  // Prevent user from setting a negative compression buffer size:
890  // 1) It doesn't make sense
891  // 2) We use a negative value to indicate that the parameters
892  // have not been set yet
893  if (compressionBufferSize < 0.0)
894  DUNE_THROW(BCRSMatrixError,"You cannot set a negative compressionBufferSize value");
895 
896  // make sure the parameters aren't changed after memory has been allocated
897  if (ready != notAllocated)
898  DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
899  avg = _avg;
900  compressionBufferSize_ = compressionBufferSize;
901  }
902 
910  {
911  // return immediately when self-assignment
912  if (&Mat==this) return *this;
913 
914  if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
915  DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");
916 
917  // make it simple: ALWAYS throw away memory for a and j_
918  // and deallocate rows only if n != Mat.n
919  deallocate(n!=Mat.n);
920 
921  // reallocate the rows if required
922  if (n>0 && n!=Mat.n) {
923  // free rows
924  for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
925  std::allocator_traits<decltype(rowAllocator_)>::destroy(rowAllocator_, riter);
926  rowAllocator_.deallocate(r,n);
927  }
928 
929  nnz_ = Mat.nonzeroes();
930 
931  // allocate a, share j_
932  j_ = Mat.j_;
933  allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);
934 
935  // build window structure
936  copyWindowStructure(Mat);
937  return *this;
938  }
939 
942  {
943 
944  if (!(ready == notAllocated || ready == built))
945  DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");
946 
947  for (size_type i=0; i<n; i++) r[i] = k;
948  return *this;
949  }
950 
951  //===== row-wise creation interface
952 
955  {
956  public:
959  : Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
960  {
961  if (Mat.build_mode == unknown && Mat.ready == building)
962  {
963  Mat.build_mode = row_wise;
964  }
965  if (i==0 && Mat.ready != building)
966  DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
967  if(Mat.build_mode!=row_wise)
968  DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
969  if(i==0 && _Mat.N()==0)
970  // empty Matrix is always built.
971  Mat.ready = built;
972  }
973 
976  {
977  // this should only be called if matrix is in creation
978  if (Mat.ready != building)
979  DUNE_THROW(BCRSMatrixError,"matrix already built up");
980 
981  // row i is defined through the pattern
982  // get memory for the row and initialize the j_ array
983  // this depends on the allocation mode
984 
985  // compute size of the row
986  size_type s = pattern.size();
987 
988  if(s>0) {
989  // update number of nonzeroes including this row
990  nnz += s;
991 
992  // alloc memory / set window
993  if (Mat.nnz_ > 0)
994  {
995  // memory is allocated in one long array
996 
997  // check if that memory is sufficient
998  if (nnz > Mat.nnz_)
999  DUNE_THROW(BCRSMatrixError,"allocated nnz too small");
1000 
1001  // set row i
1002  Mat.r[i].set(s,nullptr,current_row.getindexptr());
1003  current_row.setindexptr(current_row.getindexptr()+s);
1004  }else{
1005  // memory is allocated individually per row
1006  // allocate and set row i
1007  B* b = Mat.allocator_.allocate(s);
1008  // use placement new to call constructor that allocates
1009  // additional memory.
1010  new (b) B[s];
1011  size_type* j = Mat.sizeAllocator_.allocate(s);
1012  Mat.r[i].set(s,b,j);
1013  }
1014  }else
1015  // setup empty row
1016  Mat.r[i].set(0,nullptr,nullptr);
1017 
1018  // initialize the j array for row i from pattern
1019  std::copy(pattern.cbegin(), pattern.cend(), Mat.r[i].getindexptr());
1020 
1021  // now go to next row
1022  i++;
1023  pattern.clear();
1024 
1025  // check if this was last row
1026  if (i==Mat.n)
1027  {
1028  Mat.ready = built;
1029  if(Mat.nnz_ > 0)
1030  {
1031  // Set nnz to the exact number of nonzero blocks inserted
1032  // as some methods rely on it
1033  Mat.nnz_ = nnz;
1034  // allocate data array
1035  Mat.allocateData();
1036  Mat.setDataPointers();
1037  }
1038  }
1039  // done
1040  return *this;
1041  }
1042 
1044  bool operator!= (const CreateIterator& it) const
1045  {
1046  return (i!=it.i) || (&Mat!=&it.Mat);
1047  }
1048 
1050  bool operator== (const CreateIterator& it) const
1051  {
1052  return (i==it.i) && (&Mat==&it.Mat);
1053  }
1054 
1056  size_type index () const
1057  {
1058  return i;
1059  }
1060 
1063  {
1064  pattern.insert(j);
1065  }
1066 
1069  {
1070  return pattern.find(j) != pattern.end();
1071  }
1077  size_type size() const
1078  {
1079  return pattern.size();
1080  }
1081 
1082  private:
1083  BCRSMatrix& Mat; // the matrix we are defining
1084  size_type i; // current row to be defined
1085  size_type nnz; // count total number of nonzeros
1086  typedef std::set<size_type,std::less<size_type> > PatternType;
1087  PatternType pattern; // used to compile entries in a row
1088  row_type current_row; // row pointing to the current row to setup
1089  };
1090 
1092  friend class CreateIterator;
1093 
1096  {
1097  return CreateIterator(*this,0);
1098  }
1099 
1102  {
1103  return CreateIterator(*this,n);
1104  }
1105 
1106 
1107  //===== random creation interface
1108 
1116  {
1117  if (build_mode!=random)
1118  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1119  if (ready != building)
1120  DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1121 
1122  r[i].setsize(s);
1123  }
1124 
1127  {
1128 #ifdef DUNE_ISTL_WITH_CHECKING
1129  if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
1130  if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
1131 #endif
1132  return r[i].getsize();
1133  }
1134 
1137  {
1138  if (build_mode!=random)
1139  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1140  if (ready != building)
1141  DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1142 
1143  r[i].setsize(r[i].getsize()+s);
1144  }
1145 
1147  void endrowsizes ()
1148  {
1149  if (build_mode!=random)
1150  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1151  if (ready != building)
1152  DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1153 
1154  // compute total size, check positivity
1155  size_type total=0;
1156  for (size_type i=0; i<n; i++)
1157  {
1158  total += r[i].getsize();
1159  }
1160 
1161  if(nnz_ == 0)
1162  // allocate/check memory
1163  allocate(n,m,total,false,false);
1164  else if(nnz_ < total)
1165  DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
1166  <<"sufficient for calculated nonzeros ("<<total<<"! ");
1167 
1168  // set the window pointers correctly
1170 
1171  // initialize j_ array with m (an invalid column index)
1172  // this indicates an unused entry
1173  for (size_type k=0; k<nnz_; k++)
1174  j_.get()[k] = m;
1175  ready = rowSizesBuilt;
1176  }
1177 
1179 
1190  {
1191  if (build_mode!=random)
1192  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1193  if (ready==built)
1194  DUNE_THROW(BCRSMatrixError,"matrix already built up");
1195  if (ready==building)
1196  DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
1197  if (ready==notAllocated)
1198  DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1199 
1200  if (col >= m)
1201  DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1202 
1203  // get row range
1204  size_type* const first = r[row].getindexptr();
1205  size_type* const last = first + r[row].getsize();
1206 
1207  // find correct insertion position for new column index
1208  size_type* pos = std::lower_bound(first,last,col);
1209 
1210  // check if index is already in row
1211  if (pos!=last && *pos == col) return;
1212 
1213  // find end of already inserted column indices
1214  size_type* end = std::lower_bound(pos,last,m);
1215  if (end==last)
1216  DUNE_THROW(BCRSMatrixError,"row is too small");
1217 
1218  // insert new column index at correct position
1219  std::copy_backward(pos,end,end+1);
1220  *pos = col;
1221  }
1222 
1224 
1231  template<typename It>
1232  void setIndices(size_type row, It begin, It end)
1233  {
1234  size_type row_size = r[row].size();
1235  size_type* col_begin = r[row].getindexptr();
1236  size_type* col_end;
1237  // consistency check between allocated row size and number of passed column indices
1238  if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
1239  DUNE_THROW(BCRSMatrixError,"Given size of row " << row
1240  << " (" << row_size
1241  << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
1242  std::sort(col_begin,col_end);
1243  }
1244 
1246  void endindices ()
1247  {
1248  if (build_mode!=random)
1249  DUNE_THROW(BCRSMatrixError,"requires random build mode");
1250  if (ready==built)
1251  DUNE_THROW(BCRSMatrixError,"matrix already built up");
1252  if (ready==building)
1253  DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
1254  if (ready==notAllocated)
1255  DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1256 
1257  // check if there are undefined indices
1258  RowIterator endi=end();
1259  for (RowIterator i=begin(); i!=endi; ++i)
1260  {
1261  ColIterator endj = (*i).end();
1262  for (ColIterator j=(*i).begin(); j!=endj; ++j) {
1263  if (j.index() >= m) {
1264  dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
1265  <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
1266  nnz_ -= ((*i).end().offset() - j.offset());
1267  r[i.index()].setsize(j.offset());
1268  break;
1269  }
1270  }
1271  }
1272 
1273  allocateData();
1274  setDataPointers();
1275 
1276  // if not, set matrix to built
1277  ready = built;
1278  }
1279 
1280  //===== implicit creation interface
1281 
1283 
1295  {
1296 #ifdef DUNE_ISTL_WITH_CHECKING
1297  if (build_mode!=implicit)
1298  DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1299  if (ready==built)
1300  DUNE_THROW(BCRSMatrixError,"matrix already built up, use operator[] for entry access now");
1301  if (ready==notAllocated)
1302  DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1303  if (ready!=building)
1304  DUNE_THROW(InvalidStateException,"You may only use entry() during the 'building' stage");
1305 
1306  if (row >= n)
1307  DUNE_THROW(BCRSMatrixError,"row index exceeds matrix size");
1308  if (col >= m)
1309  DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1310 #endif
1311 
1312  size_type* begin = r[row].getindexptr();
1313  size_type* end = begin + r[row].getsize();
1314 
1315  size_type* pos = std::find(begin, end, col);
1316 
1317  //treat the case that there was a match in the array
1318  if (pos != end)
1319  if (*pos == col)
1320  {
1321  std::ptrdiff_t offset = pos - r[row].getindexptr();
1322  B* aptr = r[row].getptr() + offset;
1323 
1324  return *aptr;
1325  }
1326 
1327  //determine whether overflow has to be taken into account or not
1328  if (r[row].getsize() == avg)
1329  return overflow[std::make_pair(row,col)];
1330  else
1331  {
1332  //modify index array
1333  *end = col;
1334 
1335  //do simultaneous operations on data array a
1336  std::ptrdiff_t offset = end - r[row].getindexptr();
1337  B* apos = r[row].getptr() + offset;
1338 
1339  //increase rowsize
1340  r[row].setsize(r[row].getsize()+1);
1341 
1342  //return reference to the newly created entry
1343  return *apos;
1344  }
1345  }
1346 
1348 
1359  {
1360  if (build_mode!=implicit)
1361  DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1362  if (ready==built)
1363  DUNE_THROW(BCRSMatrixError,"matrix already built up, no more need for compression");
1364  if (ready==notAllocated)
1365  DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1366  if (ready!=building)
1367  DUNE_THROW(InvalidStateException,"You may only call compress() at the end of the 'building' stage");
1368 
1369  //calculate statistics
1370  CompressionStatistics stats;
1371  stats.overflow_total = overflow.size();
1372  stats.maximum = 0;
1373 
1374  //get insertion iterators pointing to one before start (for later use of ++it)
1375  size_type* jiit = j_.get();
1376  B* aiit = a;
1377 
1378  //get iterator to the smallest overflow element
1379  typename OverflowType::iterator oit = overflow.begin();
1380 
1381  //store a copy of index pointers on which to perform sorting
1382  std::vector<size_type*> perm;
1383 
1384  //iterate over all rows and copy elements into their position in the compressed array
1385  for (size_type i=0; i<n; i++)
1386  {
1387  //get old pointers into a and j and size without overflow changes
1388  size_type* begin = r[i].getindexptr();
1389  //B* apos = r[i].getptr();
1390  size_type size = r[i].getsize();
1391 
1392  perm.resize(size);
1393 
1394  typename std::vector<size_type*>::iterator it = perm.begin();
1395  for (size_type* iit = begin; iit < begin + size; ++iit, ++it)
1396  *it = iit;
1397 
1398  //sort permutation array
1399  std::sort(perm.begin(),perm.end(),PointerCompare<size_type>());
1400 
1401  //change row window pointer to their new positions
1402  r[i].setindexptr(jiit);
1403  r[i].setptr(aiit);
1404 
1405  for (it = perm.begin(); it != perm.end(); ++it)
1406  {
1407  //check whether there are elements in the overflow area which take precedence
1408  while ((oit!=overflow.end()) && (oit->first < std::make_pair(i,**it)))
1409  {
1410  //check whether there is enough memory to write to
1411  if (jiit > begin)
1413  "Allocated memory for BCRSMatrix exhausted during compress()!"
1414  "Please increase either the average number of entries per row or the compressionBufferSize value."
1415  );
1416  //copy an element from the overflow area to the insertion position in a and j
1417  *jiit = oit->first.second;
1418  ++jiit;
1419  *aiit = oit->second;
1420  ++aiit;
1421  ++oit;
1422  r[i].setsize(r[i].getsize()+1);
1423  }
1424 
1425  //check whether there is enough memory to write to
1426  if (jiit > begin)
1428  "Allocated memory for BCRSMatrix exhausted during compress()!"
1429  "Please increase either the average number of entries per row or the compressionBufferSize value."
1430  );
1431 
1432  //copy element from array
1433  *jiit = **it;
1434  ++jiit;
1435  B* apos = *it - j_.get() + a;
1436  *aiit = *apos;
1437  ++aiit;
1438  }
1439 
1440  //copy remaining elements from the overflow area
1441  while ((oit!=overflow.end()) && (oit->first.first == i))
1442  {
1443  //check whether there is enough memory to write to
1444  if (jiit > begin)
1446  "Allocated memory for BCRSMatrix exhausted during compress()!"
1447  "Please increase either the average number of entries per row or the compressionBufferSize value."
1448  );
1449 
1450  //copy and element from the overflow area to the insertion position in a and j
1451  *jiit = oit->first.second;
1452  ++jiit;
1453  *aiit = oit->second;
1454  ++aiit;
1455  ++oit;
1456  r[i].setsize(r[i].getsize()+1);
1457  }
1458 
1459  // update maximum row size
1460  if (r[i].getsize()>stats.maximum)
1461  stats.maximum = r[i].getsize();
1462  }
1463 
1464  // overflow area may be cleared
1465  overflow.clear();
1466 
1467  //determine average number of entries and memory usage
1468  std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
1469  nnz_ = diff;
1470  stats.avg = (double) (nnz_) / (double) n;
1471  stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
1472 
1473  //matrix is now built
1474  ready = built;
1475 
1476  return stats;
1477  }
1478 
1479  //===== vector space arithmetic
1480 
1483  {
1484 #ifdef DUNE_ISTL_WITH_CHECKING
1485  if (ready != built)
1486  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1487 #endif
1488 
1489  if (nnz_ > 0)
1490  {
1491  // process 1D array
1492  for (size_type i=0; i<nnz_; i++)
1493  a[i] *= k;
1494  }
1495  else
1496  {
1497  RowIterator endi=end();
1498  for (RowIterator i=begin(); i!=endi; ++i)
1499  {
1500  ColIterator endj = (*i).end();
1501  for (ColIterator j=(*i).begin(); j!=endj; ++j)
1502  (*j) *= k;
1503  }
1504  }
1505 
1506  return *this;
1507  }
1508 
1511  {
1512 #ifdef DUNE_ISTL_WITH_CHECKING
1513  if (ready != built)
1514  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1515 #endif
1516 
1517  if (nnz_ > 0)
1518  {
1519  // process 1D array
1520  for (size_type i=0; i<nnz_; i++)
1521  a[i] /= k;
1522  }
1523  else
1524  {
1525  RowIterator endi=end();
1526  for (RowIterator i=begin(); i!=endi; ++i)
1527  {
1528  ColIterator endj = (*i).end();
1529  for (ColIterator j=(*i).begin(); j!=endj; ++j)
1530  (*j) /= k;
1531  }
1532  }
1533 
1534  return *this;
1535  }
1536 
1537 
1544  {
1545 #ifdef DUNE_ISTL_WITH_CHECKING
1546  if (ready != built || b.ready != built)
1547  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1548  if(N()!=b.N() || M() != b.M())
1549  DUNE_THROW(RangeError, "Matrix sizes do not match!");
1550 #endif
1551  RowIterator endi=end();
1552  ConstRowIterator j=b.begin();
1553  for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1554  i->operator+=(*j);
1555  }
1556 
1557  return *this;
1558  }
1559 
1566  {
1567 #ifdef DUNE_ISTL_WITH_CHECKING
1568  if (ready != built || b.ready != built)
1569  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1570  if(N()!=b.N() || M() != b.M())
1571  DUNE_THROW(RangeError, "Matrix sizes do not match!");
1572 #endif
1573  RowIterator endi=end();
1574  ConstRowIterator j=b.begin();
1575  for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1576  i->operator-=(*j);
1577  }
1578 
1579  return *this;
1580  }
1581 
1591  {
1592 #ifdef DUNE_ISTL_WITH_CHECKING
1593  if (ready != built || b.ready != built)
1594  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1595  if(N()!=b.N() || M() != b.M())
1596  DUNE_THROW(RangeError, "Matrix sizes do not match!");
1597 #endif
1598  RowIterator endi=end();
1599  ConstRowIterator j=b.begin();
1600  for(RowIterator i=begin(); i!=endi; ++i, ++j)
1601  i->axpy(alpha, *j);
1602 
1603  return *this;
1604  }
1605 
1606  //===== linear maps
1607 
1609  template<class X, class Y>
1610  void mv (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,
1616  "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
1617  if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,
1618  "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
1619 #endif
1620  ConstRowIterator endi=end();
1621  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1622  {
1623  y[i.index()]=0;
1624  ConstColIterator endj = (*i).end();
1625  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1626  {
1627  auto&& xj = Impl::asVector(x[j.index()]);
1628  auto&& yi = Impl::asVector(y[i.index()]);
1629  Impl::asMatrix(*j).umv(xj, yi);
1630  }
1631  }
1632  }
1633 
1635  template<class X, class Y>
1636  void umv (const X& x, Y& y) const
1637  {
1638 #ifdef DUNE_ISTL_WITH_CHECKING
1639  if (ready != built)
1640  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1641  if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1642  if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1643 #endif
1644  ConstRowIterator endi=end();
1645  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1646  {
1647  ConstColIterator endj = (*i).end();
1648  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1649  {
1650  auto&& xj = Impl::asVector(x[j.index()]);
1651  auto&& yi = Impl::asVector(y[i.index()]);
1652  Impl::asMatrix(*j).umv(xj,yi);
1653  }
1654  }
1655  }
1656 
1658  template<class X, class Y>
1659  void mmv (const X& x, Y& y) const
1660  {
1661 #ifdef DUNE_ISTL_WITH_CHECKING
1662  if (ready != built)
1663  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1664  if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1665  if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1666 #endif
1667  ConstRowIterator endi=end();
1668  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1669  {
1670  ConstColIterator endj = (*i).end();
1671  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1672  {
1673  auto&& xj = Impl::asVector(x[j.index()]);
1674  auto&& yi = Impl::asVector(y[i.index()]);
1675  Impl::asMatrix(*j).mmv(xj,yi);
1676  }
1677  }
1678  }
1679 
1681  template<class X, class Y, class F>
1682  void usmv (F&& alpha, const X& x, Y& y) const
1683  {
1684 #ifdef DUNE_ISTL_WITH_CHECKING
1685  if (ready != built)
1686  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1687  if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1688  if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1689 #endif
1690  ConstRowIterator endi=end();
1691  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1692  {
1693  ConstColIterator endj = (*i).end();
1694  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1695  {
1696  auto&& xj = Impl::asVector(x[j.index()]);
1697  auto&& yi = Impl::asVector(y[i.index()]);
1698  Impl::asMatrix(*j).usmv(alpha,xj,yi);
1699  }
1700  }
1701  }
1702 
1704  template<class X, class Y>
1705  void mtv (const X& x, Y& y) const
1706  {
1707 #ifdef DUNE_ISTL_WITH_CHECKING
1708  if (ready != built)
1709  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1710  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1711  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1712 #endif
1713  for(size_type i=0; i<y.N(); ++i)
1714  y[i]=0;
1715  umtv(x,y);
1716  }
1717 
1719  template<class X, class Y>
1720  void umtv (const X& x, Y& y) const
1721  {
1722 #ifdef DUNE_ISTL_WITH_CHECKING
1723  if (ready != built)
1724  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1725  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1726  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1727 #endif
1728  ConstRowIterator endi=end();
1729  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1730  {
1731  ConstColIterator endj = (*i).end();
1732  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1733  {
1734  auto&& xi = Impl::asVector(x[i.index()]);
1735  auto&& yj = Impl::asVector(y[j.index()]);
1736  Impl::asMatrix(*j).umtv(xi,yj);
1737  }
1738  }
1739  }
1740 
1742  template<class X, class Y>
1743  void mmtv (const X& x, Y& y) const
1744  {
1745 #ifdef DUNE_ISTL_WITH_CHECKING
1746  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1747  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1748 #endif
1749  ConstRowIterator endi=end();
1750  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1751  {
1752  ConstColIterator endj = (*i).end();
1753  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1754  {
1755  auto&& xi = Impl::asVector(x[i.index()]);
1756  auto&& yj = Impl::asVector(y[j.index()]);
1757  Impl::asMatrix(*j).mmtv(xi,yj);
1758  }
1759  }
1760  }
1761 
1763  template<class X, class Y>
1764  void usmtv (const field_type& alpha, const X& x, Y& y) const
1765  {
1766 #ifdef DUNE_ISTL_WITH_CHECKING
1767  if (ready != built)
1768  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1769  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1770  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1771 #endif
1772  ConstRowIterator endi=end();
1773  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1774  {
1775  ConstColIterator endj = (*i).end();
1776  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1777  {
1778  auto&& xi = Impl::asVector(x[i.index()]);
1779  auto&& yj = Impl::asVector(y[j.index()]);
1780  Impl::asMatrix(*j).usmtv(alpha,xi,yj);
1781  }
1782  }
1783  }
1784 
1786  template<class X, class Y>
1787  void umhv (const X& x, Y& y) const
1788  {
1789 #ifdef DUNE_ISTL_WITH_CHECKING
1790  if (ready != built)
1791  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1792  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1793  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1794 #endif
1795  ConstRowIterator endi=end();
1796  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1797  {
1798  ConstColIterator endj = (*i).end();
1799  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1800  {
1801  auto&& xi = Impl::asVector(x[i.index()]);
1802  auto&& yj = Impl::asVector(y[j.index()]);
1803  Impl::asMatrix(*j).umhv(xi,yj);
1804  }
1805  }
1806  }
1807 
1809  template<class X, class Y>
1810  void mmhv (const X& x, Y& y) const
1811  {
1812 #ifdef DUNE_ISTL_WITH_CHECKING
1813  if (ready != built)
1814  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1815  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1816  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1817 #endif
1818  ConstRowIterator endi=end();
1819  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1820  {
1821  ConstColIterator endj = (*i).end();
1822  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1823  {
1824  auto&& xi = Impl::asVector(x[i.index()]);
1825  auto&& yj = Impl::asVector(y[j.index()]);
1826  Impl::asMatrix(*j).mmhv(xi,yj);
1827  }
1828  }
1829  }
1830 
1832  template<class X, class Y>
1833  void usmhv (const field_type& alpha, const X& x, Y& y) const
1834  {
1835 #ifdef DUNE_ISTL_WITH_CHECKING
1836  if (ready != built)
1837  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1838  if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1839  if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1840 #endif
1841  ConstRowIterator endi=end();
1842  for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1843  {
1844  ConstColIterator endj = (*i).end();
1845  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1846  {
1847  auto&& xi = Impl::asVector(x[i.index()]);
1848  auto&& yj = Impl::asVector(y[j.index()]);
1849  Impl::asMatrix(*j).usmhv(alpha,xi,yj);
1850  }
1851  }
1852  }
1853 
1854 
1855  //===== norms
1856 
1858  typename FieldTraits<field_type>::real_type frobenius_norm2 () const
1859  {
1860 #ifdef DUNE_ISTL_WITH_CHECKING
1861  if (ready != built)
1862  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1863 #endif
1864 
1865  typename FieldTraits<field_type>::real_type sum=0;
1866 
1867  for (auto&& row : (*this))
1868  for (auto&& entry : row)
1869  sum += Impl::asMatrix(entry).frobenius_norm2();
1870 
1871  return sum;
1872  }
1873 
1875  typename FieldTraits<field_type>::real_type frobenius_norm () const
1876  {
1877  return sqrt(frobenius_norm2());
1878  }
1879 
1881  template <typename ft = field_type,
1882  typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
1883  typename FieldTraits<ft>::real_type infinity_norm() const {
1884  if (ready != built)
1885  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1886 
1887  using real_type = typename FieldTraits<ft>::real_type;
1888  using std::max;
1889 
1890  real_type norm = 0;
1891  for (auto const &x : *this) {
1892  real_type sum = 0;
1893  for (auto const &y : x)
1894  sum += Impl::asMatrix(y).infinity_norm();
1895  norm = max(sum, norm);
1896  }
1897  return norm;
1898  }
1899 
1901  template <typename ft = field_type,
1902  typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
1903  typename FieldTraits<ft>::real_type infinity_norm_real() const {
1904  if (ready != built)
1905  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1906 
1907  using real_type = typename FieldTraits<ft>::real_type;
1908  using std::max;
1909 
1910  real_type norm = 0;
1911  for (auto const &x : *this) {
1912  real_type sum = 0;
1913  for (auto const &y : x)
1914  sum += Impl::asMatrix(y).infinity_norm_real();
1915  norm = max(sum, norm);
1916  }
1917  return norm;
1918  }
1919 
1921  template <typename ft = field_type,
1922  typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
1923  typename FieldTraits<ft>::real_type infinity_norm() const {
1924  if (ready != built)
1925  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1926 
1927  using real_type = typename FieldTraits<ft>::real_type;
1928  using std::max;
1929 
1930  real_type norm = 0;
1931  real_type isNaN = 1;
1932  for (auto const &x : *this) {
1933  real_type sum = 0;
1934  for (auto const &y : x)
1935  sum += Impl::asMatrix(y).infinity_norm();
1936  norm = max(sum, norm);
1937  isNaN += sum;
1938  }
1939 
1940  return norm * (isNaN / isNaN);
1941  }
1942 
1944  template <typename ft = field_type,
1945  typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
1946  typename FieldTraits<ft>::real_type infinity_norm_real() const {
1947  if (ready != built)
1948  DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1949 
1950  using real_type = typename FieldTraits<ft>::real_type;
1951  using std::max;
1952 
1953  real_type norm = 0;
1954  real_type isNaN = 1;
1955 
1956  for (auto const &x : *this) {
1957  real_type sum = 0;
1958  for (auto const &y : x)
1959  sum += Impl::asMatrix(y).infinity_norm_real();
1960  norm = max(sum, norm);
1961  isNaN += sum;
1962  }
1963 
1964  return norm * (isNaN / isNaN);
1965  }
1966 
1967  //===== sizes
1968 
1970  size_type N () const
1971  {
1972  return n;
1973  }
1974 
1976  size_type M () const
1977  {
1978  return m;
1979  }
1980 
1983  {
1984  // in case of row-wise allocation
1985  if( nnz_ <= 0 )
1986  nnz_ = std::accumulate( begin(), end(), size_type( 0 ), [] ( size_type s, const row_type &row ) { return s+row.getsize(); } );
1987  return nnz_;
1988  }
1989 
1992  {
1993  return ready;
1994  }
1995 
1998  {
1999  return build_mode;
2000  }
2001 
2002  //===== query
2003 
2005  bool exists (size_type i, size_type j) const
2006  {
2007 #ifdef DUNE_ISTL_WITH_CHECKING
2008  if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
2009  if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
2010 #endif
2011  return (r[i].size() && r[i].find(j) != r[i].end());
2012  }
2013 
2014 
2015  protected:
2016  // state information
2017  BuildMode build_mode; // row wise or whole matrix
2018  BuildStage ready; // indicate the stage the matrix building is in
2019 
2020  // The allocator used for memory management
2021  typename std::allocator_traits<A>::template rebind_alloc<B> allocator_;
2022 
2023  typename std::allocator_traits<A>::template rebind_alloc<row_type> rowAllocator_;
2024 
2025  typename std::allocator_traits<A>::template rebind_alloc<size_type> sizeAllocator_;
2026 
2027  // size of the matrix
2028  size_type n; // number of rows
2029  size_type m; // number of columns
2030  mutable size_type nnz_; // number of nonzeroes contained in the matrix
2031  size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
2032  // zero means that memory is allocated separately for each row.
2033 
2034  // the rows are dynamically allocated
2035  row_type* r; // [n] the individual rows having pointers into a,j arrays
2036 
2037  // dynamically allocated memory
2038  B* a; // [allocationSize] non-zero entries of the matrix in row-wise ordering
2039  // If a single array of column indices is used, it can be shared
2040  // between different matrices with the same sparsity pattern
2041  std::shared_ptr<size_type> j_; // [allocationSize] column indices of entries
2042 
2043  // additional data is needed in implicit buildmode
2046 
2047  typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
2049 
2051  {
2052  row_type current_row(a,j_.get(),0); // Pointers to current row data
2053  for (size_type i=0; i<n; i++, ++row) {
2054  // set row i
2055  size_type s = row->getsize();
2056 
2057  if (s>0) {
2058  // setup pointers and size
2059  r[i].set(s,current_row.getptr(), current_row.getindexptr());
2060  // update pointer for next row
2061  current_row.setptr(current_row.getptr()+s);
2062  current_row.setindexptr(current_row.getindexptr()+s);
2063  } else{
2064  // empty row
2065  r[i].set(0,nullptr,nullptr);
2066  }
2067  }
2068  }
2069 
2071 
2076  {
2077  size_type* jptr = j_.get();
2078  for (size_type i=0; i<n; ++i, ++row) {
2079  // set row i
2080  size_type s = row->getsize();
2081 
2082  if (s>0) {
2083  // setup pointers and size
2084  r[i].setsize(s);
2085  r[i].setindexptr(jptr);
2086  } else{
2087  // empty row
2088  r[i].set(0,nullptr,nullptr);
2089  }
2090 
2091  // advance position in global array
2092  jptr += s;
2093  }
2094  }
2095 
2097 
2102  {
2103  B* aptr = a;
2104  for (size_type i=0; i<n; ++i) {
2105  // set row i
2106  if (r[i].getsize() > 0) {
2107  // setup pointers and size
2108  r[i].setptr(aptr);
2109  } else{
2110  // empty row
2111  r[i].set(0,nullptr,nullptr);
2112  }
2113 
2114  // advance position in global array
2115  aptr += r[i].getsize();
2116  }
2117  }
2118 
2121  {
2122  setWindowPointers(Mat.begin());
2123 
2124  // copy data
2125  for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
2126 
2127  // finish off
2128  build_mode = row_wise; // dummy
2129  ready = built;
2130  }
2131 
2137  void deallocate(bool deallocateRows=true)
2138  {
2139 
2140  if (notAllocated)
2141  return;
2142 
2143  if (allocationSize_>0)
2144  {
2145  // a,j_ have been allocated as one long vector
2146  j_.reset();
2147  if (a)
2148  {
2149  for(B *aiter=a+(allocationSize_-1), *aend=a-1; aiter!=aend; --aiter)
2150  std::allocator_traits<decltype(allocator_)>::destroy(allocator_, aiter);
2151  allocator_.deallocate(a,allocationSize_);
2152  a = nullptr;
2153  }
2154  }
2155  else if (r)
2156  {
2157  // check if memory for rows have been allocated individually
2158  for (size_type i=0; i<n; i++)
2159  if (r[i].getsize()>0)
2160  {
2161  for (B *col=r[i].getptr()+(r[i].getsize()-1),
2162  *colend = r[i].getptr()-1; col!=colend; --col) {
2163  std::allocator_traits<decltype(allocator_)>::destroy(allocator_, col);
2164  }
2165  sizeAllocator_.deallocate(r[i].getindexptr(),1);
2166  allocator_.deallocate(r[i].getptr(),1);
2167  // clear out row data in case we don't want to deallocate the rows
2168  // otherwise we might run into a double free problem here later
2169  r[i].set(0,nullptr,nullptr);
2170  }
2171  }
2172 
2173  // deallocate the rows
2174  if (n>0 && deallocateRows && r) {
2175  for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
2176  std::allocator_traits<decltype(rowAllocator_)>::destroy(rowAllocator_, riter);
2177  rowAllocator_.deallocate(r,n);
2178  r = nullptr;
2179  }
2180 
2181  // Mark matrix as not built at all.
2183 
2184  }
2185 
2188  {
2189  typename std::allocator_traits<A>::template rebind_alloc<size_type>& sizeAllocator_;
2190 
2191  public:
2192  Deallocator(typename std::allocator_traits<A>::template rebind_alloc<size_type>& sizeAllocator)
2193  : sizeAllocator_(sizeAllocator)
2194  {}
2195 
2196  void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); }
2197  };
2198 
2199 
2217  void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
2218  {
2219  // Store size
2220  n = rows;
2221  m = columns;
2222  nnz_ = allocationSize;
2223  allocationSize_ = allocationSize;
2224 
2225  // allocate rows
2226  if(allocateRows) {
2227  if (n>0) {
2228  if (r)
2229  DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
2230  r = rowAllocator_.allocate(rows);
2231  // initialize row entries
2232  for(row_type* ri=r; ri!=r+rows; ++ri)
2233  std::allocator_traits<decltype(rowAllocator_)>::construct(rowAllocator_, ri, row_type());
2234  }else{
2235  r = 0;
2236  }
2237  }
2238 
2239  // allocate a and j_ array
2240  if (allocate_data)
2241  allocateData();
2242  if (allocationSize_>0) {
2243  // allocate column indices only if not yet present (enable sharing)
2244  if (!j_.get())
2246  }else{
2247  j_.reset();
2248  }
2249 
2250  // Mark the matrix as not built.
2251  ready = building;
2252  }
2253 
2255  {
2256  if (a)
2257  DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
2258  if (allocationSize_>0) {
2259  a = allocator_.allocate(allocationSize_);
2260  // use placement new to call constructor that allocates
2261  // additional memory.
2262  new (a) B[allocationSize_];
2263  } else {
2264  a = nullptr;
2265  }
2266  }
2267 
2274  {
2275  if (build_mode != implicit)
2276  DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
2277  if (ready != notAllocated)
2278  DUNE_THROW(InvalidStateException,"memory has already been allocated");
2279 
2280  // check to make sure the user has actually set the parameters
2281  if (compressionBufferSize_ < 0)
2282  DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
2283  //calculate size of overflow region, add buffer for row sort!
2284  size_type osize = (size_type) (_n*avg)*compressionBufferSize_ + 4*avg;
2285  allocationSize_ = _n*avg + osize;
2286 
2287  allocate(_n, _m, allocationSize_,true,true);
2288 
2289  //set row pointers correctly
2290  size_type* jptr = j_.get() + osize;
2291  B* aptr = a + osize;
2292  for (size_type i=0; i<n; i++)
2293  {
2294  r[i].set(0,aptr,jptr);
2295  jptr = jptr + avg;
2296  aptr = aptr + avg;
2297  }
2298 
2299  ready = building;
2300  }
2301  };
2302 
2303 
2306 } // end namespace
2307 
2308 #endif
Helper functions for determining the vector/matrix block level.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Some handy generic functions for ISTL matrices.
Col col
Definition: matrixmatrix.hh:349
Definition: allocator.hh:9
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get([[maybe_unused]] const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:291
Definition: matrixutils.hh:209
Statistics about compression achieved in implicit mode.
Definition: bcrsmatrix.hh:86
size_type overflow_total
total number of elements written to the overflow area during construction.
Definition: bcrsmatrix.hh:92
size_type maximum
maximum number of non-zeroes per row.
Definition: bcrsmatrix.hh:90
double avg
average number of non-zeroes per row.
Definition: bcrsmatrix.hh:88
double mem_ratio
fraction of wasted memory resulting from non-used overflow area.
Definition: bcrsmatrix.hh:97
A wrapper for uniform access to the BCRSMatrix during and after the build stage in implicit build mod...
Definition: bcrsmatrix.hh:115
Matrix::block_type block_type
The block_type of the underlying matrix.
Definition: bcrsmatrix.hh:123
ImplicitMatrixBuilder(Matrix &m)
Creates an ImplicitMatrixBuilder for matrix m.
Definition: bcrsmatrix.hh:168
M_ Matrix
The underlying matrix.
Definition: bcrsmatrix.hh:120
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:192
size_type M() const
The number of columns in the matrix.
Definition: bcrsmatrix.hh:215
Matrix::size_type size_type
The size_type of the underlying matrix.
Definition: bcrsmatrix.hh:126
row_object operator[](size_type i) const
Returns a proxy for entries in row i.
Definition: bcrsmatrix.hh:203
size_type N() const
The number of rows in the matrix.
Definition: bcrsmatrix.hh:209
Proxy row object for entry access.
Definition: bcrsmatrix.hh:135
block_type & operator[](size_type j) const
Returns entry in column j.
Definition: bcrsmatrix.hh:140
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: bcrsmatrix.hh:1875
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:486
std::allocator_traits< A >::template rebind_alloc< row_type > rowAllocator_
Definition: bcrsmatrix.hh:2023
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: bcrsmatrix.hh:2005
BuildStage buildStage() const
The current build stage of the matrix.
Definition: bcrsmatrix.hh:1991
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:673
friend class CreateIterator
allow CreateIterator to access internal data
Definition: bcrsmatrix.hh:1092
void copyWindowStructure(const BCRSMatrix &Mat)
Copy the window structure from another matrix.
Definition: bcrsmatrix.hh:2120
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: bcrsmatrix.hh:1810
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: bcrsmatrix.hh:1833
void umtv(const X &x, Y &y) const
y += A^T x
Definition: bcrsmatrix.hh:1720
double compressionBufferSize_
Definition: bcrsmatrix.hh:2045
size_type m
Definition: bcrsmatrix.hh:2029
RealRowIterator< const row_type > const_iterator
The const iterator over the matrix rows.
Definition: bcrsmatrix.hh:705
static constexpr unsigned int blocklevel
increment block level counter
Definition: bcrsmatrix.hh:505
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:2217
~BCRSMatrix()
destructor
Definition: bcrsmatrix.hh:822
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:909
void allocateData()
Definition: bcrsmatrix.hh:2254
void deallocate(bool deallocateRows=true)
deallocate memory of the matrix.
Definition: bcrsmatrix.hh:2137
Iterator RowIterator
rename the iterators for easier access
Definition: bcrsmatrix.hh:699
BCRSMatrix()
an empty matrix
Definition: bcrsmatrix.hh:747
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1147
void incrementrowsize(size_type i, size_type s=1)
increment size of row i by s (1 by default)
Definition: bcrsmatrix.hh:1136
void mtv(const X &x, Y &y) const
y = A^T x
Definition: bcrsmatrix.hh:1705
void umhv(const X &x, Y &y) const
y += A^H x
Definition: bcrsmatrix.hh:1787
size_type nonzeroes() const
number of blocks that are stored (the number of blocks that possibly are nonzero)
Definition: bcrsmatrix.hh:1982
size_type allocationSize_
Definition: bcrsmatrix.hh:2031
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: bcrsmatrix.hh:1903
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition: bcrsmatrix.hh:736
BuildStage ready
Definition: bcrsmatrix.hh:2018
BuildMode build_mode
Definition: bcrsmatrix.hh:2017
void setrowsize(size_type i, size_type s)
Set number of indices in row i to s.
Definition: bcrsmatrix.hh:1115
RealRowIterator< row_type > Iterator
Definition: bcrsmatrix.hh:670
size_type nnz_
Definition: bcrsmatrix.hh:2030
std::allocator_traits< A >::template rebind_alloc< size_type > sizeAllocator_
Definition: bcrsmatrix.hh:2025
RealRowIterator< row_type > iterator
The iterator over the (mutable matrix rows.
Definition: bcrsmatrix.hh:669
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: bcrsmatrix.hh:1764
ConstIterator beforeBegin() const
Definition: bcrsmatrix.hh:730
RealRowIterator< const row_type > ConstIterator
Definition: bcrsmatrix.hh:706
Iterator beforeBegin()
Definition: bcrsmatrix.hh:693
B * a
Definition: bcrsmatrix.hh:2038
BuildMode
we support two modes
Definition: bcrsmatrix.hh:508
@ implicit
Build entries randomly with an educated guess for the number of entries per row.
Definition: bcrsmatrix.hh:537
@ unknown
Build mode not set!
Definition: bcrsmatrix.hh:541
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:528
@ row_wise
Build in a row-wise manner.
Definition: bcrsmatrix.hh:519
BCRSMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: bcrsmatrix.hh:1510
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: bcrsmatrix.hh:1858
BCRSMatrix(size_type _n, size_type _m, size_type _nnz, BuildMode bm)
matrix with known number of nonzeroes
Definition: bcrsmatrix.hh:754
::Dune::CompressionStatistics< size_type > CompressionStatistics
The type for the statistics object returned by compress()
Definition: bcrsmatrix.hh:501
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition: bcrsmatrix.hh:803
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:679
row_type * r
Definition: bcrsmatrix.hh:2035
void setIndices(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition: bcrsmatrix.hh:1232
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1189
std::map< std::pair< size_type, size_type >, B > OverflowType
Definition: bcrsmatrix.hh:2047
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:702
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: bcrsmatrix.hh:1883
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:498
OverflowType overflow
Definition: bcrsmatrix.hh:2048
BCRSMatrix(size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
construct matrix with a known average number of entries per row
Definition: bcrsmatrix.hh:782
row_type & operator[](size_type i)
random access to the rows
Definition: bcrsmatrix.hh:547
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1101
BCRSMatrix & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition: bcrsmatrix.hh:1590
BCRSMatrix & operator+=(const BCRSMatrix &b)
Add the entries of another matrix to this one.
Definition: bcrsmatrix.hh:1543
Iterator beforeEnd()
Definition: bcrsmatrix.hh:686
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:739
void usmv(F &&alpha, const X &x, Y &y) const
y += alpha A x
Definition: bcrsmatrix.hh:1682
size_type getrowsize(size_type i) const
get current number of indices in row i
Definition: bcrsmatrix.hh:1126
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1976
size_type n
Definition: bcrsmatrix.hh:2028
Imp::CompressedBlockVectorWindow< B, A > row_type
implement row_type with compressed vector
Definition: bcrsmatrix.hh:495
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1095
BuildStage
Definition: bcrsmatrix.hh:467
@ rowSizesBuilt
The row sizes of the matrix are known.
Definition: bcrsmatrix.hh:478
@ built
The matrix structure is fully built.
Definition: bcrsmatrix.hh:480
@ notbuilt
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:469
@ notAllocated
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:471
@ building
Matrix is currently being built, some memory has been allocated, build mode and size are fixed.
Definition: bcrsmatrix.hh:473
BuildMode buildMode() const
The currently selected build mode of the matrix.
Definition: bcrsmatrix.hh:1997
void mmv(const X &x, Y &y) const
y -= A x
Definition: bcrsmatrix.hh:1659
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1610
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:489
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: bcrsmatrix.hh:1743
size_type avg
Definition: bcrsmatrix.hh:2044
void umv(const X &x, Y &y) const
y += A x
Definition: bcrsmatrix.hh:1636
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:2273
void setImplicitBuildModeParameters(size_type _avg, double compressionBufferSize)
Set parameters needed for creation in implicit build mode.
Definition: bcrsmatrix.hh:887
BCRSMatrix(size_type _n, size_type _m, BuildMode bm)
matrix with unknown number of nonzeroes
Definition: bcrsmatrix.hh:763
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1246
BCRSMatrix & operator-=(const BCRSMatrix &b)
Subtract the entries of another matrix from this one.
Definition: bcrsmatrix.hh:1565
BCRSMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: bcrsmatrix.hh:1482
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition: bcrsmatrix.hh:1358
B & entry(size_type row, size_type col)
Returns reference to entry (row,col) of the matrix.
Definition: bcrsmatrix.hh:1294
void setDataPointers()
Set data pointers for all rows.
Definition: bcrsmatrix.hh:2101
std::allocator_traits< A >::template rebind_alloc< B > allocator_
Definition: bcrsmatrix.hh:2021
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1970
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:831
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:859
std::shared_ptr< size_type > j_
Definition: bcrsmatrix.hh:2041
void setWindowPointers(ConstRowIterator row)
Definition: bcrsmatrix.hh:2050
void setColumnPointers(ConstRowIterator row)
Copy row sizes from iterator range starting at row and set column index pointers for all rows.
Definition: bcrsmatrix.hh:2075
ConstIterator end() const
Get const iterator to one beyond last row.
Definition: bcrsmatrix.hh:716
ConstIterator begin() const
Get const iterator to first row.
Definition: bcrsmatrix.hh:710
A allocator_type
export the allocator type
Definition: bcrsmatrix.hh:492
ConstIterator beforeEnd() const
Definition: bcrsmatrix.hh:723
Iterator access to matrix rows
Definition: bcrsmatrix.hh:577
RealRowIterator()
empty constructor, use with care!
Definition: bcrsmatrix.hh:594
bool equals(const RealRowIterator< ValueType > &other) const
equality
Definition: bcrsmatrix.hh:622
std::remove_const< T >::type ValueType
The unqualified value type.
Definition: bcrsmatrix.hh:581
RealRowIterator(const RealRowIterator< ValueType > &it)
Definition: bcrsmatrix.hh:598
bool equals(const RealRowIterator< const ValueType > &other) const
equality
Definition: bcrsmatrix.hh:629
RealRowIterator(row_type *_p, size_type _i)
constructor
Definition: bcrsmatrix.hh:589
std::ptrdiff_t distanceTo(const RealRowIterator< const ValueType > &other) const
Definition: bcrsmatrix.hh:615
size_type index() const
return index
Definition: bcrsmatrix.hh:604
std::ptrdiff_t distanceTo(const RealRowIterator< ValueType > &other) const
Definition: bcrsmatrix.hh:609
Iterator class for sequential creation of blocks
Definition: bcrsmatrix.hh:955
bool operator==(const CreateIterator &it) const
equality
Definition: bcrsmatrix.hh:1050
size_type index() const
The number of the row that the iterator currently points to.
Definition: bcrsmatrix.hh:1056
bool operator!=(const CreateIterator &it) const
inequality
Definition: bcrsmatrix.hh:1044
CreateIterator(BCRSMatrix &_Mat, size_type _i)
constructor
Definition: bcrsmatrix.hh:958
CreateIterator & operator++()
prefix increment
Definition: bcrsmatrix.hh:975
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1062
bool contains(size_type j)
return true if column index is in row
Definition: bcrsmatrix.hh:1068
size_type size() const
Get the current row size.
Definition: bcrsmatrix.hh:1077
Class used by shared_ptr to deallocate memory using the proper allocator.
Definition: bcrsmatrix.hh:2188
Deallocator(typename std::allocator_traits< A >::template rebind_alloc< size_type > &sizeAllocator)
Definition: bcrsmatrix.hh:2192
void operator()(size_type *p)
Definition: bcrsmatrix.hh:2196
Error specific to BCRSMatrix.
Definition: istlexception.hh:22
Thrown when the compression buffer used by the implicit BCRSMatrix construction is exhausted.
Definition: istlexception.hh:35
A generic dynamic dense matrix.
Definition: matrix.hh:559
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:575
T block_type
Export the type representing the components.
Definition: matrix.hh:566
Definition: matrixutils.hh:536