Dune Core Modules (2.6.0)

bcrsmatrix.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4#ifndef DUNE_ISTL_BCRSMATRIX_HH
5#define DUNE_ISTL_BCRSMATRIX_HH
6
7#include <cmath>
8#include <complex>
9#include <set>
10#include <iostream>
11#include <algorithm>
12#include <numeric>
13#include <vector>
14#include <map>
15
16#include "istlexception.hh"
17#include "bvector.hh"
18#include "matrixutils.hh"
23
28namespace Dune {
29
69 template<typename M>
70 struct MatrixDimension;
71
73
79 template<typename size_type>
81 {
83 double avg;
85 size_type maximum;
87 size_type overflow_total;
89
92 double mem_ratio;
93 };
94
96
108 template<class M_>
110 {
111
112 public:
113
115 typedef M_ Matrix;
116
119
121 typedef typename Matrix::size_type size_type;
122
124
130 {
131
132 public:
133
136 {
137 return _m.entry(_i,j);
138 }
139
140#ifndef DOXYGEN
141
143 : _m(m)
144 , _i(i)
145 {}
146
147#endif
148
149 private:
150
151 Matrix& _m;
152 size_type _i;
153
154 };
155
157
164 : _m(m)
165 {
166 if (m.buildMode() != Matrix::implicit)
167 DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix in implicit build mode");
168 if (m.buildStage() != Matrix::building)
169 DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix with set size that has not been compressed() yet");
170 }
171
173
187 ImplicitMatrixBuilder(Matrix& m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
188 : _m(m)
189 {
190 if (m.buildStage() != Matrix::notAllocated)
191 DUNE_THROW(BCRSMatrixError,"You can only set up a matrix for this ImplicitBuilder if it has no memory allocated yet");
192 m.setBuildMode(Matrix::implicit);
193 m.setImplicitBuildModeParameters(avg_cols_per_row,overflow_fraction);
194 m.setSize(rows,cols);
195 }
196
199 {
200 return row_object(_m,i);
201 }
202
204 size_type N() const
205 {
206 return _m.N();
207 }
208
210 size_type M() const
211 {
212 return _m.M();
213 }
214
215 private:
216
217 Matrix& _m;
218
219 };
220
421 template<class B, class A=std::allocator<B> >
423 {
424 friend struct MatrixDimension<BCRSMatrix>;
425 public:
439 built=3
440 };
441
442 //===== type definitions and constants
443
445 typedef typename B::field_type field_type;
446
448 typedef B block_type;
449
451 typedef A allocator_type;
452
454 typedef Imp::CompressedBlockVectorWindow<B,A> row_type;
455
457 typedef typename A::size_type size_type;
458
460 typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
461
463 enum {
465 blocklevel = B::blocklevel+1
466 };
467
502 unknown
503 };
504
505 //===== random access interface to rows of the matrix
506
509 {
510#ifdef DUNE_ISTL_WITH_CHECKING
511 if (build_mode == implicit && ready != built)
512 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
513 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
514 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
515#endif
516 return r[i];
517 }
518
521 {
522#ifdef DUNE_ISTL_WITH_CHECKING
523 if (build_mode == implicit && ready != built)
524 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
525 if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
526 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
527#endif
528 return r[i];
529 }
530
531
532 //===== iterator interface to rows of the matrix
533
535 template<class T>
537 : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
538 {
539
540 public:
542 typedef typename std::remove_const<T>::type ValueType;
543
546 friend class RealRowIterator<const ValueType>;
547 friend class RealRowIterator<ValueType>;
548
551 : p(_p), i(_i)
552 {}
553
556 : p(0), i(0)
557 {}
558
560 : p(it.p), i(it.i)
561 {}
562
563
566 {
567 return i;
568 }
569
570 std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
571 {
572 assert(other.p==p);
573 return (other.i-i);
574 }
575
576 std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
577 {
578 assert(other.p==p);
579 return (other.i-i);
580 }
581
583 bool equals (const RealRowIterator<ValueType>& other) const
584 {
585 assert(other.p==p);
586 return i==other.i;
587 }
588
590 bool equals (const RealRowIterator<const ValueType>& other) const
591 {
592 assert(other.p==p);
593 return i==other.i;
594 }
595
596 private:
598 void increment()
599 {
600 ++i;
601 }
602
604 void decrement()
605 {
606 --i;
607 }
608
609 void advance(std::ptrdiff_t diff)
610 {
611 i+=diff;
612 }
613
614 T& elementAt(std::ptrdiff_t diff) const
615 {
616 return p[i+diff];
617 }
618
620 row_type& dereference () const
621 {
622 return p[i];
623 }
624
625 row_type* p;
626 size_type i;
627 };
628
632
635 {
636 return Iterator(r,0);
637 }
638
641 {
642 return Iterator(r,n);
643 }
644
648 {
649 return Iterator(r,n-1);
650 }
651
655 {
656 return Iterator(r,-1);
657 }
658
661
663 typedef typename row_type::Iterator ColIterator;
664
668
669
672 {
673 return ConstIterator(r,0);
674 }
675
678 {
679 return ConstIterator(r,n);
680 }
681
685 {
686 return ConstIterator(r,n-1);
687 }
688
692 {
693 return ConstIterator(r,-1);
694 }
695
698
700 typedef typename row_type::ConstIterator ConstColIterator;
701
702 //===== constructors & resizers
703
704 // we use a negative overflowsize to indicate that the implicit
705 // mode parameters have not been set yet
706
709 : build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
710 allocationSize_(0), r(0), a(0),
711 avg(0), overflowsize(-1.0)
712 {}
713
716 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
717 allocationSize_(0), r(0), a(0),
718 avg(0), overflowsize(-1.0)
719 {
720 allocate(_n, _m, _nnz,true,false);
721 }
722
725 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
726 allocationSize_(0), r(0), a(0),
727 avg(0), overflowsize(-1.0)
728 {
729 allocate(_n, _m,0,true,false);
730 }
731
733
743 BCRSMatrix (size_type _n, size_type _m, size_type _avg, double _overflowsize, BuildMode bm)
744 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
745 allocationSize_(0), r(0), a(0),
746 avg(_avg), overflowsize(_overflowsize)
747 {
748 if (bm != implicit)
749 DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
750 // Prevent user from setting a negative overflowsize:
751 // 1) It doesn't make sense
752 // 2) We use a negative overflow value to indicate that the parameters
753 // have not been set yet
754 if (_overflowsize < 0.0)
755 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
756 implicit_allocate(_n,_m);
757 }
758
765 : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
766 allocationSize_(0), r(0), a(0),
767 avg(Mat.avg), overflowsize(Mat.overflowsize)
768 {
769 if (!(Mat.ready == notAllocated || Mat.ready == built))
770 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
771
772 // deep copy in global array
773 size_type _nnz = Mat.nonzeroes();
774
775 j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
776 allocate(Mat.n, Mat.m, _nnz, true, true);
777
778 // build window structure
780 }
781
784 {
785 deallocate();
786 }
787
793 {
794 if (ready == notAllocated)
795 {
796 build_mode = bm;
797 return;
798 }
799 if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
800 build_mode = bm;
801 else
802 DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
803 }
804
820 void setSize(size_type rows, size_type columns, size_type nnz=0)
821 {
822 // deallocate already setup memory
823 deallocate();
824
825 if (build_mode == implicit)
826 {
827 if (nnz>0)
828 DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");
829
830 // implicit allocates differently
831 implicit_allocate(rows,columns);
832 }
833 else
834 {
835 // allocate matrix memory
836 allocate(rows, columns, nnz, true, false);
837 }
838 }
839
848 void setImplicitBuildModeParameters(size_type _avg, double _overflow)
849 {
850 // Prevent user from setting a negative overflowsize:
851 // 1) It doesn't make sense
852 // 2) We use a negative overflow value to indicate that the parameters
853 // have not been set yet
854 if (_overflow < 0.0)
855 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
856
857 // make sure the parameters aren't changed after memory has been allocated
858 if (ready != notAllocated)
859 DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
860 avg = _avg;
861 overflowsize = _overflow;
862 }
863
871 {
872 // return immediately when self-assignment
873 if (&Mat==this) return *this;
874
875 if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
876 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");
877
878 // make it simple: ALWAYS throw away memory for a and j_
879 // and deallocate rows only if n != Mat.n
880 deallocate(n!=Mat.n);
881
882 // reallocate the rows if required
883 if (n>0 && n!=Mat.n) {
884 // free rows
885 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
886 rowAllocator_.destroy(riter);
887 rowAllocator_.deallocate(r,n);
888 }
889
890 nnz_ = Mat.nonzeroes();
891
892 // allocate a, share j_
893 j_ = Mat.j_;
894 allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);
895
896 // build window structure
898 return *this;
899 }
900
903 {
904
905 if (!(ready == notAllocated || ready == built))
906 DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");
907
908 for (size_type i=0; i<n; i++) r[i] = k;
909 return *this;
910 }
911
912 //===== row-wise creation interface
913
916 {
917 public:
920 : Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
921 {
922 if (Mat.build_mode == unknown && Mat.ready == building)
923 {
924 Mat.build_mode = row_wise;
925 }
926 if (i==0 && Mat.ready != building)
927 DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
928 if(Mat.build_mode!=row_wise)
929 DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
930 if(i==0 && _Mat.N()==0)
931 // empty Matrix is always built.
932 Mat.ready = built;
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,nullptr,nullptr);
978
979 // initialize the j array for row i from pattern
980 std::copy(pattern.cbegin(), pattern.cend(), Mat.r[i].getindexptr());
981
982 // now go to next row
983 i++;
984 pattern.clear();
985
986 // check if this was last row
987 if (i==Mat.n)
988 {
989 Mat.ready = built;
990 if(Mat.nnz_ > 0)
991 {
992 // Set nnz to the exact number of nonzero blocks inserted
993 // as some methods rely on it
994 Mat.nnz_ = nnz;
995 // allocate data array
996 Mat.allocateData();
997 Mat.setDataPointers();
998 }
999 }
1000 // done
1001 return *this;
1002 }
1003
1005 bool operator!= (const CreateIterator& it) const
1006 {
1007 return (i!=it.i) || (&Mat!=&it.Mat);
1008 }
1009
1011 bool operator== (const CreateIterator& it) const
1012 {
1013 return (i==it.i) && (&Mat==&it.Mat);
1014 }
1015
1018 {
1019 return i;
1020 }
1021
1024 {
1025 pattern.insert(j);
1026 }
1027
1030 {
1031 return pattern.find(j) != pattern.end();
1032 }
1039 {
1040 return pattern.size();
1041 }
1042
1043 private:
1044 BCRSMatrix& Mat; // the matrix we are defining
1045 size_type i; // current row to be defined
1046 size_type nnz; // count total number of nonzeros
1047 typedef std::set<size_type,std::less<size_type> > PatternType;
1048 PatternType pattern; // used to compile entries in a row
1049 row_type current_row; // row pointing to the current row to setup
1050 };
1051
1053 friend class CreateIterator;
1054
1057 {
1058 return CreateIterator(*this,0);
1059 }
1060
1063 {
1064 return CreateIterator(*this,n);
1065 }
1066
1067
1068 //===== random creation interface
1069
1077 {
1078 if (build_mode!=random)
1079 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1080 if (ready != building)
1081 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1082
1083 r[i].setsize(s);
1084 }
1085
1088 {
1089#ifdef DUNE_ISTL_WITH_CHECKING
1090 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
1091 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
1092#endif
1093 return r[i].getsize();
1094 }
1095
1098 {
1099 if (build_mode!=random)
1100 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1101 if (ready != building)
1102 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1103
1104 r[i].setsize(r[i].getsize()+s);
1105 }
1106
1109 {
1110 if (build_mode!=random)
1111 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1112 if (ready != building)
1113 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1114
1115 // compute total size, check positivity
1116 size_type total=0;
1117 for (size_type i=0; i<n; i++)
1118 {
1119 total += r[i].getsize();
1120 }
1121
1122 if(nnz_ == 0)
1123 // allocate/check memory
1124 allocate(n,m,total,false,false);
1125 else if(nnz_ < total)
1126 DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
1127 <<"sufficient for calculated nonzeros ("<<total<<"! ");
1128
1129 // set the window pointers correctly
1131
1132 // initialize j_ array with m (an invalid column index)
1133 // this indicates an unused entry
1134 for (size_type k=0; k<nnz_; k++)
1135 j_.get()[k] = m;
1136 ready = rowSizesBuilt;
1137 }
1138
1140
1151 {
1152 if (build_mode!=random)
1153 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1154 if (ready==built)
1155 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1156 if (ready==building)
1157 DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
1158 if (ready==notAllocated)
1159 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1160
1161 if (col >= m)
1162 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1163
1164 // get row range
1165 size_type* const first = r[row].getindexptr();
1166 size_type* const last = first + r[row].getsize();
1167
1168 // find correct insertion position for new column index
1169 size_type* pos = std::lower_bound(first,last,col);
1170
1171 // check if index is already in row
1172 if (pos!=last && *pos == col) return;
1173
1174 // find end of already inserted column indices
1175 size_type* end = std::lower_bound(pos,last,m);
1176 if (end==last)
1177 DUNE_THROW(BCRSMatrixError,"row is too small");
1178
1179 // insert new column index at correct position
1180 std::copy_backward(pos,end,end+1);
1181 *pos = col;
1182 }
1183
1185
1192 template<typename It>
1193 void setIndices(size_type row, It begin, It end)
1194 {
1195 size_type row_size = r[row].size();
1196 size_type* col_begin = r[row].getindexptr();
1197 size_type* col_end;
1198 // consistency check between allocated row size and number of passed column indices
1199 if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
1200 DUNE_THROW(BCRSMatrixError,"Given size of row " << row
1201 << " (" << row_size
1202 << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
1203 std::sort(col_begin,col_end);
1204 }
1205
1208 {
1209 if (build_mode!=random)
1210 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1211 if (ready==built)
1212 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1213 if (ready==building)
1214 DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
1215 if (ready==notAllocated)
1216 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1217
1218 // check if there are undefined indices
1219 RowIterator endi=end();
1220 for (RowIterator i=begin(); i!=endi; ++i)
1221 {
1222 ColIterator endj = (*i).end();
1223 for (ColIterator j=(*i).begin(); j!=endj; ++j) {
1224 if (j.index() >= m) {
1225 dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
1226 <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
1227 nnz_ -= ((*i).end().offset() - j.offset());
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<class X, class Y, class F>
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 template <typename ft = field_type,
1807 typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
1808 typename FieldTraits<ft>::real_type infinity_norm() const {
1809 if (ready != built)
1810 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1811
1812 using real_type = typename FieldTraits<ft>::real_type;
1813 using std::max;
1814
1815 real_type norm = 0;
1816 for (auto const &x : *this) {
1817 real_type sum = 0;
1818 for (auto const &y : x)
1819 sum += y.infinity_norm();
1820 norm = max(sum, norm);
1821 }
1822 return norm;
1823 }
1824
1826 template <typename ft = field_type,
1827 typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
1828 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1829 if (ready != built)
1830 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1831
1832 using real_type = typename FieldTraits<ft>::real_type;
1833 using std::max;
1834
1835 real_type norm = 0;
1836 for (auto const &x : *this) {
1837 real_type sum = 0;
1838 for (auto const &y : x)
1839 sum += y.infinity_norm_real();
1840 norm = max(sum, norm);
1841 }
1842 return norm;
1843 }
1844
1846 template <typename ft = field_type,
1847 typename std::enable_if<has_nan<ft>::value, int>::type = 0>
1848 typename FieldTraits<ft>::real_type infinity_norm() const {
1849 if (ready != built)
1850 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1851
1852 using real_type = typename FieldTraits<ft>::real_type;
1853 using std::max;
1854
1855 real_type norm = 0;
1856 real_type isNaN = 1;
1857 for (auto const &x : *this) {
1858 real_type sum = 0;
1859 for (auto const &y : x)
1860 sum += y.infinity_norm();
1861 norm = max(sum, norm);
1862 isNaN += sum;
1863 }
1864 isNaN /= isNaN;
1865 return norm * isNaN;
1866 }
1867
1869 template <typename ft = field_type,
1870 typename std::enable_if<has_nan<ft>::value, int>::type = 0>
1871 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1872 if (ready != built)
1873 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1874
1875 using real_type = typename FieldTraits<ft>::real_type;
1876 using std::max;
1877
1878 real_type norm = 0;
1879 real_type isNaN = 1;
1880 for (auto const &x : *this) {
1881 real_type sum = 0;
1882 for (auto const &y : x)
1883 sum += y.infinity_norm_real();
1884 norm = max(sum, norm);
1885 isNaN += sum;
1886 }
1887 isNaN /= isNaN;
1888 return norm * isNaN;
1889 }
1890
1891 //===== sizes
1892
1894 size_type N () const
1895 {
1896 return n;
1897 }
1898
1900 size_type M () const
1901 {
1902 return m;
1903 }
1904
1907 {
1908 // in case of row-wise allocation
1909 if( nnz_ <= 0 )
1910 nnz_ = std::accumulate( begin(), end(), size_type( 0 ), [] ( size_type s, const row_type &row ) { return s+row.getsize(); } );
1911 return nnz_;
1912 }
1913
1916 {
1917 return ready;
1918 }
1919
1922 {
1923 return build_mode;
1924 }
1925
1926 //===== query
1927
1929 bool exists (size_type i, size_type j) const
1930 {
1931#ifdef DUNE_ISTL_WITH_CHECKING
1932 if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
1933 if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
1934#endif
1935 return (r[i].size() && r[i].find(j) != r[i].end());
1936 }
1937
1938
1939 protected:
1940 // state information
1941 BuildMode build_mode; // row wise or whole matrix
1942 BuildStage ready; // indicate the stage the matrix building is in
1943
1944 // The allocator used for memory management
1945 typename A::template rebind<B>::other allocator_;
1946
1947 typename A::template rebind<row_type>::other rowAllocator_;
1948
1949 typename A::template rebind<size_type>::other sizeAllocator_;
1950
1951 // size of the matrix
1952 size_type n; // number of rows
1953 size_type m; // number of columns
1954 mutable size_type nnz_; // number of nonzeroes contained in the matrix
1955 size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
1956 // zero means that memory is allocated separately for each row.
1957
1958 // the rows are dynamically allocated
1959 row_type* r; // [n] the individual rows having pointers into a,j arrays
1960
1961 // dynamically allocated memory
1962 B* a; // [allocationSize] non-zero entries of the matrix in row-wise ordering
1963 // If a single array of column indices is used, it can be shared
1964 // between different matrices with the same sparsity pattern
1965 std::shared_ptr<size_type> j_; // [allocationSize] column indices of entries
1966
1967 // additional data is needed in implicit buildmode
1968 size_type avg;
1969 double overflowsize;
1970
1971 typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
1972 OverflowType overflow;
1973
1974 void setWindowPointers(ConstRowIterator row)
1975 {
1976 row_type current_row(a,j_.get(),0); // Pointers to current row data
1977 for (size_type i=0; i<n; i++, ++row) {
1978 // set row i
1979 size_type s = row->getsize();
1980
1981 if (s>0) {
1982 // setup pointers and size
1983 r[i].set(s,current_row.getptr(), current_row.getindexptr());
1984 // update pointer for next row
1985 current_row.setptr(current_row.getptr()+s);
1986 current_row.setindexptr(current_row.getindexptr()+s);
1987 } else{
1988 // empty row
1989 r[i].set(0,nullptr,nullptr);
1990 }
1991 }
1992 }
1993
1995
2000 {
2001 size_type* jptr = j_.get();
2002 for (size_type i=0; i<n; ++i, ++row) {
2003 // set row i
2004 size_type s = row->getsize();
2005
2006 if (s>0) {
2007 // setup pointers and size
2008 r[i].setsize(s);
2009 r[i].setindexptr(jptr);
2010 } else{
2011 // empty row
2012 r[i].set(0,nullptr,nullptr);
2013 }
2014
2015 // advance position in global array
2016 jptr += s;
2017 }
2018 }
2019
2021
2026 {
2027 B* aptr = a;
2028 for (size_type i=0; i<n; ++i) {
2029 // set row i
2030 if (r[i].getsize() > 0) {
2031 // setup pointers and size
2032 r[i].setptr(aptr);
2033 } else{
2034 // empty row
2035 r[i].set(0,nullptr,nullptr);
2036 }
2037
2038 // advance position in global array
2039 aptr += r[i].getsize();
2040 }
2041 }
2042
2045 {
2046 setWindowPointers(Mat.begin());
2047
2048 // copy data
2049 for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
2050
2051 // finish off
2052 build_mode = row_wise; // dummy
2053 ready = built;
2054 }
2055
2061 void deallocate(bool deallocateRows=true)
2062 {
2063
2064 if (notAllocated)
2065 return;
2066
2067 if (allocationSize_>0)
2068 {
2069 // a,j_ have been allocated as one long vector
2070 j_.reset();
2071 if (a)
2072 {
2073 for(B *aiter=a+(allocationSize_-1), *aend=a-1; aiter!=aend; --aiter)
2074 allocator_.destroy(aiter);
2075 allocator_.deallocate(a,allocationSize_);
2076 a = nullptr;
2077 }
2078 }
2079 else if (r)
2080 {
2081 // check if memory for rows have been allocated individually
2082 for (size_type i=0; i<n; i++)
2083 if (r[i].getsize()>0)
2084 {
2085 for (B *col=r[i].getptr()+(r[i].getsize()-1),
2086 *colend = r[i].getptr()-1; col!=colend; --col) {
2087 allocator_.destroy(col);
2088 }
2089 sizeAllocator_.deallocate(r[i].getindexptr(),1);
2090 allocator_.deallocate(r[i].getptr(),1);
2091 // clear out row data in case we don't want to deallocate the rows
2092 // otherwise we might run into a double free problem here later
2093 r[i].set(0,nullptr,nullptr);
2094 }
2095 }
2096
2097 // deallocate the rows
2098 if (n>0 && deallocateRows && r) {
2099 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
2100 rowAllocator_.destroy(riter);
2101 rowAllocator_.deallocate(r,n);
2102 r = nullptr;
2103 }
2104
2105 // Mark matrix as not built at all.
2106 ready=notAllocated;
2107
2108 }
2109
2112 {
2113 typename A::template rebind<size_type>::other& sizeAllocator_;
2114
2115 public:
2116 Deallocator(typename A::template rebind<size_type>::other& sizeAllocator)
2117 : sizeAllocator_(sizeAllocator)
2118 {}
2119
2120 void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); }
2121 };
2122
2123
2141 void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
2142 {
2143 // Store size
2144 n = rows;
2145 m = columns;
2146 nnz_ = allocationSize;
2147 allocationSize_ = allocationSize;
2148
2149 // allocate rows
2150 if(allocateRows) {
2151 if (n>0) {
2152 if (r)
2153 DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
2154 r = rowAllocator_.allocate(rows);
2155 // initialize row entries
2156 for(row_type* ri=r; ri!=r+rows; ++ri)
2157 rowAllocator_.construct(ri, row_type());
2158 }else{
2159 r = 0;
2160 }
2161 }
2162
2163 // allocate a and j_ array
2164 if (allocate_data)
2165 allocateData();
2166 if (allocationSize_>0) {
2167 // allocate column indices only if not yet present (enable sharing)
2168 if (!j_.get())
2169 j_.reset(sizeAllocator_.allocate(allocationSize_),Deallocator(sizeAllocator_));
2170 }else{
2171 j_.reset();
2172 }
2173
2174 // Mark the matrix as not built.
2175 ready = building;
2176 }
2177
2178 void allocateData()
2179 {
2180 if (a)
2181 DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
2182 if (allocationSize_>0) {
2183 a = allocator_.allocate(allocationSize_);
2184 // use placement new to call constructor that allocates
2185 // additional memory.
2186 new (a) B[allocationSize_];
2187 } else {
2188 a = nullptr;
2189 }
2190 }
2191
2198 {
2199 if (build_mode != implicit)
2200 DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
2201 if (ready != notAllocated)
2202 DUNE_THROW(InvalidStateException,"memory has already been allocated");
2203
2204 // check to make sure the user has actually set the parameters
2205 if (overflowsize < 0)
2206 DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
2207 //calculate size of overflow region, add buffer for row sort!
2208 size_type osize = (size_type) (_n*avg)*overflowsize + 4*avg;
2209 allocationSize_ = _n*avg + osize;
2210
2211 allocate(_n, _m, allocationSize_,true,true);
2212
2213 //set row pointers correctly
2214 size_type* jptr = j_.get() + osize;
2215 B* aptr = a + osize;
2216 for (size_type i=0; i<n; i++)
2217 {
2218 r[i].set(0,aptr,jptr);
2219 jptr = jptr + avg;
2220 aptr = aptr + avg;
2221 }
2222
2223 ready = building;
2224 }
2225 };
2226
2227
2230} // end namespace
2231
2232#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:916
bool operator==(const CreateIterator &it) const
equality
Definition: bcrsmatrix.hh:1011
CreateIterator & operator++()
prefix increment
Definition: bcrsmatrix.hh:936
size_type index() const
The number of the row that the iterator currently points to.
Definition: bcrsmatrix.hh:1017
bool operator!=(const CreateIterator &it) const
inequality
Definition: bcrsmatrix.hh:1005
CreateIterator(BCRSMatrix &_Mat, size_type _i)
constructor
Definition: bcrsmatrix.hh:919
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1023
bool contains(size_type j)
return true if column index is in row
Definition: bcrsmatrix.hh:1029
size_type size() const
Get the current row size.
Definition: bcrsmatrix.hh:1038
Class used by shared_ptr to deallocate memory using the proper allocator.
Definition: bcrsmatrix.hh:2112
Iterator access to matrix rows
Definition: bcrsmatrix.hh:538
RealRowIterator()
empty constructor, use with care!
Definition: bcrsmatrix.hh:555
bool equals(const RealRowIterator< ValueType > &other) const
equality
Definition: bcrsmatrix.hh:583
std::remove_const< T >::type ValueType
The unqualified value type.
Definition: bcrsmatrix.hh:542
bool equals(const RealRowIterator< const ValueType > &other) const
equality
Definition: bcrsmatrix.hh:590
RealRowIterator(row_type *_p, size_type _i)
constructor
Definition: bcrsmatrix.hh:550
size_type index() const
return index
Definition: bcrsmatrix.hh:565
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:423
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: bcrsmatrix.hh:1929
BuildStage buildStage() const
The current build stage of the matrix.
Definition: bcrsmatrix.hh:1915
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:634
friend class CreateIterator
allow CreateIterator to access internal data
Definition: bcrsmatrix.hh:1053
void copyWindowStructure(const BCRSMatrix &Mat)
Copy the window structure from another matrix.
Definition: bcrsmatrix.hh:2044
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:666
BCRSMatrix(size_type _n, size_type _m, size_type _avg, double _overflowsize, BuildMode bm)
construct matrix with a known average number of entries per row
Definition: bcrsmatrix.hh:743
void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
Allocate memory for the matrix structure.
Definition: bcrsmatrix.hh:2141
B::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:445
BCRSMatrix & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition: bcrsmatrix.hh:1551
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: bcrsmatrix.hh:1828
~BCRSMatrix()
destructor
Definition: bcrsmatrix.hh:783
void deallocate(bool deallocateRows=true)
deallocate memory of the matrix.
Definition: bcrsmatrix.hh:2061
Iterator RowIterator
rename the iterators for easier access
Definition: bcrsmatrix.hh:660
row_type & operator[](size_type i)
random access to the rows
Definition: bcrsmatrix.hh:508
BCRSMatrix()
an empty matrix
Definition: bcrsmatrix.hh:708
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1108
void incrementrowsize(size_type i, size_type s=1)
increment size of row i by s (1 by default)
Definition: bcrsmatrix.hh:1097
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:1906
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition: bcrsmatrix.hh:697
void setrowsize(size_type i, size_type s)
Set number of indices in row i to s.
Definition: bcrsmatrix.hh:1076
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:630
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:691
Iterator beforeBegin()
Definition: bcrsmatrix.hh:654
BuildMode
we support two modes
Definition: bcrsmatrix.hh:469
@ implicit
Build entries randomly with an educated guess for the number of entries per row.
Definition: bcrsmatrix.hh:498
@ unknown
Build mode not set!
Definition: bcrsmatrix.hh:502
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:489
@ row_wise
Build in a row-wise manner.
Definition: bcrsmatrix.hh:480
BCRSMatrix(size_type _n, size_type _m, size_type _nnz, BuildMode bm)
matrix with known number of nonzeroes
Definition: bcrsmatrix.hh:715
::Dune::CompressionStatistics< size_type > CompressionStatistics
The type for the statistics object returned by compress()
Definition: bcrsmatrix.hh:460
BCRSMatrix & operator-=(const BCRSMatrix &b)
Subtract the entries of another matrix from this one.
Definition: bcrsmatrix.hh:1526
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition: bcrsmatrix.hh:764
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:640
void setIndices(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition: bcrsmatrix.hh:1193
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1150
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:663
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: bcrsmatrix.hh:1800
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:457
BCRSMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: bcrsmatrix.hh:1471
BCRSMatrix & operator+=(const BCRSMatrix &b)
Add the entries of another matrix to this one.
Definition: bcrsmatrix.hh:1504
@ blocklevel
The number of blocklevels the matrix contains.
Definition: bcrsmatrix.hh:465
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1062
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:647
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:700
void usmv(F &&alpha, const X &x, Y &y) const
y += alpha A x
Definition: bcrsmatrix.hh:1631
size_type getrowsize(size_type i) const
get current number of indices in row i
Definition: bcrsmatrix.hh:1087
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:1900
Imp::CompressedBlockVectorWindow< B, A > row_type
implement row_type with compressed vector
Definition: bcrsmatrix.hh:454
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1056
BuildStage
Definition: bcrsmatrix.hh:426
@ rowSizesBuilt
The row sizes of the matrix are known.
Definition: bcrsmatrix.hh:437
@ built
The matrix structure is fully built.
Definition: bcrsmatrix.hh:439
@ notbuilt
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:428
@ notAllocated
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:430
@ building
Matrix is currently being built, some memory has been allocated, build mode and size are fixed.
Definition: bcrsmatrix.hh:432
BuildMode buildMode() const
The currently selected build mode of the matrix.
Definition: bcrsmatrix.hh:1921
void mmv(const X &x, Y &y) const
y -= A x
Definition: bcrsmatrix.hh:1612
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: bcrsmatrix.hh:1808
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:448
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:2197
BCRSMatrix(size_type _n, size_type _m, BuildMode bm)
matrix with unknown number of nonzeroes
Definition: bcrsmatrix.hh:724
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1207
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition: bcrsmatrix.hh:1319
void setDataPointers()
Set data pointers for all rows.
Definition: bcrsmatrix.hh:2025
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1894
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:792
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:820
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:870
void setColumnPointers(ConstRowIterator row)
Copy row sizes from iterator range starting at row and set column index pointers for all rows.
Definition: bcrsmatrix.hh:1999
ConstIterator end() const
Get const iterator to one beyond last row.
Definition: bcrsmatrix.hh:677
ConstIterator begin() const
Get const iterator to first row.
Definition: bcrsmatrix.hh:671
void setImplicitBuildModeParameters(size_type _avg, double _overflow)
Set parameters needed for creation in implicit build mode.
Definition: bcrsmatrix.hh:848
A allocator_type
export the allocator type
Definition: bcrsmatrix.hh:451
ConstIterator beforeEnd() const
Definition: bcrsmatrix.hh:684
Proxy row object for entry access.
Definition: bcrsmatrix.hh:130
block_type & operator[](size_type j) const
Returns entry in column j.
Definition: bcrsmatrix.hh:135
A wrapper for uniform access to the BCRSMatrix during and after the build stage in implicit build mod...
Definition: bcrsmatrix.hh:110
Matrix::block_type block_type
The block_type of the underlying matrix.
Definition: bcrsmatrix.hh:118
ImplicitMatrixBuilder(Matrix &m)
Creates an ImplicitMatrixBuilder for matrix m.
Definition: bcrsmatrix.hh:163
M_ Matrix
The underlying matrix.
Definition: bcrsmatrix.hh:115
ImplicitMatrixBuilder(Matrix &m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
Sets up matrix m for implicit construction using the given parameters and creates an ImplicitBmatrixu...
Definition: bcrsmatrix.hh:187
size_type M() const
The number of columns in the matrix.
Definition: bcrsmatrix.hh:210
Matrix::size_type size_type
The size_type of the underlying matrix.
Definition: bcrsmatrix.hh:121
row_object operator[](size_type i) const
Returns a proxy for entries in row i.
Definition: bcrsmatrix.hh:198
size_type N() const
The number of rows in the matrix.
Definition: bcrsmatrix.hh:204
The overflow error used during implicit BCRSMatrix construction was exhausted.
Definition: istlexception.hh:34
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:279
A generic dynamic dense matrix.
Definition: matrix.hh:555
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:571
T block_type
Export the type representing the components.
Definition: matrix.hh:562
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:433
Default exception class for range errors.
Definition: exceptions.hh:252
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:331
constexpr decltype(auto) elementAt(Container &&c, Index &&i)
Get element at given position from container.
Definition: hybridutilities.hh:133
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:10
Standard Dune debug streams.
Statistics about compression achieved in implicit mode.
Definition: bcrsmatrix.hh:81
size_type overflow_total
total number of elements written to the overflow area during construction.
Definition: bcrsmatrix.hh:87
size_type maximum
maximum number of non-zeroes per row.
Definition: bcrsmatrix.hh:85
double avg
average number of non-zeroes per row.
Definition: bcrsmatrix.hh:83
double mem_ratio
fraction of wasted memory resulting from non-used overflow area.
Definition: bcrsmatrix.hh:92
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)