Dune Core Modules (2.4.1)

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"
24
29namespace Dune {
30
70 template<typename M>
71 struct MatrixDimension;
72
74
80 template<typename size_type>
82 {
84 double avg;
86 size_type maximum;
88 size_type overflow_total;
90
93 double mem_ratio;
94 };
95
97
109 template<class M_>
111 {
112
113 public:
114
116 typedef M_ Matrix;
117
120
122 typedef typename Matrix::size_type size_type;
123
125
131 {
132
133 public:
134
137 {
138 return _m.entry(_i,j);
139 }
140
141#ifndef DOXYGEN
142
144 : _m(m)
145 , _i(i)
146 {}
147
148#endif
149
150 private:
151
152 Matrix& _m;
153 size_type _i;
154
155 };
156
158
165 : _m(m)
166 {
167 if (m.buildMode() != Matrix::implicit)
168 DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix in implicit build mode");
169 if (m.buildStage() != Matrix::building)
170 DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix with set size that has not been compressed() yet");
171 }
172
174
188 ImplicitMatrixBuilder(Matrix& m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
189 : _m(m)
190 {
191 if (m.buildStage() != Matrix::notAllocated)
192 DUNE_THROW(BCRSMatrixError,"You can only set up a matrix for this ImplicitBuilder if it has no memory allocated yet");
193 m.setBuildMode(Matrix::implicit);
194 m.setImplicitBuildModeParameters(avg_cols_per_row,overflow_fraction);
195 m.setSize(rows,cols);
196 }
197
200 {
201 return row_object(_m,i);
202 }
203
205 size_type N() const
206 {
207 return _m.N();
208 }
209
211 size_type M() const
212 {
213 return _m.M();
214 }
215
216 private:
217
218 Matrix& _m;
219
220 };
221
411 template<class B, class A=std::allocator<B> >
413 {
414 friend struct MatrixDimension<BCRSMatrix>;
415 public:
429 built=3
430 };
431
432 //===== type definitions and constants
433
435 typedef typename B::field_type field_type;
436
438 typedef B block_type;
439
441 typedef A allocator_type;
442
445
447 typedef typename A::size_type size_type;
448
450 typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
451
453 enum {
455 blocklevel = B::blocklevel+1
456 };
457
492 unknown
493 };
494
495 //===== random access interface to rows of the matrix
496
499 {
500#ifdef DUNE_ISTL_WITH_CHECKING
501 if (build_mode == implicit && ready != built)
502 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
503 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
504 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
505#endif
506 return r[i];
507 }
508
511 {
512#ifdef DUNE_ISTL_WITH_CHECKING
513 if (build_mode == implicit && ready != built)
514 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
515 if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
516 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
517#endif
518 return r[i];
519 }
520
521
522 //===== iterator interface to rows of the matrix
523
525 template<class T>
527 : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
528 {
529
530 public:
532 typedef typename remove_const<T>::type ValueType;
533
536 friend class RealRowIterator<const ValueType>;
537 friend class RealRowIterator<ValueType>;
538
541 : p(_p), i(_i)
542 {}
543
546 : p(0), i(0)
547 {}
548
550 : p(it.p), i(it.i)
551 {}
552
553
556 {
557 return i;
558 }
559
560 std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
561 {
562 assert(other.p==p);
563 return (other.i-i);
564 }
565
566 std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
567 {
568 assert(other.p==p);
569 return (other.i-i);
570 }
571
573 bool equals (const RealRowIterator<ValueType>& other) const
574 {
575 assert(other.p==p);
576 return i==other.i;
577 }
578
580 bool equals (const RealRowIterator<const ValueType>& other) const
581 {
582 assert(other.p==p);
583 return i==other.i;
584 }
585
586 private:
588 void increment()
589 {
590 ++i;
591 }
592
594 void decrement()
595 {
596 --i;
597 }
598
599 void advance(std::ptrdiff_t diff)
600 {
601 i+=diff;
602 }
603
604 T& elementAt(std::ptrdiff_t diff) const
605 {
606 return p[i+diff];
607 }
608
610 row_type& dereference () const
611 {
612 return p[i];
613 }
614
615 row_type* p;
616 size_type i;
617 };
618
622
625 {
626 return Iterator(r,0);
627 }
628
631 {
632 return Iterator(r,n);
633 }
634
638 {
639 return Iterator(r,n-1);
640 }
641
645 {
646 return Iterator(r,-1);
647 }
648
651
654
658
659
662 {
663 return ConstIterator(r,0);
664 }
665
668 {
669 return ConstIterator(r,n);
670 }
671
675 {
676 return ConstIterator(r,n-1);
677 }
678
682 {
683 return ConstIterator(r,-1);
684 }
685
688
691
692 //===== constructors & resizers
693
694 // we use a negative overflowsize to indicate that the implicit
695 // mode parameters have not been set yet
696
699 : build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
700 allocationSize(0), r(0), a(0),
701 avg(0), overflowsize(-1.0)
702 {}
703
706 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
707 allocationSize(0), r(0), a(0),
708 avg(0), overflowsize(-1.0)
709 {
710 allocate(_n, _m, _nnz,true,false);
711 }
712
715 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
716 allocationSize(0), r(0), a(0),
717 avg(0), overflowsize(-1.0)
718 {
719 allocate(_n, _m,0,true,false);
720 }
721
723
733 BCRSMatrix (size_type _n, size_type _m, size_type _avg, double _overflowsize, BuildMode bm)
734 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
735 allocationSize(0), r(0), a(0),
736 avg(_avg), overflowsize(_overflowsize)
737 {
738 if (bm != implicit)
739 DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
740 // Prevent user from setting a negative overflowsize:
741 // 1) It doesn't make sense
742 // 2) We use a negative overflow value to indicate that the parameters
743 // have not been set yet
744 if (_overflowsize < 0.0)
745 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
746 implicit_allocate(_n,_m);
747 }
748
755 : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
756 allocationSize(0), r(0), a(0),
757 avg(Mat.avg), overflowsize(Mat.overflowsize)
758 {
759 if (!(Mat.ready == notAllocated || Mat.ready == built))
760 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
761
762 // deep copy in global array
763 size_type _nnz = Mat.nnz_;
764
765 // in case of row-wise allocation
766 if (_nnz<=0)
767 {
768 _nnz = 0;
769 for (size_type i=0; i<Mat.n; i++)
770 _nnz += Mat.r[i].getsize();
771 }
772
773 j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
774 allocate(Mat.n, Mat.m, _nnz, true, true);
775
776 // build window structure
778 }
779
782 {
783 deallocate();
784 }
785
791 {
792 if (ready == notAllocated)
793 {
794 build_mode = bm;
795 return;
796 }
797 if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
798 build_mode = bm;
799 else
800 DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
801 }
802
818 void setSize(size_type rows, size_type columns, size_type nnz=0)
819 {
820 // deallocate already setup memory
821 deallocate();
822
823 if (build_mode == implicit)
824 {
825 if (nnz>0)
826 DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");
827
828 // implicit allocates differently
829 implicit_allocate(rows,columns);
830 }
831 else
832 {
833 // allocate matrix memory
834 allocate(rows, columns, nnz, true, false);
835 }
836 }
837
846 void setImplicitBuildModeParameters(size_type _avg, double _overflow)
847 {
848 // Prevent user from setting a negative overflowsize:
849 // 1) It doesn't make sense
850 // 2) We use a negative overflow value to indicate that the parameters
851 // have not been set yet
852 if (_overflow < 0.0)
853 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
854
855 // make sure the parameters aren't changed after memory has been allocated
856 if (ready != notAllocated)
857 DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
858 avg = _avg;
859 overflowsize = _overflow;
860 }
861
869 {
870 // return immediately when self-assignment
871 if (&Mat==this) return *this;
872
873 if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
874 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");
875
876 // make it simple: ALWAYS throw away memory for a and j_
877 // and deallocate rows only if n != Mat.n
878 deallocate(n!=Mat.n);
879
880 // reallocate the rows if required
881 if (n>0 && n!=Mat.n) {
882 // free rows
883 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
884 rowAllocator_.destroy(riter);
885 rowAllocator_.deallocate(r,n);
886 }
887
888 nnz_ = Mat.nnz_;
889 if (nnz_ <= 0)
890 {
891 for (size_type i=0; i<Mat.n; i++)
892 nnz_ += Mat.r[i].getsize();
893 }
894
895 // allocate a, share j_
896 j_ = Mat.j_;
897 allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);
898
899 // build window structure
901 return *this;
902 }
903
906 {
907
908 if (!(ready == notAllocated || ready == built))
909 DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");
910
911 for (size_type i=0; i<n; i++) r[i] = k;
912 return *this;
913 }
914
915 //===== row-wise creation interface
916
919 {
920 public:
923 : Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
924 {
925 if (Mat.build_mode == unknown && Mat.ready == building)
926 {
927 Mat.build_mode = row_wise;
928 }
929 if (i==0 && Mat.ready != building)
930 DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
931 if(Mat.build_mode!=row_wise)
932 DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
933 }
934
937 {
938 // this should only be called if matrix is in creation
939 if (Mat.ready != building)
940 DUNE_THROW(BCRSMatrixError,"matrix already built up");
941
942 // row i is defined through the pattern
943 // get memory for the row and initialize the j_ array
944 // this depends on the allocation mode
945
946 // compute size of the row
947 size_type s = pattern.size();
948
949 if(s>0) {
950 // update number of nonzeroes including this row
951 nnz += s;
952
953 // alloc memory / set window
954 if (Mat.nnz_ > 0)
955 {
956 // memory is allocated in one long array
957
958 // check if that memory is sufficient
959 if (nnz > Mat.nnz_)
960 DUNE_THROW(BCRSMatrixError,"allocated nnz too small");
961
962 // set row i
963 Mat.r[i].set(s,nullptr,current_row.getindexptr());
964 current_row.setindexptr(current_row.getindexptr()+s);
965 }else{
966 // memory is allocated individually per row
967 // allocate and set row i
968 B* b = Mat.allocator_.allocate(s);
969 // use placement new to call constructor that allocates
970 // additional memory.
971 new (b) B[s];
972 size_type* j = Mat.sizeAllocator_.allocate(s);
973 Mat.r[i].set(s,b,j);
974 }
975 }else
976 // setup empty row
977 Mat.r[i].set(0,0,0);
978
979 // initialize the j array for row i from pattern
980 size_type k=0;
981 size_type *j = Mat.r[i].getindexptr();
982 for (typename PatternType::const_iterator it=pattern.begin(); it!=pattern.end(); ++it)
983 j[k++] = *it;
984
985 // now go to next row
986 i++;
987 pattern.clear();
988
989 // check if this was last row
990 if (i==Mat.n)
991 {
992 Mat.ready = built;
993 if(Mat.nnz_ > 0)
994 {
995 // Set nnz to the exact number of nonzero blocks inserted
996 // as some methods rely on it
997 Mat.nnz_ = nnz;
998 // allocate data array
999 Mat.allocateData();
1000 Mat.setDataPointers();
1001 }
1002 }
1003 // done
1004 return *this;
1005 }
1006
1008 bool operator!= (const CreateIterator& it) const
1009 {
1010 return (i!=it.i) || (&Mat!=&it.Mat);
1011 }
1012
1014 bool operator== (const CreateIterator& it) const
1015 {
1016 return (i==it.i) && (&Mat==&it.Mat);
1017 }
1018
1021 {
1022 return i;
1023 }
1024
1027 {
1028 pattern.insert(j);
1029 }
1030
1033 {
1034 if (pattern.find(j)!=pattern.end())
1035 return true;
1036 else
1037 return false;
1038 }
1045 {
1046 return pattern.size();
1047 }
1048
1049 private:
1050 BCRSMatrix& Mat; // the matrix we are defining
1051 size_type i; // current row to be defined
1052 size_type nnz; // count total number of nonzeros
1053 typedef std::set<size_type,std::less<size_type> > PatternType;
1054 PatternType pattern; // used to compile entries in a row
1055 row_type current_row; // row pointing to the current row to setup
1056 };
1057
1059 friend class CreateIterator;
1060
1063 {
1064 return CreateIterator(*this,0);
1065 }
1066
1069 {
1070 return CreateIterator(*this,n);
1071 }
1072
1073
1074 //===== random creation interface
1075
1078 {
1079 if (build_mode!=random)
1080 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1081 if (ready != building)
1082 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1083
1084 r[i].setsize(s);
1085 }
1086
1089 {
1090#ifdef DUNE_ISTL_WITH_CHECKING
1091 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
1092 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
1093#endif
1094 return r[i].getsize();
1095 }
1096
1099 {
1100 if (build_mode!=random)
1101 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1102 if (ready != building)
1103 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1104
1105 r[i].setsize(r[i].getsize()+s);
1106 }
1107
1110 {
1111 if (build_mode!=random)
1112 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1113 if (ready != building)
1114 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1115
1116 // compute total size, check positivity
1117 size_type total=0;
1118 for (size_type i=0; i<n; i++)
1119 {
1120 total += r[i].getsize();
1121 }
1122
1123 if(nnz_ == 0)
1124 // allocate/check memory
1125 allocate(n,m,total,false,false);
1126 else if(nnz_ < total)
1127 DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
1128 <<"sufficient for calculated nonzeros ("<<total<<"! ");
1129
1130 // set the window pointers correctly
1132
1133 // initialize j_ array with m (an invalid column index)
1134 // this indicates an unused entry
1135 for (size_type k=0; k<nnz_; k++)
1136 j_.get()[k] = m;
1137 ready = rowSizesBuilt;
1138 }
1139
1141
1152 {
1153 if (build_mode!=random)
1154 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1155 if (ready==built)
1156 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1157 if (ready==building)
1158 DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
1159 if (ready==notAllocated)
1160 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1161
1162 if (col >= m)
1163 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1164
1165 // get row range
1166 size_type* const first = r[row].getindexptr();
1167 size_type* const last = first + r[row].getsize();
1168
1169 // find correct insertion position for new column index
1170 size_type* pos = std::lower_bound(first,last,col);
1171
1172 // check if index is already in row
1173 if (pos!=last && *pos == col) return;
1174
1175 // find end of already inserted column indices
1176 size_type* end = std::lower_bound(pos,last,m);
1177 if (end==last)
1178 DUNE_THROW(BCRSMatrixError,"row is too small");
1179
1180 // insert new column index at correct position
1181 std::copy_backward(pos,end,end+1);
1182 *pos = col;
1183 }
1184
1186
1193 template<typename It>
1194 void setIndices(size_type row, It begin, It end)
1195 {
1196 size_type row_size = r[row].size();
1197 size_type* col_begin = r[row].getindexptr();
1198 size_type* col_end;
1199 // consistency check between allocated row size and number of passed column indices
1200 if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
1201 DUNE_THROW(BCRSMatrixError,"Given size of row " << row
1202 << " (" << row_size
1203 << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
1204 std::sort(col_begin,col_end);
1205 }
1206
1209 {
1210 if (build_mode!=random)
1211 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1212 if (ready==built)
1213 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1214 if (ready==building)
1215 DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
1216 if (ready==notAllocated)
1217 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1218
1219 // check if there are undefined indices
1220 RowIterator endi=end();
1221 for (RowIterator i=begin(); i!=endi; ++i)
1222 {
1223 ColIterator endj = (*i).end();
1224 for (ColIterator j=(*i).begin(); j!=endj; ++j) {
1225 if (j.index() >= m) {
1226 dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
1227 <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
1228 r[i.index()].setsize(j.offset());
1229 break;
1230 }
1231 }
1232 }
1233
1234 allocateData();
1236
1237 // if not, set matrix to built
1238 ready = built;
1239 }
1240
1241 //===== implicit creation interface
1242
1244
1256 {
1257#ifdef DUNE_ISTL_WITH_CHECKING
1258 if (build_mode!=implicit)
1259 DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1260 if (ready==built)
1261 DUNE_THROW(BCRSMatrixError,"matrix already built up, use operator[] for entry access now");
1262 if (ready==notAllocated)
1263 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1264 if (ready!=building)
1265 DUNE_THROW(InvalidStateException,"You may only use entry() during the 'building' stage");
1266
1267 if (row >= n)
1268 DUNE_THROW(BCRSMatrixError,"row index exceeds matrix size");
1269 if (col >= m)
1270 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1271#endif
1272
1273 size_type* begin = r[row].getindexptr();
1274 size_type* end = begin + r[row].getsize();
1275
1276 size_type* pos = std::find(begin, end, col);
1277
1278 //treat the case that there was a match in the array
1279 if (pos != end)
1280 if (*pos == col)
1281 {
1282 std::ptrdiff_t offset = pos - r[row].getindexptr();
1283 B* aptr = r[row].getptr() + offset;
1284
1285 return *aptr;
1286 }
1287
1288 //determine whether overflow has to be taken into account or not
1289 if (r[row].getsize() == avg)
1290 return overflow[std::make_pair(row,col)];
1291 else
1292 {
1293 //modify index array
1294 *end = col;
1295
1296 //do simulatenous operations on data array a
1297 std::ptrdiff_t offset = end - r[row].getindexptr();
1298 B* apos = r[row].getptr() + offset;
1299
1300 //increase rowsize
1301 r[row].setsize(r[row].getsize()+1);
1302
1303 //return reference to the newly created entry
1304 return *apos;
1305 }
1306 }
1307
1309
1320 {
1321 if (build_mode!=implicit)
1322 DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1323 if (ready==built)
1324 DUNE_THROW(BCRSMatrixError,"matrix already built up, no more need for compression");
1325 if (ready==notAllocated)
1326 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1327 if (ready!=building)
1328 DUNE_THROW(InvalidStateException,"You may only call compress() at the end of the 'building' stage");
1329
1330 //calculate statistics
1332 stats.overflow_total = overflow.size();
1333 stats.maximum = 0;
1334
1335 //get insertion iterators pointing to one before start (for later use of ++it)
1336 size_type* jiit = j_.get();
1337 B* aiit = a;
1338
1339 //get iterator to the smallest overflow element
1340 typename OverflowType::iterator oit = overflow.begin();
1341
1342 //store a copy of index pointers on which to perform sortation
1343 std::vector<size_type*> perm;
1344
1345 //iterate over all rows and copy elements into their position in the compressed array
1346 for (size_type i=0; i<n; i++)
1347 {
1348 //get old pointers into a and j and size without overflow changes
1349 size_type* begin = r[i].getindexptr();
1350 //B* apos = r[i].getptr();
1351 size_type size = r[i].getsize();
1352
1353 perm.resize(size);
1354
1355 typename std::vector<size_type*>::iterator it = perm.begin();
1356 for (size_type* iit = begin; iit < begin + size; ++iit, ++it)
1357 *it = iit;
1358
1359 //sort permutation array
1360 std::sort(perm.begin(),perm.end(),PointerCompare<size_type>());
1361
1362 //change row window pointer to their new positions
1363 r[i].setindexptr(jiit);
1364 r[i].setptr(aiit);
1365
1366 for (it = perm.begin(); it != perm.end(); ++it)
1367 {
1368 //check whether there are elements in the overflow area which take precedence
1369 while ((oit!=overflow.end()) && (oit->first < std::make_pair(i,**it)))
1370 {
1371 //check whether there is enough memory to write to
1372 if (jiit > begin)
1374 "Allocated memory for BCRSMatrix exhausted during compress()!"
1375 "Please increase either the average number of entries per row or the overflow fraction."
1376 );
1377 //copy an element from the overflow area to the insertion position in a and j
1378 *jiit = oit->first.second;
1379 ++jiit;
1380 *aiit = oit->second;
1381 ++aiit;
1382 ++oit;
1383 r[i].setsize(r[i].getsize()+1);
1384 }
1385
1386 //check whether there is enough memory to write to
1387 if (jiit > begin)
1389 "Allocated memory for BCRSMatrix exhausted during compress()!"
1390 "Please increase either the average number of entries per row or the overflow fraction."
1391 );
1392
1393 //copy element from array
1394 *jiit = **it;
1395 ++jiit;
1396 B* apos = *it - j_.get() + a;
1397 *aiit = *apos;
1398 ++aiit;
1399 }
1400
1401 //copy remaining elements from the overflow area
1402 while ((oit!=overflow.end()) && (oit->first.first == i))
1403 {
1404 //check whether there is enough memory to write to
1405 if (jiit > begin)
1407 "Allocated memory for BCRSMatrix exhausted during compress()!"
1408 "Please increase either the average number of entries per row or the overflow fraction."
1409 );
1410
1411 //copy and element from the overflow area to the insertion position in a and j
1412 *jiit = oit->first.second;
1413 ++jiit;
1414 *aiit = oit->second;
1415 ++aiit;
1416 ++oit;
1417 r[i].setsize(r[i].getsize()+1);
1418 }
1419
1420 // update maximum row size
1421 if (r[i].getsize()>stats.maximum)
1422 stats.maximum = r[i].getsize();
1423 }
1424
1425 // overflow area may be cleared
1426 overflow.clear();
1427
1428 //determine average number of entries and memory usage
1429 std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
1430 nnz_ = diff;
1431 stats.avg = (double) (nnz_) / (double) n;
1432 stats.mem_ratio = (double) (nnz_) / (double) allocationSize;
1433
1434 //matrix is now built
1435 ready = built;
1436
1437 return stats;
1438 }
1439
1440 //===== vector space arithmetic
1441
1444 {
1445#ifdef DUNE_ISTL_WITH_CHECKING
1446 if (ready != built)
1447 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1448#endif
1449
1450 if (nnz_ > 0)
1451 {
1452 // process 1D array
1453 for (size_type i=0; i<nnz_; i++)
1454 a[i] *= k;
1455 }
1456 else
1457 {
1458 RowIterator endi=end();
1459 for (RowIterator i=begin(); i!=endi; ++i)
1460 {
1461 ColIterator endj = (*i).end();
1462 for (ColIterator j=(*i).begin(); j!=endj; ++j)
1463 (*j) *= k;
1464 }
1465 }
1466
1467 return *this;
1468 }
1469
1472 {
1473#ifdef DUNE_ISTL_WITH_CHECKING
1474 if (ready != built)
1475 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1476#endif
1477
1478 if (nnz_ > 0)
1479 {
1480 // process 1D array
1481 for (size_type i=0; i<nnz_; i++)
1482 a[i] /= k;
1483 }
1484 else
1485 {
1486 RowIterator endi=end();
1487 for (RowIterator i=begin(); i!=endi; ++i)
1488 {
1489 ColIterator endj = (*i).end();
1490 for (ColIterator j=(*i).begin(); j!=endj; ++j)
1491 (*j) /= k;
1492 }
1493 }
1494
1495 return *this;
1496 }
1497
1498
1505 {
1506#ifdef DUNE_ISTL_WITH_CHECKING
1507 if (ready != built || b.ready != built)
1508 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1509 if(N()!=b.N() || M() != b.M())
1510 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1511#endif
1512 RowIterator endi=end();
1513 ConstRowIterator j=b.begin();
1514 for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1515 i->operator+=(*j);
1516 }
1517
1518 return *this;
1519 }
1520
1527 {
1528#ifdef DUNE_ISTL_WITH_CHECKING
1529 if (ready != built || b.ready != built)
1530 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1531 if(N()!=b.N() || M() != b.M())
1532 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1533#endif
1534 RowIterator endi=end();
1535 ConstRowIterator j=b.begin();
1536 for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1537 i->operator-=(*j);
1538 }
1539
1540 return *this;
1541 }
1542
1552 {
1553#ifdef DUNE_ISTL_WITH_CHECKING
1554 if (ready != built || b.ready != built)
1555 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1556 if(N()!=b.N() || M() != b.M())
1557 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1558#endif
1559 RowIterator endi=end();
1560 ConstRowIterator j=b.begin();
1561 for(RowIterator i=begin(); i!=endi; ++i, ++j)
1562 i->axpy(alpha, *j);
1563
1564 return *this;
1565 }
1566
1567 //===== linear maps
1568
1570 template<class X, class Y>
1571 void mv (const X& x, Y& y) const
1572 {
1573#ifdef DUNE_ISTL_WITH_CHECKING
1574 if (ready != built)
1575 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1576 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,
1577 "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
1578 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,
1579 "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
1580#endif
1581 ConstRowIterator endi=end();
1582 for (ConstRowIterator i=begin(); i!=endi; ++i)
1583 {
1584 y[i.index()]=0;
1585 ConstColIterator endj = (*i).end();
1586 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1587 (*j).umv(x[j.index()],y[i.index()]);
1588 }
1589 }
1590
1592 template<class X, class Y>
1593 void umv (const X& x, Y& y) const
1594 {
1595#ifdef DUNE_ISTL_WITH_CHECKING
1596 if (ready != built)
1597 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1598 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1599 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1600#endif
1601 ConstRowIterator endi=end();
1602 for (ConstRowIterator i=begin(); i!=endi; ++i)
1603 {
1604 ConstColIterator endj = (*i).end();
1605 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1606 (*j).umv(x[j.index()],y[i.index()]);
1607 }
1608 }
1609
1611 template<class X, class Y>
1612 void mmv (const X& x, Y& y) const
1613 {
1614#ifdef DUNE_ISTL_WITH_CHECKING
1615 if (ready != built)
1616 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1617 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1618 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1619#endif
1620 ConstRowIterator endi=end();
1621 for (ConstRowIterator i=begin(); i!=endi; ++i)
1622 {
1623 ConstColIterator endj = (*i).end();
1624 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1625 (*j).mmv(x[j.index()],y[i.index()]);
1626 }
1627 }
1628
1630 template<typename F, class X, class Y>
1631 void usmv (F&& alpha, const X& x, Y& y) const
1632 {
1633#ifdef DUNE_ISTL_WITH_CHECKING
1634 if (ready != built)
1635 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1636 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1637 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1638#endif
1639 ConstRowIterator endi=end();
1640 for (ConstRowIterator i=begin(); i!=endi; ++i)
1641 {
1642 ConstColIterator endj = (*i).end();
1643 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1644 (*j).usmv(alpha,x[j.index()],y[i.index()]);
1645 }
1646 }
1647
1649 template<class X, class Y>
1650 void mtv (const X& x, Y& y) const
1651 {
1652#ifdef DUNE_ISTL_WITH_CHECKING
1653 if (ready != built)
1654 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1655 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1656 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1657#endif
1658 for(size_type i=0; i<y.N(); ++i)
1659 y[i]=0;
1660 umtv(x,y);
1661 }
1662
1664 template<class X, class Y>
1665 void umtv (const X& x, Y& y) const
1666 {
1667#ifdef DUNE_ISTL_WITH_CHECKING
1668 if (ready != built)
1669 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1670 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1671 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1672#endif
1673 ConstRowIterator endi=end();
1674 for (ConstRowIterator i=begin(); i!=endi; ++i)
1675 {
1676 ConstColIterator endj = (*i).end();
1677 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1678 (*j).umtv(x[i.index()],y[j.index()]);
1679 }
1680 }
1681
1683 template<class X, class Y>
1684 void mmtv (const X& x, Y& y) const
1685 {
1686#ifdef DUNE_ISTL_WITH_CHECKING
1687 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1688 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1689#endif
1690 ConstRowIterator endi=end();
1691 for (ConstRowIterator i=begin(); i!=endi; ++i)
1692 {
1693 ConstColIterator endj = (*i).end();
1694 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1695 (*j).mmtv(x[i.index()],y[j.index()]);
1696 }
1697 }
1698
1700 template<class X, class Y>
1701 void usmtv (const field_type& alpha, const X& x, Y& y) const
1702 {
1703#ifdef DUNE_ISTL_WITH_CHECKING
1704 if (ready != built)
1705 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1706 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1707 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1708#endif
1709 ConstRowIterator endi=end();
1710 for (ConstRowIterator i=begin(); i!=endi; ++i)
1711 {
1712 ConstColIterator endj = (*i).end();
1713 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1714 (*j).usmtv(alpha,x[i.index()],y[j.index()]);
1715 }
1716 }
1717
1719 template<class X, class Y>
1720 void umhv (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=begin(); i!=endi; ++i)
1730 {
1731 ConstColIterator endj = (*i).end();
1732 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1733 (*j).umhv(x[i.index()],y[j.index()]);
1734 }
1735 }
1736
1738 template<class X, class Y>
1739 void mmhv (const X& x, Y& y) const
1740 {
1741#ifdef DUNE_ISTL_WITH_CHECKING
1742 if (ready != built)
1743 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1744 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1745 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1746#endif
1747 ConstRowIterator endi=end();
1748 for (ConstRowIterator i=begin(); i!=endi; ++i)
1749 {
1750 ConstColIterator endj = (*i).end();
1751 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1752 (*j).mmhv(x[i.index()],y[j.index()]);
1753 }
1754 }
1755
1757 template<class X, class Y>
1758 void usmhv (const field_type& alpha, const X& x, Y& y) const
1759 {
1760#ifdef DUNE_ISTL_WITH_CHECKING
1761 if (ready != built)
1762 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1763 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1764 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1765#endif
1766 ConstRowIterator endi=end();
1767 for (ConstRowIterator i=begin(); i!=endi; ++i)
1768 {
1769 ConstColIterator endj = (*i).end();
1770 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1771 (*j).usmhv(alpha,x[i.index()],y[j.index()]);
1772 }
1773 }
1774
1775
1776 //===== norms
1777
1779 typename FieldTraits<field_type>::real_type frobenius_norm2 () const
1780 {
1781#ifdef DUNE_ISTL_WITH_CHECKING
1782 if (ready != built)
1783 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1784#endif
1785
1786 double sum=0;
1787
1788 ConstRowIterator endi=end();
1789 for (ConstRowIterator i=begin(); i!=endi; ++i)
1790 {
1791 ConstColIterator endj = (*i).end();
1792 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1793 sum += (*j).frobenius_norm2();
1794 }
1795
1796 return sum;
1797 }
1798
1800 typename FieldTraits<field_type>::real_type frobenius_norm () const
1801 {
1802 return sqrt(frobenius_norm2());
1803 }
1804
1806 typename FieldTraits<field_type>::real_type infinity_norm () const
1807 {
1808 if (ready != built)
1809 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1810
1811 double max=0;
1812 ConstRowIterator endi=end();
1813 for (ConstRowIterator i=begin(); i!=endi; ++i)
1814 {
1815 double sum=0;
1816 ConstColIterator endj = (*i).end();
1817 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1818 sum += (*j).infinity_norm();
1819 max = std::max(max,sum);
1820 }
1821 return max;
1822 }
1823
1825 typename FieldTraits<field_type>::real_type infinity_norm_real () const
1826 {
1827#ifdef DUNE_ISTL_WITH_CHECKING
1828 if (ready != built)
1829 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1830#endif
1831
1832 double max=0;
1833 ConstRowIterator endi=end();
1834 for (ConstRowIterator i=begin(); i!=endi; ++i)
1835 {
1836 double sum=0;
1837 ConstColIterator endj = (*i).end();
1838 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1839 sum += (*j).infinity_norm_real();
1840 max = std::max(max,sum);
1841 }
1842 return max;
1843 }
1844
1845
1846 //===== sizes
1847
1849 size_type N () const
1850 {
1851 return n;
1852 }
1853
1855 size_type M () const
1856 {
1857 return m;
1858 }
1859
1862 {
1863 return nnz_;
1864 }
1865
1868 {
1869 return ready;
1870 }
1871
1874 {
1875 return build_mode;
1876 }
1877
1878 //===== query
1879
1881 bool exists (size_type i, size_type j) const
1882 {
1883#ifdef DUNE_ISTL_WITH_CHECKING
1884 if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
1885 if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
1886#endif
1887 if (r[i].size() && r[i].find(j)!=r[i].end())
1888 return true;
1889 else
1890 return false;
1891 }
1892
1893
1894 protected:
1895 // state information
1896 BuildMode build_mode; // row wise or whole matrix
1897 BuildStage ready; // indicate the stage the matrix building is in
1898
1899 // The allocator used for memory management
1900 typename A::template rebind<B>::other allocator_;
1901
1902 typename A::template rebind<row_type>::other rowAllocator_;
1903
1904 typename A::template rebind<size_type>::other sizeAllocator_;
1905
1906 // size of the matrix
1907 size_type n; // number of rows
1908 size_type m; // number of columns
1909 size_type nnz_; // number of nonzeroes contained in the matrix
1910 size_type allocationSize; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize
1911 // zero means that memory is allocated separately for each row.
1912
1913 // the rows are dynamically allocated
1914 row_type* r; // [n] the individual rows having pointers into a,j arrays
1915
1916 // dynamically allocated memory
1917 B* a; // [allocationSize] non-zero entries of the matrix in row-wise ordering
1918 // If a single array of column indices is used, it can be shared
1919 // between different matrices with the same sparsity pattern
1920 std::shared_ptr<size_type> j_; // [allocationSize] column indices of entries
1921
1922 // additional data is needed in implicit buildmode
1923 size_type avg;
1924 double overflowsize;
1925
1926 typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
1927 OverflowType overflow;
1928
1929 void setWindowPointers(ConstRowIterator row)
1930 {
1931 row_type current_row(a,j_.get(),0); // Pointers to current row data
1932 for (size_type i=0; i<n; i++, ++row) {
1933 // set row i
1934 size_type s = row->getsize();
1935
1936 if (s>0) {
1937 // setup pointers and size
1938 r[i].set(s,current_row.getptr(), current_row.getindexptr());
1939 // update pointer for next row
1940 current_row.setptr(current_row.getptr()+s);
1941 current_row.setindexptr(current_row.getindexptr()+s);
1942 } else{
1943 // empty row
1944 r[i].set(0,0,0);
1945 }
1946 }
1947 }
1948
1950
1955 {
1956 size_type* jptr = j_.get();
1957 for (size_type i=0; i<n; ++i, ++row) {
1958 // set row i
1959 size_type s = row->getsize();
1960
1961 if (s>0) {
1962 // setup pointers and size
1963 r[i].setsize(s);
1964 r[i].setindexptr(jptr);
1965 } else{
1966 // empty row
1967 r[i].set(0,0,0);
1968 }
1969
1970 // advance position in global array
1971 jptr += s;
1972 }
1973 }
1974
1976
1981 {
1982 B* aptr = a;
1983 for (size_type i=0; i<n; ++i) {
1984 // set row i
1985 if (r[i].getsize() > 0) {
1986 // setup pointers and size
1987 r[i].setptr(aptr);
1988 } else{
1989 // empty row
1990 r[i].set(0,0,0);
1991 }
1992
1993 // advance position in global array
1994 aptr += r[i].getsize();
1995 }
1996 }
1997
2000 {
2001 setWindowPointers(Mat.begin());
2002
2003 // copy data
2004 for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
2005
2006 // finish off
2007 build_mode = row_wise; // dummy
2008 ready = built;
2009 }
2010
2016 void deallocate(bool deallocateRows=true)
2017 {
2018
2019 if (notAllocated)
2020 return;
2021
2022 if (allocationSize>0)
2023 {
2024 // a,j_ have been allocated as one long vector
2025 j_.reset();
2026 if (a)
2027 {
2028 for(B *aiter=a+(allocationSize-1), *aend=a-1; aiter!=aend; --aiter)
2029 allocator_.destroy(aiter);
2030 allocator_.deallocate(a,allocationSize);
2031 a = nullptr;
2032 }
2033 }
2034 else if (r)
2035 {
2036 // check if memory for rows have been allocated individually
2037 for (size_type i=0; i<n; i++)
2038 if (r[i].getsize()>0)
2039 {
2040 for (B *col=r[i].getptr()+(r[i].getsize()-1),
2041 *colend = r[i].getptr()-1; col!=colend; --col) {
2042 allocator_.destroy(col);
2043 }
2044 sizeAllocator_.deallocate(r[i].getindexptr(),1);
2045 allocator_.deallocate(r[i].getptr(),1);
2046 // clear out row data in case we don't want to deallocate the rows
2047 // otherwise we might run into a double free problem here later
2048 r[i].set(0,nullptr,nullptr);
2049 }
2050 }
2051
2052 // deallocate the rows
2053 if (n>0 && deallocateRows && r) {
2054 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
2055 rowAllocator_.destroy(riter);
2056 rowAllocator_.deallocate(r,n);
2057 r = nullptr;
2058 }
2059
2060 // Mark matrix as not built at all.
2061 ready=notAllocated;
2062
2063 }
2064
2067 {
2068 typename A::template rebind<size_type>::other& sizeAllocator_;
2069
2070 public:
2071 Deallocator(typename A::template rebind<size_type>::other& sizeAllocator)
2072 : sizeAllocator_(sizeAllocator)
2073 {}
2074
2075 void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); }
2076 };
2077
2078
2096 void allocate(size_type rows, size_type columns, size_type allocationSize_, bool allocateRows, bool allocate_data)
2097 {
2098 // Store size
2099 n = rows;
2100 m = columns;
2101 nnz_ = allocationSize_;
2102 allocationSize = allocationSize_;
2103
2104 // allocate rows
2105 if(allocateRows) {
2106 if (n>0) {
2107 if (r)
2108 DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
2109 r = rowAllocator_.allocate(rows);
2110 }else{
2111 r = 0;
2112 }
2113 }
2114
2115 // allocate a and j_ array
2116 if (allocate_data)
2117 allocateData();
2118 if (allocationSize>0) {
2119 // allocate column indices only if not yet present (enable sharing)
2120 if (!j_.get())
2121 j_.reset(sizeAllocator_.allocate(allocationSize),Deallocator(sizeAllocator_));
2122 }else{
2123 j_.reset();
2124 for(row_type* ri=r; ri!=r+rows; ++ri)
2125 rowAllocator_.construct(ri, row_type());
2126 }
2127
2128 // Mark the matrix as not built.
2129 ready = building;
2130 }
2131
2132 void allocateData()
2133 {
2134 if (a)
2135 DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
2136 if (allocationSize>0) {
2137 a = allocator_.allocate(allocationSize);
2138 // use placement new to call constructor that allocates
2139 // additional memory.
2140 new (a) B[allocationSize];
2141 } else {
2142 a = nullptr;
2143 }
2144 }
2145
2152 {
2153 if (build_mode != implicit)
2154 DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
2155 if (ready != notAllocated)
2156 DUNE_THROW(InvalidStateException,"memory has already been allocated");
2157
2158 // check to make sure the user has actually set the parameters
2159 if (overflowsize < 0)
2160 DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
2161 //calculate size of overflow region, add buffer for row sort!
2162 size_type osize = (size_type) (_n*avg)*overflowsize + 4*avg;
2163 allocationSize = _n*avg + osize;
2164
2165 allocate(_n, _m, allocationSize,true,true);
2166
2167 //set row pointers correctly
2168 size_type* jptr = j_.get() + osize;
2169 B* aptr = a + osize;
2170 for (size_type i=0; i<n; i++)
2171 {
2172 r[i].set(0,aptr,jptr);
2173 jptr = jptr + avg;
2174 aptr = aptr + avg;
2175 }
2176
2177 ready = building;
2178 }
2179 };
2180
2181
2184} // end namespace
2185
2186#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:919
bool operator==(const CreateIterator &it) const
equality
Definition: bcrsmatrix.hh:1014
CreateIterator & operator++()
prefix increment
Definition: bcrsmatrix.hh:936
size_type index() const
dereferencing
Definition: bcrsmatrix.hh:1020
bool operator!=(const CreateIterator &it) const
inequality
Definition: bcrsmatrix.hh:1008
CreateIterator(BCRSMatrix &_Mat, size_type _i)
constructor
Definition: bcrsmatrix.hh:922
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1026
bool contains(size_type j)
return true if column index is in row
Definition: bcrsmatrix.hh:1032
size_type size() const
Get the current row size.
Definition: bcrsmatrix.hh:1044
Class used by shared_ptr to deallocate memory using the proper allocator.
Definition: bcrsmatrix.hh:2067
Iterator access to matrix rows
Definition: bcrsmatrix.hh:528
RealRowIterator()
empty constructor, use with care!
Definition: bcrsmatrix.hh:545
bool equals(const RealRowIterator< ValueType > &other) const
equality
Definition: bcrsmatrix.hh:573
bool equals(const RealRowIterator< const ValueType > &other) const
equality
Definition: bcrsmatrix.hh:580
RealRowIterator(row_type *_p, size_type _i)
constructor
Definition: bcrsmatrix.hh:540
size_type index() const
return index
Definition: bcrsmatrix.hh:555
remove_const< T >::type ValueType
The unqualified value type.
Definition: bcrsmatrix.hh:532
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:413
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: bcrsmatrix.hh:1881
BuildStage buildStage() const
The current build stage of the matrix.
Definition: bcrsmatrix.hh:1867
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:624
friend class CreateIterator
allow CreateIterator to access internal data
Definition: bcrsmatrix.hh:1059
void copyWindowStructure(const BCRSMatrix &Mat)
Copy the window structure from another matrix.
Definition: bcrsmatrix.hh:1999
B & entry(size_type row, size_type col)
Returns reference to entry (row,col) of the matrix.
Definition: bcrsmatrix.hh:1255
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: bcrsmatrix.hh:1739
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: bcrsmatrix.hh:1758
void umtv(const X &x, Y &y) const
y += A^T x
Definition: bcrsmatrix.hh:1665
RealRowIterator< const row_type > const_iterator
The const iterator over the matrix rows.
Definition: bcrsmatrix.hh:656
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:733
B::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:435
BCRSMatrix & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition: bcrsmatrix.hh:1551
~BCRSMatrix()
destructor
Definition: bcrsmatrix.hh:781
void deallocate(bool deallocateRows=true)
deallocate memory of the matrix.
Definition: bcrsmatrix.hh:2016
Iterator RowIterator
rename the iterators for easier access
Definition: bcrsmatrix.hh:650
row_type & operator[](size_type i)
random access to the rows
Definition: bcrsmatrix.hh:498
BCRSMatrix()
an empty matrix
Definition: bcrsmatrix.hh:698
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1109
void incrementrowsize(size_type i, size_type s=1)
increment size of row i by s (1 by default)
Definition: bcrsmatrix.hh:1098
void mtv(const X &x, Y &y) const
y = A^T x
Definition: bcrsmatrix.hh:1650
void umhv(const X &x, Y &y) const
y += A^H x
Definition: bcrsmatrix.hh:1720
size_type nonzeroes() const
number of blocks that are stored (the number of blocks that possibly are nonzero)
Definition: bcrsmatrix.hh:1861
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition: bcrsmatrix.hh:687
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:2096
void setrowsize(size_type i, size_type s)
set number of indices in row i to s
Definition: bcrsmatrix.hh:1077
BCRSMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: bcrsmatrix.hh:1443
RealRowIterator< row_type > iterator
The iterator over the (mutable matrix rows.
Definition: bcrsmatrix.hh:620
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: bcrsmatrix.hh:1701
ConstIterator beforeBegin() const
Definition: bcrsmatrix.hh:681
Iterator beforeBegin()
Definition: bcrsmatrix.hh:644
BuildMode
we support two modes
Definition: bcrsmatrix.hh:459
@ implicit
Build entries randomly with an educated guess on entries per row.
Definition: bcrsmatrix.hh:488
@ unknown
Build mode not set!
Definition: bcrsmatrix.hh:492
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:479
@ row_wise
Build in a row-wise manner.
Definition: bcrsmatrix.hh:470
BCRSMatrix(size_type _n, size_type _m, size_type _nnz, BuildMode bm)
matrix with known number of nonzeroes
Definition: bcrsmatrix.hh:705
::Dune::CompressionStatistics< size_type > CompressionStatistics
The type for the statistics object returned by compress()
Definition: bcrsmatrix.hh:450
BCRSMatrix & operator-=(const BCRSMatrix &b)
Substract the entries of another matrix to this one.
Definition: bcrsmatrix.hh:1526
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition: bcrsmatrix.hh:754
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:630
void setIndices(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition: bcrsmatrix.hh:1194
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1151
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:653
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: bcrsmatrix.hh:1800
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:447
BCRSMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: bcrsmatrix.hh:1471
BCRSMatrix & operator+=(const BCRSMatrix &b)
Add the entries of another matrix to this one.
Definition: bcrsmatrix.hh:1504
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1068
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: bcrsmatrix.hh:1779
Iterator beforeEnd()
Definition: bcrsmatrix.hh:637
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:690
void usmv(F &&alpha, const X &x, Y &y) const
y += alpha A x
Definition: bcrsmatrix.hh:1631
CompressedBlockVectorWindow< B, A > row_type
implement row_type with compressed vector
Definition: bcrsmatrix.hh:444
size_type getrowsize(size_type i) const
get current number of indices in row i
Definition: bcrsmatrix.hh:1088
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1855
@ blocklevel
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:455
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1062
BuildStage
Definition: bcrsmatrix.hh:416
@ rowSizesBuilt
The row sizes of the matrix are known.
Definition: bcrsmatrix.hh:427
@ built
The matrix structure is fully built.
Definition: bcrsmatrix.hh:429
@ notbuilt
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:418
@ notAllocated
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:420
@ building
Matrix is currently being built, some memory has been allocated, build mode and size are fixed.
Definition: bcrsmatrix.hh:422
BuildMode buildMode() const
The currently selected build mode of the matrix.
Definition: bcrsmatrix.hh:1873
void mmv(const X &x, Y &y) const
y -= A x
Definition: bcrsmatrix.hh:1612
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: bcrsmatrix.hh:1825
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1571
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:438
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: bcrsmatrix.hh:1684
void umv(const X &x, Y &y) const
y += A x
Definition: bcrsmatrix.hh:1593
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:2151
BCRSMatrix(size_type _n, size_type _m, BuildMode bm)
matrix with unknown number of nonzeroes
Definition: bcrsmatrix.hh:714
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1208
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition: bcrsmatrix.hh:1319
void setDataPointers()
Set data pointers for all rows.
Definition: bcrsmatrix.hh:1980
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1849
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:790
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:818
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: bcrsmatrix.hh:1806
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:868
void setColumnPointers(ConstRowIterator row)
Copy row sizes from iterator range starting at row and set column index pointers for all rows.
Definition: bcrsmatrix.hh:1954
ConstIterator end() const
Get const iterator to one beyond last row.
Definition: bcrsmatrix.hh:667
ConstIterator begin() const
Get const iterator to first row.
Definition: bcrsmatrix.hh:661
void setImplicitBuildModeParameters(size_type _avg, double _overflow)
Set parameters needed for creation in implicit build mode.
Definition: bcrsmatrix.hh:846
A allocator_type
export the allocator type
Definition: bcrsmatrix.hh:441
ConstIterator beforeEnd() const
Definition: bcrsmatrix.hh:674
compressed_block_vector_unmanaged< B, A >::ConstIterator ConstIterator
make iterators available as types
Definition: bvector.hh:952
void set(size_type _n, B *_p, size_type *_j)
set size and pointer
Definition: bvector.hh:1024
size_type * getindexptr()
get pointer
Definition: bvector.hh:1056
void setsize(size_type _n)
set size only
Definition: bvector.hh:1032
compressed_block_vector_unmanaged< B, A >::Iterator Iterator
make iterators available as types
Definition: bvector.hh:949
void setptr(B *_p)
set pointer only
Definition: bvector.hh:1038
B * getptr()
get pointer
Definition: bvector.hh:1050
size_type getsize() const
get size
Definition: bvector.hh:1073
void setindexptr(size_type *_j)
set pointer only
Definition: bvector.hh:1044
Proxy row object for entry access.
Definition: bcrsmatrix.hh:131
block_type & operator[](size_type j) const
Returns entry in column j.
Definition: bcrsmatrix.hh:136
A wrapper for uniform access to the BCRSMatrix during and after the build stage in implicit build mod...
Definition: bcrsmatrix.hh:111
Matrix::block_type block_type
The block_type of the underlying matrix.
Definition: bcrsmatrix.hh:119
ImplicitMatrixBuilder(Matrix &m)
Creates an ImplicitMatrixBuilder for matrix m.
Definition: bcrsmatrix.hh:164
M_ Matrix
The underlying matrix.
Definition: bcrsmatrix.hh:116
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:188
size_type M() const
The number of columns in the matrix.
Definition: bcrsmatrix.hh:211
Matrix::size_type size_type
The size_type of the underlying matrix.
Definition: bcrsmatrix.hh:122
row_object operator[](size_type i) const
Returns a proxy for entries in row i.
Definition: bcrsmatrix.hh:199
size_type N() const
The number of rows in the matrix.
Definition: bcrsmatrix.hh:205
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:306
A generic dynamic dense matrix.
Definition: matrix.hh:25
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:41
T block_type
Export the type representing the components.
Definition: matrix.hh:32
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:430
Default exception class for range errors.
Definition: exceptions.hh:279
size_type size() const
number of blocks in the array (are of size 1 here)
Definition: basearray.hh:752
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
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:10
Fallback implementation of the nullptr object in C++0x.
Standard Dune debug streams.
Statistics about compression achieved in implicit mode.
Definition: bcrsmatrix.hh:82
size_type overflow_total
total number of elements written to the overflow area during construction.
Definition: bcrsmatrix.hh:88
size_type maximum
maximum number of non-zeroes per row.
Definition: bcrsmatrix.hh:86
double avg
average number of non-zeroes per row.
Definition: bcrsmatrix.hh:84
double mem_ratio
fraction of wasted memory resulting from non-used overflow area.
Definition: bcrsmatrix.hh:93
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)