dune-istl  2.2.1
bcrsmatrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_BCRSMATRIX_HH
5 #define DUNE_BCRSMATRIX_HH
6 
7 #include<cmath>
8 #include<complex>
9 #include<set>
10 #include<iostream>
11 #include<algorithm>
12 #include<numeric>
13 #include<vector>
14 
15 #include "istlexception.hh"
16 #include "bvector.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>
22 
27 namespace Dune {
28 
68  template<typename M>
70 
178  template<class B, class A=std::allocator<B> >
180  {
181  friend struct MatrixDimension<BCRSMatrix>;
182 
183  private:
184  enum BuildStage{
186  notbuilt=0,
191  rowSizesBuilt=1,
193  built=2
194  };
195 
196  public:
197 
198  //===== type definitions and constants
199 
201  typedef typename B::field_type field_type;
202 
204  typedef B block_type;
205 
207  typedef A allocator_type;
208 
211 
213  typedef typename A::size_type size_type;
214 
216  enum {
218  blocklevel = B::blocklevel+1
219  };
220 
222  enum BuildMode {
247  };
248 
249 
250  //===== random access interface to rows of the matrix
251 
254  {
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");
259 #endif
260  return r[i];
261  }
262 
265  {
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");
269 #endif
270  return r[i];
271  }
272 
273 
274  //===== iterator interface to rows of the matrix
275 
277  template<class T>
279  : public RandomAccessIteratorFacade<RealRowIterator<T>, T>
280  {
281 
282  public:
285 
286  friend class RandomAccessIteratorFacade<RealRowIterator<const ValueType>, const ValueType>;
287  friend class RandomAccessIteratorFacade<RealRowIterator<ValueType>, ValueType>;
288  friend class RealRowIterator<const ValueType>;
289  friend class RealRowIterator<ValueType>;
290 
293  : p(_p), i(_i)
294  {}
295 
298  : p(0), i(0)
299  {}
300 
302  : p(it.p), i(it.i)
303  {}
304 
305 
307  size_type index () const
308  {
309  return i;
310  }
311 
312  std::ptrdiff_t distanceTo(const RealRowIterator<ValueType>& other) const
313  {
314  assert(other.p==p);
315  return (other.i-i);
316  }
317 
318  std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
319  {
320  assert(other.p==p);
321  return (other.i-i);
322  }
323 
325  bool equals (const RealRowIterator<ValueType>& other) const
326  {
327  assert(other.p==p);
328  return i==other.i;
329  }
330 
332  bool equals (const RealRowIterator<const ValueType>& other) const
333  {
334  assert(other.p==p);
335  return i==other.i;
336  }
337 
338  private:
340  void increment()
341  {
342  ++i;
343  }
344 
346  void decrement()
347  {
348  --i;
349  }
350 
351  void advance(std::ptrdiff_t diff)
352  {
353  i+=diff;
354  }
355 
356  T& elementAt(std::ptrdiff_t diff) const
357  {
358  return p[i+diff];
359  }
360 
362  row_type& dereference () const
363  {
364  return p[i];
365  }
366 
367  row_type* p;
368  size_type i;
369  };
370 
374 
377  {
378  return Iterator(r,0);
379  }
380 
383  {
384  return Iterator(r,n);
385  }
386 
390  {
391  return Iterator(r,n-1);
392  }
393 
397  {
398  return Iterator(r,-1);
399  }
400 
403 
405  typedef typename row_type::Iterator ColIterator;
406 
410 
411 
414  {
415  return ConstIterator(r,0);
416  }
417 
420  {
421  return ConstIterator(r,n);
422  }
423 
427  {
428  return ConstIterator(r,n-1);
429  }
430 
434  {
435  return ConstIterator(r,-1);
436  }
437 
440 
443 
444  //===== constructors & resizers
445 
448  : build_mode(unknown), ready(notbuilt), n(0), m(0), nnz(0),
449  r(0), a(0)
450  {}
451 
454  : build_mode(bm), ready(notbuilt)
455  {
456  allocate(_n, _m, _nnz);
457  }
458 
461  : build_mode(bm), ready(notbuilt)
462  {
463  allocate(_n, _m);
464  }
465 
471  BCRSMatrix (const BCRSMatrix& Mat)
472  : n(Mat.n), nnz(0)
473  {
474  // deep copy in global array
475  size_type _nnz = Mat.nnz;
476 
477  // in case of row-wise allocation
478  if (_nnz<=0)
479  {
480  _nnz = 0;
481  for (size_type i=0; i<n; i++)
482  _nnz += Mat.r[i].getsize();
483  }
484 
485  j = Mat.j; // enable column index sharing, release array in case of row-wise allocation
486  allocate(Mat.n, Mat.m, _nnz);
487 
488  // build window structure
489  copyWindowStructure(Mat);
490  }
491 
494  {
495  deallocate();
496  }
497 
503  {
504  if(ready==notbuilt)
505  build_mode = bm;
506  else
507  DUNE_THROW(InvalidStateException, "Matrix structure is already built (ready="<<ready<<").");
508  }
509 
525  void setSize(size_type rows, size_type columns, size_type nnz=0)
526  {
527  // deallocate already setup memory
528  deallocate();
529 
530  // allocate matrix memory
531  allocate(rows, columns, nnz);
532  }
533 
541  {
542  // return immediately when self-assignment
543  if (&Mat==this) return *this;
544 
545  // make it simple: ALWAYS throw away memory for a and j
546  deallocate(false);
547 
548  // reallocate the rows if required
549  if (n>0 && n!=Mat.n) {
550  // free rows
551  for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
552  rowAllocator_.destroy(riter);
553  rowAllocator_.deallocate(r,n);
554  }
555 
556  nnz=Mat.nnz;
557  if (nnz<=0)
558  {
559  for (size_type i=0; i<Mat.n; i++)
560  nnz += Mat.r[i].getsize();
561  }
562 
563  // allocate a, share j
564  j = Mat.j;
565  allocate(Mat.n, Mat.m, nnz, n!=Mat.n);
566 
567  // build window structure
568  copyWindowStructure(Mat);
569  return *this;
570  }
571 
574  {
575  for (size_type i=0; i<n; i++) r[i] = k;
576  return *this;
577  }
578 
579  //===== row-wise creation interface
580 
583  {
584  public:
587  : Mat(_Mat), i(_i), nnz(0), current_row(Mat.a, Mat.j.get(), 0)
588  {
589  if (i==0 && Mat.ready)
590  DUNE_THROW(ISTLError,"creation only allowed for uninitialized matrix");
591  if(Mat.build_mode!=row_wise)
592  {
593  if(Mat.build_mode==unknown)
594  Mat.build_mode=row_wise;
595  else
596  DUNE_THROW(ISTLError,"creation only allowed if row wise allocation was requested in the constructor");
597  }
598  }
599 
602  {
603  // this should only be called if matrix is in creation
604  if (Mat.ready)
605  DUNE_THROW(ISTLError,"matrix already built up");
606 
607  // row i is defined through the pattern
608  // get memory for the row and initialize the j array
609  // this depends on the allocation mode
610 
611  // compute size of the row
612  size_type s = pattern.size();
613 
614  if(s>0){
615  // update number of nonzeroes including this row
616  nnz += s;
617 
618  // alloc memory / set window
619  if (Mat.nnz>0)
620  {
621  // memory is allocated in one long array
622 
623  // check if that memory is sufficient
624  if (nnz>Mat.nnz)
625  DUNE_THROW(ISTLError,"allocated nnz too small");
626 
627  // set row i
628  Mat.r[i].set(s,current_row.getptr(),current_row.getindexptr());
629  current_row.setptr(current_row.getptr()+s);
630  current_row.setindexptr(current_row.getindexptr()+s);
631  }else{
632  // memory is allocated individually per row
633  // allocate and set row i
634  B* a = Mat.allocator_.allocate(s);
635  new (a) B[s];
636  size_type* j = Mat.sizeAllocator_.allocate(s);
637  new (j) size_type[s];
638  Mat.r[i].set(s,a,j);
639  }
640  }else
641  // setup empty row
642  Mat.r[i].set(0,0,0);
643 
644  // initialize the j array for row i from pattern
645  size_type k=0;
646  size_type *j = Mat.r[i].getindexptr();
647  for (typename PatternType::const_iterator it=pattern.begin(); it!=pattern.end(); ++it)
648  j[k++] = *it;
649 
650  // now go to next row
651  i++;
652  pattern.clear();
653 
654  // check if this was last row
655  if (i==Mat.n)
656  {
657  Mat.ready = built;
658  if(Mat.nnz>0)
659  // Set nnz to the exact number of nonzero blocks inserted
660  // as some methods rely on it
661  Mat.nnz=nnz;
662  }
663  // done
664  return *this;
665  }
666 
668  bool operator!= (const CreateIterator& it) const
669  {
670  return (i!=it.i) || (&Mat!=&it.Mat);
671  }
672 
674  bool operator== (const CreateIterator& it) const
675  {
676  return (i==it.i) && (&Mat==&it.Mat);
677  }
678 
680  size_type index () const
681  {
682  return i;
683  }
684 
686  void insert (size_type j)
687  {
688  pattern.insert(j);
689  }
690 
693  {
694  if (pattern.find(j)!=pattern.end())
695  return true;
696  else
697  return false;
698  }
704  size_type size() const
705  {
706  return pattern.size();
707  }
708 
709  private:
710  BCRSMatrix& Mat; // the matrix we are defining
711  size_type i; // current row to be defined
712  size_type nnz; // count total number of nonzeros
713  typedef std::set<size_type,std::less<size_type> > PatternType;
714  PatternType pattern; // used to compile entries in a row
715  row_type current_row; // row pointing to the current row to setup
716  };
717 
719  friend class CreateIterator;
720 
723  {
724  return CreateIterator(*this,0);
725  }
726 
729  {
730  return CreateIterator(*this,n);
731  }
732 
733 
734  //===== random creation interface
735 
738  {
739  if (build_mode!=random)
740  DUNE_THROW(ISTLError,"requires random build mode");
741  if (ready)
742  DUNE_THROW(ISTLError,"matrix row sizes already built up");
743 
744  r[i].setsize(s);
745  }
746 
749  {
750 #ifdef DUNE_ISTL_WITH_CHECKING
751  if (r==0) DUNE_THROW(ISTLError,"row not initialized yet");
752  if (i>=n) DUNE_THROW(ISTLError,"index out of range");
753 #endif
754  return r[i].getsize();
755  }
756 
759  {
760  if (build_mode!=random)
761  DUNE_THROW(ISTLError,"requires random build mode");
762  if (ready)
763  DUNE_THROW(ISTLError,"matrix row sizes already built up");
764 
765  r[i].setsize(r[i].getsize()+s);
766  }
767 
769  void endrowsizes ()
770  {
771  if (build_mode!=random)
772  DUNE_THROW(ISTLError,"requires random build mode");
773  if (ready)
774  DUNE_THROW(ISTLError,"matrix row sizes already built up");
775 
776  // compute total size, check positivity
777  size_type total=0;
778  for (size_type i=0; i<n; i++)
779  {
780  if (r[i].getsize()<0)
781  DUNE_THROW(ISTLError,"rowsize must be nonnegative");
782  total += r[i].getsize();
783  }
784 
785  if(nnz==0)
786  // allocate/check memory
787  allocate(n,m,total,false);
788  else if(nnz<total)
789  DUNE_THROW(ISTLError,"Specified number of nonzeros ("<<nnz<<") not "
790  <<"sufficient for calculated nonzeros ("<<total<<"! ");
791 
792  // set the window pointers correctly
793  setWindowPointers(begin());
794 
795  // initialize j array with m (an invalid column index)
796  // this indicates an unused entry
797  for (size_type k=0; k<nnz; k++)
798  j.get()[k] = m;
799  ready = rowSizesBuilt;
800  }
801 
803 
814  {
815  if (build_mode!=random)
816  DUNE_THROW(ISTLError,"requires random build mode");
817  if (ready==built)
818  DUNE_THROW(ISTLError,"matrix already built up");
819  if (ready==notbuilt)
820  DUNE_THROW(ISTLError,"matrix row sizes not built up yet");
821 
822  if (col >= m)
823  DUNE_THROW(ISTLError,"column index exceeds matrix size");
824 
825  // get row range
826  size_type* const first = r[row].getindexptr();
827  size_type* const last = first + r[row].getsize();
828 
829  // find correct insertion position for new column index
830  size_type* pos = std::lower_bound(first,last,col);
831 
832  // check if index is already in row
833  if (pos!=last && *pos == col) return;
834 
835  // find end of already inserted column indices
836  size_type* end = std::lower_bound(pos,last,m);
837  if (end==last)
838  DUNE_THROW(ISTLError,"row is too small");
839 
840  // insert new column index at correct position
841  std::copy_backward(pos,end,end+1);
842  *pos = col;
843 
844  }
845 
847  void endindices ()
848  {
849  if (build_mode!=random)
850  DUNE_THROW(ISTLError,"requires random build mode");
851  if (ready==built)
852  DUNE_THROW(ISTLError,"matrix already built up");
853  if (ready==notbuilt)
854  DUNE_THROW(ISTLError,"row sizes are not built up yet");
855 
856  // check if there are undefined indices
857  RowIterator endi=end();
858  for (RowIterator i=begin(); i!=endi; ++i)
859  {
860  ColIterator endj = (*i).end();
861  for (ColIterator j=(*i).begin(); j!=endj; ++j){
862  if (j.index()<0)
863  {
864  std::cout << "j[" << j.offset() << "]=" << j.index() << std::endl;
865  DUNE_THROW(ISTLError,"undefined index detected");
866  }
867  if (j.index()>=m){
868  dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
869  <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
870  r[i.index()].setsize(j.offset());
871  break;
872  }
873  }
874  }
875 
876  // if not, set matrix to built
877  ready = built;
878  }
879 
880  //===== vector space arithmetic
881 
884  {
885  if (nnz>0)
886  {
887  // process 1D array
888  for (size_type i=0; i<nnz; i++)
889  a[i] *= k;
890  }
891  else
892  {
893  RowIterator endi=end();
894  for (RowIterator i=begin(); i!=endi; ++i)
895  {
896  ColIterator endj = (*i).end();
897  for (ColIterator j=(*i).begin(); j!=endj; ++j)
898  (*j) *= k;
899  }
900  }
901 
902  return *this;
903  }
904 
907  {
908  if (nnz>0)
909  {
910  // process 1D array
911  for (size_type i=0; i<nnz; i++)
912  a[i] /= k;
913  }
914  else
915  {
916  RowIterator endi=end();
917  for (RowIterator i=begin(); i!=endi; ++i)
918  {
919  ColIterator endj = (*i).end();
920  for (ColIterator j=(*i).begin(); j!=endj; ++j)
921  (*j) /= k;
922  }
923  }
924 
925  return *this;
926  }
927 
928 
935  {
936 #ifdef DUNE_ISTL_WITH_CHECKING
937  if(N()!=b.N() || M() != b.M())
938  DUNE_THROW(RangeError, "Matrix sizes do not match!");
939 #endif
940  RowIterator endi=end();
941  ConstRowIterator j=b.begin();
942  for (RowIterator i=begin(); i!=endi; ++i, ++j){
943  i->operator+=(*j);
944  }
945 
946  return *this;
947  }
948 
955  {
956 #ifdef DUNE_ISTL_WITH_CHECKING
957  if(N()!=b.N() || M() != b.M())
958  DUNE_THROW(RangeError, "Matrix sizes do not match!");
959 #endif
960  RowIterator endi=end();
961  ConstRowIterator j=b.begin();
962  for (RowIterator i=begin(); i!=endi; ++i, ++j){
963  i->operator-=(*j);
964  }
965 
966  return *this;
967  }
968 
978  {
979 #ifdef DUNE_ISTL_WITH_CHECKING
980  if(N()!=b.N() || M() != b.M())
981  DUNE_THROW(RangeError, "Matrix sizes do not match!");
982 #endif
983  RowIterator endi=end();
984  ConstRowIterator j=b.begin();
985  for(RowIterator i=begin(); i!=endi; ++i, ++j)
986  i->axpy(alpha, *j);
987 
988  return *this;
989  }
990 
991  //===== linear maps
992 
994  template<class X, class Y>
995  void mv (const X& x, Y& y) const
996  {
997 #ifdef DUNE_ISTL_WITH_CHECKING
998  if (x.N()!=M()) DUNE_THROW(ISTLError,
999  "Size mismatch: M: " << N() << "x" << M() << " x: " << x.N());
1000  if (y.N()!=N()) DUNE_THROW(ISTLError,
1001  "Size mismatch: M: " << N() << "x" << M() << " y: " << y.N());
1002 #endif
1003  ConstRowIterator endi=end();
1004  for (ConstRowIterator i=begin(); i!=endi; ++i)
1005  {
1006  y[i.index()]=0;
1007  ConstColIterator endj = (*i).end();
1008  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1009  (*j).umv(x[j.index()],y[i.index()]);
1010  }
1011  }
1012 
1014  template<class X, class Y>
1015  void umv (const X& x, Y& y) const
1016  {
1017 #ifdef DUNE_ISTL_WITH_CHECKING
1018  if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1019  if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1020 #endif
1021  ConstRowIterator endi=end();
1022  for (ConstRowIterator i=begin(); i!=endi; ++i)
1023  {
1024  ConstColIterator endj = (*i).end();
1025  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1026  (*j).umv(x[j.index()],y[i.index()]);
1027  }
1028  }
1029 
1031  template<class X, class Y>
1032  void mmv (const X& x, Y& y) const
1033  {
1034 #ifdef DUNE_ISTL_WITH_CHECKING
1035  if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1036  if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1037 #endif
1038  ConstRowIterator endi=end();
1039  for (ConstRowIterator i=begin(); i!=endi; ++i)
1040  {
1041  ConstColIterator endj = (*i).end();
1042  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1043  (*j).mmv(x[j.index()],y[i.index()]);
1044  }
1045  }
1046 
1048  template<class X, class Y>
1049  void usmv (const field_type& alpha, const X& x, Y& y) const
1050  {
1051 #ifdef DUNE_ISTL_WITH_CHECKING
1052  if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1053  if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1054 #endif
1055  ConstRowIterator endi=end();
1056  for (ConstRowIterator i=begin(); i!=endi; ++i)
1057  {
1058  ConstColIterator endj = (*i).end();
1059  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1060  (*j).usmv(alpha,x[j.index()],y[i.index()]);
1061  }
1062  }
1063 
1065  template<class X, class Y>
1066  void mtv (const X& x, Y& y) const
1067  {
1068 #ifdef DUNE_ISTL_WITH_CHECKING
1069  if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1070  if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1071 #endif
1072  for(size_type i=0; i<y.N(); ++i)
1073  y[i]=0;
1074  umtv(x,y);
1075  }
1076 
1078  template<class X, class Y>
1079  void umtv (const X& x, Y& y) const
1080  {
1081 #ifdef DUNE_ISTL_WITH_CHECKING
1082  if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1083  if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1084 #endif
1085  ConstRowIterator endi=end();
1086  for (ConstRowIterator i=begin(); i!=endi; ++i)
1087  {
1088  ConstColIterator endj = (*i).end();
1089  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1090  (*j).umtv(x[i.index()],y[j.index()]);
1091  }
1092  }
1093 
1095  template<class X, class Y>
1096  void mmtv (const X& x, Y& y) const
1097  {
1098 #ifdef DUNE_ISTL_WITH_CHECKING
1099  if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1100  if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1101 #endif
1102  ConstRowIterator endi=end();
1103  for (ConstRowIterator i=begin(); i!=endi; ++i)
1104  {
1105  ConstColIterator endj = (*i).end();
1106  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1107  (*j).mmtv(x[i.index()],y[j.index()]);
1108  }
1109  }
1110 
1112  template<class X, class Y>
1113  void usmtv (const field_type& alpha, const X& x, Y& y) const
1114  {
1115 #ifdef DUNE_ISTL_WITH_CHECKING
1116  if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1117  if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1118 #endif
1119  ConstRowIterator endi=end();
1120  for (ConstRowIterator i=begin(); i!=endi; ++i)
1121  {
1122  ConstColIterator endj = (*i).end();
1123  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1124  (*j).usmtv(alpha,x[i.index()],y[j.index()]);
1125  }
1126  }
1127 
1129  template<class X, class Y>
1130  void umhv (const X& x, Y& y) const
1131  {
1132 #ifdef DUNE_ISTL_WITH_CHECKING
1133  if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1134  if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1135 #endif
1136  ConstRowIterator endi=end();
1137  for (ConstRowIterator i=begin(); i!=endi; ++i)
1138  {
1139  ConstColIterator endj = (*i).end();
1140  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1141  (*j).umhv(x[i.index()],y[j.index()]);
1142  }
1143  }
1144 
1146  template<class X, class Y>
1147  void mmhv (const X& x, Y& y) const
1148  {
1149 #ifdef DUNE_ISTL_WITH_CHECKING
1150  if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1151  if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1152 #endif
1153  ConstRowIterator endi=end();
1154  for (ConstRowIterator i=begin(); i!=endi; ++i)
1155  {
1156  ConstColIterator endj = (*i).end();
1157  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1158  (*j).mmhv(x[i.index()],y[j.index()]);
1159  }
1160  }
1161 
1163  template<class X, class Y>
1164  void usmhv (const field_type& alpha, const X& x, Y& y) const
1165  {
1166 #ifdef DUNE_ISTL_WITH_CHECKING
1167  if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
1168  if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
1169 #endif
1170  ConstRowIterator endi=end();
1171  for (ConstRowIterator i=begin(); i!=endi; ++i)
1172  {
1173  ConstColIterator endj = (*i).end();
1174  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1175  (*j).usmhv(alpha,x[i.index()],y[j.index()]);
1176  }
1177  }
1178 
1179 
1180  //===== norms
1181 
1183  double frobenius_norm2 () const
1184  {
1185  double sum=0;
1186 
1187  ConstRowIterator endi=end();
1188  for (ConstRowIterator i=begin(); i!=endi; ++i)
1189  {
1190  ConstColIterator endj = (*i).end();
1191  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1192  sum += (*j).frobenius_norm2();
1193  }
1194 
1195  return sum;
1196  }
1197 
1199  double frobenius_norm () const
1200  {
1201  return sqrt(frobenius_norm2());
1202  }
1203 
1205  double infinity_norm () const
1206  {
1207  double max=0;
1208  ConstRowIterator endi=end();
1209  for (ConstRowIterator i=begin(); i!=endi; ++i)
1210  {
1211  double sum=0;
1212  ConstColIterator endj = (*i).end();
1213  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1214  sum += (*j).infinity_norm();
1215  max = std::max(max,sum);
1216  }
1217  return max;
1218  }
1219 
1221  double infinity_norm_real () const
1222  {
1223  double max=0;
1224  ConstRowIterator endi=end();
1225  for (ConstRowIterator i=begin(); i!=endi; ++i)
1226  {
1227  double sum=0;
1228  ConstColIterator endj = (*i).end();
1229  for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
1230  sum += (*j).infinity_norm_real();
1231  max = std::max(max,sum);
1232  }
1233  return max;
1234  }
1235 
1236 
1237  //===== sizes
1238 
1240  size_type N () const
1241  {
1242  return n;
1243  }
1244 
1246  size_type M () const
1247  {
1248  return m;
1249  }
1250 
1253  {
1254  return nnz;
1255  }
1256 
1257 
1258  //===== query
1259 
1261  bool exists (size_type i, size_type j) const
1262  {
1263 #ifdef DUNE_ISTL_WITH_CHECKING
1264  if (i<0 || i>=n) DUNE_THROW(ISTLError,"row index out of range");
1265  if (j<0 || j>=m) DUNE_THROW(ISTLError,"column index out of range");
1266 #endif
1267  if (r[i].size() && r[i].find(j)!=r[i].end())
1268  return true;
1269  else
1270  return false;
1271  }
1272 
1273 
1274  private:
1275  // state information
1276  BuildMode build_mode; // row wise or whole matrix
1277  BuildStage ready; // indicate the stage the matrix building is in
1278 
1279  // The allocator used for memory management
1280  typename A::template rebind<B>::other allocator_;
1281 
1282  typename A::template rebind<row_type>::other rowAllocator_;
1283 
1284  typename A::template rebind<size_type>::other sizeAllocator_;
1285 
1286  // size of the matrix
1287  size_type n; // number of rows
1288  size_type m; // number of columns
1289  size_type nnz; // number of nonzeros allocated in the a and j array below
1290  // zero means that memory is allocated separately for each row.
1291 
1292  // the rows are dynamically allocated
1293  row_type* r; // [n] the individual rows having pointers into a,j arrays
1294 
1295  // dynamically allocated memory
1296  B* a; // [nnz] non-zero entries of the matrix in row-wise ordering
1297  // If a single array of column indices is used, it can be shared
1298  // between different matrices with the same sparsity pattern
1299  Dune::shared_ptr<size_type> j; // [nnz] column indices of entries
1300 
1301 
1302  void setWindowPointers(ConstRowIterator row)
1303  {
1304  row_type current_row(a,j.get(),0); // Pointers to current row data
1305  for (size_type i=0; i<n; i++, ++row){
1306  // set row i
1307  size_type s = row->getsize();
1308 
1309  if (s>0){
1310  // setup pointers and size
1311  r[i].set(s,current_row.getptr(), current_row.getindexptr());
1312  // update pointer for next row
1313  current_row.setptr(current_row.getptr()+s);
1314  current_row.setindexptr(current_row.getindexptr()+s);
1315  } else{
1316  // empty row
1317  r[i].set(0,0,0);
1318  }
1319  }
1320  }
1321 
1323  void copyWindowStructure(const BCRSMatrix& Mat)
1324  {
1325  setWindowPointers(Mat.begin());
1326 
1327  // copy data
1328  for (size_type i=0; i<n; i++) r[i] = Mat.r[i];
1329 
1330  // finish off
1331  build_mode = row_wise; // dummy
1332  ready = built;
1333  }
1334 
1340  void deallocate(bool deallocateRows=true)
1341  {
1342 
1343  if (nnz>0)
1344  {
1345  // a,j have been allocated as one long vector
1346  j.reset();
1347  for(B *aiter=a+(nnz-1), *aend=a-1; aiter!=aend; --aiter)
1348  allocator_.destroy(aiter);
1349  allocator_.deallocate(a,n);
1350  }
1351  else
1352  {
1353  // check if memory for rows have been allocated individually
1354  for (size_type i=0; i<n; i++)
1355  if (r[i].getsize()>0)
1356  {
1357  for (B *col=r[i].getptr()+(r[i].getsize()-1),
1358  *colend = r[i].getptr()-1; col!=colend; --col) {
1359  allocator_.destroy(col);
1360  }
1361  sizeAllocator_.deallocate(r[i].getindexptr(),1);
1362  allocator_.deallocate(r[i].getptr(),1);
1363  }
1364  }
1365 
1366  // deallocate the rows
1367  if (n>0 && deallocateRows) {
1368  for(row_type *riter=r+(n-1), *rend=r-1; riter!=rend; --riter)
1369  rowAllocator_.destroy(riter);
1370  rowAllocator_.deallocate(r,n);
1371  }
1372 
1373  // Mark matrix as not built at all.
1374  ready=notbuilt;
1375 
1376  }
1377 
1379  class Deallocator
1380  {
1381  typename A::template rebind<size_type>::other& sizeAllocator_;
1382 
1383  public:
1384  Deallocator(typename A::template rebind<size_type>::other& sizeAllocator)
1385  : sizeAllocator_(sizeAllocator)
1386  {}
1387 
1388  void operator()(size_type* p) { sizeAllocator_.deallocate(p,1); }
1389  };
1390 
1391 
1409  void allocate(size_type rows, size_type columns, size_type nnz_=0, bool allocateRows=true)
1410  {
1411  // Store size
1412  n = rows;
1413  m = columns;
1414  nnz = nnz_;
1415 
1416  // allocate rows
1417  if(allocateRows){
1418  if (n>0){
1419  r = rowAllocator_.allocate(rows);
1420  new (r) row_type[rows];
1421  }else{
1422  r = 0;
1423  }
1424  }
1425 
1426 
1427  // allocate a and j array
1428  if (nnz>0){
1429  a = allocator_.allocate(nnz);
1430  // allocate column indices only if not yet present (enable sharing)
1431  if (!j.get())
1432  j.reset(sizeAllocator_.allocate(nnz),Deallocator(sizeAllocator_));
1433  }else{
1434  a = 0;
1435  j.reset();
1436  for(row_type* ri=r; ri!=r+rows;++ri)
1437  rowAllocator_.construct(ri, row_type());
1438  }
1439 
1440  // Mark the matrix as not built.
1441  ready = notbuilt;
1442  }
1443 
1444  };
1445 
1446 
1449 } // end namespace
1450 
1451 #endif