dune-istl
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=4 sw=2 sts=2: 00003 00004 #ifndef DUNE_BCRSMATRIX_HH 00005 #define DUNE_BCRSMATRIX_HH 00006 00007 #include<cmath> 00008 #include<complex> 00009 #include<set> 00010 #include<iostream> 00011 #include<algorithm> 00012 #include<numeric> 00013 #include<vector> 00014 00015 #include "istlexception.hh" 00016 #include "bvector.hh" 00017 #include <dune/common/shared_ptr.hh> 00018 #include <dune/common/stdstreams.hh> 00019 #include <dune/common/iteratorfacades.hh> 00020 #include <dune/common/typetraits.hh> 00021 #include <dune/common/static_assert.hh> 00022 00027 namespace Dune { 00028 00068 template<typename M> 00069 struct MatrixDimension; 00070 00178 template<class B, class A=std::allocator<B> > 00179 class BCRSMatrix 00180 { 00181 friend struct MatrixDimension<BCRSMatrix>; 00182 00183 private: 00184 enum BuildStage{ 00186 notbuilt=0, 00191 rowSizesBuilt=1, 00193 built=2 00194 }; 00195 00196 public: 00197 00198 //===== type definitions and constants 00199 00201 typedef typename B::field_type field_type; 00202 00204 typedef B block_type; 00205 00207 typedef A allocator_type; 00208 00210 typedef CompressedBlockVectorWindow<B,A> row_type; 00211 00213 typedef typename A::size_type size_type; 00214 00216 enum { 00218 blocklevel = B::blocklevel+1 00219 }; 00220 00222 enum BuildMode { 00233 row_wise, 00242 random, 00246 unknown 00247 }; 00248 00249 00250 //===== random access interface to rows of the matrix 00251 00253 row_type& operator[] (size_type i) 00254 { 00255 #ifdef DUNE_ISTL_WITH_CHECKING 00256 if (r==0) DUNE_THROW(ISTLError,"row not initialized yet"); 00257 if (i>=n) DUNE_THROW(ISTLError,"index out of range"); 00258 if (r[i].getptr()==0) DUNE_THROW(ISTLError,"row not initialized yet"); 00259 #endif 00260 return r[i]; 00261 } 00262 00264 const row_type& operator[] (size_type i) const 00265 { 00266 #ifdef DUNE_ISTL_WITH_CHECKING 00267 if (built!=ready) DUNE_THROW(ISTLError,"row not initialized yet"); 00268 if (i>=n) DUNE_THROW(ISTLError,"index out of range"); 00269 #endif 00270 return r[i]; 00271 } 00272 00273 00274 //===== iterator interface to rows of the matrix 00275 00277 template<class T> 00278 class RealRowIterator 00279 : public RandomAccessIteratorFacade<RealRowIterator<T>, T> 00280 { 00281 00282 public: 00284 typedef typename remove_const<T>::type ValueType; 00285 00286 friend class RandomAccessIteratorFacade<RealRowIterator<const ValueType>, const ValueType>; 00287 friend class RandomAccessIteratorFacade<RealRowIterator<ValueType>, ValueType>; 00288 friend class RealRowIterator<const ValueType>; 00289 friend class RealRowIterator<ValueType>; 00290 00292 RealRowIterator (row_type* _p, size_type _i) 00293 : p(_p), i(_i) 00294 {} 00295 00297 RealRowIterator () 00298 : p(0), i(0) 00299 {} 00300 00301 RealRowIterator(const RealRowIterator<ValueType>& it) 00302 : p(it.p), i(it.i) 00303 {} 00304 00305 00307 size_type index () const 00308 { 00309 return i; 00310 } 00311 00312 std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const 00313 { 00314 assert(other.p==p); 00315 return (other.i-i); 00316 } 00317 00318 std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const 00319 { 00320 assert(other.p==p); 00321 return (other.i-i); 00322 } 00323 00325 bool equals (const RealRowIterator<ValueType>& other) const 00326 { 00327 assert(other.p==p); 00328 return i==other.i; 00329 } 00330 00332 bool equals (const RealRowIterator<const ValueType>& other) const 00333 { 00334 assert(other.p==p); 00335 return i==other.i; 00336 } 00337 00338 private: 00340 void increment() 00341 { 00342 ++i; 00343 } 00344 00346 void decrement() 00347 { 00348 --i; 00349 } 00350 00351 void advance(std::ptrdiff_t diff) 00352 { 00353 i+=diff; 00354 } 00355 00356 T& elementAt(std::ptrdiff_t diff) const 00357 { 00358 return p[i+diff]; 00359 } 00360 00362 row_type& dereference () const 00363 { 00364 return p[i]; 00365 } 00366 00367 row_type* p; 00368 size_type i; 00369 }; 00370 00372 typedef RealRowIterator<row_type> iterator; 00373 typedef RealRowIterator<row_type> Iterator; 00374 00376 Iterator begin () 00377 { 00378 return Iterator(r,0); 00379 } 00380 00382 Iterator end () 00383 { 00384 return Iterator(r,n); 00385 } 00386 00389 Iterator beforeEnd () 00390 { 00391 return Iterator(r,n-1); 00392 } 00393 00396 Iterator beforeBegin () 00397 { 00398 return Iterator(r,-1); 00399 } 00400 00402 typedef Iterator RowIterator; 00403 00405 typedef typename row_type::Iterator ColIterator; 00406 00408 typedef RealRowIterator<const row_type> const_iterator; 00409 typedef RealRowIterator<const row_type> ConstIterator; 00410 00411 00413 ConstIterator begin () const 00414 { 00415 return ConstIterator(r,0); 00416 } 00417 00419 ConstIterator end () const 00420 { 00421 return ConstIterator(r,n); 00422 } 00423 00426 ConstIterator beforeEnd() const 00427 { 00428 return ConstIterator(r,n-1); 00429 } 00430 00433 ConstIterator beforeBegin () const 00434 { 00435 return ConstIterator(r,-1); 00436 } 00437 00439 typedef ConstIterator ConstRowIterator; 00440 00442 typedef typename row_type::ConstIterator ConstColIterator; 00443 00444 //===== constructors & resizers 00445 00447 BCRSMatrix () 00448 : build_mode(unknown), ready(notbuilt), n(0), m(0), nnz(0), 00449 r(0), a(0) 00450 {} 00451 00453 BCRSMatrix (size_type _n, size_type _m, size_type _nnz, BuildMode bm) 00454 : build_mode(bm), ready(notbuilt) 00455 { 00456 allocate(_n, _m, _nnz); 00457 } 00458 00460 BCRSMatrix (size_type _n, size_type _m, BuildMode bm) 00461 : build_mode(bm), ready(notbuilt) 00462 { 00463 allocate(_n, _m); 00464 } 00465 00471 BCRSMatrix (const BCRSMatrix& Mat) 00472 : n(Mat.n), nnz(0) 00473 { 00474 // deep copy in global array 00475 size_type _nnz = Mat.nnz; 00476 00477 // in case of row-wise allocation 00478 if (_nnz<=0) 00479 { 00480 _nnz = 0; 00481 for (size_type i=0; i<n; i++) 00482 _nnz += Mat.r[i].getsize(); 00483 } 00484 00485 j = Mat.j; // enable column index sharing, release array in case of row-wise allocation 00486 allocate(Mat.n, Mat.m, _nnz); 00487 00488 // build window structure 00489 copyWindowStructure(Mat); 00490 } 00491 00493 ~BCRSMatrix () 00494 { 00495 deallocate(); 00496 } 00497 00502 void setBuildMode(BuildMode bm) 00503 { 00504 if(ready==notbuilt) 00505 build_mode = bm; 00506 else 00507 DUNE_THROW(InvalidStateException, "Matrix structure is already built (ready="<<ready<<")."); 00508 } 00509 00525 void setSize(size_type rows, size_type columns, size_type nnz=0) 00526 { 00527 // deallocate already setup memory 00528 deallocate(); 00529 00530 // allocate matrix memory 00531 allocate(rows, columns, nnz); 00532 } 00533 00540 BCRSMatrix& operator= (const BCRSMatrix& Mat) 00541 { 00542 // return immediately when self-assignment 00543 if (&Mat==this) return *this; 00544 00545 // make it simple: ALWAYS throw away memory for a and j 00546 deallocate(false); 00547 00548 // reallocate the rows if required 00549 if (n>0 && n!=Mat.n) { 00550 // free rows 00551 int i=n; 00552 while (i) 00553 r[--i].~row_type(); 00554 rowAllocator_.deallocate(r,n); 00555 } 00556 00557 nnz=Mat.nnz; 00558 if (nnz<=0) 00559 { 00560 for (size_type i=0; i<Mat.n; i++) 00561 nnz += Mat.r[i].getsize(); 00562 } 00563 00564 // allocate a, share j 00565 j = Mat.j; 00566 allocate(Mat.n, Mat.m, nnz, n!=Mat.n); 00567 00568 // build window structure 00569 copyWindowStructure(Mat); 00570 return *this; 00571 } 00572 00574 BCRSMatrix& operator= (const field_type& k) 00575 { 00576 for (size_type i=0; i<n; i++) r[i] = k; 00577 return *this; 00578 } 00579 00580 //===== row-wise creation interface 00581 00583 class CreateIterator 00584 { 00585 public: 00587 CreateIterator (BCRSMatrix& _Mat, size_type _i) 00588 : Mat(_Mat), i(_i), nnz(0), current_row(Mat.a, Mat.j.get(), 0) 00589 { 00590 if (i==0 && Mat.ready) 00591 DUNE_THROW(ISTLError,"creation only allowed for uninitialized matrix"); 00592 if(Mat.build_mode!=row_wise) 00593 { 00594 if(Mat.build_mode==unknown) 00595 Mat.build_mode=row_wise; 00596 else 00597 DUNE_THROW(ISTLError,"creation only allowed if row wise allocation was requested in the constructor"); 00598 } 00599 } 00600 00602 CreateIterator& operator++() 00603 { 00604 // this should only be called if matrix is in creation 00605 if (Mat.ready) 00606 DUNE_THROW(ISTLError,"matrix already built up"); 00607 00608 // row i is defined through the pattern 00609 // get memory for the row and initialize the j array 00610 // this depends on the allocation mode 00611 00612 // compute size of the row 00613 size_type s = pattern.size(); 00614 00615 if(s>0){ 00616 // update number of nonzeroes including this row 00617 nnz += s; 00618 00619 // alloc memory / set window 00620 if (Mat.nnz>0) 00621 { 00622 // memory is allocated in one long array 00623 00624 // check if that memory is sufficient 00625 if (nnz>Mat.nnz) 00626 DUNE_THROW(ISTLError,"allocated nnz too small"); 00627 00628 // set row i 00629 Mat.r[i].set(s,current_row.getptr(),current_row.getindexptr()); 00630 current_row.setptr(current_row.getptr()+s); 00631 current_row.setindexptr(current_row.getindexptr()+s); 00632 }else{ 00633 // memory is allocated individually per row 00634 // allocate and set row i 00635 B* a = Mat.allocator_.allocate(s); 00636 new (a) B[s]; 00637 size_type* j = Mat.sizeAllocator_.allocate(s); 00638 new (j) size_type[s]; 00639 Mat.r[i].set(s,a,j); 00640 } 00641 }else 00642 // setup empty row 00643 Mat.r[i].set(0,0,0); 00644 00645 // initialize the j array for row i from pattern 00646 size_type k=0; 00647 size_type *j = Mat.r[i].getindexptr(); 00648 for (typename PatternType::const_iterator it=pattern.begin(); it!=pattern.end(); ++it) 00649 j[k++] = *it; 00650 00651 // now go to next row 00652 i++; 00653 pattern.clear(); 00654 00655 // check if this was last row 00656 if (i==Mat.n) 00657 { 00658 Mat.ready = built; 00659 if(Mat.nnz>0) 00660 // Set nnz to the exact number of nonzero blocks inserted 00661 // as some methods rely on it 00662 Mat.nnz=nnz; 00663 } 00664 // done 00665 return *this; 00666 } 00667 00669 bool operator!= (const CreateIterator& it) const 00670 { 00671 return (i!=it.i) || (&Mat!=&it.Mat); 00672 } 00673 00675 bool operator== (const CreateIterator& it) const 00676 { 00677 return (i==it.i) && (&Mat==&it.Mat); 00678 } 00679 00681 size_type index () const 00682 { 00683 return i; 00684 } 00685 00687 void insert (size_type j) 00688 { 00689 pattern.insert(j); 00690 } 00691 00693 bool contains (size_type j) 00694 { 00695 if (pattern.find(j)!=pattern.end()) 00696 return true; 00697 else 00698 return false; 00699 } 00705 size_type size() const 00706 { 00707 return pattern.size(); 00708 } 00709 00710 private: 00711 BCRSMatrix& Mat; // the matrix we are defining 00712 size_type i; // current row to be defined 00713 size_type nnz; // count total number of nonzeros 00714 typedef std::set<size_type,std::less<size_type> > PatternType; 00715 PatternType pattern; // used to compile entries in a row 00716 row_type current_row; // row pointing to the current row to setup 00717 }; 00718 00720 friend class CreateIterator; 00721 00723 CreateIterator createbegin () 00724 { 00725 return CreateIterator(*this,0); 00726 } 00727 00729 CreateIterator createend () 00730 { 00731 return CreateIterator(*this,n); 00732 } 00733 00734 00735 //===== random creation interface 00736 00738 void setrowsize (size_type i, size_type s) 00739 { 00740 if (build_mode!=random) 00741 DUNE_THROW(ISTLError,"requires random build mode"); 00742 if (ready) 00743 DUNE_THROW(ISTLError,"matrix row sizes already built up"); 00744 00745 r[i].setsize(s); 00746 } 00747 00749 size_type getrowsize (size_type i) const 00750 { 00751 #ifdef DUNE_ISTL_WITH_CHECKING 00752 if (r==0) DUNE_THROW(ISTLError,"row not initialized yet"); 00753 if (i>=n) DUNE_THROW(ISTLError,"index out of range"); 00754 #endif 00755 return r[i].getsize(); 00756 } 00757 00759 void incrementrowsize (size_type i, size_type s = 1) 00760 { 00761 if (build_mode!=random) 00762 DUNE_THROW(ISTLError,"requires random build mode"); 00763 if (ready) 00764 DUNE_THROW(ISTLError,"matrix row sizes already built up"); 00765 00766 r[i].setsize(r[i].getsize()+s); 00767 } 00768 00770 void endrowsizes () 00771 { 00772 if (build_mode!=random) 00773 DUNE_THROW(ISTLError,"requires random build mode"); 00774 if (ready) 00775 DUNE_THROW(ISTLError,"matrix row sizes already built up"); 00776 00777 // compute total size, check positivity 00778 size_type total=0; 00779 for (size_type i=0; i<n; i++) 00780 { 00781 if (r[i].getsize()<0) 00782 DUNE_THROW(ISTLError,"rowsize must be nonnegative"); 00783 total += r[i].getsize(); 00784 } 00785 00786 if(nnz==0) 00787 // allocate/check memory 00788 allocate(n,m,total,false); 00789 else if(nnz<total) 00790 DUNE_THROW(ISTLError,"Specified number of nonzeros ("<<nnz<<") not " 00791 <<"sufficient for calculated nonzeros ("<<total<<"! "); 00792 00793 // set the window pointers correctly 00794 setWindowPointers(begin()); 00795 00796 // initialize j array with m (an invalid column index) 00797 // this indicates an unused entry 00798 for (size_type k=0; k<nnz; k++) 00799 j.get()[k] = m; 00800 ready = rowSizesBuilt; 00801 } 00802 00804 00814 void addindex (size_type row, size_type col) 00815 { 00816 if (build_mode!=random) 00817 DUNE_THROW(ISTLError,"requires random build mode"); 00818 if (ready==built) 00819 DUNE_THROW(ISTLError,"matrix already built up"); 00820 if (ready==notbuilt) 00821 DUNE_THROW(ISTLError,"matrix row sizes not built up yet"); 00822 00823 if (col >= m) 00824 DUNE_THROW(ISTLError,"column index exceeds matrix size"); 00825 00826 // get row range 00827 size_type* const first = r[row].getindexptr(); 00828 size_type* const last = first + r[row].getsize(); 00829 00830 // find correct insertion position for new column index 00831 size_type* pos = std::lower_bound(first,last,col); 00832 00833 // check if index is already in row 00834 if (pos!=last && *pos == col) return; 00835 00836 // find end of already inserted column indices 00837 size_type* end = std::lower_bound(pos,last,m); 00838 if (end==last) 00839 DUNE_THROW(ISTLError,"row is too small"); 00840 00841 // insert new column index at correct position 00842 std::copy_backward(pos,end,end+1); 00843 *pos = col; 00844 00845 } 00846 00848 void endindices () 00849 { 00850 if (build_mode!=random) 00851 DUNE_THROW(ISTLError,"requires random build mode"); 00852 if (ready==built) 00853 DUNE_THROW(ISTLError,"matrix already built up"); 00854 if (ready==notbuilt) 00855 DUNE_THROW(ISTLError,"row sizes are not built up yet"); 00856 00857 // check if there are undefined indices 00858 RowIterator endi=end(); 00859 for (RowIterator i=begin(); i!=endi; ++i) 00860 { 00861 ColIterator endj = (*i).end(); 00862 for (ColIterator j=(*i).begin(); j!=endj; ++j){ 00863 if (j.index()<0) 00864 { 00865 std::cout << "j[" << j.offset() << "]=" << j.index() << std::endl; 00866 DUNE_THROW(ISTLError,"undefined index detected"); 00867 } 00868 if (j.index()>=m){ 00869 dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset() 00870 <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl; 00871 r[i.index()].setsize(j.offset()); 00872 break; 00873 } 00874 } 00875 } 00876 00877 // if not, set matrix to built 00878 ready = built; 00879 } 00880 00881 //===== vector space arithmetic 00882 00884 BCRSMatrix& operator*= (const field_type& k) 00885 { 00886 if (nnz>0) 00887 { 00888 // process 1D array 00889 for (size_type i=0; i<nnz; i++) 00890 a[i] *= k; 00891 } 00892 else 00893 { 00894 RowIterator endi=end(); 00895 for (RowIterator i=begin(); i!=endi; ++i) 00896 { 00897 ColIterator endj = (*i).end(); 00898 for (ColIterator j=(*i).begin(); j!=endj; ++j) 00899 (*j) *= k; 00900 } 00901 } 00902 00903 return *this; 00904 } 00905 00907 BCRSMatrix& operator/= (const field_type& k) 00908 { 00909 if (nnz>0) 00910 { 00911 // process 1D array 00912 for (size_type i=0; i<nnz; i++) 00913 a[i] /= k; 00914 } 00915 else 00916 { 00917 RowIterator endi=end(); 00918 for (RowIterator i=begin(); i!=endi; ++i) 00919 { 00920 ColIterator endj = (*i).end(); 00921 for (ColIterator j=(*i).begin(); j!=endj; ++j) 00922 (*j) /= k; 00923 } 00924 } 00925 00926 return *this; 00927 } 00928 00929 00935 BCRSMatrix& operator+= (const BCRSMatrix& b) 00936 { 00937 #ifdef DUNE_ISTL_WITH_CHECKING 00938 if(N()!=b.N() || M() != b.M()) 00939 DUNE_THROW(RangeError, "Matrix sizes do not match!"); 00940 #endif 00941 RowIterator endi=end(); 00942 ConstRowIterator j=b.begin(); 00943 for (RowIterator i=begin(); i!=endi; ++i, ++j){ 00944 i->operator+=(*j); 00945 } 00946 00947 return *this; 00948 } 00949 00955 BCRSMatrix& operator-= (const BCRSMatrix& b) 00956 { 00957 #ifdef DUNE_ISTL_WITH_CHECKING 00958 if(N()!=b.N() || M() != b.M()) 00959 DUNE_THROW(RangeError, "Matrix sizes do not match!"); 00960 #endif 00961 RowIterator endi=end(); 00962 ConstRowIterator j=b.begin(); 00963 for (RowIterator i=begin(); i!=endi; ++i, ++j){ 00964 i->operator-=(*j); 00965 } 00966 00967 return *this; 00968 } 00969 00978 BCRSMatrix& axpy(field_type alpha, const BCRSMatrix& b) 00979 { 00980 #ifdef DUNE_ISTL_WITH_CHECKING 00981 if(N()!=b.N() || M() != b.M()) 00982 DUNE_THROW(RangeError, "Matrix sizes do not match!"); 00983 #endif 00984 RowIterator endi=end(); 00985 ConstRowIterator j=b.begin(); 00986 for(RowIterator i=begin(); i!=endi; ++i, ++j) 00987 i->axpy(alpha, *j); 00988 00989 return *this; 00990 } 00991 00992 //===== linear maps 00993 00995 template<class X, class Y> 00996 void mv (const X& x, Y& y) const 00997 { 00998 #ifdef DUNE_ISTL_WITH_CHECKING 00999 if (x.N()!=M()) DUNE_THROW(ISTLError, 01000 "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N()); 01001 if (y.N()!=N()) DUNE_THROW(ISTLError, 01002 "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N()); 01003 #endif 01004 ConstRowIterator endi=end(); 01005 for (ConstRowIterator i=begin(); i!=endi; ++i) 01006 { 01007 y[i.index()]=0; 01008 ConstColIterator endj = (*i).end(); 01009 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01010 (*j).umv(x[j.index()],y[i.index()]); 01011 } 01012 } 01013 01015 template<class X, class Y> 01016 void umv (const X& x, Y& y) const 01017 { 01018 #ifdef DUNE_ISTL_WITH_CHECKING 01019 if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01020 if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01021 #endif 01022 ConstRowIterator endi=end(); 01023 for (ConstRowIterator i=begin(); i!=endi; ++i) 01024 { 01025 ConstColIterator endj = (*i).end(); 01026 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01027 (*j).umv(x[j.index()],y[i.index()]); 01028 } 01029 } 01030 01032 template<class X, class Y> 01033 void mmv (const X& x, Y& y) const 01034 { 01035 #ifdef DUNE_ISTL_WITH_CHECKING 01036 if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01037 if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01038 #endif 01039 ConstRowIterator endi=end(); 01040 for (ConstRowIterator i=begin(); i!=endi; ++i) 01041 { 01042 ConstColIterator endj = (*i).end(); 01043 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01044 (*j).mmv(x[j.index()],y[i.index()]); 01045 } 01046 } 01047 01049 template<class X, class Y> 01050 void usmv (const field_type& alpha, const X& x, Y& y) const 01051 { 01052 #ifdef DUNE_ISTL_WITH_CHECKING 01053 if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01054 if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01055 #endif 01056 ConstRowIterator endi=end(); 01057 for (ConstRowIterator i=begin(); i!=endi; ++i) 01058 { 01059 ConstColIterator endj = (*i).end(); 01060 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01061 (*j).usmv(alpha,x[j.index()],y[i.index()]); 01062 } 01063 } 01064 01066 template<class X, class Y> 01067 void mtv (const X& x, Y& y) const 01068 { 01069 #ifdef DUNE_ISTL_WITH_CHECKING 01070 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01071 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01072 #endif 01073 for(size_type i=0; i<y.N(); ++i) 01074 y[i]=0; 01075 umtv(x,y); 01076 } 01077 01079 template<class X, class Y> 01080 void umtv (const X& x, Y& y) const 01081 { 01082 #ifdef DUNE_ISTL_WITH_CHECKING 01083 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01084 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01085 #endif 01086 ConstRowIterator endi=end(); 01087 for (ConstRowIterator i=begin(); i!=endi; ++i) 01088 { 01089 ConstColIterator endj = (*i).end(); 01090 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01091 (*j).umtv(x[i.index()],y[j.index()]); 01092 } 01093 } 01094 01096 template<class X, class Y> 01097 void mmtv (const X& x, Y& y) const 01098 { 01099 #ifdef DUNE_ISTL_WITH_CHECKING 01100 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01101 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01102 #endif 01103 ConstRowIterator endi=end(); 01104 for (ConstRowIterator i=begin(); i!=endi; ++i) 01105 { 01106 ConstColIterator endj = (*i).end(); 01107 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01108 (*j).mmtv(x[i.index()],y[j.index()]); 01109 } 01110 } 01111 01113 template<class X, class Y> 01114 void usmtv (const field_type& alpha, const X& x, Y& y) const 01115 { 01116 #ifdef DUNE_ISTL_WITH_CHECKING 01117 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01118 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01119 #endif 01120 ConstRowIterator endi=end(); 01121 for (ConstRowIterator i=begin(); i!=endi; ++i) 01122 { 01123 ConstColIterator endj = (*i).end(); 01124 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01125 (*j).usmtv(alpha,x[i.index()],y[j.index()]); 01126 } 01127 } 01128 01130 template<class X, class Y> 01131 void umhv (const X& x, Y& y) const 01132 { 01133 #ifdef DUNE_ISTL_WITH_CHECKING 01134 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01135 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01136 #endif 01137 ConstRowIterator endi=end(); 01138 for (ConstRowIterator i=begin(); i!=endi; ++i) 01139 { 01140 ConstColIterator endj = (*i).end(); 01141 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01142 (*j).umhv(x[i.index()],y[j.index()]); 01143 } 01144 } 01145 01147 template<class X, class Y> 01148 void mmhv (const X& x, Y& y) const 01149 { 01150 #ifdef DUNE_ISTL_WITH_CHECKING 01151 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01152 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01153 #endif 01154 ConstRowIterator endi=end(); 01155 for (ConstRowIterator i=begin(); i!=endi; ++i) 01156 { 01157 ConstColIterator endj = (*i).end(); 01158 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01159 (*j).mmhv(x[i.index()],y[j.index()]); 01160 } 01161 } 01162 01164 template<class X, class Y> 01165 void usmhv (const field_type& alpha, const X& x, Y& y) const 01166 { 01167 #ifdef DUNE_ISTL_WITH_CHECKING 01168 if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range"); 01169 if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range"); 01170 #endif 01171 ConstRowIterator endi=end(); 01172 for (ConstRowIterator i=begin(); i!=endi; ++i) 01173 { 01174 ConstColIterator endj = (*i).end(); 01175 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01176 (*j).usmhv(alpha,x[i.index()],y[j.index()]); 01177 } 01178 } 01179 01180 01181 //===== norms 01182 01184 double frobenius_norm2 () const 01185 { 01186 double sum=0; 01187 01188 ConstRowIterator endi=end(); 01189 for (ConstRowIterator i=begin(); i!=endi; ++i) 01190 { 01191 ConstColIterator endj = (*i).end(); 01192 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01193 sum += (*j).frobenius_norm2(); 01194 } 01195 01196 return sum; 01197 } 01198 01200 double frobenius_norm () const 01201 { 01202 return sqrt(frobenius_norm2()); 01203 } 01204 01206 double infinity_norm () const 01207 { 01208 double max=0; 01209 ConstRowIterator endi=end(); 01210 for (ConstRowIterator i=begin(); i!=endi; ++i) 01211 { 01212 double sum=0; 01213 ConstColIterator endj = (*i).end(); 01214 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01215 sum += (*j).infinity_norm(); 01216 max = std::max(max,sum); 01217 } 01218 return max; 01219 } 01220 01222 double infinity_norm_real () const 01223 { 01224 double max=0; 01225 ConstRowIterator endi=end(); 01226 for (ConstRowIterator i=begin(); i!=endi; ++i) 01227 { 01228 double sum=0; 01229 ConstColIterator endj = (*i).end(); 01230 for (ConstColIterator j=(*i).begin(); j!=endj; ++j) 01231 sum += (*j).infinity_norm_real(); 01232 max = std::max(max,sum); 01233 } 01234 return max; 01235 } 01236 01237 01238 //===== sizes 01239 01241 size_type N () const 01242 { 01243 return n; 01244 } 01245 01247 size_type M () const 01248 { 01249 return m; 01250 } 01251 01253 size_type nonzeroes () const 01254 { 01255 return nnz; 01256 } 01257 01258 01259 //===== query 01260 01262 bool exists (size_type i, size_type j) const 01263 { 01264 #ifdef DUNE_ISTL_WITH_CHECKING 01265 if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range"); 01266 if (j<0 || i>=m) DUNE_THROW(ISTLError,"index out of range"); 01267 #endif 01268 if (r[i].size() && r[i].find(j)!=r[i].end()) 01269 return true; 01270 else 01271 return false; 01272 } 01273 01274 01275 private: 01276 // state information 01277 BuildMode build_mode; // row wise or whole matrix 01278 BuildStage ready; // indicate the stage the matrix building is in 01279 01280 // The allocator used for memory management 01281 typename A::template rebind<B>::other allocator_; 01282 01283 typename A::template rebind<row_type>::other rowAllocator_; 01284 01285 typename A::template rebind<size_type>::other sizeAllocator_; 01286 01287 // size of the matrix 01288 size_type n; // number of rows 01289 size_type m; // number of columns 01290 size_type nnz; // number of nonzeros allocated in the a and j array below 01291 // zero means that memory is allocated separately for each row. 01292 01293 // the rows are dynamically allocated 01294 row_type* r; // [n] the individual rows having pointers into a,j arrays 01295 01296 // dynamically allocated memory 01297 B* a; // [nnz] non-zero entries of the matrix in row-wise ordering 01298 // If a single array of column indices is used, it can be shared 01299 // between different matrices with the same sparsity pattern 01300 Dune::shared_ptr<size_type> j; // [nnz] column indices of entries 01301 01302 01303 void setWindowPointers(ConstRowIterator row) 01304 { 01305 row_type current_row(a,j.get(),0); // Pointers to current row data 01306 for (size_type i=0; i<n; i++, ++row){ 01307 // set row i 01308 size_type s = row->getsize(); 01309 01310 if (s>0){ 01311 // setup pointers and size 01312 r[i].set(s,current_row.getptr(), current_row.getindexptr()); 01313 // update pointer for next row 01314 current_row.setptr(current_row.getptr()+s); 01315 current_row.setindexptr(current_row.getindexptr()+s); 01316 } else{ 01317 // empty row 01318 r[i].set(0,0,0); 01319 } 01320 } 01321 } 01322 01324 void copyWindowStructure(const BCRSMatrix& Mat) 01325 { 01326 setWindowPointers(Mat.begin()); 01327 01328 // copy data 01329 for (size_type i=0; i<n; i++) r[i] = Mat.r[i]; 01330 01331 // finish off 01332 build_mode = row_wise; // dummy 01333 ready = built; 01334 } 01335 01341 void deallocate(bool deallocateRows=true) 01342 { 01343 01344 if (nnz>0) 01345 { 01346 // a,j have been allocated as one long vector 01347 j.reset(); 01348 int i=nnz; 01349 while (i) 01350 a[--i].~B(); 01351 allocator_.deallocate(a,n); 01352 } 01353 else 01354 { 01355 // check if memory for rows have been allocated individually 01356 for (size_type i=0; i<n; i++) 01357 if (r[i].getsize()>0) 01358 { 01359 int j=r[i].getsize(); 01360 while (j) { 01361 r[i].getindexptr()[--j].~size_type(); 01362 r[i].getptr()[j].~B(); 01363 } 01364 sizeAllocator_.deallocate(r[i].getindexptr(),1); 01365 allocator_.deallocate(r[i].getptr(),1); 01366 } 01367 } 01368 01369 // deallocate the rows 01370 if (n>0 && deallocateRows) { 01371 int i=n; 01372 while (i) 01373 r[--i].~row_type(); 01374 rowAllocator_.deallocate(r,n); 01375 } 01376 01377 // Mark matrix as not built at all. 01378 ready=notbuilt; 01379 01380 } 01381 01383 class Deallocator 01384 { 01385 typename A::template rebind<size_type>::other& sizeAllocator_; 01386 01387 public: 01388 Deallocator(typename A::template rebind<size_type>::other& sizeAllocator) 01389 : sizeAllocator_(sizeAllocator) 01390 {} 01391 01392 void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); } 01393 }; 01394 01395 01413 void allocate(size_type rows, size_type columns, size_type nnz_=0, bool allocateRows=true) 01414 { 01415 // Store size 01416 n = rows; 01417 m = columns; 01418 nnz = nnz_; 01419 01420 // allocate rows 01421 if(allocateRows){ 01422 if (n>0){ 01423 r = rowAllocator_.allocate(rows); 01424 new (r) row_type[rows]; 01425 }else{ 01426 r = 0; 01427 } 01428 } 01429 01430 01431 // allocate a and j array 01432 if (nnz>0){ 01433 a = allocator_.allocate(nnz); 01434 // allocate column indices only if not yet present (enable sharing) 01435 if (!j.get()) 01436 j.reset(sizeAllocator_.allocate(nnz),Deallocator(sizeAllocator_)); 01437 }else{ 01438 a = 0; 01439 j.reset(); 01440 } 01441 // Mark the matrix as not built. 01442 ready = notbuilt; 01443 } 01444 01445 }; 01446 01447 01450 } // end namespace 01451 01452 #endif