Dune Core Modules (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"
26
28
33namespace 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
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;
157 size_type _i;
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:
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
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
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
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
804 : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
805 allocationSize_(0), r(0), a(0),
806 avg(Mat.avg), compressionBufferSize_(Mat.compressionBufferSize_)
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
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
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
1057 {
1058 return i;
1059 }
1060
1063 {
1064 pattern.insert(j);
1065 }
1066
1069 {
1070 return pattern.find(j) != pattern.end();
1071 }
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
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
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();
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
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
2044 size_type avg;
2045 double compressionBufferSize_;
2046
2047 typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
2048 OverflowType overflow;
2049
2050 void setWindowPointers(ConstRowIterator row)
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.
2182 ready=notAllocated;
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())
2245 j_.reset(sizeAllocator_.allocate(allocationSize_),Deallocator(sizeAllocator_));
2246 }else{
2247 j_.reset();
2248 }
2249
2250 // Mark the matrix as not built.
2251 ready = building;
2252 }
2253
2254 void allocateData()
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...
Error specific to BCRSMatrix.
Definition: istlexception.hh:22
Iterator class for sequential creation of blocks
Definition: bcrsmatrix.hh:955
bool operator==(const CreateIterator &it) const
equality
Definition: bcrsmatrix.hh:1050
CreateIterator & operator++()
prefix increment
Definition: bcrsmatrix.hh:975
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
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
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
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
size_type index() const
return index
Definition: bcrsmatrix.hh:604
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:486
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
B & entry(size_type row, size_type col)
Returns reference to entry (row,col) of the matrix.
Definition: bcrsmatrix.hh:1294
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
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 & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition: bcrsmatrix.hh:1590
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: bcrsmatrix.hh:1903
~BCRSMatrix()
destructor
Definition: bcrsmatrix.hh:822
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
row_type & operator[](size_type i)
random access to the rows
Definition: bcrsmatrix.hh:547
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
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition: bcrsmatrix.hh:736
void setrowsize(size_type i, size_type s)
Set number of indices in row i to s.
Definition: bcrsmatrix.hh:1115
BCRSMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: bcrsmatrix.hh:1482
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
Iterator beforeBegin()
Definition: bcrsmatrix.hh:693
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(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 & operator-=(const BCRSMatrix &b)
Subtract the entries of another matrix from this one.
Definition: bcrsmatrix.hh:1565
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition: bcrsmatrix.hh:803
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:679
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
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:702
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: bcrsmatrix.hh:1875
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:498
BCRSMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: bcrsmatrix.hh:1510
BCRSMatrix & operator+=(const BCRSMatrix &b)
Add the entries of another matrix to this one.
Definition: bcrsmatrix.hh:1543
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
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1101
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: bcrsmatrix.hh:1858
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
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
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: bcrsmatrix.hh:1883
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
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
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition: bcrsmatrix.hh:1358
void setDataPointers()
Set data pointers for all rows.
Definition: bcrsmatrix.hh:2101
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
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:909
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
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 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
Thrown when the compression buffer used by the implicit BCRSMatrix construction is exhausted.
Definition: istlexception.hh:35
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: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
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:432
Default exception class for range errors.
Definition: exceptions.hh:252
Traits for type conversions and type information.
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:133
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:289
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
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: alignedallocator.hh:11
Implements a scalar matrix view wrapper around an existing scalar.
Implements a scalar vector view wrapper around an existing scalar.
Standard Dune debug streams.
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)