Dune Core Modules (2.5.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
16#include "istlexception.hh"
17#include "bvector.hh"
18#include "matrixutils.hh"
23
28namespace Dune {
29
69 template<typename M>
70 struct MatrixDimension;
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
119
121 typedef typename Matrix::size_type size_type;
122
124
130 {
131
132 public:
133
136 {
137 return _m.entry(_i,j);
138 }
139
140#ifndef DOXYGEN
141
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
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
421 template<class B, class A=std::allocator<B> >
423 {
424 friend struct MatrixDimension<BCRSMatrix>;
425 public:
439 built=3
440 };
441
442 //===== type definitions and constants
443
445 typedef typename B::field_type field_type;
446
448 typedef B block_type;
449
451 typedef A allocator_type;
452
455
457 typedef typename A::size_type size_type;
458
460 typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
461
463 enum {
465 blocklevel = B::blocklevel+1
466 };
467
502 unknown
503 };
504
505 //===== random access interface to rows of the matrix
506
509 {
510#ifdef DUNE_ISTL_WITH_CHECKING
511 if (build_mode == implicit && ready != built)
512 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
513 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
514 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
515#endif
516 return r[i];
517 }
518
521 {
522#ifdef DUNE_ISTL_WITH_CHECKING
523 if (build_mode == implicit && ready != built)
524 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
525 if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
526 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
527#endif
528 return r[i];
529 }
530
531
532 //===== iterator interface to rows of the matrix
533
535 template<class T>
537 : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
538 {
539
540 public:
542 typedef typename std::remove_const<T>::type ValueType;
543
546 friend class RealRowIterator<const ValueType>;
547 friend class RealRowIterator<ValueType>;
548
551 : p(_p), i(_i)
552 {}
553
556 : p(0), i(0)
557 {}
558
560 : p(it.p), i(it.i)
561 {}
562
563
566 {
567 return i;
568 }
569
570 std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
571 {
572 assert(other.p==p);
573 return (other.i-i);
574 }
575
576 std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
577 {
578 assert(other.p==p);
579 return (other.i-i);
580 }
581
583 bool equals (const RealRowIterator<ValueType>& other) const
584 {
585 assert(other.p==p);
586 return i==other.i;
587 }
588
590 bool equals (const RealRowIterator<const ValueType>& other) const
591 {
592 assert(other.p==p);
593 return i==other.i;
594 }
595
596 private:
598 void increment()
599 {
600 ++i;
601 }
602
604 void decrement()
605 {
606 --i;
607 }
608
609 void advance(std::ptrdiff_t diff)
610 {
611 i+=diff;
612 }
613
614 T& elementAt(std::ptrdiff_t diff) const
615 {
616 return p[i+diff];
617 }
618
620 row_type& dereference () const
621 {
622 return p[i];
623 }
624
625 row_type* p;
626 size_type i;
627 };
628
632
635 {
636 return Iterator(r,0);
637 }
638
641 {
642 return Iterator(r,n);
643 }
644
648 {
649 return Iterator(r,n-1);
650 }
651
655 {
656 return Iterator(r,-1);
657 }
658
661
664
668
669
672 {
673 return ConstIterator(r,0);
674 }
675
678 {
679 return ConstIterator(r,n);
680 }
681
685 {
686 return ConstIterator(r,n-1);
687 }
688
692 {
693 return ConstIterator(r,-1);
694 }
695
698
701
702 //===== constructors & resizers
703
704 // we use a negative overflowsize to indicate that the implicit
705 // mode parameters have not been set yet
706
709 : build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
710 allocationSize_(0), r(0), a(0),
711 avg(0), overflowsize(-1.0)
712 {}
713
716 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
717 allocationSize_(0), r(0), a(0),
718 avg(0), overflowsize(-1.0)
719 {
720 allocate(_n, _m, _nnz,true,false);
721 }
722
725 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
726 allocationSize_(0), r(0), a(0),
727 avg(0), overflowsize(-1.0)
728 {
729 allocate(_n, _m,0,true,false);
730 }
731
733
743 BCRSMatrix (size_type _n, size_type _m, size_type _avg, double _overflowsize, BuildMode bm)
744 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
745 allocationSize_(0), r(0), a(0),
746 avg(_avg), overflowsize(_overflowsize)
747 {
748 if (bm != implicit)
749 DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
750 // Prevent user from setting a negative overflowsize:
751 // 1) It doesn't make sense
752 // 2) We use a negative overflow value to indicate that the parameters
753 // have not been set yet
754 if (_overflowsize < 0.0)
755 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
756 implicit_allocate(_n,_m);
757 }
758
765 : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
766 allocationSize_(0), r(0), a(0),
767 avg(Mat.avg), overflowsize(Mat.overflowsize)
768 {
769 if (!(Mat.ready == notAllocated || Mat.ready == built))
770 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
771
772 // deep copy in global array
773 size_type _nnz = Mat.nnz_;
774
775 // in case of row-wise allocation
776 if (_nnz<=0)
777 {
778 _nnz = 0;
779 for (size_type i=0; i<Mat.n; i++)
780 _nnz += Mat.r[i].getsize();
781 }
782
783 j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
784 allocate(Mat.n, Mat.m, _nnz, true, true);
785
786 // build window structure
788 }
789
792 {
793 deallocate();
794 }
795
801 {
802 if (ready == notAllocated)
803 {
804 build_mode = bm;
805 return;
806 }
807 if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
808 build_mode = bm;
809 else
810 DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
811 }
812
828 void setSize(size_type rows, size_type columns, size_type nnz=0)
829 {
830 // deallocate already setup memory
831 deallocate();
832
833 if (build_mode == implicit)
834 {
835 if (nnz>0)
836 DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");
837
838 // implicit allocates differently
839 implicit_allocate(rows,columns);
840 }
841 else
842 {
843 // allocate matrix memory
844 allocate(rows, columns, nnz, true, false);
845 }
846 }
847
856 void setImplicitBuildModeParameters(size_type _avg, double _overflow)
857 {
858 // Prevent user from setting a negative overflowsize:
859 // 1) It doesn't make sense
860 // 2) We use a negative overflow value to indicate that the parameters
861 // have not been set yet
862 if (_overflow < 0.0)
863 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
864
865 // make sure the parameters aren't changed after memory has been allocated
866 if (ready != notAllocated)
867 DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
868 avg = _avg;
869 overflowsize = _overflow;
870 }
871
879 {
880 // return immediately when self-assignment
881 if (&Mat==this) return *this;
882
883 if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
884 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");
885
886 // make it simple: ALWAYS throw away memory for a and j_
887 // and deallocate rows only if n != Mat.n
888 deallocate(n!=Mat.n);
889
890 // reallocate the rows if required
891 if (n>0 && n!=Mat.n) {
892 // free rows
893 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
894 rowAllocator_.destroy(riter);
895 rowAllocator_.deallocate(r,n);
896 }
897
898 nnz_ = Mat.nnz_;
899 if (nnz_ <= 0)
900 {
901 for (size_type i=0; i<Mat.n; i++)
902 nnz_ += Mat.r[i].getsize();
903 }
904
905 // allocate a, share j_
906 j_ = Mat.j_;
907 allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);
908
909 // build window structure
911 return *this;
912 }
913
916 {
917
918 if (!(ready == notAllocated || ready == built))
919 DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");
920
921 for (size_type i=0; i<n; i++) r[i] = k;
922 return *this;
923 }
924
925 //===== row-wise creation interface
926
929 {
930 public:
933 : Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
934 {
935 if (Mat.build_mode == unknown && Mat.ready == building)
936 {
937 Mat.build_mode = row_wise;
938 }
939 if (i==0 && Mat.ready != building)
940 DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
941 if(Mat.build_mode!=row_wise)
942 DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
943 }
944
947 {
948 // this should only be called if matrix is in creation
949 if (Mat.ready != building)
950 DUNE_THROW(BCRSMatrixError,"matrix already built up");
951
952 // row i is defined through the pattern
953 // get memory for the row and initialize the j_ array
954 // this depends on the allocation mode
955
956 // compute size of the row
957 size_type s = pattern.size();
958
959 if(s>0) {
960 // update number of nonzeroes including this row
961 nnz += s;
962
963 // alloc memory / set window
964 if (Mat.nnz_ > 0)
965 {
966 // memory is allocated in one long array
967
968 // check if that memory is sufficient
969 if (nnz > Mat.nnz_)
970 DUNE_THROW(BCRSMatrixError,"allocated nnz too small");
971
972 // set row i
973 Mat.r[i].set(s,nullptr,current_row.getindexptr());
974 current_row.setindexptr(current_row.getindexptr()+s);
975 }else{
976 // memory is allocated individually per row
977 // allocate and set row i
978 B* b = Mat.allocator_.allocate(s);
979 // use placement new to call constructor that allocates
980 // additional memory.
981 new (b) B[s];
982 size_type* j = Mat.sizeAllocator_.allocate(s);
983 Mat.r[i].set(s,b,j);
984 }
985 }else
986 // setup empty row
987 Mat.r[i].set(0,nullptr,nullptr);
988
989 // initialize the j array for row i from pattern
990 std::copy(pattern.cbegin(), pattern.cend(), Mat.r[i].getindexptr());
991
992 // now go to next row
993 i++;
994 pattern.clear();
995
996 // check if this was last row
997 if (i==Mat.n)
998 {
999 Mat.ready = built;
1000 if(Mat.nnz_ > 0)
1001 {
1002 // Set nnz to the exact number of nonzero blocks inserted
1003 // as some methods rely on it
1004 Mat.nnz_ = nnz;
1005 // allocate data array
1006 Mat.allocateData();
1007 Mat.setDataPointers();
1008 }
1009 }
1010 // done
1011 return *this;
1012 }
1013
1015 bool operator!= (const CreateIterator& it) const
1016 {
1017 return (i!=it.i) || (&Mat!=&it.Mat);
1018 }
1019
1021 bool operator== (const CreateIterator& it) const
1022 {
1023 return (i==it.i) && (&Mat==&it.Mat);
1024 }
1025
1028 {
1029 return i;
1030 }
1031
1034 {
1035 pattern.insert(j);
1036 }
1037
1040 {
1041 return pattern.find(j) != pattern.end();
1042 }
1049 {
1050 return pattern.size();
1051 }
1052
1053 private:
1054 BCRSMatrix& Mat; // the matrix we are defining
1055 size_type i; // current row to be defined
1056 size_type nnz; // count total number of nonzeros
1057 typedef std::set<size_type,std::less<size_type> > PatternType;
1058 PatternType pattern; // used to compile entries in a row
1059 row_type current_row; // row pointing to the current row to setup
1060 };
1061
1063 friend class CreateIterator;
1064
1067 {
1068 return CreateIterator(*this,0);
1069 }
1070
1073 {
1074 return CreateIterator(*this,n);
1075 }
1076
1077
1078 //===== random creation interface
1079
1087 {
1088 if (build_mode!=random)
1089 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1090 if (ready != building)
1091 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1092
1093 r[i].setsize(s);
1094 }
1095
1098 {
1099#ifdef DUNE_ISTL_WITH_CHECKING
1100 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
1101 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
1102#endif
1103 return r[i].getsize();
1104 }
1105
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 r[i].setsize(r[i].getsize()+s);
1115 }
1116
1119 {
1120 if (build_mode!=random)
1121 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1122 if (ready != building)
1123 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1124
1125 // compute total size, check positivity
1126 size_type total=0;
1127 for (size_type i=0; i<n; i++)
1128 {
1129 total += r[i].getsize();
1130 }
1131
1132 if(nnz_ == 0)
1133 // allocate/check memory
1134 allocate(n,m,total,false,false);
1135 else if(nnz_ < total)
1136 DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
1137 <<"sufficient for calculated nonzeros ("<<total<<"! ");
1138
1139 // set the window pointers correctly
1141
1142 // initialize j_ array with m (an invalid column index)
1143 // this indicates an unused entry
1144 for (size_type k=0; k<nnz_; k++)
1145 j_.get()[k] = m;
1146 ready = rowSizesBuilt;
1147 }
1148
1150
1161 {
1162 if (build_mode!=random)
1163 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1164 if (ready==built)
1165 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1166 if (ready==building)
1167 DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
1168 if (ready==notAllocated)
1169 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1170
1171 if (col >= m)
1172 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1173
1174 // get row range
1175 size_type* const first = r[row].getindexptr();
1176 size_type* const last = first + r[row].getsize();
1177
1178 // find correct insertion position for new column index
1179 size_type* pos = std::lower_bound(first,last,col);
1180
1181 // check if index is already in row
1182 if (pos!=last && *pos == col) return;
1183
1184 // find end of already inserted column indices
1185 size_type* end = std::lower_bound(pos,last,m);
1186 if (end==last)
1187 DUNE_THROW(BCRSMatrixError,"row is too small");
1188
1189 // insert new column index at correct position
1190 std::copy_backward(pos,end,end+1);
1191 *pos = col;
1192 }
1193
1195
1202 template<typename It>
1203 void setIndices(size_type row, It begin, It end)
1204 {
1205 size_type row_size = r[row].size();
1206 size_type* col_begin = r[row].getindexptr();
1207 size_type* col_end;
1208 // consistency check between allocated row size and number of passed column indices
1209 if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
1210 DUNE_THROW(BCRSMatrixError,"Given size of row " << row
1211 << " (" << row_size
1212 << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
1213 std::sort(col_begin,col_end);
1214 }
1215
1218 {
1219 if (build_mode!=random)
1220 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1221 if (ready==built)
1222 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1223 if (ready==building)
1224 DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
1225 if (ready==notAllocated)
1226 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1227
1228 // check if there are undefined indices
1229 RowIterator endi=end();
1230 for (RowIterator i=begin(); i!=endi; ++i)
1231 {
1232 ColIterator endj = (*i).end();
1233 for (ColIterator j=(*i).begin(); j!=endj; ++j) {
1234 if (j.index() >= m) {
1235 dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
1236 <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
1237 r[i.index()].setsize(j.offset());
1238 break;
1239 }
1240 }
1241 }
1242
1243 allocateData();
1245
1246 // if not, set matrix to built
1247 ready = built;
1248 }
1249
1250 //===== implicit creation interface
1251
1253
1265 {
1266#ifdef DUNE_ISTL_WITH_CHECKING
1267 if (build_mode!=implicit)
1268 DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1269 if (ready==built)
1270 DUNE_THROW(BCRSMatrixError,"matrix already built up, use operator[] for entry access now");
1271 if (ready==notAllocated)
1272 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1273 if (ready!=building)
1274 DUNE_THROW(InvalidStateException,"You may only use entry() during the 'building' stage");
1275
1276 if (row >= n)
1277 DUNE_THROW(BCRSMatrixError,"row index exceeds matrix size");
1278 if (col >= m)
1279 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1280#endif
1281
1282 size_type* begin = r[row].getindexptr();
1283 size_type* end = begin + r[row].getsize();
1284
1285 size_type* pos = std::find(begin, end, col);
1286
1287 //treat the case that there was a match in the array
1288 if (pos != end)
1289 if (*pos == col)
1290 {
1291 std::ptrdiff_t offset = pos - r[row].getindexptr();
1292 B* aptr = r[row].getptr() + offset;
1293
1294 return *aptr;
1295 }
1296
1297 //determine whether overflow has to be taken into account or not
1298 if (r[row].getsize() == avg)
1299 return overflow[std::make_pair(row,col)];
1300 else
1301 {
1302 //modify index array
1303 *end = col;
1304
1305 //do simulatenous operations on data array a
1306 std::ptrdiff_t offset = end - r[row].getindexptr();
1307 B* apos = r[row].getptr() + offset;
1308
1309 //increase rowsize
1310 r[row].setsize(r[row].getsize()+1);
1311
1312 //return reference to the newly created entry
1313 return *apos;
1314 }
1315 }
1316
1318
1329 {
1330 if (build_mode!=implicit)
1331 DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1332 if (ready==built)
1333 DUNE_THROW(BCRSMatrixError,"matrix already built up, no more need for compression");
1334 if (ready==notAllocated)
1335 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1336 if (ready!=building)
1337 DUNE_THROW(InvalidStateException,"You may only call compress() at the end of the 'building' stage");
1338
1339 //calculate statistics
1341 stats.overflow_total = overflow.size();
1342 stats.maximum = 0;
1343
1344 //get insertion iterators pointing to one before start (for later use of ++it)
1345 size_type* jiit = j_.get();
1346 B* aiit = a;
1347
1348 //get iterator to the smallest overflow element
1349 typename OverflowType::iterator oit = overflow.begin();
1350
1351 //store a copy of index pointers on which to perform sortation
1352 std::vector<size_type*> perm;
1353
1354 //iterate over all rows and copy elements into their position in the compressed array
1355 for (size_type i=0; i<n; i++)
1356 {
1357 //get old pointers into a and j and size without overflow changes
1358 size_type* begin = r[i].getindexptr();
1359 //B* apos = r[i].getptr();
1360 size_type size = r[i].getsize();
1361
1362 perm.resize(size);
1363
1364 typename std::vector<size_type*>::iterator it = perm.begin();
1365 for (size_type* iit = begin; iit < begin + size; ++iit, ++it)
1366 *it = iit;
1367
1368 //sort permutation array
1369 std::sort(perm.begin(),perm.end(),PointerCompare<size_type>());
1370
1371 //change row window pointer to their new positions
1372 r[i].setindexptr(jiit);
1373 r[i].setptr(aiit);
1374
1375 for (it = perm.begin(); it != perm.end(); ++it)
1376 {
1377 //check whether there are elements in the overflow area which take precedence
1378 while ((oit!=overflow.end()) && (oit->first < std::make_pair(i,**it)))
1379 {
1380 //check whether there is enough memory to write to
1381 if (jiit > begin)
1383 "Allocated memory for BCRSMatrix exhausted during compress()!"
1384 "Please increase either the average number of entries per row or the overflow fraction."
1385 );
1386 //copy an element from the overflow area to the insertion position in a and j
1387 *jiit = oit->first.second;
1388 ++jiit;
1389 *aiit = oit->second;
1390 ++aiit;
1391 ++oit;
1392 r[i].setsize(r[i].getsize()+1);
1393 }
1394
1395 //check whether there is enough memory to write to
1396 if (jiit > begin)
1398 "Allocated memory for BCRSMatrix exhausted during compress()!"
1399 "Please increase either the average number of entries per row or the overflow fraction."
1400 );
1401
1402 //copy element from array
1403 *jiit = **it;
1404 ++jiit;
1405 B* apos = *it - j_.get() + a;
1406 *aiit = *apos;
1407 ++aiit;
1408 }
1409
1410 //copy remaining elements from the overflow area
1411 while ((oit!=overflow.end()) && (oit->first.first == i))
1412 {
1413 //check whether there is enough memory to write to
1414 if (jiit > begin)
1416 "Allocated memory for BCRSMatrix exhausted during compress()!"
1417 "Please increase either the average number of entries per row or the overflow fraction."
1418 );
1419
1420 //copy and element from the overflow area to the insertion position in a and j
1421 *jiit = oit->first.second;
1422 ++jiit;
1423 *aiit = oit->second;
1424 ++aiit;
1425 ++oit;
1426 r[i].setsize(r[i].getsize()+1);
1427 }
1428
1429 // update maximum row size
1430 if (r[i].getsize()>stats.maximum)
1431 stats.maximum = r[i].getsize();
1432 }
1433
1434 // overflow area may be cleared
1435 overflow.clear();
1436
1437 //determine average number of entries and memory usage
1438 std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
1439 nnz_ = diff;
1440 stats.avg = (double) (nnz_) / (double) n;
1441 stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
1442
1443 //matrix is now built
1444 ready = built;
1445
1446 return stats;
1447 }
1448
1449 //===== vector space arithmetic
1450
1453 {
1454#ifdef DUNE_ISTL_WITH_CHECKING
1455 if (ready != built)
1456 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1457#endif
1458
1459 if (nnz_ > 0)
1460 {
1461 // process 1D array
1462 for (size_type i=0; i<nnz_; i++)
1463 a[i] *= k;
1464 }
1465 else
1466 {
1467 RowIterator endi=end();
1468 for (RowIterator i=begin(); i!=endi; ++i)
1469 {
1470 ColIterator endj = (*i).end();
1471 for (ColIterator j=(*i).begin(); j!=endj; ++j)
1472 (*j) *= k;
1473 }
1474 }
1475
1476 return *this;
1477 }
1478
1481 {
1482#ifdef DUNE_ISTL_WITH_CHECKING
1483 if (ready != built)
1484 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1485#endif
1486
1487 if (nnz_ > 0)
1488 {
1489 // process 1D array
1490 for (size_type i=0; i<nnz_; i++)
1491 a[i] /= k;
1492 }
1493 else
1494 {
1495 RowIterator endi=end();
1496 for (RowIterator i=begin(); i!=endi; ++i)
1497 {
1498 ColIterator endj = (*i).end();
1499 for (ColIterator j=(*i).begin(); j!=endj; ++j)
1500 (*j) /= k;
1501 }
1502 }
1503
1504 return *this;
1505 }
1506
1507
1514 {
1515#ifdef DUNE_ISTL_WITH_CHECKING
1516 if (ready != built || b.ready != built)
1517 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1518 if(N()!=b.N() || M() != b.M())
1519 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1520#endif
1521 RowIterator endi=end();
1522 ConstRowIterator j=b.begin();
1523 for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1524 i->operator+=(*j);
1525 }
1526
1527 return *this;
1528 }
1529
1536 {
1537#ifdef DUNE_ISTL_WITH_CHECKING
1538 if (ready != built || b.ready != built)
1539 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1540 if(N()!=b.N() || M() != b.M())
1541 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1542#endif
1543 RowIterator endi=end();
1544 ConstRowIterator j=b.begin();
1545 for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1546 i->operator-=(*j);
1547 }
1548
1549 return *this;
1550 }
1551
1561 {
1562#ifdef DUNE_ISTL_WITH_CHECKING
1563 if (ready != built || b.ready != built)
1564 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1565 if(N()!=b.N() || M() != b.M())
1566 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1567#endif
1568 RowIterator endi=end();
1569 ConstRowIterator j=b.begin();
1570 for(RowIterator i=begin(); i!=endi; ++i, ++j)
1571 i->axpy(alpha, *j);
1572
1573 return *this;
1574 }
1575
1576 //===== linear maps
1577
1579 template<class X, class Y>
1580 void mv (const X& x, Y& y) const
1581 {
1582#ifdef DUNE_ISTL_WITH_CHECKING
1583 if (ready != built)
1584 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1585 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,
1586 "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
1587 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,
1588 "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
1589#endif
1590 ConstRowIterator endi=end();
1591 for (ConstRowIterator i=begin(); i!=endi; ++i)
1592 {
1593 y[i.index()]=0;
1594 ConstColIterator endj = (*i).end();
1595 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1596 (*j).umv(x[j.index()],y[i.index()]);
1597 }
1598 }
1599
1601 template<class X, class Y>
1602 void umv (const X& x, Y& y) const
1603 {
1604#ifdef DUNE_ISTL_WITH_CHECKING
1605 if (ready != built)
1606 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1607 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1608 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1609#endif
1610 ConstRowIterator endi=end();
1611 for (ConstRowIterator i=begin(); i!=endi; ++i)
1612 {
1613 ConstColIterator endj = (*i).end();
1614 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1615 (*j).umv(x[j.index()],y[i.index()]);
1616 }
1617 }
1618
1620 template<class X, class Y>
1621 void mmv (const X& x, Y& y) const
1622 {
1623#ifdef DUNE_ISTL_WITH_CHECKING
1624 if (ready != built)
1625 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1626 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1627 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1628#endif
1629 ConstRowIterator endi=end();
1630 for (ConstRowIterator i=begin(); i!=endi; ++i)
1631 {
1632 ConstColIterator endj = (*i).end();
1633 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1634 (*j).mmv(x[j.index()],y[i.index()]);
1635 }
1636 }
1637
1639 template<class X, class Y, class F>
1640 void usmv (F&& alpha, const X& x, Y& y) const
1641 {
1642#ifdef DUNE_ISTL_WITH_CHECKING
1643 if (ready != built)
1644 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1645 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1646 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1647#endif
1648 ConstRowIterator endi=end();
1649 for (ConstRowIterator i=begin(); i!=endi; ++i)
1650 {
1651 ConstColIterator endj = (*i).end();
1652 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1653 (*j).usmv(alpha,x[j.index()],y[i.index()]);
1654 }
1655 }
1656
1658 template<class X, class Y>
1659 void mtv (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()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1665 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1666#endif
1667 for(size_type i=0; i<y.N(); ++i)
1668 y[i]=0;
1669 umtv(x,y);
1670 }
1671
1673 template<class X, class Y>
1674 void umtv (const X& x, Y& y) const
1675 {
1676#ifdef DUNE_ISTL_WITH_CHECKING
1677 if (ready != built)
1678 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1679 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1680 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1681#endif
1682 ConstRowIterator endi=end();
1683 for (ConstRowIterator i=begin(); i!=endi; ++i)
1684 {
1685 ConstColIterator endj = (*i).end();
1686 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1687 (*j).umtv(x[i.index()],y[j.index()]);
1688 }
1689 }
1690
1692 template<class X, class Y>
1693 void mmtv (const X& x, Y& y) const
1694 {
1695#ifdef DUNE_ISTL_WITH_CHECKING
1696 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1697 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1698#endif
1699 ConstRowIterator endi=end();
1700 for (ConstRowIterator i=begin(); i!=endi; ++i)
1701 {
1702 ConstColIterator endj = (*i).end();
1703 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1704 (*j).mmtv(x[i.index()],y[j.index()]);
1705 }
1706 }
1707
1709 template<class X, class Y>
1710 void usmtv (const field_type& alpha, const X& x, Y& y) const
1711 {
1712#ifdef DUNE_ISTL_WITH_CHECKING
1713 if (ready != built)
1714 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1715 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1716 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1717#endif
1718 ConstRowIterator endi=end();
1719 for (ConstRowIterator i=begin(); i!=endi; ++i)
1720 {
1721 ConstColIterator endj = (*i).end();
1722 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1723 (*j).usmtv(alpha,x[i.index()],y[j.index()]);
1724 }
1725 }
1726
1728 template<class X, class Y>
1729 void umhv (const X& x, Y& y) const
1730 {
1731#ifdef DUNE_ISTL_WITH_CHECKING
1732 if (ready != built)
1733 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1734 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1735 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1736#endif
1737 ConstRowIterator endi=end();
1738 for (ConstRowIterator i=begin(); i!=endi; ++i)
1739 {
1740 ConstColIterator endj = (*i).end();
1741 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1742 (*j).umhv(x[i.index()],y[j.index()]);
1743 }
1744 }
1745
1747 template<class X, class Y>
1748 void mmhv (const X& x, Y& y) const
1749 {
1750#ifdef DUNE_ISTL_WITH_CHECKING
1751 if (ready != built)
1752 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1753 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1754 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1755#endif
1756 ConstRowIterator endi=end();
1757 for (ConstRowIterator i=begin(); i!=endi; ++i)
1758 {
1759 ConstColIterator endj = (*i).end();
1760 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1761 (*j).mmhv(x[i.index()],y[j.index()]);
1762 }
1763 }
1764
1766 template<class X, class Y>
1767 void usmhv (const field_type& alpha, const X& x, Y& y) const
1768 {
1769#ifdef DUNE_ISTL_WITH_CHECKING
1770 if (ready != built)
1771 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1772 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1773 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1774#endif
1775 ConstRowIterator endi=end();
1776 for (ConstRowIterator i=begin(); i!=endi; ++i)
1777 {
1778 ConstColIterator endj = (*i).end();
1779 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1780 (*j).usmhv(alpha,x[i.index()],y[j.index()]);
1781 }
1782 }
1783
1784
1785 //===== norms
1786
1788 typename FieldTraits<field_type>::real_type frobenius_norm2 () const
1789 {
1790#ifdef DUNE_ISTL_WITH_CHECKING
1791 if (ready != built)
1792 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1793#endif
1794
1795 double sum=0;
1796
1797 ConstRowIterator endi=end();
1798 for (ConstRowIterator i=begin(); i!=endi; ++i)
1799 {
1800 ConstColIterator endj = (*i).end();
1801 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1802 sum += (*j).frobenius_norm2();
1803 }
1804
1805 return sum;
1806 }
1807
1809 typename FieldTraits<field_type>::real_type frobenius_norm () const
1810 {
1811 return sqrt(frobenius_norm2());
1812 }
1813
1815 template <typename ft = field_type,
1816 typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
1817 typename FieldTraits<ft>::real_type infinity_norm() const {
1818 if (ready != built)
1819 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1820
1821 using real_type = typename FieldTraits<ft>::real_type;
1822 using std::max;
1823
1824 real_type norm = 0;
1825 for (auto const &x : *this) {
1826 real_type sum = 0;
1827 for (auto const &y : x)
1828 sum += y.infinity_norm();
1829 norm = max(sum, norm);
1830 }
1831 return norm;
1832 }
1833
1835 template <typename ft = field_type,
1836 typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
1837 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1838 if (ready != built)
1839 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1840
1841 using real_type = typename FieldTraits<ft>::real_type;
1842 using std::max;
1843
1844 real_type norm = 0;
1845 for (auto const &x : *this) {
1846 real_type sum = 0;
1847 for (auto const &y : x)
1848 sum += y.infinity_norm_real();
1849 norm = max(sum, norm);
1850 }
1851 return norm;
1852 }
1853
1855 template <typename ft = field_type,
1856 typename std::enable_if<has_nan<ft>::value, int>::type = 0>
1857 typename FieldTraits<ft>::real_type infinity_norm() const {
1858 if (ready != built)
1859 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1860
1861 using real_type = typename FieldTraits<ft>::real_type;
1862 using std::max;
1863
1864 real_type norm = 0;
1865 real_type isNaN = 1;
1866 for (auto const &x : *this) {
1867 real_type sum = 0;
1868 for (auto const &y : x)
1869 sum += y.infinity_norm();
1870 norm = max(sum, norm);
1871 isNaN += sum;
1872 }
1873 isNaN /= isNaN;
1874 return norm * isNaN;
1875 }
1876
1878 template <typename ft = field_type,
1879 typename std::enable_if<has_nan<ft>::value, int>::type = 0>
1880 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1881 if (ready != built)
1882 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1883
1884 using real_type = typename FieldTraits<ft>::real_type;
1885 using std::max;
1886
1887 real_type norm = 0;
1888 real_type isNaN = 1;
1889 for (auto const &x : *this) {
1890 real_type sum = 0;
1891 for (auto const &y : x)
1892 sum += y.infinity_norm_real();
1893 norm = max(sum, norm);
1894 isNaN += sum;
1895 }
1896 isNaN /= isNaN;
1897 return norm * isNaN;
1898 }
1899
1900 //===== sizes
1901
1903 size_type N () const
1904 {
1905 return n;
1906 }
1907
1909 size_type M () const
1910 {
1911 return m;
1912 }
1913
1916 {
1917 return nnz_;
1918 }
1919
1922 {
1923 return ready;
1924 }
1925
1928 {
1929 return build_mode;
1930 }
1931
1932 //===== query
1933
1935 bool exists (size_type i, size_type j) const
1936 {
1937#ifdef DUNE_ISTL_WITH_CHECKING
1938 if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
1939 if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
1940#endif
1941 return (r[i].size() && r[i].find(j) != r[i].end());
1942 }
1943
1944
1945 protected:
1946 // state information
1947 BuildMode build_mode; // row wise or whole matrix
1948 BuildStage ready; // indicate the stage the matrix building is in
1949
1950 // The allocator used for memory management
1951 typename A::template rebind<B>::other allocator_;
1952
1953 typename A::template rebind<row_type>::other rowAllocator_;
1954
1955 typename A::template rebind<size_type>::other sizeAllocator_;
1956
1957 // size of the matrix
1958 size_type n; // number of rows
1959 size_type m; // number of columns
1960 size_type nnz_; // number of nonzeroes contained in the matrix
1961 size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
1962 // zero means that memory is allocated separately for each row.
1963
1964 // the rows are dynamically allocated
1965 row_type* r; // [n] the individual rows having pointers into a,j arrays
1966
1967 // dynamically allocated memory
1968 B* a; // [allocationSize] non-zero entries of the matrix in row-wise ordering
1969 // If a single array of column indices is used, it can be shared
1970 // between different matrices with the same sparsity pattern
1971 std::shared_ptr<size_type> j_; // [allocationSize] column indices of entries
1972
1973 // additional data is needed in implicit buildmode
1974 size_type avg;
1975 double overflowsize;
1976
1977 typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
1978 OverflowType overflow;
1979
1980 void setWindowPointers(ConstRowIterator row)
1981 {
1982 row_type current_row(a,j_.get(),0); // Pointers to current row data
1983 for (size_type i=0; i<n; i++, ++row) {
1984 // set row i
1985 size_type s = row->getsize();
1986
1987 if (s>0) {
1988 // setup pointers and size
1989 r[i].set(s,current_row.getptr(), current_row.getindexptr());
1990 // update pointer for next row
1991 current_row.setptr(current_row.getptr()+s);
1992 current_row.setindexptr(current_row.getindexptr()+s);
1993 } else{
1994 // empty row
1995 r[i].set(0,nullptr,nullptr);
1996 }
1997 }
1998 }
1999
2001
2006 {
2007 size_type* jptr = j_.get();
2008 for (size_type i=0; i<n; ++i, ++row) {
2009 // set row i
2010 size_type s = row->getsize();
2011
2012 if (s>0) {
2013 // setup pointers and size
2014 r[i].setsize(s);
2015 r[i].setindexptr(jptr);
2016 } else{
2017 // empty row
2018 r[i].set(0,nullptr,nullptr);
2019 }
2020
2021 // advance position in global array
2022 jptr += s;
2023 }
2024 }
2025
2027
2032 {
2033 B* aptr = a;
2034 for (size_type i=0; i<n; ++i) {
2035 // set row i
2036 if (r[i].getsize() > 0) {
2037 // setup pointers and size
2038 r[i].setptr(aptr);
2039 } else{
2040 // empty row
2041 r[i].set(0,nullptr,nullptr);
2042 }
2043
2044 // advance position in global array
2045 aptr += r[i].getsize();
2046 }
2047 }
2048
2051 {
2052 setWindowPointers(Mat.begin());
2053
2054 // copy data
2055 for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
2056
2057 // finish off
2058 build_mode = row_wise; // dummy
2059 ready = built;
2060 }
2061
2067 void deallocate(bool deallocateRows=true)
2068 {
2069
2070 if (notAllocated)
2071 return;
2072
2073 if (allocationSize_>0)
2074 {
2075 // a,j_ have been allocated as one long vector
2076 j_.reset();
2077 if (a)
2078 {
2079 for(B *aiter=a+(allocationSize_-1), *aend=a-1; aiter!=aend; --aiter)
2080 allocator_.destroy(aiter);
2081 allocator_.deallocate(a,allocationSize_);
2082 a = nullptr;
2083 }
2084 }
2085 else if (r)
2086 {
2087 // check if memory for rows have been allocated individually
2088 for (size_type i=0; i<n; i++)
2089 if (r[i].getsize()>0)
2090 {
2091 for (B *col=r[i].getptr()+(r[i].getsize()-1),
2092 *colend = r[i].getptr()-1; col!=colend; --col) {
2093 allocator_.destroy(col);
2094 }
2095 sizeAllocator_.deallocate(r[i].getindexptr(),1);
2096 allocator_.deallocate(r[i].getptr(),1);
2097 // clear out row data in case we don't want to deallocate the rows
2098 // otherwise we might run into a double free problem here later
2099 r[i].set(0,nullptr,nullptr);
2100 }
2101 }
2102
2103 // deallocate the rows
2104 if (n>0 && deallocateRows && r) {
2105 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
2106 rowAllocator_.destroy(riter);
2107 rowAllocator_.deallocate(r,n);
2108 r = nullptr;
2109 }
2110
2111 // Mark matrix as not built at all.
2112 ready=notAllocated;
2113
2114 }
2115
2118 {
2119 typename A::template rebind<size_type>::other& sizeAllocator_;
2120
2121 public:
2122 Deallocator(typename A::template rebind<size_type>::other& sizeAllocator)
2123 : sizeAllocator_(sizeAllocator)
2124 {}
2125
2126 void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); }
2127 };
2128
2129
2147 void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
2148 {
2149 // Store size
2150 n = rows;
2151 m = columns;
2152 nnz_ = allocationSize;
2153 allocationSize_ = allocationSize;
2154
2155 // allocate rows
2156 if(allocateRows) {
2157 if (n>0) {
2158 if (r)
2159 DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
2160 r = rowAllocator_.allocate(rows);
2161 }else{
2162 r = 0;
2163 }
2164 }
2165
2166 // allocate a and j_ array
2167 if (allocate_data)
2168 allocateData();
2169 if (allocationSize_>0) {
2170 // allocate column indices only if not yet present (enable sharing)
2171 if (!j_.get())
2172 j_.reset(sizeAllocator_.allocate(allocationSize_),Deallocator(sizeAllocator_));
2173 }else{
2174 j_.reset();
2175 for(row_type* ri=r; ri!=r+rows; ++ri)
2176 rowAllocator_.construct(ri, row_type());
2177 }
2178
2179 // Mark the matrix as not built.
2180 ready = building;
2181 }
2182
2183 void allocateData()
2184 {
2185 if (a)
2186 DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
2187 if (allocationSize_>0) {
2188 a = allocator_.allocate(allocationSize_);
2189 // use placement new to call constructor that allocates
2190 // additional memory.
2191 new (a) B[allocationSize_];
2192 } else {
2193 a = nullptr;
2194 }
2195 }
2196
2203 {
2204 if (build_mode != implicit)
2205 DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
2206 if (ready != notAllocated)
2207 DUNE_THROW(InvalidStateException,"memory has already been allocated");
2208
2209 // check to make sure the user has actually set the parameters
2210 if (overflowsize < 0)
2211 DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
2212 //calculate size of overflow region, add buffer for row sort!
2213 size_type osize = (size_type) (_n*avg)*overflowsize + 4*avg;
2214 allocationSize_ = _n*avg + osize;
2215
2216 allocate(_n, _m, allocationSize_,true,true);
2217
2218 //set row pointers correctly
2219 size_type* jptr = j_.get() + osize;
2220 B* aptr = a + osize;
2221 for (size_type i=0; i<n; i++)
2222 {
2223 r[i].set(0,aptr,jptr);
2224 jptr = jptr + avg;
2225 aptr = aptr + avg;
2226 }
2227
2228 ready = building;
2229 }
2230 };
2231
2232
2235} // end namespace
2236
2237#endif
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Error specific to BCRSMatrix.
Definition: istlexception.hh:21
Iterator class for sequential creation of blocks
Definition: bcrsmatrix.hh:929
bool operator==(const CreateIterator &it) const
equality
Definition: bcrsmatrix.hh:1021
CreateIterator & operator++()
prefix increment
Definition: bcrsmatrix.hh:946
size_type index() const
The number of the row that the iterator currently points to.
Definition: bcrsmatrix.hh:1027
bool operator!=(const CreateIterator &it) const
inequality
Definition: bcrsmatrix.hh:1015
CreateIterator(BCRSMatrix &_Mat, size_type _i)
constructor
Definition: bcrsmatrix.hh:932
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1033
bool contains(size_type j)
return true if column index is in row
Definition: bcrsmatrix.hh:1039
size_type size() const
Get the current row size.
Definition: bcrsmatrix.hh:1048
Class used by shared_ptr to deallocate memory using the proper allocator.
Definition: bcrsmatrix.hh:2118
Iterator access to matrix rows
Definition: bcrsmatrix.hh:538
RealRowIterator()
empty constructor, use with care!
Definition: bcrsmatrix.hh:555
bool equals(const RealRowIterator< ValueType > &other) const
equality
Definition: bcrsmatrix.hh:583
std::remove_const< T >::type ValueType
The unqualified value type.
Definition: bcrsmatrix.hh:542
bool equals(const RealRowIterator< const ValueType > &other) const
equality
Definition: bcrsmatrix.hh:590
RealRowIterator(row_type *_p, size_type _i)
constructor
Definition: bcrsmatrix.hh:550
size_type index() const
return index
Definition: bcrsmatrix.hh:565
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:423
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: bcrsmatrix.hh:1935
BuildStage buildStage() const
The current build stage of the matrix.
Definition: bcrsmatrix.hh:1921
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:634
friend class CreateIterator
allow CreateIterator to access internal data
Definition: bcrsmatrix.hh:1063
void copyWindowStructure(const BCRSMatrix &Mat)
Copy the window structure from another matrix.
Definition: bcrsmatrix.hh:2050
B & entry(size_type row, size_type col)
Returns reference to entry (row,col) of the matrix.
Definition: bcrsmatrix.hh:1264
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: bcrsmatrix.hh:1748
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: bcrsmatrix.hh:1767
void umtv(const X &x, Y &y) const
y += A^T x
Definition: bcrsmatrix.hh:1674
RealRowIterator< const row_type > const_iterator
The const iterator over the matrix rows.
Definition: bcrsmatrix.hh:666
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:743
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:2147
B::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:445
BCRSMatrix & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition: bcrsmatrix.hh:1560
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: bcrsmatrix.hh:1837
~BCRSMatrix()
destructor
Definition: bcrsmatrix.hh:791
void deallocate(bool deallocateRows=true)
deallocate memory of the matrix.
Definition: bcrsmatrix.hh:2067
Iterator RowIterator
rename the iterators for easier access
Definition: bcrsmatrix.hh:660
row_type & operator[](size_type i)
random access to the rows
Definition: bcrsmatrix.hh:508
BCRSMatrix()
an empty matrix
Definition: bcrsmatrix.hh:708
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1118
void incrementrowsize(size_type i, size_type s=1)
increment size of row i by s (1 by default)
Definition: bcrsmatrix.hh:1107
void mtv(const X &x, Y &y) const
y = A^T x
Definition: bcrsmatrix.hh:1659
void umhv(const X &x, Y &y) const
y += A^H x
Definition: bcrsmatrix.hh:1729
size_type nonzeroes() const
number of blocks that are stored (the number of blocks that possibly are nonzero)
Definition: bcrsmatrix.hh:1915
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition: bcrsmatrix.hh:697
void setrowsize(size_type i, size_type s)
Set number of indices in row i to s.
Definition: bcrsmatrix.hh:1086
BCRSMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: bcrsmatrix.hh:1452
RealRowIterator< row_type > iterator
The iterator over the (mutable matrix rows.
Definition: bcrsmatrix.hh:630
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: bcrsmatrix.hh:1710
ConstIterator beforeBegin() const
Definition: bcrsmatrix.hh:691
Iterator beforeBegin()
Definition: bcrsmatrix.hh:654
@ blocklevel
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:465
BuildMode
we support two modes
Definition: bcrsmatrix.hh:469
@ implicit
Build entries randomly with an educated guess for the number of entries per row.
Definition: bcrsmatrix.hh:498
@ unknown
Build mode not set!
Definition: bcrsmatrix.hh:502
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:489
@ row_wise
Build in a row-wise manner.
Definition: bcrsmatrix.hh:480
BCRSMatrix(size_type _n, size_type _m, size_type _nnz, BuildMode bm)
matrix with known number of nonzeroes
Definition: bcrsmatrix.hh:715
::Dune::CompressionStatistics< size_type > CompressionStatistics
The type for the statistics object returned by compress()
Definition: bcrsmatrix.hh:460
BCRSMatrix & operator-=(const BCRSMatrix &b)
Subtract the entries of another matrix from this one.
Definition: bcrsmatrix.hh:1535
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition: bcrsmatrix.hh:764
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:640
void setIndices(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition: bcrsmatrix.hh:1203
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1160
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:663
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: bcrsmatrix.hh:1809
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:457
BCRSMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: bcrsmatrix.hh:1480
BCRSMatrix & operator+=(const BCRSMatrix &b)
Add the entries of another matrix to this one.
Definition: bcrsmatrix.hh:1513
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1072
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: bcrsmatrix.hh:1788
Iterator beforeEnd()
Definition: bcrsmatrix.hh:647
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:700
void usmv(F &&alpha, const X &x, Y &y) const
y += alpha A x
Definition: bcrsmatrix.hh:1640
CompressedBlockVectorWindow< B, A > row_type
implement row_type with compressed vector
Definition: bcrsmatrix.hh:454
size_type getrowsize(size_type i) const
get current number of indices in row i
Definition: bcrsmatrix.hh:1097
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1909
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1066
BuildStage
Definition: bcrsmatrix.hh:426
@ rowSizesBuilt
The row sizes of the matrix are known.
Definition: bcrsmatrix.hh:437
@ built
The matrix structure is fully built.
Definition: bcrsmatrix.hh:439
@ notbuilt
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:428
@ notAllocated
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:430
@ building
Matrix is currently being built, some memory has been allocated, build mode and size are fixed.
Definition: bcrsmatrix.hh:432
BuildMode buildMode() const
The currently selected build mode of the matrix.
Definition: bcrsmatrix.hh:1927
void mmv(const X &x, Y &y) const
y -= A x
Definition: bcrsmatrix.hh:1621
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: bcrsmatrix.hh:1817
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1580
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:448
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: bcrsmatrix.hh:1693
void umv(const X &x, Y &y) const
y += A x
Definition: bcrsmatrix.hh:1602
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:2202
BCRSMatrix(size_type _n, size_type _m, BuildMode bm)
matrix with unknown number of nonzeroes
Definition: bcrsmatrix.hh:724
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1217
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition: bcrsmatrix.hh:1328
void setDataPointers()
Set data pointers for all rows.
Definition: bcrsmatrix.hh:2031
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1903
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:800
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:828
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:878
void setColumnPointers(ConstRowIterator row)
Copy row sizes from iterator range starting at row and set column index pointers for all rows.
Definition: bcrsmatrix.hh:2005
ConstIterator end() const
Get const iterator to one beyond last row.
Definition: bcrsmatrix.hh:677
ConstIterator begin() const
Get const iterator to first row.
Definition: bcrsmatrix.hh:671
void setImplicitBuildModeParameters(size_type _avg, double _overflow)
Set parameters needed for creation in implicit build mode.
Definition: bcrsmatrix.hh:856
A allocator_type
export the allocator type
Definition: bcrsmatrix.hh:451
ConstIterator beforeEnd() const
Definition: bcrsmatrix.hh:684
compressed_block_vector_unmanaged< B, A >::ConstIterator ConstIterator
make iterators available as types
Definition: bvector.hh:1052
void set(size_type _n, B *_p, size_type *_j)
set size and pointer
Definition: bvector.hh:1104
size_type * getindexptr()
get pointer
Definition: bvector.hh:1136
void setsize(size_type _n)
set size only
Definition: bvector.hh:1112
compressed_block_vector_unmanaged< B, A >::Iterator Iterator
make iterators available as types
Definition: bvector.hh:1049
void setptr(B *_p)
set pointer only
Definition: bvector.hh:1118
B * getptr()
get pointer
Definition: bvector.hh:1130
size_type getsize() const
get size
Definition: bvector.hh:1153
void setindexptr(size_type *_j)
set pointer only
Definition: bvector.hh:1124
Proxy row object for entry access.
Definition: bcrsmatrix.hh:130
block_type & operator[](size_type j) const
Returns entry in column j.
Definition: bcrsmatrix.hh:135
A wrapper for uniform access to the BCRSMatrix during and after the build stage in implicit build mod...
Definition: bcrsmatrix.hh:110
Matrix::block_type block_type
The block_type of the underlying matrix.
Definition: bcrsmatrix.hh:118
ImplicitMatrixBuilder(Matrix &m)
Creates an ImplicitMatrixBuilder for matrix m.
Definition: bcrsmatrix.hh:163
M_ Matrix
The underlying matrix.
Definition: bcrsmatrix.hh:115
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 M() const
The number of columns in the matrix.
Definition: bcrsmatrix.hh:210
Matrix::size_type size_type
The size_type of the underlying matrix.
Definition: bcrsmatrix.hh:121
row_object operator[](size_type i) const
Returns a proxy for entries in row i.
Definition: bcrsmatrix.hh:198
size_type N() const
The number of rows in the matrix.
Definition: bcrsmatrix.hh:204
The overflow error used during implicit BCRSMatrix construction was exhausted.
Definition: istlexception.hh:34
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:279
A generic dynamic dense matrix.
Definition: matrix.hh:555
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:571
T block_type
Export the type representing the components.
Definition: matrix.hh:562
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:433
Default exception class for range errors.
Definition: exceptions.hh:252
size_type size() const
number of blocks in the array (are of size 1 here)
Definition: basearray.hh:735
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
constexpr decltype(auto) elementAt(Container &&c, Index &&i)
Get element at given position from container.
Definition: hybridutilities.hh:138
DWarnType dwarn(std::cerr)
Stream for warnings indicating problems.
Definition: stdstreams.hh:159
This file implements iterator facade classes for writing stl conformant iterators.
Some handy generic functions for ISTL matrices.
Dune namespace.
Definition: alignment.hh:11
Standard Dune debug streams.
Statistics about compression achieved in implicit mode.
Definition: bcrsmatrix.hh:81
size_type overflow_total
total number of elements written to the overflow area during construction.
Definition: bcrsmatrix.hh:87
size_type maximum
maximum number of non-zeroes per row.
Definition: bcrsmatrix.hh:85
double avg
average number of non-zeroes per row.
Definition: bcrsmatrix.hh:83
double mem_ratio
fraction of wasted memory resulting from non-used overflow area.
Definition: bcrsmatrix.hh:92
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)