dune-istl 3.0-git
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"
19#include <dune/common/stdstreams.hh>
20#include <dune/common/iteratorfacades.hh>
21#include <dune/common/typetraits.hh>
22#include <dune/common/ftraits.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;
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
410 template<class B, class A=std::allocator<B> >
412 {
413 friend struct MatrixDimension<BCRSMatrix>;
414 public:
430
431 //===== type definitions and constants
432
434 typedef typename B::field_type field_type;
435
437 typedef B block_type;
438
440 typedef A allocator_type;
441
444
446 typedef typename A::size_type size_type;
447
449 typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
450
452 enum {
454 blocklevel = B::blocklevel+1
455 };
456
493
494 //===== random access interface to rows of the matrix
495
498 {
499#ifdef DUNE_ISTL_WITH_CHECKING
500 if (build_mode == implicit && ready != built)
501 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
502 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
503 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
504#endif
505 return r[i];
506 }
507
510 {
511#ifdef DUNE_ISTL_WITH_CHECKING
512 if (build_mode == implicit && ready != built)
513 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
514 if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
515 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
516#endif
517 return r[i];
518 }
519
520
521 //===== iterator interface to rows of the matrix
522
524 template<class T>
526 : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
527 {
528
529 public:
531 typedef typename std::remove_const<T>::type ValueType;
532
533 friend class RandomAccessIteratorFacade<RealRowIterator<const ValueType>, const ValueType>;
534 friend class RandomAccessIteratorFacade<RealRowIterator<ValueType>, ValueType>;
535 friend class RealRowIterator<const ValueType>;
536 friend class RealRowIterator<ValueType>;
537
540 : p(_p), i(_i)
541 {}
542
545 : p(0), i(0)
546 {}
547
549 : p(it.p), i(it.i)
550 {}
551
552
555 {
556 return i;
557 }
558
559 std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
560 {
561 assert(other.p==p);
562 return (other.i-i);
563 }
564
566 {
567 assert(other.p==p);
568 return (other.i-i);
569 }
570
573 {
574 assert(other.p==p);
575 return i==other.i;
576 }
577
580 {
581 assert(other.p==p);
582 return i==other.i;
583 }
584
585 private:
587 void increment()
588 {
589 ++i;
590 }
591
593 void decrement()
594 {
595 --i;
596 }
597
598 void advance(std::ptrdiff_t diff)
599 {
600 i+=diff;
601 }
602
603 T& elementAt(std::ptrdiff_t diff) const
604 {
605 return p[i+diff];
606 }
607
609 row_type& dereference () const
610 {
611 return p[i];
612 }
613
614 row_type* p;
615 size_type i;
616 };
617
621
624 {
625 return Iterator(r,0);
626 }
627
630 {
631 return Iterator(r,n);
632 }
633
637 {
638 return Iterator(r,n-1);
639 }
640
644 {
645 return Iterator(r,-1);
646 }
647
650
653
657
658
661 {
662 return ConstIterator(r,0);
663 }
664
667 {
668 return ConstIterator(r,n);
669 }
670
674 {
675 return ConstIterator(r,n-1);
676 }
677
681 {
682 return ConstIterator(r,-1);
683 }
684
687
690
691 //===== constructors & resizers
692
693 // we use a negative overflowsize to indicate that the implicit
694 // mode parameters have not been set yet
695
698 : build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
699 allocationSize_(0), r(0), a(0),
700 avg(0), overflowsize(-1.0)
701 {}
702
705 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
706 allocationSize_(0), r(0), a(0),
707 avg(0), overflowsize(-1.0)
708 {
709 allocate(_n, _m, _nnz,true,false);
710 }
711
714 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
715 allocationSize_(0), r(0), a(0),
716 avg(0), overflowsize(-1.0)
717 {
718 allocate(_n, _m,0,true,false);
719 }
720
722
733 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
734 allocationSize_(0), r(0), a(0),
736 {
737 if (bm != implicit)
738 DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
739 // Prevent user from setting a negative overflowsize:
740 // 1) It doesn't make sense
741 // 2) We use a negative overflow value to indicate that the parameters
742 // have not been set yet
743 if (_overflowsize < 0.0)
744 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
745 implicit_allocate(_n,_m);
746 }
747
754 : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
755 allocationSize_(0), r(0), a(0),
757 {
758 if (!(Mat.ready == notAllocated || Mat.ready == built))
759 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
760
761 // deep copy in global array
762 size_type _nnz = Mat.nnz_;
763
764 // in case of row-wise allocation
765 if (_nnz<=0)
766 {
767 _nnz = 0;
768 for (size_type i=0; i<Mat.n; i++)
769 _nnz += Mat.r[i].getsize();
770 }
771
772 j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
773 allocate(Mat.n, Mat.m, _nnz, true, true);
774
775 // build window structure
777 }
778
781 {
782 deallocate();
783 }
784
790 {
791 if (ready == notAllocated)
792 {
793 build_mode = bm;
794 return;
795 }
796 if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
797 build_mode = bm;
798 else
799 DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
800 }
801
818 {
819 // deallocate already setup memory
820 deallocate();
821
822 if (build_mode == implicit)
823 {
824 if (nnz>0)
825 DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");
826
827 // implicit allocates differently
829 }
830 else
831 {
832 // allocate matrix memory
833 allocate(rows, columns, nnz, true, false);
834 }
835 }
836
846 {
847 // Prevent user from setting a negative overflowsize:
848 // 1) It doesn't make sense
849 // 2) We use a negative overflow value to indicate that the parameters
850 // have not been set yet
851 if (_overflow < 0.0)
852 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
853
854 // make sure the parameters aren't changed after memory has been allocated
855 if (ready != notAllocated)
856 DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
857 avg = _avg;
859 }
860
868 {
869 // return immediately when self-assignment
870 if (&Mat==this) return *this;
871
872 if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
873 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");
874
875 // make it simple: ALWAYS throw away memory for a and j_
876 // and deallocate rows only if n != Mat.n
877 deallocate(n!=Mat.n);
878
879 // reallocate the rows if required
880 if (n>0 && n!=Mat.n) {
881 // free rows
882 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
883 rowAllocator_.destroy(riter);
884 rowAllocator_.deallocate(r,n);
885 }
886
887 nnz_ = Mat.nnz_;
888 if (nnz_ <= 0)
889 {
890 for (size_type i=0; i<Mat.n; i++)
891 nnz_ += Mat.r[i].getsize();
892 }
893
894 // allocate a, share j_
895 j_ = Mat.j_;
896 allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);
897
898 // build window structure
900 return *this;
901 }
902
905 {
906
907 if (!(ready == notAllocated || ready == built))
908 DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");
909
910 for (size_type i=0; i<n; i++) r[i] = k;
911 return *this;
912 }
913
914 //===== row-wise creation interface
915
918 {
919 public:
922 : Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
923 {
924 if (Mat.build_mode == unknown && Mat.ready == building)
925 {
926 Mat.build_mode = row_wise;
927 }
928 if (i==0 && Mat.ready != building)
929 DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
930 if(Mat.build_mode!=row_wise)
931 DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
932 }
933
936 {
937 // this should only be called if matrix is in creation
938 if (Mat.ready != building)
939 DUNE_THROW(BCRSMatrixError,"matrix already built up");
940
941 // row i is defined through the pattern
942 // get memory for the row and initialize the j_ array
943 // this depends on the allocation mode
944
945 // compute size of the row
946 size_type s = pattern.size();
947
948 if(s>0) {
949 // update number of nonzeroes including this row
950 nnz += s;
951
952 // alloc memory / set window
953 if (Mat.nnz_ > 0)
954 {
955 // memory is allocated in one long array
956
957 // check if that memory is sufficient
958 if (nnz > Mat.nnz_)
959 DUNE_THROW(BCRSMatrixError,"allocated nnz too small");
960
961 // set row i
962 Mat.r[i].set(s,nullptr,current_row.getindexptr());
963 current_row.setindexptr(current_row.getindexptr()+s);
964 }else{
965 // memory is allocated individually per row
966 // allocate and set row i
967 B* b = Mat.allocator_.allocate(s);
968 // use placement new to call constructor that allocates
969 // additional memory.
970 new (b) B[s];
971 size_type* j = Mat.sizeAllocator_.allocate(s);
972 Mat.r[i].set(s,b,j);
973 }
974 }else
975 // setup empty row
976 Mat.r[i].set(0,nullptr,nullptr);
977
978 // initialize the j array for row i from pattern
979 std::copy(pattern.cbegin(), pattern.cend(), Mat.r[i].getindexptr());
980
981 // now go to next row
982 i++;
983 pattern.clear();
984
985 // check if this was last row
986 if (i==Mat.n)
987 {
988 Mat.ready = built;
989 if(Mat.nnz_ > 0)
990 {
991 // Set nnz to the exact number of nonzero blocks inserted
992 // as some methods rely on it
993 Mat.nnz_ = nnz;
994 // allocate data array
995 Mat.allocateData();
996 Mat.setDataPointers();
997 }
998 }
999 // done
1000 return *this;
1001 }
1002
1004 bool operator!= (const CreateIterator& it) const
1005 {
1006 return (i!=it.i) || (&Mat!=&it.Mat);
1007 }
1008
1010 bool operator== (const CreateIterator& it) const
1011 {
1012 return (i==it.i) && (&Mat==&it.Mat);
1013 }
1014
1017 {
1018 return i;
1019 }
1020
1023 {
1024 pattern.insert(j);
1025 }
1026
1029 {
1030 return pattern.find(j) != pattern.end();
1031 }
1038 {
1039 return pattern.size();
1040 }
1041
1042 private:
1043 BCRSMatrix& Mat; // the matrix we are defining
1044 size_type i; // current row to be defined
1045 size_type nnz; // count total number of nonzeros
1046 typedef std::set<size_type,std::less<size_type> > PatternType;
1047 PatternType pattern; // used to compile entries in a row
1048 row_type current_row; // row pointing to the current row to setup
1049 };
1050
1052 friend class CreateIterator;
1053
1056 {
1057 return CreateIterator(*this,0);
1058 }
1059
1062 {
1063 return CreateIterator(*this,n);
1064 }
1065
1066
1067 //===== random creation interface
1068
1076 {
1077 if (build_mode!=random)
1078 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1079 if (ready != building)
1080 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1081
1082 r[i].setsize(s);
1083 }
1084
1087 {
1088#ifdef DUNE_ISTL_WITH_CHECKING
1089 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
1090 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
1091#endif
1092 return r[i].getsize();
1093 }
1094
1097 {
1098 if (build_mode!=random)
1099 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1100 if (ready != building)
1101 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1102
1103 r[i].setsize(r[i].getsize()+s);
1104 }
1105
1108 {
1109 if (build_mode!=random)
1110 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1111 if (ready != building)
1112 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1113
1114 // compute total size, check positivity
1115 size_type total=0;
1116 for (size_type i=0; i<n; i++)
1117 {
1118 total += r[i].getsize();
1119 }
1120
1121 if(nnz_ == 0)
1122 // allocate/check memory
1123 allocate(n,m,total,false,false);
1124 else if(nnz_ < total)
1125 DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
1126 <<"sufficient for calculated nonzeros ("<<total<<"! ");
1127
1128 // set the window pointers correctly
1130
1131 // initialize j_ array with m (an invalid column index)
1132 // this indicates an unused entry
1133 for (size_type k=0; k<nnz_; k++)
1134 j_.get()[k] = m;
1136 }
1137
1139
1150 {
1151 if (build_mode!=random)
1152 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1153 if (ready==built)
1154 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1155 if (ready==building)
1156 DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
1157 if (ready==notAllocated)
1158 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1159
1160 if (col >= m)
1161 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1162
1163 // get row range
1164 size_type* const first = r[row].getindexptr();
1165 size_type* const last = first + r[row].getsize();
1166
1167 // find correct insertion position for new column index
1168 size_type* pos = std::lower_bound(first,last,col);
1169
1170 // check if index is already in row
1171 if (pos!=last && *pos == col) return;
1172
1173 // find end of already inserted column indices
1174 size_type* end = std::lower_bound(pos,last,m);
1175 if (end==last)
1176 DUNE_THROW(BCRSMatrixError,"row is too small");
1177
1178 // insert new column index at correct position
1179 std::copy_backward(pos,end,end+1);
1180 *pos = col;
1181 }
1182
1184
1191 template<typename It>
1193 {
1194 size_type row_size = r[row].size();
1195 size_type* col_begin = r[row].getindexptr();
1197 // consistency check between allocated row size and number of passed column indices
1198 if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
1199 DUNE_THROW(BCRSMatrixError,"Given size of row " << row
1200 << " (" << row_size
1201 << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
1202 std::sort(col_begin,col_end);
1203 }
1204
1207 {
1208 if (build_mode!=random)
1209 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1210 if (ready==built)
1211 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1212 if (ready==building)
1213 DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
1214 if (ready==notAllocated)
1215 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1216
1217 // check if there are undefined indices
1219 for (RowIterator i=begin(); i!=endi; ++i)
1220 {
1221 ColIterator endj = (*i).end();
1222 for (ColIterator j=(*i).begin(); j!=endj; ++j) {
1223 if (j.index() >= m) {
1224 dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
1225 <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
1226 r[i.index()].setsize(j.offset());
1227 break;
1228 }
1229 }
1230 }
1231
1232 allocateData();
1234
1235 // if not, set matrix to built
1236 ready = built;
1237 }
1238
1239 //===== implicit creation interface
1240
1242
1254 {
1255#ifdef DUNE_ISTL_WITH_CHECKING
1256 if (build_mode!=implicit)
1257 DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1258 if (ready==built)
1259 DUNE_THROW(BCRSMatrixError,"matrix already built up, use operator[] for entry access now");
1260 if (ready==notAllocated)
1261 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1262 if (ready!=building)
1263 DUNE_THROW(InvalidStateException,"You may only use entry() during the 'building' stage");
1264
1265 if (row >= n)
1266 DUNE_THROW(BCRSMatrixError,"row index exceeds matrix size");
1267 if (col >= m)
1268 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1269#endif
1270
1271 size_type* begin = r[row].getindexptr();
1272 size_type* end = begin + r[row].getsize();
1273
1274 size_type* pos = std::find(begin, end, col);
1275
1276 //treat the case that there was a match in the array
1277 if (pos != end)
1278 if (*pos == col)
1279 {
1280 std::ptrdiff_t offset = pos - r[row].getindexptr();
1281 B* aptr = r[row].getptr() + offset;
1282
1283 return *aptr;
1284 }
1285
1286 //determine whether overflow has to be taken into account or not
1287 if (r[row].getsize() == avg)
1288 return overflow[std::make_pair(row,col)];
1289 else
1290 {
1291 //modify index array
1292 *end = col;
1293
1294 //do simulatenous operations on data array a
1295 std::ptrdiff_t offset = end - r[row].getindexptr();
1296 B* apos = r[row].getptr() + offset;
1297
1298 //increase rowsize
1299 r[row].setsize(r[row].getsize()+1);
1300
1301 //return reference to the newly created entry
1302 return *apos;
1303 }
1304 }
1305
1307
1318 {
1319 if (build_mode!=implicit)
1320 DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1321 if (ready==built)
1322 DUNE_THROW(BCRSMatrixError,"matrix already built up, no more need for compression");
1323 if (ready==notAllocated)
1324 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1325 if (ready!=building)
1326 DUNE_THROW(InvalidStateException,"You may only call compress() at the end of the 'building' stage");
1327
1328 //calculate statistics
1330 stats.overflow_total = overflow.size();
1331 stats.maximum = 0;
1332
1333 //get insertion iterators pointing to one before start (for later use of ++it)
1334 size_type* jiit = j_.get();
1335 B* aiit = a;
1336
1337 //get iterator to the smallest overflow element
1338 typename OverflowType::iterator oit = overflow.begin();
1339
1340 //store a copy of index pointers on which to perform sortation
1341 std::vector<size_type*> perm;
1342
1343 //iterate over all rows and copy elements into their position in the compressed array
1344 for (size_type i=0; i<n; i++)
1345 {
1346 //get old pointers into a and j and size without overflow changes
1347 size_type* begin = r[i].getindexptr();
1348 //B* apos = r[i].getptr();
1349 size_type size = r[i].getsize();
1350
1351 perm.resize(size);
1352
1353 typename std::vector<size_type*>::iterator it = perm.begin();
1354 for (size_type* iit = begin; iit < begin + size; ++iit, ++it)
1355 *it = iit;
1356
1357 //sort permutation array
1358 std::sort(perm.begin(),perm.end(),PointerCompare<size_type>());
1359
1360 //change row window pointer to their new positions
1361 r[i].setindexptr(jiit);
1362 r[i].setptr(aiit);
1363
1364 for (it = perm.begin(); it != perm.end(); ++it)
1365 {
1366 //check whether there are elements in the overflow area which take precedence
1367 while ((oit!=overflow.end()) && (oit->first < std::make_pair(i,**it)))
1368 {
1369 //check whether there is enough memory to write to
1370 if (jiit > begin)
1372 "Allocated memory for BCRSMatrix exhausted during compress()!"
1373 "Please increase either the average number of entries per row or the overflow fraction."
1374 );
1375 //copy an element from the overflow area to the insertion position in a and j
1376 *jiit = oit->first.second;
1377 ++jiit;
1378 *aiit = oit->second;
1379 ++aiit;
1380 ++oit;
1381 r[i].setsize(r[i].getsize()+1);
1382 }
1383
1384 //check whether there is enough memory to write to
1385 if (jiit > begin)
1387 "Allocated memory for BCRSMatrix exhausted during compress()!"
1388 "Please increase either the average number of entries per row or the overflow fraction."
1389 );
1390
1391 //copy element from array
1392 *jiit = **it;
1393 ++jiit;
1394 B* apos = *it - j_.get() + a;
1395 *aiit = *apos;
1396 ++aiit;
1397 }
1398
1399 //copy remaining elements from the overflow area
1400 while ((oit!=overflow.end()) && (oit->first.first == i))
1401 {
1402 //check whether there is enough memory to write to
1403 if (jiit > begin)
1405 "Allocated memory for BCRSMatrix exhausted during compress()!"
1406 "Please increase either the average number of entries per row or the overflow fraction."
1407 );
1408
1409 //copy and element from the overflow area to the insertion position in a and j
1410 *jiit = oit->first.second;
1411 ++jiit;
1412 *aiit = oit->second;
1413 ++aiit;
1414 ++oit;
1415 r[i].setsize(r[i].getsize()+1);
1416 }
1417
1418 // update maximum row size
1419 if (r[i].getsize()>stats.maximum)
1420 stats.maximum = r[i].getsize();
1421 }
1422
1423 // overflow area may be cleared
1424 overflow.clear();
1425
1426 //determine average number of entries and memory usage
1427 std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
1428 nnz_ = diff;
1429 stats.avg = (double) (nnz_) / (double) n;
1430 stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
1431
1432 //matrix is now built
1433 ready = built;
1434
1435 return stats;
1436 }
1437
1438 //===== vector space arithmetic
1439
1442 {
1443#ifdef DUNE_ISTL_WITH_CHECKING
1444 if (ready != built)
1445 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1446#endif
1447
1448 if (nnz_ > 0)
1449 {
1450 // process 1D array
1451 for (size_type i=0; i<nnz_; i++)
1452 a[i] *= k;
1453 }
1454 else
1455 {
1457 for (RowIterator i=begin(); i!=endi; ++i)
1458 {
1459 ColIterator endj = (*i).end();
1460 for (ColIterator j=(*i).begin(); j!=endj; ++j)
1461 (*j) *= k;
1462 }
1463 }
1464
1465 return *this;
1466 }
1467
1470 {
1471#ifdef DUNE_ISTL_WITH_CHECKING
1472 if (ready != built)
1473 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1474#endif
1475
1476 if (nnz_ > 0)
1477 {
1478 // process 1D array
1479 for (size_type i=0; i<nnz_; i++)
1480 a[i] /= k;
1481 }
1482 else
1483 {
1485 for (RowIterator i=begin(); i!=endi; ++i)
1486 {
1487 ColIterator endj = (*i).end();
1488 for (ColIterator j=(*i).begin(); j!=endj; ++j)
1489 (*j) /= k;
1490 }
1491 }
1492
1493 return *this;
1494 }
1495
1496
1503 {
1504#ifdef DUNE_ISTL_WITH_CHECKING
1505 if (ready != built || b.ready != built)
1506 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1507 if(N()!=b.N() || M() != b.M())
1508 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1509#endif
1511 ConstRowIterator j=b.begin();
1512 for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1513 i->operator+=(*j);
1514 }
1515
1516 return *this;
1517 }
1518
1525 {
1526#ifdef DUNE_ISTL_WITH_CHECKING
1527 if (ready != built || b.ready != built)
1528 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1529 if(N()!=b.N() || M() != b.M())
1530 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1531#endif
1533 ConstRowIterator j=b.begin();
1534 for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1535 i->operator-=(*j);
1536 }
1537
1538 return *this;
1539 }
1540
1550 {
1551#ifdef DUNE_ISTL_WITH_CHECKING
1552 if (ready != built || b.ready != built)
1553 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1554 if(N()!=b.N() || M() != b.M())
1555 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1556#endif
1558 ConstRowIterator j=b.begin();
1559 for(RowIterator i=begin(); i!=endi; ++i, ++j)
1560 i->axpy(alpha, *j);
1561
1562 return *this;
1563 }
1564
1565 //===== linear maps
1566
1568 template<class X, class Y>
1569 void mv (const X& x, Y& y) const
1570 {
1571#ifdef DUNE_ISTL_WITH_CHECKING
1572 if (ready != built)
1573 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1574 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,
1575 "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
1576 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,
1577 "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
1578#endif
1580 for (ConstRowIterator i=begin(); i!=endi; ++i)
1581 {
1582 y[i.index()]=0;
1583 ConstColIterator endj = (*i).end();
1584 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1585 (*j).umv(x[j.index()],y[i.index()]);
1586 }
1587 }
1588
1590 template<class X, class Y>
1591 void umv (const X& x, Y& y) const
1592 {
1593#ifdef DUNE_ISTL_WITH_CHECKING
1594 if (ready != built)
1595 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1596 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1597 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1598#endif
1600 for (ConstRowIterator i=begin(); i!=endi; ++i)
1601 {
1602 ConstColIterator endj = (*i).end();
1603 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1604 (*j).umv(x[j.index()],y[i.index()]);
1605 }
1606 }
1607
1609 template<class X, class Y>
1610 void mmv (const X& x, Y& y) const
1611 {
1612#ifdef DUNE_ISTL_WITH_CHECKING
1613 if (ready != built)
1614 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1615 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1616 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1617#endif
1619 for (ConstRowIterator i=begin(); i!=endi; ++i)
1620 {
1621 ConstColIterator endj = (*i).end();
1622 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1623 (*j).mmv(x[j.index()],y[i.index()]);
1624 }
1625 }
1626
1628 template<class X, class Y, class F>
1629 void usmv (F&& alpha, const X& x, Y& y) const
1630 {
1631#ifdef DUNE_ISTL_WITH_CHECKING
1632 if (ready != built)
1633 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1634 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1635 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1636#endif
1638 for (ConstRowIterator i=begin(); i!=endi; ++i)
1639 {
1640 ConstColIterator endj = (*i).end();
1641 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1642 (*j).usmv(alpha,x[j.index()],y[i.index()]);
1643 }
1644 }
1645
1647 template<class X, class Y>
1648 void mtv (const X& x, Y& y) const
1649 {
1650#ifdef DUNE_ISTL_WITH_CHECKING
1651 if (ready != built)
1652 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1653 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1654 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1655#endif
1656 for(size_type i=0; i<y.N(); ++i)
1657 y[i]=0;
1658 umtv(x,y);
1659 }
1660
1662 template<class X, class Y>
1663 void umtv (const X& x, Y& y) const
1664 {
1665#ifdef DUNE_ISTL_WITH_CHECKING
1666 if (ready != built)
1667 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1668 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1669 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1670#endif
1672 for (ConstRowIterator i=begin(); i!=endi; ++i)
1673 {
1674 ConstColIterator endj = (*i).end();
1675 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1676 (*j).umtv(x[i.index()],y[j.index()]);
1677 }
1678 }
1679
1681 template<class X, class Y>
1682 void mmtv (const X& x, Y& y) const
1683 {
1684#ifdef DUNE_ISTL_WITH_CHECKING
1685 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1686 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1687#endif
1689 for (ConstRowIterator i=begin(); i!=endi; ++i)
1690 {
1691 ConstColIterator endj = (*i).end();
1692 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1693 (*j).mmtv(x[i.index()],y[j.index()]);
1694 }
1695 }
1696
1698 template<class X, class Y>
1699 void usmtv (const field_type& alpha, const X& x, Y& y) const
1700 {
1701#ifdef DUNE_ISTL_WITH_CHECKING
1702 if (ready != built)
1703 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1704 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1705 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1706#endif
1708 for (ConstRowIterator i=begin(); i!=endi; ++i)
1709 {
1710 ConstColIterator endj = (*i).end();
1711 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1712 (*j).usmtv(alpha,x[i.index()],y[j.index()]);
1713 }
1714 }
1715
1717 template<class X, class Y>
1718 void umhv (const X& x, Y& y) const
1719 {
1720#ifdef DUNE_ISTL_WITH_CHECKING
1721 if (ready != built)
1722 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1723 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1724 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1725#endif
1727 for (ConstRowIterator i=begin(); i!=endi; ++i)
1728 {
1729 ConstColIterator endj = (*i).end();
1730 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1731 (*j).umhv(x[i.index()],y[j.index()]);
1732 }
1733 }
1734
1736 template<class X, class Y>
1737 void mmhv (const X& x, Y& y) const
1738 {
1739#ifdef DUNE_ISTL_WITH_CHECKING
1740 if (ready != built)
1741 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1742 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1743 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1744#endif
1746 for (ConstRowIterator i=begin(); i!=endi; ++i)
1747 {
1748 ConstColIterator endj = (*i).end();
1749 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1750 (*j).mmhv(x[i.index()],y[j.index()]);
1751 }
1752 }
1753
1755 template<class X, class Y>
1756 void usmhv (const field_type& alpha, const X& x, Y& y) const
1757 {
1758#ifdef DUNE_ISTL_WITH_CHECKING
1759 if (ready != built)
1760 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1761 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1762 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1763#endif
1765 for (ConstRowIterator i=begin(); i!=endi; ++i)
1766 {
1767 ConstColIterator endj = (*i).end();
1768 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1769 (*j).usmhv(alpha,x[i.index()],y[j.index()]);
1770 }
1771 }
1772
1773
1774 //===== norms
1775
1778 {
1779#ifdef DUNE_ISTL_WITH_CHECKING
1780 if (ready != built)
1781 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1782#endif
1783
1784 double sum=0;
1785
1787 for (ConstRowIterator i=begin(); i!=endi; ++i)
1788 {
1789 ConstColIterator endj = (*i).end();
1790 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1791 sum += (*j).frobenius_norm2();
1792 }
1793
1794 return sum;
1795 }
1796
1799 {
1800 return sqrt(frobenius_norm2());
1801 }
1802
1804 template <typename ft = field_type,
1805 typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
1807 if (ready != built)
1808 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1809
1810 using real_type = typename FieldTraits<ft>::real_type;
1811 using std::max;
1812
1813 real_type norm = 0;
1814 for (auto const &x : *this) {
1815 real_type sum = 0;
1816 for (auto const &y : x)
1817 sum += y.infinity_norm();
1818 norm = max(sum, norm);
1819 }
1820 return norm;
1821 }
1822
1824 template <typename ft = field_type,
1825 typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
1827 if (ready != built)
1828 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1829
1830 using real_type = typename FieldTraits<ft>::real_type;
1831 using std::max;
1832
1833 real_type norm = 0;
1834 for (auto const &x : *this) {
1835 real_type sum = 0;
1836 for (auto const &y : x)
1837 sum += y.infinity_norm_real();
1838 norm = max(sum, norm);
1839 }
1840 return norm;
1841 }
1842
1844 template <typename ft = field_type,
1845 typename std::enable_if<has_nan<ft>::value, int>::type = 0>
1847 if (ready != built)
1848 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1849
1850 using real_type = typename FieldTraits<ft>::real_type;
1851 using std::max;
1852
1853 real_type norm = 0;
1854 real_type isNaN = 1;
1855 for (auto const &x : *this) {
1856 real_type sum = 0;
1857 for (auto const &y : x)
1858 sum += y.infinity_norm();
1859 norm = max(sum, norm);
1860 isNaN += sum;
1861 }
1862 isNaN /= isNaN;
1863 return norm * isNaN;
1864 }
1865
1867 template <typename ft = field_type,
1868 typename std::enable_if<has_nan<ft>::value, int>::type = 0>
1870 if (ready != built)
1871 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1872
1873 using real_type = typename FieldTraits<ft>::real_type;
1874 using std::max;
1875
1876 real_type norm = 0;
1877 real_type isNaN = 1;
1878 for (auto const &x : *this) {
1879 real_type sum = 0;
1880 for (auto const &y : x)
1881 sum += y.infinity_norm_real();
1882 norm = max(sum, norm);
1883 isNaN += sum;
1884 }
1885 isNaN /= isNaN;
1886 return norm * isNaN;
1887 }
1888
1889 //===== sizes
1890
1892 size_type N () const
1893 {
1894 return n;
1895 }
1896
1898 size_type M () const
1899 {
1900 return m;
1901 }
1902
1905 {
1906 return nnz_;
1907 }
1908
1911 {
1912 return ready;
1913 }
1914
1917 {
1918 return build_mode;
1919 }
1920
1921 //===== query
1922
1924 bool exists (size_type i, size_type j) const
1925 {
1926#ifdef DUNE_ISTL_WITH_CHECKING
1927 if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
1928 if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
1929#endif
1930 return (r[i].size() && r[i].find(j) != r[i].end());
1931 }
1932
1933
1934 protected:
1935 // state information
1936 BuildMode build_mode; // row wise or whole matrix
1937 BuildStage ready; // indicate the stage the matrix building is in
1938
1939 // The allocator used for memory management
1940 typename A::template rebind<B>::other allocator_;
1941
1943
1945
1946 // size of the matrix
1947 size_type n; // number of rows
1948 size_type m; // number of columns
1949 size_type nnz_; // number of nonzeroes contained in the matrix
1950 size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
1951 // zero means that memory is allocated separately for each row.
1952
1953 // the rows are dynamically allocated
1954 row_type* r; // [n] the individual rows having pointers into a,j arrays
1955
1956 // dynamically allocated memory
1957 B* a; // [allocationSize] non-zero entries of the matrix in row-wise ordering
1958 // If a single array of column indices is used, it can be shared
1959 // between different matrices with the same sparsity pattern
1960 std::shared_ptr<size_type> j_; // [allocationSize] column indices of entries
1961
1962 // additional data is needed in implicit buildmode
1965
1966 typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
1968
1970 {
1971 row_type current_row(a,j_.get(),0); // Pointers to current row data
1972 for (size_type i=0; i<n; i++, ++row) {
1973 // set row i
1974 size_type s = row->getsize();
1975
1976 if (s>0) {
1977 // setup pointers and size
1978 r[i].set(s,current_row.getptr(), current_row.getindexptr());
1979 // update pointer for next row
1980 current_row.setptr(current_row.getptr()+s);
1981 current_row.setindexptr(current_row.getindexptr()+s);
1982 } else{
1983 // empty row
1984 r[i].set(0,nullptr,nullptr);
1985 }
1986 }
1987 }
1988
1990
1995 {
1996 size_type* jptr = j_.get();
1997 for (size_type i=0; i<n; ++i, ++row) {
1998 // set row i
1999 size_type s = row->getsize();
2000
2001 if (s>0) {
2002 // setup pointers and size
2003 r[i].setsize(s);
2004 r[i].setindexptr(jptr);
2005 } else{
2006 // empty row
2007 r[i].set(0,nullptr,nullptr);
2008 }
2009
2010 // advance position in global array
2011 jptr += s;
2012 }
2013 }
2014
2016
2021 {
2022 B* aptr = a;
2023 for (size_type i=0; i<n; ++i) {
2024 // set row i
2025 if (r[i].getsize() > 0) {
2026 // setup pointers and size
2027 r[i].setptr(aptr);
2028 } else{
2029 // empty row
2030 r[i].set(0,nullptr,nullptr);
2031 }
2032
2033 // advance position in global array
2034 aptr += r[i].getsize();
2035 }
2036 }
2037
2040 {
2041 setWindowPointers(Mat.begin());
2042
2043 // copy data
2044 for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
2045
2046 // finish off
2047 build_mode = row_wise; // dummy
2048 ready = built;
2049 }
2050
2057 {
2058
2059 if (notAllocated)
2060 return;
2061
2062 if (allocationSize_>0)
2063 {
2064 // a,j_ have been allocated as one long vector
2065 j_.reset();
2066 if (a)
2067 {
2068 for(B *aiter=a+(allocationSize_-1), *aend=a-1; aiter!=aend; --aiter)
2069 allocator_.destroy(aiter);
2070 allocator_.deallocate(a,allocationSize_);
2071 a = nullptr;
2072 }
2073 }
2074 else if (r)
2075 {
2076 // check if memory for rows have been allocated individually
2077 for (size_type i=0; i<n; i++)
2078 if (r[i].getsize()>0)
2079 {
2080 for (B *col=r[i].getptr()+(r[i].getsize()-1),
2081 *colend = r[i].getptr()-1; col!=colend; --col) {
2082 allocator_.destroy(col);
2083 }
2084 sizeAllocator_.deallocate(r[i].getindexptr(),1);
2085 allocator_.deallocate(r[i].getptr(),1);
2086 // clear out row data in case we don't want to deallocate the rows
2087 // otherwise we might run into a double free problem here later
2088 r[i].set(0,nullptr,nullptr);
2089 }
2090 }
2091
2092 // deallocate the rows
2093 if (n>0 && deallocateRows && r) {
2094 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
2095 rowAllocator_.destroy(riter);
2096 rowAllocator_.deallocate(r,n);
2097 r = nullptr;
2098 }
2099
2100 // Mark matrix as not built at all.
2102
2103 }
2104
2107 {
2108 typename A::template rebind<size_type>::other& sizeAllocator_;
2109
2110 public:
2112 : sizeAllocator_(sizeAllocator)
2113 {}
2114
2115 void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); }
2116 };
2117
2118
2137 {
2138 // Store size
2139 n = rows;
2140 m = columns;
2143
2144 // allocate rows
2145 if(allocateRows) {
2146 if (n>0) {
2147 if (r)
2148 DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
2149 r = rowAllocator_.allocate(rows);
2150 }else{
2151 r = 0;
2152 }
2153 }
2154
2155 // allocate a and j_ array
2156 if (allocate_data)
2157 allocateData();
2158 if (allocationSize_>0) {
2159 // allocate column indices only if not yet present (enable sharing)
2160 if (!j_.get())
2162 }else{
2163 j_.reset();
2164 for(row_type* ri=r; ri!=r+rows; ++ri)
2165 rowAllocator_.construct(ri, row_type());
2166 }
2167
2168 // Mark the matrix as not built.
2169 ready = building;
2170 }
2171
2173 {
2174 if (a)
2175 DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
2176 if (allocationSize_>0) {
2177 a = allocator_.allocate(allocationSize_);
2178 // use placement new to call constructor that allocates
2179 // additional memory.
2180 new (a) B[allocationSize_];
2181 } else {
2182 a = nullptr;
2183 }
2184 }
2185
2192 {
2193 if (build_mode != implicit)
2194 DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
2195 if (ready != notAllocated)
2196 DUNE_THROW(InvalidStateException,"memory has already been allocated");
2197
2198 // check to make sure the user has actually set the parameters
2199 if (overflowsize < 0)
2200 DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
2201 //calculate size of overflow region, add buffer for row sort!
2203 allocationSize_ = _n*avg + osize;
2204
2205 allocate(_n, _m, allocationSize_,true,true);
2206
2207 //set row pointers correctly
2208 size_type* jptr = j_.get() + osize;
2209 B* aptr = a + osize;
2210 for (size_type i=0; i<n; i++)
2211 {
2212 r[i].set(0,aptr,jptr);
2213 jptr = jptr + avg;
2214 aptr = aptr + avg;
2215 }
2216
2217 ready = building;
2218 }
2219 };
2220
2221
2224} // end namespace
2225
2226#endif
Some handy generic functions for ISTL matrices.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Col col
Definition matrixmatrix.hh:347
Definition basearray.hh:19
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition dependency.hh:292
size_type size() const
number of blocks in the array (are of size 1 here)
Definition basearray.hh:756
Definition matrixutils.hh:231
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
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
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 sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:412
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition bcrsmatrix.hh:1924
BuildStage buildStage() const
The current build stage of the matrix.
Definition bcrsmatrix.hh:1910
Iterator begin()
Get iterator to first row.
Definition bcrsmatrix.hh:623
void copyWindowStructure(const BCRSMatrix &Mat)
Copy the window structure from another matrix.
Definition bcrsmatrix.hh:2039
B & entry(size_type row, size_type col)
Returns reference to entry (row,col) of the matrix.
Definition bcrsmatrix.hh:1253
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition bcrsmatrix.hh:1737
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition bcrsmatrix.hh:1756
void umtv(const X &x, Y &y) const
y += A^T x
Definition bcrsmatrix.hh:1663
size_type m
Definition bcrsmatrix.hh:1948
RealRowIterator< const row_type > const_iterator
The const iterator over the matrix rows.
Definition bcrsmatrix.hh:655
double overflowsize
Definition bcrsmatrix.hh:1964
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:732
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:2136
B::field_type field_type
export the type representing the field
Definition bcrsmatrix.hh:434
BCRSMatrix & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition bcrsmatrix.hh:1549
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition bcrsmatrix.hh:1826
~BCRSMatrix()
destructor
Definition bcrsmatrix.hh:780
void allocateData()
Definition bcrsmatrix.hh:2172
void deallocate(bool deallocateRows=true)
deallocate memory of the matrix.
Definition bcrsmatrix.hh:2056
Iterator RowIterator
rename the iterators for easier access
Definition bcrsmatrix.hh:649
row_type & operator[](size_type i)
random access to the rows
Definition bcrsmatrix.hh:497
BCRSMatrix()
an empty matrix
Definition bcrsmatrix.hh:697
void endrowsizes()
indicate that size of all rows is defined
Definition bcrsmatrix.hh:1107
void incrementrowsize(size_type i, size_type s=1)
increment size of row i by s (1 by default)
Definition bcrsmatrix.hh:1096
void mtv(const X &x, Y &y) const
y = A^T x
Definition bcrsmatrix.hh:1648
void umhv(const X &x, Y &y) const
y += A^H x
Definition bcrsmatrix.hh:1718
size_type nonzeroes() const
number of blocks that are stored (the number of blocks that possibly are nonzero)
Definition bcrsmatrix.hh:1904
size_type allocationSize_
Definition bcrsmatrix.hh:1950
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition bcrsmatrix.hh:686
BuildStage ready
Definition bcrsmatrix.hh:1937
BuildMode build_mode
Definition bcrsmatrix.hh:1936
void setrowsize(size_type i, size_type s)
Set number of indices in row i to s.
Definition bcrsmatrix.hh:1075
RealRowIterator< row_type > Iterator
Definition bcrsmatrix.hh:620
size_type nnz_
Definition bcrsmatrix.hh:1949
BCRSMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition bcrsmatrix.hh:1441
RealRowIterator< row_type > iterator
The iterator over the (mutable matrix rows.
Definition bcrsmatrix.hh:619
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition bcrsmatrix.hh:1699
ConstIterator beforeBegin() const
Definition bcrsmatrix.hh:680
RealRowIterator< const row_type > ConstIterator
Definition bcrsmatrix.hh:656
Iterator beforeBegin()
Definition bcrsmatrix.hh:643
B * a
Definition bcrsmatrix.hh:1957
BuildMode
we support two modes
Definition bcrsmatrix.hh:458
@ implicit
Build entries randomly with an educated guess for the number of entries per row.
Definition bcrsmatrix.hh:487
@ unknown
Build mode not set!
Definition bcrsmatrix.hh:491
@ random
Build entries randomly.
Definition bcrsmatrix.hh:478
@ row_wise
Build in a row-wise manner.
Definition bcrsmatrix.hh:469
BCRSMatrix(size_type _n, size_type _m, size_type _nnz, BuildMode bm)
matrix with known number of nonzeroes
Definition bcrsmatrix.hh:704
::Dune::CompressionStatistics< size_type > CompressionStatistics
The type for the statistics object returned by compress()
Definition bcrsmatrix.hh:449
BCRSMatrix & operator-=(const BCRSMatrix &b)
Subtract the entries of another matrix from this one.
Definition bcrsmatrix.hh:1524
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition bcrsmatrix.hh:753
Iterator end()
Get iterator to one beyond last row.
Definition bcrsmatrix.hh:629
row_type * r
Definition bcrsmatrix.hh:1954
void setIndices(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition bcrsmatrix.hh:1192
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition bcrsmatrix.hh:1149
std::map< std::pair< size_type, size_type >, B > OverflowType
Definition bcrsmatrix.hh:1966
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition bcrsmatrix.hh:652
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition bcrsmatrix.hh:1798
A::size_type size_type
The type for the index access and the size.
Definition bcrsmatrix.hh:446
BCRSMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition bcrsmatrix.hh:1469
OverflowType overflow
Definition bcrsmatrix.hh:1967
BCRSMatrix & operator+=(const BCRSMatrix &b)
Add the entries of another matrix to this one.
Definition bcrsmatrix.hh:1502
CreateIterator createend()
get create iterator pointing to one after the last block
Definition bcrsmatrix.hh:1061
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition bcrsmatrix.hh:1777
Iterator beforeEnd()
Definition bcrsmatrix.hh:636
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition bcrsmatrix.hh:689
A::template rebind< B >::other allocator_
Definition bcrsmatrix.hh:1940
void usmv(F &&alpha, const X &x, Y &y) const
y += alpha A x
Definition bcrsmatrix.hh:1629
CompressedBlockVectorWindow< B, A > row_type
implement row_type with compressed vector
Definition bcrsmatrix.hh:443
size_type getrowsize(size_type i) const
get current number of indices in row i
Definition bcrsmatrix.hh:1086
size_type M() const
number of columns (counted in blocks)
Definition bcrsmatrix.hh:1898
size_type n
Definition bcrsmatrix.hh:1947
CreateIterator createbegin()
get initial create iterator
Definition bcrsmatrix.hh:1055
BuildStage
Definition bcrsmatrix.hh:415
@ rowSizesBuilt
The row sizes of the matrix are known.
Definition bcrsmatrix.hh:426
@ built
The matrix structure is fully built.
Definition bcrsmatrix.hh:428
@ notbuilt
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition bcrsmatrix.hh:417
@ notAllocated
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition bcrsmatrix.hh:419
@ building
Matrix is currently being built, some memory has been allocated, build mode and size are fixed.
Definition bcrsmatrix.hh:421
BuildMode buildMode() const
The currently selected build mode of the matrix.
Definition bcrsmatrix.hh:1916
A::template rebind< row_type >::other rowAllocator_
Definition bcrsmatrix.hh:1942
void mmv(const X &x, Y &y) const
y -= A x
Definition bcrsmatrix.hh:1610
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition bcrsmatrix.hh:1806
void mv(const X &x, Y &y) const
y = A x
Definition bcrsmatrix.hh:1569
B block_type
export the type representing the components
Definition bcrsmatrix.hh:437
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition bcrsmatrix.hh:1682
size_type avg
Definition bcrsmatrix.hh:1963
A::template rebind< size_type >::other sizeAllocator_
Definition bcrsmatrix.hh:1944
void umv(const X &x, Y &y) const
y += A x
Definition bcrsmatrix.hh:1591
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:2191
BCRSMatrix(size_type _n, size_type _m, BuildMode bm)
matrix with unknown number of nonzeroes
Definition bcrsmatrix.hh:713
void endindices()
indicate that all indices are defined, check consistency
Definition bcrsmatrix.hh:1206
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition bcrsmatrix.hh:1317
@ blocklevel
The number of blocklevels the matrix contains.
Definition bcrsmatrix.hh:454
void setDataPointers()
Set data pointers for all rows.
Definition bcrsmatrix.hh:2020
size_type N() const
number of rows (counted in blocks)
Definition bcrsmatrix.hh:1892
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition bcrsmatrix.hh:789
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition bcrsmatrix.hh:817
std::shared_ptr< size_type > j_
Definition bcrsmatrix.hh:1960
void setWindowPointers(ConstRowIterator row)
Definition bcrsmatrix.hh:1969
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition bcrsmatrix.hh:867
void setColumnPointers(ConstRowIterator row)
Copy row sizes from iterator range starting at row and set column index pointers for all rows.
Definition bcrsmatrix.hh:1994
ConstIterator end() const
Get const iterator to one beyond last row.
Definition bcrsmatrix.hh:666
ConstIterator begin() const
Get const iterator to first row.
Definition bcrsmatrix.hh:660
void setImplicitBuildModeParameters(size_type _avg, double _overflow)
Set parameters needed for creation in implicit build mode.
Definition bcrsmatrix.hh:845
A allocator_type
export the allocator type
Definition bcrsmatrix.hh:440
ConstIterator beforeEnd() const
Definition bcrsmatrix.hh:673
Iterator access to matrix rows
Definition bcrsmatrix.hh:527
RealRowIterator()
empty constructor, use with care!
Definition bcrsmatrix.hh:544
bool equals(const RealRowIterator< ValueType > &other) const
equality
Definition bcrsmatrix.hh:572
std::remove_const< T >::type ValueType
The unqualified value type.
Definition bcrsmatrix.hh:531
RealRowIterator(const RealRowIterator< ValueType > &it)
Definition bcrsmatrix.hh:548
bool equals(const RealRowIterator< const ValueType > &other) const
equality
Definition bcrsmatrix.hh:579
RealRowIterator(row_type *_p, size_type _i)
constructor
Definition bcrsmatrix.hh:539
std::ptrdiff_t distanceTo(const RealRowIterator< const ValueType > &other) const
Definition bcrsmatrix.hh:565
size_type index() const
return index
Definition bcrsmatrix.hh:554
std::ptrdiff_t distanceTo(const RealRowIterator< ValueType > &other) const
Definition bcrsmatrix.hh:559
Iterator class for sequential creation of blocks
Definition bcrsmatrix.hh:918
bool operator==(const CreateIterator &it) const
equality
Definition bcrsmatrix.hh:1010
CreateIterator & operator++()
prefix increment
Definition bcrsmatrix.hh:935
size_type index() const
The number of the row that the iterator currently points to.
Definition bcrsmatrix.hh:1016
bool operator!=(const CreateIterator &it) const
inequality
Definition bcrsmatrix.hh:1004
CreateIterator(BCRSMatrix &_Mat, size_type _i)
constructor
Definition bcrsmatrix.hh:921
void insert(size_type j)
put column index in row
Definition bcrsmatrix.hh:1022
bool contains(size_type j)
return true if column index is in row
Definition bcrsmatrix.hh:1028
size_type size() const
Get the current row size.
Definition bcrsmatrix.hh:1037
Class used by shared_ptr to deallocate memory using the proper allocator.
Definition bcrsmatrix.hh:2107
Deallocator(typename A::template rebind< size_type >::other &sizeAllocator)
Definition bcrsmatrix.hh:2111
void operator()(size_type *p)
Definition bcrsmatrix.hh:2115
Definition bvector.hh:1050
void set(size_type _n, B *_p, size_type *_j)
set size and pointer
Definition bvector.hh:1149
size_type * getindexptr()
get pointer
Definition bvector.hh:1181
void setsize(size_type _n)
set size only
Definition bvector.hh:1157
void setptr(B *_p)
set pointer only
Definition bvector.hh:1163
B * getptr()
get pointer
Definition bvector.hh:1175
size_type getsize() const
get size
Definition bvector.hh:1198
void setindexptr(size_type *_j)
set pointer only
Definition bvector.hh:1169
Error specific to BCRSMatrix.
Definition istlexception.hh:21
The overflow error used during implicit BCRSMatrix construction was exhausted.
Definition istlexception.hh:34
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