dune-istl  2.2.0
bcrsmatrix.hh
Go to the documentation of this file.
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