DUNE PDELab (git)

bcrsmatrix.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5
6#ifndef DUNE_ISTL_BCRSMATRIX_HH
7#define DUNE_ISTL_BCRSMATRIX_HH
8
9#include <cmath>
10#include <complex>
11#include <set>
12#include <iostream>
13#include <algorithm>
14#include <numeric>
15#include <vector>
16#include <map>
17#include <memory>
18
19#include "istlexception.hh"
20#include "bvector.hh"
21#include "matrixutils.hh"
28
30
35namespace Dune {
36
76 template<typename M>
77 struct MatrixDimension;
78
80
86 template<typename size_type>
88 {
90 double avg;
92 size_type maximum;
94 size_type overflow_total;
96
99 double mem_ratio;
100 };
101
103
115 template<class M_>
117 {
118
119 public:
120
122 typedef M_ Matrix;
123
126
128 typedef typename Matrix::size_type size_type;
129
131
137 {
138
139 public:
140
143 {
144 return _m.entry(_i,j);
145 }
146
147#ifndef DOXYGEN
148
150 : _m(m)
151 , _i(i)
152 {}
153
154#endif
155
156 private:
157
158 Matrix& _m;
159 size_type _i;
160
161 };
162
164
171 : _m(m)
172 {
173 if (m.buildMode() != Matrix::implicit)
174 DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix in implicit build mode");
175 if (m.buildStage() != Matrix::building)
176 DUNE_THROW(BCRSMatrixError,"You can only create an ImplicitBuilder for a matrix with set size that has not been compressed() yet");
177 }
178
180
194 ImplicitMatrixBuilder(Matrix& m, size_type rows, size_type cols, size_type avg_cols_per_row, double overflow_fraction)
195 : _m(m)
196 {
197 if (m.buildStage() != Matrix::notAllocated)
198 DUNE_THROW(BCRSMatrixError,"You can only set up a matrix for this ImplicitBuilder if it has no memory allocated yet");
199 m.setBuildMode(Matrix::implicit);
200 m.setImplicitBuildModeParameters(avg_cols_per_row,overflow_fraction);
201 m.setSize(rows,cols);
202 }
203
206 {
207 return row_object(_m,i);
208 }
209
211 size_type N() const
212 {
213 return _m.N();
214 }
215
217 size_type M() const
218 {
219 return _m.M();
220 }
221
222 private:
223
224 Matrix& _m;
225
226 };
227
464 template<class B, class A=std::allocator<B> >
466 {
467 friend struct MatrixDimension<BCRSMatrix>;
468 public:
482 built=3
483 };
484
485 //===== type definitions and constants
486
488 using field_type = typename Imp::BlockTraits<B>::field_type;
489
491 typedef B block_type;
492
494 typedef A allocator_type;
495
497 typedef typename A::size_type size_type;
498
500 typedef Imp::CompressedBlockVectorWindow<B,size_type> row_type;
501
503 typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
504
539 unknown
540 };
541
542 //===== random access interface to rows of the matrix
543
546 {
547#ifdef DUNE_ISTL_WITH_CHECKING
548 if (build_mode == implicit && ready != built)
549 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
550 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
551 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
552#endif
553 return r[i];
554 }
555
558 {
559#ifdef DUNE_ISTL_WITH_CHECKING
560 if (build_mode == implicit && ready != built)
561 DUNE_THROW(BCRSMatrixError,"You cannot use operator[] in implicit build mode before calling compress()");
562 if (built!=ready) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
563 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
564#endif
565 return r[i];
566 }
567
568
569 //===== iterator interface to rows of the matrix
570
572 template<class T>
574 : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
575 {
576
577 public:
579 typedef typename std::remove_const<T>::type ValueType;
580
583 friend class RealRowIterator<const ValueType>;
584 friend class RealRowIterator<ValueType>;
585
588 : p(_p), i(_i)
589 {}
590
593 : p(0), i(0)
594 {}
595
597 : p(it.p), i(it.i)
598 {}
599
600
603 {
604 return i;
605 }
606
607 std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
608 {
609 assert(other.p==p);
610 return (other.i-i);
611 }
612
613 std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
614 {
615 assert(other.p==p);
616 return (other.i-i);
617 }
618
620 bool equals (const RealRowIterator<ValueType>& other) const
621 {
622 assert(other.p==p);
623 return i==other.i;
624 }
625
627 bool equals (const RealRowIterator<const ValueType>& other) const
628 {
629 assert(other.p==p);
630 return i==other.i;
631 }
632
633 private:
635 void increment()
636 {
637 ++i;
638 }
639
641 void decrement()
642 {
643 --i;
644 }
645
646 void advance(std::ptrdiff_t diff)
647 {
648 i+=diff;
649 }
650
651 T& elementAt(std::ptrdiff_t diff) const
652 {
653 return p[i+diff];
654 }
655
657 row_type& dereference () const
658 {
659 return p[i];
660 }
661
662 row_type* p;
663 size_type i;
664 };
665
669
672 {
673 return Iterator(r,0);
674 }
675
678 {
679 return Iterator(r,n);
680 }
681
685 {
686 return Iterator(r,n-1);
687 }
688
692 {
693 return Iterator(r,-1);
694 }
695
698
700 typedef typename row_type::Iterator ColIterator;
701
705
706
709 {
710 return ConstIterator(r,0);
711 }
712
715 {
716 return ConstIterator(r,n);
717 }
718
722 {
723 return ConstIterator(r,n-1);
724 }
725
729 {
730 return ConstIterator(r,-1);
731 }
732
735
737 typedef typename row_type::ConstIterator ConstColIterator;
738
739 //===== constructors & resizers
740
741 // we use a negative compressionBufferSize to indicate that the implicit
742 // mode parameters have not been set yet
743
746 : build_mode(unknown), ready(notAllocated), n(0), m(0), nnz_(0),
747 allocationSize_(0), r(0), a(0),
748 avg(0), compressionBufferSize_(-1.0)
749 {}
750
753 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
754 allocationSize_(0), r(0), a(0),
755 avg(0), compressionBufferSize_(-1.0)
756 {
757 allocate(_n, _m, _nnz,true,false);
758 }
759
762 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
763 allocationSize_(0), r(0), a(0),
764 avg(0), compressionBufferSize_(-1.0)
765 {
766 allocate(_n, _m,0,true,false);
767 }
768
770
781 BCRSMatrix (size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
782 : build_mode(bm), ready(notAllocated), n(0), m(0), nnz_(0),
783 allocationSize_(0), r(0), a(0),
784 avg(_avg), compressionBufferSize_(compressionBufferSize)
785 {
786 if (bm != implicit)
787 DUNE_THROW(BCRSMatrixError,"Only call this constructor when using the implicit build mode");
788 // Prevent user from setting a negative compression buffer size:
789 // 1) It doesn't make sense
790 // 2) We use a negative value to indicate that the parameters
791 // have not been set yet
792 if (compressionBufferSize_ < 0.0)
793 DUNE_THROW(BCRSMatrixError,"You cannot set a negative overflow fraction");
794 implicit_allocate(_n,_m);
795 }
796
803 : build_mode(Mat.build_mode), ready(notAllocated), n(0), m(0), nnz_(0),
804 allocationSize_(0), r(0), a(0),
805 avg(Mat.avg), compressionBufferSize_(Mat.compressionBufferSize_)
806 {
807 if (!(Mat.ready == notAllocated || Mat.ready == built))
808 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
809
810 // deep copy in global array
811 size_type _nnz = Mat.nonzeroes();
812
813 j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
814 allocate(Mat.n, Mat.m, _nnz, true, true);
815
816 // build window structure
818 }
819
822 {
823 deallocate();
824 }
825
831 {
832 if (ready == notAllocated)
833 {
834 build_mode = bm;
835 return;
836 }
837 if (ready == building && (build_mode == unknown || build_mode == random || build_mode == row_wise) && (bm == row_wise || bm == random))
838 build_mode = bm;
839 else
840 DUNE_THROW(InvalidStateException, "Matrix structure cannot be changed at this stage anymore (ready == "<<ready<<").");
841 }
842
858 void setSize(size_type rows, size_type columns, size_type nnz=0)
859 {
860 // deallocate already setup memory
861 deallocate();
862
863 if (build_mode == implicit)
864 {
865 if (nnz>0)
866 DUNE_THROW(Dune::BCRSMatrixError,"number of non-zeroes may not be set in implicit mode, use setImplicitBuildModeParameters() instead");
867
868 // implicit allocates differently
869 implicit_allocate(rows,columns);
870 }
871 else
872 {
873 // allocate matrix memory
874 allocate(rows, columns, nnz, true, false);
875 }
876 }
877
886 void setImplicitBuildModeParameters(size_type _avg, double compressionBufferSize)
887 {
888 // Prevent user from setting a negative compression buffer size:
889 // 1) It doesn't make sense
890 // 2) We use a negative value to indicate that the parameters
891 // have not been set yet
892 if (compressionBufferSize < 0.0)
893 DUNE_THROW(BCRSMatrixError,"You cannot set a negative compressionBufferSize value");
894
895 // make sure the parameters aren't changed after memory has been allocated
896 if (ready != notAllocated)
897 DUNE_THROW(InvalidStateException,"You cannot modify build mode parameters at this stage anymore");
898 avg = _avg;
899 compressionBufferSize_ = compressionBufferSize;
900 }
901
909 {
910 // return immediately when self-assignment
911 if (&Mat==this) return *this;
912
913 if (!((ready == notAllocated || ready == built) && (Mat.ready == notAllocated || Mat.ready == built)))
914 DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copied when both target and source are empty or fully built)");
915
916 // make it simple: ALWAYS throw away memory for a and j_
917 // and deallocate rows only if n != Mat.n
918 deallocate(n!=Mat.n);
919
920 // reallocate the rows if required
921 if (n>0 && n!=Mat.n) {
922 // free rows
923 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
924 std::allocator_traits<decltype(rowAllocator_)>::destroy(rowAllocator_, riter);
925 rowAllocator_.deallocate(r,n);
926 }
927
928 nnz_ = Mat.nonzeroes();
929
930 // allocate a, share j_
931 j_ = Mat.j_;
932 allocate(Mat.n, Mat.m, nnz_, n!=Mat.n, true);
933
934 // build window structure
936 return *this;
937 }
938
941 {
942
943 if (!(ready == notAllocated || ready == built))
944 DUNE_THROW(InvalidStateException,"Scalar assignment only works on fully built BCRSMatrix)");
945
946 for (size_type i=0; i<n; i++) r[i] = k;
947 return *this;
948 }
949
950 //===== row-wise creation interface
951
954 {
955 public:
958 : Mat(_Mat), i(_i), nnz(0), current_row(nullptr, Mat.j_.get(), 0)
959 {
960 if (Mat.build_mode == unknown && Mat.ready == building)
961 {
962 Mat.build_mode = row_wise;
963 }
964 if (i==0 && Mat.ready != building)
965 DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
966 if(Mat.build_mode!=row_wise)
967 DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
968 if(i==0 && _Mat.N()==0)
969 // empty Matrix is always built.
970 Mat.ready = built;
971 }
972
975 {
976 // this should only be called if matrix is in creation
977 if (Mat.ready != building)
978 DUNE_THROW(BCRSMatrixError,"matrix already built up");
979
980 // row i is defined through the pattern
981 // get memory for the row and initialize the j_ array
982 // this depends on the allocation mode
983
984 // compute size of the row
985 size_type s = pattern.size();
986
987 if(s>0) {
988 // update number of nonzeroes including this row
989 nnz += s;
990
991 // alloc memory / set window
992 if (Mat.nnz_ > 0)
993 {
994 // memory is allocated in one long array
995
996 // check if that memory is sufficient
997 if (nnz > Mat.nnz_)
998 DUNE_THROW(BCRSMatrixError,"allocated nnz too small");
999
1000 // set row i
1001 Mat.r[i].set(s,nullptr,current_row.getindexptr());
1002 current_row.setindexptr(current_row.getindexptr()+s);
1003 }else{
1004 // memory is allocated individually per row
1005 // allocate and set row i
1006 B* b = Mat.allocator_.allocate(s);
1007 // use placement new to call constructor that allocates
1008 // additional memory.
1009 new (b) B[s];
1010 size_type* j = Mat.sizeAllocator_.allocate(s);
1011 Mat.r[i].set(s,b,j);
1012 }
1013 }else
1014 // setup empty row
1015 Mat.r[i].set(0,nullptr,nullptr);
1016
1017 // initialize the j array for row i from pattern
1018 std::copy(pattern.cbegin(), pattern.cend(), Mat.r[i].getindexptr());
1019
1020 // now go to next row
1021 i++;
1022 pattern.clear();
1023
1024 // check if this was last row
1025 if (i==Mat.n)
1026 {
1027 Mat.ready = built;
1028 if(Mat.nnz_ > 0)
1029 {
1030 // Set nnz to the exact number of nonzero blocks inserted
1031 // as some methods rely on it
1032 Mat.nnz_ = nnz;
1033 // allocate data array
1034 Mat.allocateData();
1035 Mat.setDataPointers();
1036 }
1037 }
1038 // done
1039 return *this;
1040 }
1041
1043 bool operator!= (const CreateIterator& it) const
1044 {
1045 return (i!=it.i) || (&Mat!=&it.Mat);
1046 }
1047
1049 bool operator== (const CreateIterator& it) const
1050 {
1051 return (i==it.i) && (&Mat==&it.Mat);
1052 }
1053
1056 {
1057 return i;
1058 }
1059
1062 {
1063 pattern.insert(j);
1064 }
1065
1068 {
1069 return pattern.find(j) != pattern.end();
1070 }
1077 {
1078 return pattern.size();
1079 }
1080
1081 private:
1082 BCRSMatrix& Mat; // the matrix we are defining
1083 size_type i; // current row to be defined
1084 size_type nnz; // count total number of nonzeros
1085 typedef std::set<size_type,std::less<size_type> > PatternType;
1086 PatternType pattern; // used to compile entries in a row
1087 row_type current_row; // row pointing to the current row to setup
1088 };
1089
1091 friend class CreateIterator;
1092
1095 {
1096 return CreateIterator(*this,0);
1097 }
1098
1101 {
1102 return CreateIterator(*this,n);
1103 }
1104
1105
1106 //===== random creation interface
1107
1115 {
1116 if (build_mode!=random)
1117 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1118 if (ready != building)
1119 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1120
1121 r[i].setsize(s);
1122 }
1123
1126 {
1127#ifdef DUNE_ISTL_WITH_CHECKING
1128 if (r==0) DUNE_THROW(BCRSMatrixError,"row not initialized yet");
1129 if (i>=n) DUNE_THROW(BCRSMatrixError,"index out of range");
1130#endif
1131 return r[i].getsize();
1132 }
1133
1136 {
1137 if (build_mode!=random)
1138 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1139 if (ready != building)
1140 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1141
1142 r[i].setsize(r[i].getsize()+s);
1143 }
1144
1147 {
1148 if (build_mode!=random)
1149 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1150 if (ready != building)
1151 DUNE_THROW(BCRSMatrixError,"matrix row sizes already built up");
1152
1153 // compute total size, check positivity
1154 size_type total=0;
1155 for (size_type i=0; i<n; i++)
1156 {
1157 total += r[i].getsize();
1158 }
1159
1160 if(nnz_ == 0)
1161 // allocate/check memory
1162 allocate(n,m,total,false,false);
1163 else if(nnz_ < total)
1164 DUNE_THROW(BCRSMatrixError,"Specified number of nonzeros ("<<nnz_<<") not "
1165 <<"sufficient for calculated nonzeros ("<<total<<"! ");
1166
1167 // set the window pointers correctly
1169
1170 // initialize j_ array with m (an invalid column index)
1171 // this indicates an unused entry
1172 for (size_type k=0; k<nnz_; k++)
1173 j_.get()[k] = m;
1174 ready = rowSizesBuilt;
1175 }
1176
1178
1189 {
1190 if (build_mode!=random)
1191 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1192 if (ready==built)
1193 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1194 if (ready==building)
1195 DUNE_THROW(BCRSMatrixError,"matrix row sizes not built up yet");
1196 if (ready==notAllocated)
1197 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1198
1199 if (col >= m)
1200 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1201
1202 // get row range
1203 size_type* const first = r[row].getindexptr();
1204 size_type* const last = first + r[row].getsize();
1205
1206 // find correct insertion position for new column index
1207 size_type* pos = std::lower_bound(first,last,col);
1208
1209 // check if index is already in row
1210 if (pos!=last && *pos == col) return;
1211
1212 // find end of already inserted column indices
1213 size_type* end = std::lower_bound(pos,last,m);
1214 if (end==last)
1215 DUNE_THROW(BCRSMatrixError,"row is too small");
1216
1217 // insert new column index at correct position
1218 std::copy_backward(pos,end,end+1);
1219 *pos = col;
1220 }
1221
1223
1231 template<typename It>
1233 {
1234 size_type row_size = r[row].size();
1235 size_type* col_begin = r[row].getindexptr();
1236 size_type* col_end;
1237 // consistency check between allocated row size and number of passed column indices
1238 if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
1239 DUNE_THROW(BCRSMatrixError,"Given size of row " << row
1240 << " (" << row_size
1241 << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
1242 }
1243
1244
1246
1254 template<typename It>
1255 void setIndices(size_type row, It begin, It end)
1256 {
1257 size_type row_size = r[row].size();
1258 size_type* col_begin = r[row].getindexptr();
1259 size_type* col_end;
1260 // consistency check between allocated row size and number of passed column indices
1261 if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
1262 DUNE_THROW(BCRSMatrixError,"Given size of row " << row
1263 << " (" << row_size
1264 << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
1265 std::sort(col_begin,col_end);
1266 }
1267
1270 {
1271 if (build_mode!=random)
1272 DUNE_THROW(BCRSMatrixError,"requires random build mode");
1273 if (ready==built)
1274 DUNE_THROW(BCRSMatrixError,"matrix already built up");
1275 if (ready==building)
1276 DUNE_THROW(BCRSMatrixError,"row sizes are not built up yet");
1277 if (ready==notAllocated)
1278 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1279
1280 // check if there are undefined indices
1281 RowIterator endi=end();
1282 for (RowIterator i=begin(); i!=endi; ++i)
1283 {
1284 ColIterator endj = (*i).end();
1285 for (ColIterator j=(*i).begin(); j!=endj; ++j) {
1286 if (j.index() >= m) {
1287 dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
1288 <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
1289 nnz_ -= ((*i).end().offset() - j.offset());
1290 r[i.index()].setsize(j.offset());
1291 break;
1292 }
1293 }
1294 }
1295
1296 allocateData();
1298
1299 // if not, set matrix to built
1300 ready = built;
1301 }
1302
1303 //===== implicit creation interface
1304
1306
1318 {
1319#ifdef DUNE_ISTL_WITH_CHECKING
1320 if (build_mode!=implicit)
1321 DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1322 if (ready==built)
1323 DUNE_THROW(BCRSMatrixError,"matrix already built up, use operator[] for entry access now");
1324 if (ready==notAllocated)
1325 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1326 if (ready!=building)
1327 DUNE_THROW(InvalidStateException,"You may only use entry() during the 'building' stage");
1328
1329 if (row >= n)
1330 DUNE_THROW(BCRSMatrixError,"row index exceeds matrix size");
1331 if (col >= m)
1332 DUNE_THROW(BCRSMatrixError,"column index exceeds matrix size");
1333#endif
1334
1335 size_type* begin = r[row].getindexptr();
1336 size_type* end = begin + r[row].getsize();
1337
1338 size_type* pos = std::find(begin, end, col);
1339
1340 //treat the case that there was a match in the array
1341 if (pos != end)
1342 if (*pos == col)
1343 {
1344 std::ptrdiff_t offset = pos - r[row].getindexptr();
1345 B* aptr = r[row].getptr() + offset;
1346
1347 return *aptr;
1348 }
1349
1350 //determine whether overflow has to be taken into account or not
1351 if (r[row].getsize() == avg)
1352 return overflow[std::make_pair(row,col)];
1353 else
1354 {
1355 //modify index array
1356 *end = col;
1357
1358 //do simultaneous operations on data array a
1359 std::ptrdiff_t offset = end - r[row].getindexptr();
1360 B* apos = r[row].getptr() + offset;
1361
1362 //increase rowsize
1363 r[row].setsize(r[row].getsize()+1);
1364
1365 //return reference to the newly created entry
1366 return *apos;
1367 }
1368 }
1369
1371
1382 {
1383 if (build_mode!=implicit)
1384 DUNE_THROW(BCRSMatrixError,"requires implicit build mode");
1385 if (ready==built)
1386 DUNE_THROW(BCRSMatrixError,"matrix already built up, no more need for compression");
1387 if (ready==notAllocated)
1388 DUNE_THROW(BCRSMatrixError,"matrix size not set and no memory allocated yet");
1389 if (ready!=building)
1390 DUNE_THROW(InvalidStateException,"You may only call compress() at the end of the 'building' stage");
1391
1392 //calculate statistics
1394 stats.overflow_total = overflow.size();
1395 stats.maximum = 0;
1396
1397 //get insertion iterators pointing to one before start (for later use of ++it)
1398 size_type* jiit = j_.get();
1399 B* aiit = a;
1400
1401 //get iterator to the smallest overflow element
1402 typename OverflowType::iterator oit = overflow.begin();
1403
1404 //store a copy of index pointers on which to perform sorting
1405 std::vector<size_type*> perm;
1406
1407 //iterate over all rows and copy elements into their position in the compressed array
1408 for (size_type i=0; i<n; i++)
1409 {
1410 //get old pointers into a and j and size without overflow changes
1411 size_type* begin = r[i].getindexptr();
1412 //B* apos = r[i].getptr();
1413 size_type size = r[i].getsize();
1414
1415 perm.resize(size);
1416
1417 typename std::vector<size_type*>::iterator it = perm.begin();
1418 for (size_type* iit = begin; iit < begin + size; ++iit, ++it)
1419 *it = iit;
1420
1421 //sort permutation array
1422 std::sort(perm.begin(),perm.end(),PointerCompare<size_type>());
1423
1424 //change row window pointer to their new positions
1425 r[i].setindexptr(jiit);
1426 r[i].setptr(aiit);
1427
1428 for (it = perm.begin(); it != perm.end(); ++it)
1429 {
1430 //check whether there are elements in the overflow area which take precedence
1431 while ((oit!=overflow.end()) && (oit->first < std::make_pair(i,**it)))
1432 {
1433 //check whether there is enough memory to write to
1434 if (jiit > begin)
1436 "Allocated memory for BCRSMatrix exhausted during compress()!"
1437 "Please increase either the average number of entries per row or the compressionBufferSize value."
1438 );
1439 //copy an element from the overflow area to the insertion position in a and j
1440 *jiit = oit->first.second;
1441 ++jiit;
1442 *aiit = oit->second;
1443 ++aiit;
1444 ++oit;
1445 r[i].setsize(r[i].getsize()+1);
1446 }
1447
1448 //check whether there is enough memory to write to
1449 if (jiit > begin)
1451 "Allocated memory for BCRSMatrix exhausted during compress()!"
1452 "Please increase either the average number of entries per row or the compressionBufferSize value."
1453 );
1454
1455 //copy element from array
1456 *jiit = **it;
1457 ++jiit;
1458 B* apos = *it - j_.get() + a;
1459 *aiit = *apos;
1460 ++aiit;
1461 }
1462
1463 //copy remaining elements from the overflow area
1464 while ((oit!=overflow.end()) && (oit->first.first == i))
1465 {
1466 //check whether there is enough memory to write to
1467 if (jiit > begin)
1469 "Allocated memory for BCRSMatrix exhausted during compress()!"
1470 "Please increase either the average number of entries per row or the compressionBufferSize value."
1471 );
1472
1473 //copy and element from the overflow area to the insertion position in a and j
1474 *jiit = oit->first.second;
1475 ++jiit;
1476 *aiit = oit->second;
1477 ++aiit;
1478 ++oit;
1479 r[i].setsize(r[i].getsize()+1);
1480 }
1481
1482 // update maximum row size
1483 if (r[i].getsize()>stats.maximum)
1484 stats.maximum = r[i].getsize();
1485 }
1486
1487 // overflow area may be cleared
1488 overflow.clear();
1489
1490 //determine average number of entries and memory usage
1491 if ( n == 0)
1492 {
1493 stats.avg = 0;
1494 stats.mem_ratio = 1;
1495 }
1496 else
1497 {
1498 std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
1499 nnz_ = diff;
1500 stats.avg = (double) (nnz_) / (double) n;
1501 stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
1502 }
1503
1504 //matrix is now built
1505 ready = built;
1506
1507 return stats;
1508 }
1509
1510 //===== vector space arithmetic
1511
1514 {
1515#ifdef DUNE_ISTL_WITH_CHECKING
1516 if (ready != built)
1517 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1518#endif
1519
1520 if (nnz_ > 0)
1521 {
1522 // process 1D array
1523 for (size_type i=0; i<nnz_; i++)
1524 a[i] *= k;
1525 }
1526 else
1527 {
1528 RowIterator endi=end();
1529 for (RowIterator i=begin(); i!=endi; ++i)
1530 {
1531 ColIterator endj = (*i).end();
1532 for (ColIterator j=(*i).begin(); j!=endj; ++j)
1533 (*j) *= k;
1534 }
1535 }
1536
1537 return *this;
1538 }
1539
1542 {
1543#ifdef DUNE_ISTL_WITH_CHECKING
1544 if (ready != built)
1545 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1546#endif
1547
1548 if (nnz_ > 0)
1549 {
1550 // process 1D array
1551 for (size_type i=0; i<nnz_; i++)
1552 a[i] /= k;
1553 }
1554 else
1555 {
1556 RowIterator endi=end();
1557 for (RowIterator i=begin(); i!=endi; ++i)
1558 {
1559 ColIterator endj = (*i).end();
1560 for (ColIterator j=(*i).begin(); j!=endj; ++j)
1561 (*j) /= k;
1562 }
1563 }
1564
1565 return *this;
1566 }
1567
1568
1575 {
1576#ifdef DUNE_ISTL_WITH_CHECKING
1577 if (ready != built || b.ready != built)
1578 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1579 if(N()!=b.N() || M() != b.M())
1580 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1581#endif
1582 RowIterator endi=end();
1583 ConstRowIterator j=b.begin();
1584 for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1585 i->operator+=(*j);
1586 }
1587
1588 return *this;
1589 }
1590
1597 {
1598#ifdef DUNE_ISTL_WITH_CHECKING
1599 if (ready != built || b.ready != built)
1600 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1601 if(N()!=b.N() || M() != b.M())
1602 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1603#endif
1604 RowIterator endi=end();
1605 ConstRowIterator j=b.begin();
1606 for (RowIterator i=begin(); i!=endi; ++i, ++j) {
1607 i->operator-=(*j);
1608 }
1609
1610 return *this;
1611 }
1612
1622 {
1623#ifdef DUNE_ISTL_WITH_CHECKING
1624 if (ready != built || b.ready != built)
1625 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1626 if(N()!=b.N() || M() != b.M())
1627 DUNE_THROW(RangeError, "Matrix sizes do not match!");
1628#endif
1629 RowIterator endi=end();
1630 ConstRowIterator j=b.begin();
1631 for(RowIterator i=begin(); i!=endi; ++i, ++j)
1632 i->axpy(alpha, *j);
1633
1634 return *this;
1635 }
1636
1637 //===== linear maps
1638
1640 template<class X, class Y>
1641 void mv (const X& x, Y& y) const
1642 {
1643#ifdef DUNE_ISTL_WITH_CHECKING
1644 if (ready != built)
1645 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1646 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,
1647 "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
1648 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,
1649 "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
1650#endif
1651 ConstRowIterator endi=end();
1652 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1653 {
1654 y[i.index()]=0;
1655 ConstColIterator endj = (*i).end();
1656 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1657 {
1658 auto&& xj = Impl::asVector(x[j.index()]);
1659 auto&& yi = Impl::asVector(y[i.index()]);
1660 Impl::asMatrix(*j).umv(xj, yi);
1661 }
1662 }
1663 }
1664
1666 template<class X, class Y>
1667 void umv (const X& x, Y& y) const
1668 {
1669#ifdef DUNE_ISTL_WITH_CHECKING
1670 if (ready != built)
1671 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1672 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1673 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1674#endif
1675 ConstRowIterator endi=end();
1676 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1677 {
1678 ConstColIterator endj = (*i).end();
1679 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1680 {
1681 auto&& xj = Impl::asVector(x[j.index()]);
1682 auto&& yi = Impl::asVector(y[i.index()]);
1683 Impl::asMatrix(*j).umv(xj,yi);
1684 }
1685 }
1686 }
1687
1689 template<class X, class Y>
1690 void mmv (const X& x, Y& y) const
1691 {
1692#ifdef DUNE_ISTL_WITH_CHECKING
1693 if (ready != built)
1694 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1695 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1696 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1697#endif
1698 ConstRowIterator endi=end();
1699 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1700 {
1701 ConstColIterator endj = (*i).end();
1702 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1703 {
1704 auto&& xj = Impl::asVector(x[j.index()]);
1705 auto&& yi = Impl::asVector(y[i.index()]);
1706 Impl::asMatrix(*j).mmv(xj,yi);
1707 }
1708 }
1709 }
1710
1712 template<class X, class Y, class F>
1713 void usmv (F&& alpha, const X& x, Y& y) const
1714 {
1715#ifdef DUNE_ISTL_WITH_CHECKING
1716 if (ready != built)
1717 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1718 if (x.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1719 if (y.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1720#endif
1721 ConstRowIterator endi=end();
1722 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1723 {
1724 ConstColIterator endj = (*i).end();
1725 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1726 {
1727 auto&& xj = Impl::asVector(x[j.index()]);
1728 auto&& yi = Impl::asVector(y[i.index()]);
1729 Impl::asMatrix(*j).usmv(alpha,xj,yi);
1730 }
1731 }
1732 }
1733
1735 template<class X, class Y>
1736 void mtv (const X& x, Y& y) const
1737 {
1738#ifdef DUNE_ISTL_WITH_CHECKING
1739 if (ready != built)
1740 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1741 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1742 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1743#endif
1744 for(size_type i=0; i<y.N(); ++i)
1745 y[i]=0;
1746 umtv(x,y);
1747 }
1748
1750 template<class X, class Y>
1751 void umtv (const X& x, Y& y) const
1752 {
1753#ifdef DUNE_ISTL_WITH_CHECKING
1754 if (ready != built)
1755 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1756 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1757 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1758#endif
1759 ConstRowIterator endi=end();
1760 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1761 {
1762 ConstColIterator endj = (*i).end();
1763 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1764 {
1765 auto&& xi = Impl::asVector(x[i.index()]);
1766 auto&& yj = Impl::asVector(y[j.index()]);
1767 Impl::asMatrix(*j).umtv(xi,yj);
1768 }
1769 }
1770 }
1771
1773 template<class X, class Y>
1774 void mmtv (const X& x, Y& y) const
1775 {
1776#ifdef DUNE_ISTL_WITH_CHECKING
1777 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1778 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1779#endif
1780 ConstRowIterator endi=end();
1781 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1782 {
1783 ConstColIterator endj = (*i).end();
1784 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1785 {
1786 auto&& xi = Impl::asVector(x[i.index()]);
1787 auto&& yj = Impl::asVector(y[j.index()]);
1788 Impl::asMatrix(*j).mmtv(xi,yj);
1789 }
1790 }
1791 }
1792
1794 template<class X, class Y>
1795 void usmtv (const field_type& alpha, const X& x, Y& y) const
1796 {
1797#ifdef DUNE_ISTL_WITH_CHECKING
1798 if (ready != built)
1799 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1800 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1801 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1802#endif
1803 ConstRowIterator endi=end();
1804 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1805 {
1806 ConstColIterator endj = (*i).end();
1807 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1808 {
1809 auto&& xi = Impl::asVector(x[i.index()]);
1810 auto&& yj = Impl::asVector(y[j.index()]);
1811 Impl::asMatrix(*j).usmtv(alpha,xi,yj);
1812 }
1813 }
1814 }
1815
1817 template<class X, class Y>
1818 void umhv (const X& x, Y& y) const
1819 {
1820#ifdef DUNE_ISTL_WITH_CHECKING
1821 if (ready != built)
1822 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1823 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1824 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1825#endif
1826 ConstRowIterator endi=end();
1827 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1828 {
1829 ConstColIterator endj = (*i).end();
1830 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1831 {
1832 auto&& xi = Impl::asVector(x[i.index()]);
1833 auto&& yj = Impl::asVector(y[j.index()]);
1834 Impl::asMatrix(*j).umhv(xi,yj);
1835 }
1836 }
1837 }
1838
1840 template<class X, class Y>
1841 void mmhv (const X& x, Y& y) const
1842 {
1843#ifdef DUNE_ISTL_WITH_CHECKING
1844 if (ready != built)
1845 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1846 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1847 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1848#endif
1849 ConstRowIterator endi=end();
1850 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1851 {
1852 ConstColIterator endj = (*i).end();
1853 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1854 {
1855 auto&& xi = Impl::asVector(x[i.index()]);
1856 auto&& yj = Impl::asVector(y[j.index()]);
1857 Impl::asMatrix(*j).mmhv(xi,yj);
1858 }
1859 }
1860 }
1861
1863 template<class X, class Y>
1864 void usmhv (const field_type& alpha, const X& x, Y& y) const
1865 {
1866#ifdef DUNE_ISTL_WITH_CHECKING
1867 if (ready != built)
1868 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1869 if (x.N()!=N()) DUNE_THROW(BCRSMatrixError,"index out of range");
1870 if (y.N()!=M()) DUNE_THROW(BCRSMatrixError,"index out of range");
1871#endif
1872 ConstRowIterator endi=end();
1873 for (ConstRowIterator i=this->begin(); i!=endi; ++i)
1874 {
1875 ConstColIterator endj = (*i).end();
1876 for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1877 {
1878 auto&& xi = Impl::asVector(x[i.index()]);
1879 auto&& yj = Impl::asVector(y[j.index()]);
1880 Impl::asMatrix(*j).usmhv(alpha,xi,yj);
1881 }
1882 }
1883 }
1884
1885
1886 //===== norms
1887
1889 typename FieldTraits<field_type>::real_type frobenius_norm2 () const
1890 {
1891#ifdef DUNE_ISTL_WITH_CHECKING
1892 if (ready != built)
1893 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1894#endif
1895
1896 typename FieldTraits<field_type>::real_type sum=0;
1897
1898 for (auto&& row : (*this))
1899 for (auto&& entry : row)
1900 sum += Impl::asMatrix(entry).frobenius_norm2();
1901
1902 return sum;
1903 }
1904
1906 typename FieldTraits<field_type>::real_type frobenius_norm () const
1907 {
1908 return sqrt(frobenius_norm2());
1909 }
1910
1912 template <typename ft = field_type,
1913 typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
1914 typename FieldTraits<ft>::real_type infinity_norm() const {
1915 if (ready != built)
1916 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1917
1918 using real_type = typename FieldTraits<ft>::real_type;
1919 using std::max;
1920
1921 real_type norm = 0;
1922 for (auto const &x : *this) {
1923 real_type sum = 0;
1924 for (auto const &y : x)
1925 sum += Impl::asMatrix(y).infinity_norm();
1926 norm = max(sum, norm);
1927 }
1928 return norm;
1929 }
1930
1932 template <typename ft = field_type,
1933 typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
1934 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1935 if (ready != built)
1936 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1937
1938 using real_type = typename FieldTraits<ft>::real_type;
1939 using std::max;
1940
1941 real_type norm = 0;
1942 for (auto const &x : *this) {
1943 real_type sum = 0;
1944 for (auto const &y : x)
1945 sum += Impl::asMatrix(y).infinity_norm_real();
1946 norm = max(sum, norm);
1947 }
1948 return norm;
1949 }
1950
1952 template <typename ft = field_type,
1953 typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
1954 typename FieldTraits<ft>::real_type infinity_norm() const {
1955 if (ready != built)
1956 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1957
1958 using real_type = typename FieldTraits<ft>::real_type;
1959 using std::max;
1960
1961 real_type norm = 0;
1962 real_type isNaN = 1;
1963 for (auto const &x : *this) {
1964 real_type sum = 0;
1965 for (auto const &y : x)
1966 sum += Impl::asMatrix(y).infinity_norm();
1967 norm = max(sum, norm);
1968 isNaN += sum;
1969 }
1970
1971 return norm * (isNaN / isNaN);
1972 }
1973
1975 template <typename ft = field_type,
1976 typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
1977 typename FieldTraits<ft>::real_type infinity_norm_real() const {
1978 if (ready != built)
1979 DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
1980
1981 using real_type = typename FieldTraits<ft>::real_type;
1982 using std::max;
1983
1984 real_type norm = 0;
1985 real_type isNaN = 1;
1986
1987 for (auto const &x : *this) {
1988 real_type sum = 0;
1989 for (auto const &y : x)
1990 sum += Impl::asMatrix(y).infinity_norm_real();
1991 norm = max(sum, norm);
1992 isNaN += sum;
1993 }
1994
1995 return norm * (isNaN / isNaN);
1996 }
1997
1998 //===== sizes
1999
2001 size_type N () const
2002 {
2003 return n;
2004 }
2005
2007 size_type M () const
2008 {
2009 return m;
2010 }
2011
2014 {
2015 // in case of row-wise allocation
2016 if( nnz_ <= 0 )
2017 nnz_ = std::accumulate( begin(), end(), size_type( 0 ), [] ( size_type s, const row_type &row ) { return s+row.getsize(); } );
2018 return nnz_;
2019 }
2020
2023 {
2024 return ready;
2025 }
2026
2029 {
2030 return build_mode;
2031 }
2032
2033 //===== query
2034
2036 bool exists (size_type i, size_type j) const
2037 {
2038#ifdef DUNE_ISTL_WITH_CHECKING
2039 if (i<0 || i>=n) DUNE_THROW(BCRSMatrixError,"row index out of range");
2040 if (j<0 || j>=m) DUNE_THROW(BCRSMatrixError,"column index out of range");
2041#endif
2042 return (r[i].size() && r[i].find(j) != r[i].end());
2043 }
2044
2045
2046 protected:
2047 // state information
2048 BuildMode build_mode; // row wise or whole matrix
2049 BuildStage ready; // indicate the stage the matrix building is in
2050
2051 // The allocator used for memory management
2052 typename std::allocator_traits<A>::template rebind_alloc<B> allocator_;
2053
2054 typename std::allocator_traits<A>::template rebind_alloc<row_type> rowAllocator_;
2055
2056 typename std::allocator_traits<A>::template rebind_alloc<size_type> sizeAllocator_;
2057
2058 // size of the matrix
2059 size_type n; // number of rows
2060 size_type m; // number of columns
2061 mutable size_type nnz_; // number of nonzeroes contained in the matrix
2062 size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
2063 // zero means that memory is allocated separately for each row.
2064
2065 // the rows are dynamically allocated
2066 row_type* r; // [n] the individual rows having pointers into a,j arrays
2067
2068 // dynamically allocated memory
2069 B* a; // [allocationSize] non-zero entries of the matrix in row-wise ordering
2070 // If a single array of column indices is used, it can be shared
2071 // between different matrices with the same sparsity pattern
2072 std::shared_ptr<size_type> j_; // [allocationSize] column indices of entries
2073
2074 // additional data is needed in implicit buildmode
2075 size_type avg;
2076 double compressionBufferSize_;
2077
2078 typedef std::map<std::pair<size_type,size_type>, B> OverflowType;
2079 OverflowType overflow;
2080
2081 void setWindowPointers(ConstRowIterator row)
2082 {
2083 row_type current_row(a,j_.get(),0); // Pointers to current row data
2084 for (size_type i=0; i<n; i++, ++row) {
2085 // set row i
2086 size_type s = row->getsize();
2087
2088 if (s>0) {
2089 // setup pointers and size
2090 r[i].set(s,current_row.getptr(), current_row.getindexptr());
2091 // update pointer for next row
2092 current_row.setptr(current_row.getptr()+s);
2093 current_row.setindexptr(current_row.getindexptr()+s);
2094 } else{
2095 // empty row
2096 r[i].set(0,nullptr,nullptr);
2097 }
2098 }
2099 }
2100
2102
2107 {
2108 size_type* jptr = j_.get();
2109 for (size_type i=0; i<n; ++i, ++row) {
2110 // set row i
2111 size_type s = row->getsize();
2112
2113 if (s>0) {
2114 // setup pointers and size
2115 r[i].setsize(s);
2116 r[i].setindexptr(jptr);
2117 } else{
2118 // empty row
2119 r[i].set(0,nullptr,nullptr);
2120 }
2121
2122 // advance position in global array
2123 jptr += s;
2124 }
2125 }
2126
2128
2133 {
2134 B* aptr = a;
2135 for (size_type i=0; i<n; ++i) {
2136 // set row i
2137 if (r[i].getsize() > 0) {
2138 // setup pointers and size
2139 r[i].setptr(aptr);
2140 } else{
2141 // empty row
2142 r[i].set(0,nullptr,nullptr);
2143 }
2144
2145 // advance position in global array
2146 aptr += r[i].getsize();
2147 }
2148 }
2149
2152 {
2153 setWindowPointers(Mat.begin());
2154
2155 // copy data
2156 for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
2157
2158 // finish off
2159 build_mode = row_wise; // dummy
2160 ready = built;
2161 }
2162
2168 void deallocate(bool deallocateRows=true)
2169 {
2170
2171 if (notAllocated)
2172 return;
2173
2174 if (allocationSize_>0)
2175 {
2176 // a,j_ have been allocated as one long vector
2177 j_.reset();
2178 if (a)
2179 {
2180 for(B *aiter=a+(allocationSize_-1), *aend=a-1; aiter!=aend; --aiter)
2181 std::allocator_traits<decltype(allocator_)>::destroy(allocator_, aiter);
2182 allocator_.deallocate(a,allocationSize_);
2183 a = nullptr;
2184 }
2185 }
2186 else if (r)
2187 {
2188 // check if memory for rows have been allocated individually
2189 for (size_type i=0; i<n; i++)
2190 if (r[i].getsize()>0)
2191 {
2192 for (B *col=r[i].getptr()+(r[i].getsize()-1),
2193 *colend = r[i].getptr()-1; col!=colend; --col) {
2194 std::allocator_traits<decltype(allocator_)>::destroy(allocator_, col);
2195 }
2196 sizeAllocator_.deallocate(r[i].getindexptr(),1);
2197 allocator_.deallocate(r[i].getptr(),1);
2198 // clear out row data in case we don't want to deallocate the rows
2199 // otherwise we might run into a double free problem here later
2200 r[i].set(0,nullptr,nullptr);
2201 }
2202 }
2203
2204 // deallocate the rows
2205 if (n>0 && deallocateRows && r) {
2206 for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
2207 std::allocator_traits<decltype(rowAllocator_)>::destroy(rowAllocator_, riter);
2208 rowAllocator_.deallocate(r,n);
2209 r = nullptr;
2210 }
2211
2212 // Mark matrix as not built at all.
2213 ready=notAllocated;
2214
2215 }
2216
2235 void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
2236 {
2237 // Store size
2238 n = rows;
2239 m = columns;
2240 nnz_ = allocationSize;
2241 allocationSize_ = allocationSize;
2242
2243 // allocate rows
2244 if(allocateRows) {
2245 if (n>0) {
2246 if (r)
2247 DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
2248 r = rowAllocator_.allocate(rows);
2249 // initialize row entries
2250 for(row_type* ri=r; ri!=r+rows; ++ri)
2251 std::allocator_traits<decltype(rowAllocator_)>::construct(rowAllocator_, ri, row_type());
2252 }else{
2253 r = 0;
2254 }
2255 }
2256
2257 // allocate a and j_ array
2258 if (allocate_data)
2259 allocateData();
2260 // allocate column indices only if not yet present (enable sharing)
2261 if (allocationSize_>0) {
2262 // we copy allocator and size to the deleter since _j may outlive this class
2263 if (!j_.get())
2264 j_.reset(sizeAllocator_.allocate(allocationSize_),
2265 [alloc = sizeAllocator_, size = allocationSize_](auto ptr) mutable {
2266 alloc.deallocate(ptr, size);
2267 });
2268 }else{
2269 j_.reset();
2270 }
2271
2272 // Mark the matrix as not built.
2273 ready = building;
2274 }
2275
2276 void allocateData()
2277 {
2278 if (a)
2279 DUNE_THROW(InvalidStateException,"Cannot allocate data array (already allocated)");
2280 if (allocationSize_>0) {
2281 a = allocator_.allocate(allocationSize_);
2282 // use placement new to call constructor that allocates
2283 // additional memory.
2284 new (a) B[allocationSize_];
2285 } else {
2286 a = nullptr;
2287 }
2288 }
2289
2296 {
2297 if (build_mode != implicit)
2298 DUNE_THROW(InvalidStateException,"implicit_allocate() may only be called in implicit build mode");
2299 if (ready != notAllocated)
2300 DUNE_THROW(InvalidStateException,"memory has already been allocated");
2301
2302 // check to make sure the user has actually set the parameters
2303 if (compressionBufferSize_ < 0)
2304 DUNE_THROW(InvalidStateException,"You have to set the implicit build mode parameters before starting to build the matrix");
2305 //calculate size of overflow region, add buffer for row sort!
2306 size_type osize = (size_type) (_n*avg)*compressionBufferSize_ + 4*avg;
2307 allocationSize_ = _n*avg + osize;
2308
2309 allocate(_n, _m, allocationSize_,true,true);
2310
2311 //set row pointers correctly
2312 size_type* jptr = j_.get() + osize;
2313 B* aptr = a + osize;
2314 for (size_type i=0; i<n; i++)
2315 {
2316 r[i].set(0,aptr,jptr);
2317 jptr = jptr + avg;
2318 aptr = aptr + avg;
2319 }
2320
2321 ready = building;
2322 }
2323 };
2324
2325
2326 template<class B, class A>
2327 struct FieldTraits< BCRSMatrix<B, A> >
2328 {
2329 using field_type = typename BCRSMatrix<B, A>::field_type;
2330 using real_type = typename FieldTraits<field_type>::real_type;
2331 };
2332
2335} // end namespace
2336
2337#endif
Helper functions for determining the vector/matrix block level.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Error specific to BCRSMatrix.
Definition: istlexception.hh:24
Iterator class for sequential creation of blocks
Definition: bcrsmatrix.hh:954
bool operator==(const CreateIterator &it) const
equality
Definition: bcrsmatrix.hh:1049
CreateIterator & operator++()
prefix increment
Definition: bcrsmatrix.hh:974
size_type index() const
The number of the row that the iterator currently points to.
Definition: bcrsmatrix.hh:1055
bool operator!=(const CreateIterator &it) const
inequality
Definition: bcrsmatrix.hh:1043
CreateIterator(BCRSMatrix &_Mat, size_type _i)
constructor
Definition: bcrsmatrix.hh:957
void insert(size_type j)
put column index in row
Definition: bcrsmatrix.hh:1061
bool contains(size_type j)
return true if column index is in row
Definition: bcrsmatrix.hh:1067
size_type size() const
Get the current row size.
Definition: bcrsmatrix.hh:1076
Iterator access to matrix rows
Definition: bcrsmatrix.hh:575
RealRowIterator()
empty constructor, use with care!
Definition: bcrsmatrix.hh:592
bool equals(const RealRowIterator< ValueType > &other) const
equality
Definition: bcrsmatrix.hh:620
std::remove_const< T >::type ValueType
The unqualified value type.
Definition: bcrsmatrix.hh:579
bool equals(const RealRowIterator< const ValueType > &other) const
equality
Definition: bcrsmatrix.hh:627
RealRowIterator(row_type *_p, size_type _i)
constructor
Definition: bcrsmatrix.hh:587
size_type index() const
return index
Definition: bcrsmatrix.hh:602
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: bcrsmatrix.hh:488
bool exists(size_type i, size_type j) const
return true if (i,j) is in pattern
Definition: bcrsmatrix.hh:2036
BuildStage buildStage() const
The current build stage of the matrix.
Definition: bcrsmatrix.hh:2022
Iterator begin()
Get iterator to first row.
Definition: bcrsmatrix.hh:671
friend class CreateIterator
allow CreateIterator to access internal data
Definition: bcrsmatrix.hh:1091
void copyWindowStructure(const BCRSMatrix &Mat)
Copy the window structure from another matrix.
Definition: bcrsmatrix.hh:2151
B & entry(size_type row, size_type col)
Returns reference to entry (row,col) of the matrix.
Definition: bcrsmatrix.hh:1317
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: bcrsmatrix.hh:1841
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: bcrsmatrix.hh:1864
void umtv(const X &x, Y &y) const
y += A^T x
Definition: bcrsmatrix.hh:1751
RealRowIterator< const row_type > const_iterator
The const iterator over the matrix rows.
Definition: bcrsmatrix.hh:703
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:2235
BCRSMatrix & axpy(field_type alpha, const BCRSMatrix &b)
Add the scaled entries of another matrix to this one.
Definition: bcrsmatrix.hh:1621
FieldTraits< ft >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: bcrsmatrix.hh:1934
~BCRSMatrix()
destructor
Definition: bcrsmatrix.hh:821
void deallocate(bool deallocateRows=true)
deallocate memory of the matrix.
Definition: bcrsmatrix.hh:2168
Iterator RowIterator
rename the iterators for easier access
Definition: bcrsmatrix.hh:697
row_type & operator[](size_type i)
random access to the rows
Definition: bcrsmatrix.hh:545
BCRSMatrix()
an empty matrix
Definition: bcrsmatrix.hh:745
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1146
void incrementrowsize(size_type i, size_type s=1)
increment size of row i by s (1 by default)
Definition: bcrsmatrix.hh:1135
void mtv(const X &x, Y &y) const
y = A^T x
Definition: bcrsmatrix.hh:1736
void umhv(const X &x, Y &y) const
y += A^H x
Definition: bcrsmatrix.hh:1818
size_type nonzeroes() const
number of blocks that are stored (the number of blocks that possibly are nonzero)
Definition: bcrsmatrix.hh:2013
ConstIterator ConstRowIterator
rename the const row iterator for easier access
Definition: bcrsmatrix.hh:734
void setIndicesNoSort(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition: bcrsmatrix.hh:1232
void setrowsize(size_type i, size_type s)
Set number of indices in row i to s.
Definition: bcrsmatrix.hh:1114
BCRSMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: bcrsmatrix.hh:1513
RealRowIterator< row_type > iterator
The iterator over the (mutable matrix rows.
Definition: bcrsmatrix.hh:667
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: bcrsmatrix.hh:1795
ConstIterator beforeBegin() const
Definition: bcrsmatrix.hh:728
Iterator beforeBegin()
Definition: bcrsmatrix.hh:691
BuildMode
we support two modes
Definition: bcrsmatrix.hh:506
@ implicit
Build entries randomly with an educated guess for the number of entries per row.
Definition: bcrsmatrix.hh:535
@ unknown
Build mode not set!
Definition: bcrsmatrix.hh:539
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:526
@ row_wise
Build in a row-wise manner.
Definition: bcrsmatrix.hh:517
BCRSMatrix(size_type _n, size_type _m, size_type _nnz, BuildMode bm)
matrix with known number of nonzeroes
Definition: bcrsmatrix.hh:752
Imp::CompressedBlockVectorWindow< B, size_type > row_type
implement row_type with compressed vector
Definition: bcrsmatrix.hh:500
::Dune::CompressionStatistics< size_type > CompressionStatistics
The type for the statistics object returned by compress()
Definition: bcrsmatrix.hh:503
BCRSMatrix & operator-=(const BCRSMatrix &b)
Subtract the entries of another matrix from this one.
Definition: bcrsmatrix.hh:1596
BCRSMatrix(const BCRSMatrix &Mat)
copy constructor
Definition: bcrsmatrix.hh:802
Iterator end()
Get iterator to one beyond last row.
Definition: bcrsmatrix.hh:677
void setIndices(size_type row, It begin, It end)
Set all column indices for row from the given iterator range.
Definition: bcrsmatrix.hh:1255
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1188
row_type::Iterator ColIterator
Iterator for the entries of each row.
Definition: bcrsmatrix.hh:700
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: bcrsmatrix.hh:1906
A::size_type size_type
The type for the index access and the size.
Definition: bcrsmatrix.hh:497
BCRSMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: bcrsmatrix.hh:1541
BCRSMatrix & operator+=(const BCRSMatrix &b)
Add the entries of another matrix to this one.
Definition: bcrsmatrix.hh:1574
BCRSMatrix(size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
construct matrix with a known average number of entries per row
Definition: bcrsmatrix.hh:781
CreateIterator createend()
get create iterator pointing to one after the last block
Definition: bcrsmatrix.hh:1100
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: bcrsmatrix.hh:1889
Iterator beforeEnd()
Definition: bcrsmatrix.hh:684
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:737
void usmv(F &&alpha, const X &x, Y &y) const
y += alpha A x
Definition: bcrsmatrix.hh:1713
size_type getrowsize(size_type i) const
get current number of indices in row i
Definition: bcrsmatrix.hh:1125
size_type M() const
number of columns (counted in blocks)
Definition: bcrsmatrix.hh:2007
CreateIterator createbegin()
get initial create iterator
Definition: bcrsmatrix.hh:1094
BuildStage
Definition: bcrsmatrix.hh:469
@ rowSizesBuilt
The row sizes of the matrix are known.
Definition: bcrsmatrix.hh:480
@ built
The matrix structure is fully built.
Definition: bcrsmatrix.hh:482
@ notbuilt
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:471
@ notAllocated
Matrix is not built at all, no memory has been allocated, build mode and size can still be set.
Definition: bcrsmatrix.hh:473
@ building
Matrix is currently being built, some memory has been allocated, build mode and size are fixed.
Definition: bcrsmatrix.hh:475
BuildMode buildMode() const
The currently selected build mode of the matrix.
Definition: bcrsmatrix.hh:2028
void mmv(const X &x, Y &y) const
y -= A x
Definition: bcrsmatrix.hh:1690
FieldTraits< ft >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: bcrsmatrix.hh:1914
void mv(const X &x, Y &y) const
y = A x
Definition: bcrsmatrix.hh:1641
B block_type
export the type representing the components
Definition: bcrsmatrix.hh:491
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: bcrsmatrix.hh:1774
void umv(const X &x, Y &y) const
y += A x
Definition: bcrsmatrix.hh:1667
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:2295
void setImplicitBuildModeParameters(size_type _avg, double compressionBufferSize)
Set parameters needed for creation in implicit build mode.
Definition: bcrsmatrix.hh:886
BCRSMatrix(size_type _n, size_type _m, BuildMode bm)
matrix with unknown number of nonzeroes
Definition: bcrsmatrix.hh:761
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1269
CompressionStatistics compress()
Finishes the buildstage in implicit mode.
Definition: bcrsmatrix.hh:1381
void setDataPointers()
Set data pointers for all rows.
Definition: bcrsmatrix.hh:2132
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:2001
void setBuildMode(BuildMode bm)
Sets the build mode of the matrix.
Definition: bcrsmatrix.hh:830
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:858
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:908
void setColumnPointers(ConstRowIterator row)
Copy row sizes from iterator range starting at row and set column index pointers for all rows.
Definition: bcrsmatrix.hh:2106
ConstIterator end() const
Get const iterator to one beyond last row.
Definition: bcrsmatrix.hh:714
ConstIterator begin() const
Get const iterator to first row.
Definition: bcrsmatrix.hh:708
A allocator_type
export the allocator type
Definition: bcrsmatrix.hh:494
ConstIterator beforeEnd() const
Definition: bcrsmatrix.hh:721
Proxy row object for entry access.
Definition: bcrsmatrix.hh:137
block_type & operator[](size_type j) const
Returns entry in column j.
Definition: bcrsmatrix.hh:142
A wrapper for uniform access to the BCRSMatrix during and after the build stage in implicit build mod...
Definition: bcrsmatrix.hh:117
Matrix::block_type block_type
The block_type of the underlying matrix.
Definition: bcrsmatrix.hh:125
ImplicitMatrixBuilder(Matrix &m)
Creates an ImplicitMatrixBuilder for matrix m.
Definition: bcrsmatrix.hh:170
M_ Matrix
The underlying matrix.
Definition: bcrsmatrix.hh:122
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:194
size_type M() const
The number of columns in the matrix.
Definition: bcrsmatrix.hh:217
Matrix::size_type size_type
The size_type of the underlying matrix.
Definition: bcrsmatrix.hh:128
row_object operator[](size_type i) const
Returns a proxy for entries in row i.
Definition: bcrsmatrix.hh:205
size_type N() const
The number of rows in the matrix.
Definition: bcrsmatrix.hh:211
Thrown when the compression buffer used by the implicit BCRSMatrix construction is exhausted.
Definition: istlexception.hh:37
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
A generic dynamic dense matrix.
Definition: matrix.hh:561
A::size_type size_type
Type for indices and sizes.
Definition: matrix.hh:577
T block_type
Export the type representing the components.
Definition: matrix.hh:568
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:435
Default exception class for range errors.
Definition: exceptions.hh:254
Traits for type conversions and type information.
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr decltype(auto) elementAt(Container &&c, Index &&i)
Get element at given position from container.
Definition: hybridutilities.hh:126
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
DWarnType dwarn(std::cerr)
Stream for warnings indicating problems.
Definition: stdstreams.hh:162
This file implements iterator facade classes for writing stl conformant iterators.
Some handy generic functions for ISTL matrices.
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is NaN.
Definition: fvector.hh:627
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
Implements a scalar matrix view wrapper around an existing scalar.
Implements a scalar vector view wrapper around an existing scalar.
Standard Dune debug streams.
Statistics about compression achieved in implicit mode.
Definition: bcrsmatrix.hh:88
size_type overflow_total
total number of elements written to the overflow area during construction.
Definition: bcrsmatrix.hh:94
size_type maximum
maximum number of non-zeroes per row.
Definition: bcrsmatrix.hh:92
double avg
average number of non-zeroes per row.
Definition: bcrsmatrix.hh:90
double mem_ratio
fraction of wasted memory resulting from non-used overflow area.
Definition: bcrsmatrix.hh:99
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)