4 #ifndef DUNE_BCRSMATRIX_HH
5 #define DUNE_BCRSMATRIX_HH
17 #include <dune/common/shared_ptr.hh>
18 #include <dune/common/stdstreams.hh>
19 #include <dune/common/iteratorfacades.hh>
20 #include <dune/common/typetraits.hh>
21 #include <dune/common/static_assert.hh>
69 struct MatrixDimension;
178 template<
class B,
class A=std::allocator<B> >
255 #ifdef DUNE_ISTL_WITH_CHECKING
256 if (r==0) DUNE_THROW(
ISTLError,
"row not initialized yet");
257 if (i>=n) DUNE_THROW(
ISTLError,
"index out of range");
258 if (r[i].getptr()==0) DUNE_THROW(
ISTLError,
"row not initialized yet");
266 #ifdef DUNE_ISTL_WITH_CHECKING
267 if (built!=ready) DUNE_THROW(
ISTLError,
"row not initialized yet");
268 if (i>=n) DUNE_THROW(
ISTLError,
"index out of range");
279 :
public RandomAccessIteratorFacade<RealRowIterator<T>, T>
351 void advance(std::ptrdiff_t diff)
356 T& elementAt(std::ptrdiff_t diff)
const
448 : build_mode(
unknown), ready(notbuilt), n(0), m(0), nnz(0),
454 : build_mode(bm), ready(notbuilt)
456 allocate(_n, _m, _nnz);
461 : build_mode(bm), ready(notbuilt)
486 allocate(Mat.n, Mat.m, _nnz);
489 copyWindowStructure(Mat);
507 DUNE_THROW(InvalidStateException,
"Matrix structure is already built (ready="<<ready<<
").");
531 allocate(rows, columns, nnz);
543 if (&Mat==
this)
return *
this;
549 if (n>0 && n!=Mat.n) {
554 rowAllocator_.deallocate(r,n);
566 allocate(Mat.n, Mat.m, nnz, n!=Mat.n);
569 copyWindowStructure(Mat);
588 : Mat(_Mat), i(_i), nnz(0), current_row(Mat.a, Mat.j.
get(), 0)
590 if (i==0 && Mat.ready)
591 DUNE_THROW(
ISTLError,
"creation only allowed for uninitialized matrix");
597 DUNE_THROW(
ISTLError,
"creation only allowed if row wise allocation was requested in the constructor");
606 DUNE_THROW(
ISTLError,
"matrix already built up");
626 DUNE_THROW(
ISTLError,
"allocated nnz too small");
635 B* a = Mat.allocator_.allocate(s);
637 size_type* j = Mat.sizeAllocator_.allocate(s);
648 for (
typename PatternType::const_iterator it=pattern.begin(); it!=pattern.end(); ++it)
671 return (i!=it.i) || (&Mat!=&it.Mat);
677 return (i==it.i) && (&Mat==&it.Mat);
695 if (pattern.find(j)!=pattern.end())
707 return pattern.size();
714 typedef std::set<size_type,std::less<size_type> > PatternType;
741 DUNE_THROW(
ISTLError,
"requires random build mode");
743 DUNE_THROW(
ISTLError,
"matrix row sizes already built up");
751 #ifdef DUNE_ISTL_WITH_CHECKING
752 if (r==0) DUNE_THROW(
ISTLError,
"row not initialized yet");
753 if (i>=n) DUNE_THROW(
ISTLError,
"index out of range");
762 DUNE_THROW(
ISTLError,
"requires random build mode");
764 DUNE_THROW(
ISTLError,
"matrix row sizes already built up");
766 r[i].
setsize(r[i].getsize()+s);
773 DUNE_THROW(
ISTLError,
"requires random build mode");
775 DUNE_THROW(
ISTLError,
"matrix row sizes already built up");
781 if (r[i].getsize()<0)
782 DUNE_THROW(
ISTLError,
"rowsize must be nonnegative");
788 allocate(n,m,total,
false);
790 DUNE_THROW(
ISTLError,
"Specified number of nonzeros ("<<nnz<<
") not "
791 <<
"sufficient for calculated nonzeros ("<<total<<
"! ");
794 setWindowPointers(
begin());
800 ready = rowSizesBuilt;
817 DUNE_THROW(
ISTLError,
"requires random build mode");
819 DUNE_THROW(
ISTLError,
"matrix already built up");
821 DUNE_THROW(
ISTLError,
"matrix row sizes not built up yet");
824 DUNE_THROW(
ISTLError,
"column index exceeds matrix size");
831 size_type* pos = std::lower_bound(first,last,col);
834 if (pos!=last && *pos == col)
return;
839 DUNE_THROW(
ISTLError,
"row is too small");
842 std::copy_backward(pos,end,end+1);
851 DUNE_THROW(
ISTLError,
"requires random build mode");
853 DUNE_THROW(
ISTLError,
"matrix already built up");
855 DUNE_THROW(
ISTLError,
"row sizes are not built up yet");
865 std::cout <<
"j[" << j.offset() <<
"]=" << j.index() << std::endl;
866 DUNE_THROW(
ISTLError,
"undefined index detected");
869 dwarn <<
"WARNING: size of row "<< i.index()<<
" is "<<j.offset()<<
". But was specified as being "<< (*i).end().offset()
870 <<
". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
871 r[i.index()].
setsize(j.offset());
937 #ifdef DUNE_ISTL_WITH_CHECKING
938 if(
N()!=b.
N() ||
M() != b.
M())
939 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
957 #ifdef DUNE_ISTL_WITH_CHECKING
958 if(
N()!=b.
N() ||
M() != b.
M())
959 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
980 #ifdef DUNE_ISTL_WITH_CHECKING
981 if(
N()!=b.
N() ||
M() != b.
M())
982 DUNE_THROW(RangeError,
"Matrix sizes do not match!");
995 template<
class X,
class Y>
996 void mv (
const X& x, Y& y)
const
998 #ifdef DUNE_ISTL_WITH_CHECKING
1000 "Size mismatch: M: " <<
N() <<
"x" <<
M() <<
" x: " << x.N());
1002 "Size mismatch: M: " <<
N() <<
"x" <<
M() <<
" y: " << y.N());
1010 (*j).umv(x[j.index()],y[i.index()]);
1015 template<
class X,
class Y>
1016 void umv (
const X& x, Y& y)
const
1018 #ifdef DUNE_ISTL_WITH_CHECKING
1019 if (x.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1020 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1027 (*j).umv(x[j.index()],y[i.index()]);
1032 template<
class X,
class Y>
1033 void mmv (
const X& x, Y& y)
const
1035 #ifdef DUNE_ISTL_WITH_CHECKING
1036 if (x.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1037 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1044 (*j).mmv(x[j.index()],y[i.index()]);
1049 template<
class X,
class Y>
1052 #ifdef DUNE_ISTL_WITH_CHECKING
1053 if (x.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1054 if (y.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1061 (*j).usmv(alpha,x[j.index()],y[i.index()]);
1066 template<
class X,
class Y>
1067 void mtv (
const X& x, Y& y)
const
1069 #ifdef DUNE_ISTL_WITH_CHECKING
1070 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1071 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1079 template<
class X,
class Y>
1082 #ifdef DUNE_ISTL_WITH_CHECKING
1083 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1084 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1091 (*j).umtv(x[i.index()],y[j.index()]);
1096 template<
class X,
class Y>
1099 #ifdef DUNE_ISTL_WITH_CHECKING
1100 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1101 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1108 (*j).mmtv(x[i.index()],y[j.index()]);
1113 template<
class X,
class Y>
1116 #ifdef DUNE_ISTL_WITH_CHECKING
1117 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1118 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1125 (*j).usmtv(alpha,x[i.index()],y[j.index()]);
1130 template<
class X,
class Y>
1133 #ifdef DUNE_ISTL_WITH_CHECKING
1134 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1135 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1142 (*j).umhv(x[i.index()],y[j.index()]);
1147 template<
class X,
class Y>
1150 #ifdef DUNE_ISTL_WITH_CHECKING
1151 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1152 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1159 (*j).mmhv(x[i.index()],y[j.index()]);
1164 template<
class X,
class Y>
1167 #ifdef DUNE_ISTL_WITH_CHECKING
1168 if (x.N()!=
N()) DUNE_THROW(
ISTLError,
"index out of range");
1169 if (y.N()!=
M()) DUNE_THROW(
ISTLError,
"index out of range");
1176 (*j).usmhv(alpha,x[i.index()],y[j.index()]);
1193 sum += (*j).frobenius_norm2();
1215 sum += (*j).infinity_norm();
1216 max = std::max(max,sum);
1231 sum += (*j).infinity_norm_real();
1232 max = std::max(max,sum);
1264 #ifdef DUNE_ISTL_WITH_CHECKING
1265 if (i<0 || i>=n) DUNE_THROW(
ISTLError,
"index out of range");
1266 if (j<0 || i>=m) DUNE_THROW(
ISTLError,
"index out of range");
1268 if (r[i].size() && r[i].find(j)!=r[i].
end())
1281 typename A::template rebind<B>::other allocator_;
1283 typename A::template rebind<row_type>::other rowAllocator_;
1285 typename A::template rebind<size_type>::other sizeAllocator_;
1300 Dune::shared_ptr<size_type> j;
1312 r[i].
set(s,current_row.getptr(), current_row.getindexptr());
1314 current_row.setptr(current_row.getptr()+s);
1315 current_row.setindexptr(current_row.getindexptr()+s);
1324 void copyWindowStructure(
const BCRSMatrix& Mat)
1326 setWindowPointers(Mat.begin());
1329 for (
size_type i=0; i<n; i++) r[i] = Mat.r[i];
1341 void deallocate(
bool deallocateRows=
true)
1351 allocator_.deallocate(a,n);
1357 if (r[i].getsize()>0)
1364 sizeAllocator_.deallocate(r[i].getindexptr(),1);
1365 allocator_.deallocate(r[i].getptr(),1);
1370 if (n>0 && deallocateRows) {
1374 rowAllocator_.deallocate(r,n);
1385 typename A::template rebind<size_type>::other& sizeAllocator_;
1388 Deallocator(
typename A::template rebind<size_type>::other& sizeAllocator)
1389 : sizeAllocator_(sizeAllocator)
1392 void operator()(
size_type* p) { sizeAllocator_.deallocate(p,1); }
1423 r = rowAllocator_.allocate(rows);
1433 a = allocator_.allocate(nnz);
1436 j.reset(sizeAllocator_.allocate(nnz),Deallocator(sizeAllocator_));