I started the linear algebra things. The module is called

Iterative Solver Template Library (sounds like something big ...)

I started with the vectors. It is not ready yet, don't use it.
Explanations follow...

Peter

[[Imported from SVN: r22]]

istl/bvector.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_BVECTOR_HH__
+#define __DUNE_BVECTOR_HH__
+
+#include <math.h>
+#include <complex>
+
+#include "../common/exceptions.hh"
+
+#include "spmv.hh"
+
+/*! \file __FILE__
+
+   This file implements a vector space as a tensor product of
+   a given vector space. The number of components can be given at
+   run-time.
+ */
+
+namespace Dune {
+
+  /** @defgroup SPMV Dune Sparse Matrix Vector Templates
+              @addtogroup SPMV
+              @{
+   */
+
+  /**! Construct a vector space out of a tensor product of other
+           vector spaces. The number of components is given at run-time
+
+           Vector has a window mode where it is assumed that the
+       object provides a window into a larger vector. Dynamic
+           memory management is disabled in this case. The window
+           is manipulated with the setwindow function.
+
+           Error checking: no error checking is provided normally.
+           Setting the compile time switch DUNE_ISTL_WITH_CHECKING
+           enables error checking.
+   */
+  template<class B, class A=SPMV_Allocator>
+  class BlockVector
+  {
+  public:
+
+    //===== type definitions and constants
+
+    //! export the type representing the field
+    typedef typename B::field_type field_type;
+
+    //! export the type representing the components
+    typedef B block_type;
+
+    //! export the allocator type
+    typedef A allocator_type;
+
+    //! increment block level counter
+    enum {blocklevel = B::blocklevel+1};
+
+
+    //===== constructors and such
+
+    //! makes empty vector in non-window mode
+    BlockVector ()
+    {
+      n = 0;
+      p = 0;
+      wmode = false;
+    }
+
+    //! makes empty vector and selects mode
+    BlockVector (bool _wmode)     // this is the only way to make a vector in window mode
+    {
+      n = 0;
+      p = 0;
+      wmode = _wmode;
+    }
+
+    //! make vector with _n components
+    BlockVector (int _n)
+    {
+      wmode = false;     // this is non-window mode
+      n = _n;
+      if (n>0)
+        p = A::template malloc<B>(n);
+      else
+      {
+        n = 0;
+        p = 0;
+      }
+    }
+
+    //! copy constructor
+    BlockVector (const BlockVector& a)
+    {
+      wmode = false;     // this is non-window mode
+
+      // allocate memory with same size as a
+      n = a.n;
+      if (n>0)
+        p = A::template malloc<B>(n);
+      else
+      {
+        n = 0;
+        p = 0;
+      }
+
+      // und kopiere Elemente
+      for (int i=0; i<n; i++) p[i]=a.p[i];
+    }
+
+    //! free dynamic memory
+    ~BlockVector ()
+    {
+      if (n>0 && !wmode) A::template free<B>(p);
+    }
+
+    //! reallocate vector to given size, any data is lost
+    void resize (int _n)
+    {
+      if (n==_n) return;
+
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (wmode) DUNE_THROW(ISTLError,"no resize in window mode possible");
+#endif
+
+      if (n>0) A::template free<B>(p);
+      n = _n;
+      if (n>0)
+        p = A::template malloc<B>(n);
+      else
+      {
+        n = 0;
+        p = 0;
+      }
+    }
+
+    //! assignment
+    BlockVector& operator= (const BlockVector& a)
+    {
+      if (&a!=this)     // check if this and a are different objects
+      {
+        // adjust size of array
+        if (n!=a.n)           // check if size is different
+        {
+
+#ifdef DUNE_ISTL_WITH_CHECKING
+          if (wmode) DUNE_THROW(ISTLError,"no resizing in window mode possible");
+#endif
+
+          if (n>0) A::template free<B>(p);                 // delete old memory
+          n = a.n;
+          if (n>0)
+            p = A::template malloc<B>(n);
+          else
+          {
+            n = 0;
+            p = 0;
+          }
+        }
+        // copy data
+        for (int i=0; i<n; i++) p[i]=a.p[i];
+      }
+      return *this;
+    }
+
+    //===== assignment from scalar
+    BlockVector& operator= (const field_type& k)
+    {
+      for (int i=0; i<n; i++)
+        p[i] = k;
+      return *this;
+    }
+
+    //===== access to components
+
+    //! random access to blocks
+    B& operator[] (int i)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+      return p[i];
+    }
+
+    //! same for read only access
+    const B& operator[] (int i) const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+      return p[i];
+    }
+
+    //! random access via round brackets supplied to be consistent with matrix
+    B& operator() (int i)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+      return p[i];
+    }
+
+    //! same for read only access
+    const B& operator() (int i) const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+      return p[i];
+    }
+
+    //! Iterator class for sequential access
+    class Iterator
+    {
+    public:
+      //! constructor
+      Iterator (B* _p)
+      {
+        p = _p;
+      }
+
+      //! prefix increment
+      Iterator& operator++()
+      {
+        ++p;
+        return *this;
+      }
+
+      //! equality
+      bool operator== (const Iterator& i) const
+      {
+        return p==i.p;
+      }
+
+      //! inequality
+      bool operator!= (const Iterator& i) const
+      {
+        return p!=i.p;
+      }
+
+      //! dereferencing
+      B& operator* ()
+      {
+        return *p;
+      }
+
+      //! arrow
+      B* operator-> ()
+      {
+        return p;
+      }
+
+    private:
+      B* _p;
+    };
+
+    //! begin iterator
+    Iterator begin () const
+    {
+      return Iterator(p);
+    }
+
+    //! end iterator
+    Iterator end () const
+    {
+      return Iterator(p+n);
+    }
+
+
+    //===== vector space arithmetic
+
+    //! vector space addition
+    BlockVector& operator+= (const BlockVector& y)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (n!=y.n) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+      for (int i=0; i<n; ++i) p[i] += y.p[i];
+      return *this;
+    }
+
+    //! vector space subtraction
+    BlockVector& operator-= (const BlockVector& y)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (n!=y.n) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+      for (int i=0; i<n; ++i) p[i] -= y.p[i];
+      return *this;
+    }
+
+    //! vector space multiplication with scalar
+    BlockVector& operator*= (const field_type& k)
+    {
+      for (int i=0; i<n; ++i) p[i] *= k;
+      return *this;
+    }
+
+    //! vector space division by scalar
+    BlockVector& operator/= (const field_type& k)
+    {
+      for (int i=0; i<n; ++i) p[i] /= k;
+      return *this;
+    }
+
+    //! vector space axpy operation
+    BlockVector& axpy (const field_type& a, const BlockVector& y)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (n!=y.n) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+      for (int i=0; i<n; ++i) p[i].axpy(a,y.p[i]);
+      return *this;
+    }
+
+
+    //===== Euclidean scalar product
+
+    //! scalar product
+    field_type operator* (const BlockVector& y) const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (n!=y.n) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+      field_type sum=0;
+      for (int i=0; i<n; ++i) sum += p[i]*y.p[i];
+      return sum;
+    }
+
+
+    //===== norms
+
+    //! one norm (sum over absolute values of entries)
+    double one_norm () const
+    {
+      double sum=0;
+      for (int i=0; i<n; ++i) sum += p[i].one_norm();
+      return sum;
+    }
+
+    //! simplified one norm (uses Manhattan norm for complex values)
+    double one_norm_real () const
+    {
+      double sum=0;
+      for (int i=0; i<n; ++i) sum += p[i].one_norm_real();
+      return sum;
+    }
+
+    //! two norm sqrt(sum over squared values of entries)
+    double two_norm () const
+    {
+      double sum=0;
+      for (int i=0; i<n; ++i) sum += p[i].two_norm2();
+      return sqrt(sum);
+    }
+
+    //! sqare of two norm (sum over squared values of entries), need for block recursion
+    double two_norm2 () const
+    {
+      double sum=0;
+      for (int i=0; i<n; ++i) sum += p[i].two_norm2();
+      return sum;
+    }
+
+    //! infinity norm (maximum of absolute values of entries)
+    double infinity_norm () const
+    {
+      double max=0;
+      for (int i=0; i<n; ++i) max = std::max(max,p[i].infinity_norm());
+      return max;
+    }
+
+    //! simplified infinity norm (uses Manhattan norm for complex values)
+    double infinity_norm_real () const
+    {
+      double max=0;
+      for (int i=0; i<n; ++i) max = std::max(max,p[i].infinity_norm_real());
+      return max;
+    }
+
+
+    //===== sizes
+
+    //! number of blocks in the vector (are of size 1 here)
+    int N () const
+    {
+      return n;
+    }
+
+    //! dimension of the vector space
+    int dim () const
+    {
+      int d=0;
+      for (int i=0; i<n; i++)
+        d += p[i].dim();
+      return d;
+    }
+
+
+    //===== in window mode provide a function to set the window
+
+    //! set the window in window mode -- you must know what you do!
+    void set (const int& _n, B* _p)
+    {
+      n = _n;
+      p = _p;
+    }
+
+  private:
+    bool wmode;     // true if object provides window in larger vector
+    int n;          // number of elements in array
+    B *p;           // pointer to dynamically allocated built-in array
+  };
+
+
+
+  /** @} end documentation */
+
+} // end namespace
+
+#endif

istl/fvector.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_FVECTOR_HH__
+#define __DUNE_FVECTOR_HH__
+
+#include <math.h>
+#include <complex>
+
+#include "../common/exceptions.hh"
+
+#include "spmv.hh"
+
+/*! \file __FILE__
+
+   This file implements a vector constructed from a given type
+   representing a field and a compile-time given size.
+ */
+
+namespace Dune {
+
+  /** @defgroup SPMV Dune Sparse Matrix Vector Templates
+              @addtogroup SPMV
+              @{
+   */
+
+  // forward declaration of template
+  template<class K, int n> class FieldVector;
+
+  // template meta program for assignment from scalar
+  template<int I>
+  struct fvmeta_assignscalar {
+    template<class K, int n>
+    static K& assignscalar (FieldVector<K,n>& x, const K& k)
+    {
+      x[I] = fvmeta_assignscalar<I-1>::assignscalar(x,k);
+      return x[I];
+    }
+  };
+  template<>
+  struct fvmeta_assignscalar<0> {
+    template<class K, int n>
+    static K& assignscalar (FieldVector<K,n>& x, const K& k)
+    {
+      x[0] = k;
+      return x[0];
+    }
+  };
+
+  // template meta program for operator+=
+  template<int I>
+  struct fvmeta_plusequal {
+    template<class K, int n>
+    static void plusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
+    {
+      fvmeta_plusequal<I-1>::plusequal(x,y);
+      x[I] += y[I];
+    }
+  };
+  template<>
+  struct fvmeta_plusequal<0> {
+    template<class K, int n>
+    static void plusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
+    {
+      x[0] += y[0];
+    }
+  };
+
+  // template meta program for operator-=
+  template<int I>
+  struct fvmeta_minusequal {
+    template<class K, int n>
+    static void minusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
+    {
+      fvmeta_minusequal<I-1>::minusequal(x,y);
+      x[I] -= y[I];
+    }
+  };
+  template<>
+  struct fvmeta_minusequal<0> {
+    template<class K, int n>
+    static void minusequal (FieldVector<K,n>& x, const FieldVector<K,n>& y)
+    {
+      x[0] -= y[0];
+    }
+  };
+
+  // template meta program for operator*=
+  template<int I>
+  struct fvmeta_multequal {
+    template<class K, int n>
+    static void multequal (FieldVector<K,n>& x, const K& k)
+    {
+      fvmeta_multequal<I-1>::multequal(x,k);
+      x[I] *= k;
+    }
+  };
+  template<>
+  struct fvmeta_multequal<0> {
+    template<class K, int n>
+    static void multequal (FieldVector<K,n>& x, const K& k)
+    {
+      x[0] *= k;
+    }
+  };
+
+  // template meta program for operator/=
+  template<int I>
+  struct fvmeta_divequal {
+    template<class K, int n>
+    static void divequal (FieldVector<K,n>& x, const K& k)
+    {
+      fvmeta_divequal<I-1>::divequal(x,k);
+      x[I] /= k;
+    }
+  };
+  template<>
+  struct fvmeta_divequal<0> {
+    template<class K, int n>
+    static void divequal (FieldVector<K,n>& x, const K& k)
+    {
+      x[0] /= k;
+    }
+  };
+
+  // template meta program for axpy
+  template<int I>
+  struct fvmeta_axpy {
+    template<class K, int n>
+    static void axpy (FieldVector<K,n>& x, const K& a, const FieldVector<K,n>& y)
+    {
+      fvmeta_axpy<I-1>::axpy(x,a,y);
+      x[I] += a*y[I];
+    }
+  };
+  template<>
+  struct fvmeta_axpy<0> {
+    template<class K, int n>
+    static void axpy (FieldVector<K,n>& x, const K& a, const FieldVector<K,n>& y)
+    {
+      x[0] += a*y[0];
+    }
+  };
+
+  // template meta program for dot
+  template<int I>
+  struct fvmeta_dot {
+    template<class K, int n>
+    static K dot (const FieldVector<K,n>& x, const FieldVector<K,n>& y)
+    {
+      return x[I]*y[I] + fvmeta_dot<I-1>::dot(x,y);
+    }
+  };
+  template<>
+  struct fvmeta_dot<0> {
+    template<class K, int n>
+    static K dot (const FieldVector<K,n>& x, const FieldVector<K,n>& y)
+    {
+      return x[0]*y[0];
+    }
+  };
+
+  // some abs functions needed for the norms
+  template<class K>
+  double fvmeta_abs (const K& k)
+  {
+    return (k >= 0) ? k : -k;
+  }
+
+  template<class K>
+  double fvmeta_abs (const std::complex<K>& c)
+  {
+    return sqrt(c.real()*c.real() + c.imag()*c.imag());
+  }
+
+  template<class K>
+  double fvmeta_absreal (const K& k)
+  {
+    return fvmeta_abs(k);
+  }
+
+  template<class K>
+  double fvmeta_absreal (const std::complex<K>& c)
+  {
+    return fvmeta_abs(c.real()) + fvmeta_abs(c.imag());
+  }
+
+  template<class K>
+  double fvmeta_abs2 (const K& k)
+  {
+    return k*k;
+  }
+
+  template<class K>
+  double fvmeta_abs2 (const std::complex<K>& c)
+  {
+    return c.real()*c.real() + c.imag()*c.imag();
+  }
+
+  // template meta program for one_norm
+  template<int I>
+  struct fvmeta_one_norm {
+    template<class K, int n>
+    static double one_norm (const FieldVector<K,n>& x)
+    {
+      return fvmeta_abs(x[I]) + fvmeta_one_norm<I-1>::one_norm(x);
+    }
+  };
+  template<>
+  struct fvmeta_one_norm<0> {
+    template<class K, int n>
+    static double one_norm (const FieldVector<K,n>& x)
+    {
+      return fvmeta_abs(x[0]);
+    }
+  };
+
+  // template meta program for one_norm_real
+  template<int I>
+  struct fvmeta_one_norm_real {
+    template<class K, int n>
+    static double one_norm_real (const FieldVector<K,n>& x)
+    {
+      return fvmeta_absreal(x[I]) + fvmeta_one_norm_real<I-1>::one_norm_real(x);
+    }
+  };
+  template<>
+  struct fvmeta_one_norm_real<0> {
+    template<class K, int n>
+    static double one_norm_real (const FieldVector<K,n>& x)
+    {
+      return fvmeta_absreal(x[0]);
+    }
+  };
+
+  // template meta program for two_norm squared
+  template<int I>
+  struct fvmeta_two_norm2 {
+    template<class K, int n>
+    static double two_norm2 (const FieldVector<K,n>& x)
+    {
+      return fvmeta_abs2(x[I]) + fvmeta_two_norm2<I-1>::two_norm2(x);
+    }
+  };
+  template<>
+  struct fvmeta_two_norm2<0> {
+    template<class K, int n>
+    static double two_norm2 (const FieldVector<K,n>& x)
+    {
+      return fvmeta_abs2(x[0]);
+    }
+  };
+
+  // template meta program for infinity norm
+  template<int I>
+  struct fvmeta_infinity_norm {
+    template<class K, int n>
+    static double infinity_norm (const FieldVector<K,n>& x)
+    {
+      return std::max(fvmeta_abs(x[I]),fvmeta_infinity_norm<I-1>::infinity_norm(x));
+    }
+  };
+  template<>
+  struct fvmeta_infinity_norm<0> {
+    template<class K, int n>
+    static double infinity_norm (const FieldVector<K,n>& x)
+    {
+      return fvmeta_abs(x[0]);
+    }
+  };
+
+  // template meta program for simplified infinity norm
+  template<int I>
+  struct fvmeta_infinity_norm_real {
+    template<class K, int n>
+    static double infinity_norm_real (const FieldVector<K,n>& x)
+    {
+      return std::max(fvmeta_absreal(x[I]),fvmeta_infinity_norm_real<I-1>::infinity_norm_real(x));
+    }
+  };
+  template<>
+  struct fvmeta_infinity_norm_real<0> {
+    template<class K, int n>
+    static double infinity_norm_real (const FieldVector<K,n>& x)
+    {
+      return fvmeta_absreal(x[0]);
+    }
+  };
+
+  /**! Construct a vector space out of a tensor product of fields.
+           K is the field type (use float, double, complex, etc) and n
+       is the number of components.
+
+           It is generally assumed that K is a numerical type compatible with double
+       (E.g. norms are always computed in double precision).
+
+           Implementation of all members uses template meta programs where appropriate
+   */
+  template<class K, int n>
+  class FieldVector
+  {
+  public:
+    // standard constructor and everything is sufficient ...
+
+
+    //===== type definitions and constants
+
+    //! export the type representing the field
+    typedef K field_type;
+
+    //! export the type representing the components
+    typedef K block_type;
+
+    //! We are at the leaf of the block recursion
+    enum {blocklevel = 0};
+
+
+    //===== assignment from scalar
+    FieldVector& operator= (const K& k)
+    {
+      fvmeta_assignscalar<n-1>::assignscalar(*this,k);
+      return *this;
+    }
+
+
+    //===== access to components
+
+    //! random access
+    K& operator[] (int i)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+      return p[i];
+    }
+
+    //! same for read only access
+    const K& operator[] (int i) const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+      return p[i];
+    }
+
+    //! random access via round brackets supplied to be consistent with matrix
+    K& operator() (int i)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+      return p[i];
+    }
+
+    //! same for read only access
+    const K& operator() (int i) const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+      return p[i];
+    }
+
+    //! Iterator class for sequential access
+    class Iterator
+    {
+    public:
+      //! constructor
+      Iterator (K* _p)
+      {
+        p = _p;
+      }
+
+      //! prefix increment
+      Iterator& operator++()
+      {
+        ++p;
+        return *this;
+      }
+
+      //! equality
+      bool operator== (const Iterator& i) const
+      {
+        return p==i.p;
+      }
+
+      //! inequality
+      bool operator!= (const Iterator& i) const
+      {
+        return p!=i.p;
+      }
+
+      //! dereferencing
+      K& operator* ()
+      {
+        return *p;
+      }
+
+      //! arrow
+      K* operator-> ()
+      {
+        return p;
+      }
+
+    private:
+      K* _p;
+    };
+
+    //! begin iterator
+    Iterator begin () const
+    {
+      return Iterator(p);
+    }
+
+    //! end iterator
+    Iterator end () const
+    {
+      return Iterator(p+n);
+    }
+
+
+    //===== vector space arithmetic
+
+    //! vector space addition
+    FieldVector& operator+= (const FieldVector& y)
+    {
+      fvmeta_plusequal<n-1>::plusequal(*this,y);
+      return *this;
+    }
+
+    //! vector space subtraction
+    FieldVector& operator-= (const FieldVector& y)
+    {
+      fvmeta_minusequal<n-1>::minusequal(*this,y);
+      return *this;
+    }
+
+    //! vector space multiplication with scalar
+    FieldVector& operator*= (const K& k)
+    {
+      fvmeta_multequal<n-1>::multequal(*this,k);
+      return *this;
+    }
+
+    //! vector space division by scalar
+    FieldVector& operator/= (const K& k)
+    {
+      fvmeta_divequal<n-1>::divequal(*this,k);
+      return *this;
+    }
+
+    //! vector space axpy operation
+    FieldVector& axpy (const K& a, const FieldVector& y)
+    {
+      fvmeta_axpy<n-1>::axpy(*this,a,y);
+      return *this;
+    }
+
+
+    //===== Euclidean scalar product
+
+    //! scalar product
+    const K operator* (const FieldVector& y) const
+    {
+      return fvmeta_dot<n-1>::dot(*this,y);
+    }
+
+
+    //===== norms
+
+    //! one norm (sum over absolute values of entries)
+    double one_norm () const
+    {
+      return fvmeta_one_norm<n-1>::one_norm(*this);
+    }
+
+    //! simplified one norm (uses Manhattan norm for complex values)
+    double one_norm_real () const
+    {
+      return fvmeta_one_norm_real<n-1>::one_norm_real(*this);
+    }
+
+    //! two norm sqrt(sum over squared values of entries)
+    double two_norm () const
+    {
+      return sqrt(fvmeta_two_norm2<n-1>::two_norm2(*this));
+    }
+
+    //! sqare of two norm (sum over squared values of entries), need for block recursion
+    double two_norm2 () const
+    {
+      return fvmeta_two_norm2<n-1>::two_norm2(*this);
+    }
+
+    //! infinity norm (maximum of absolute values of entries)
+    double infinity_norm () const
+    {
+      return fvmeta_infinity_norm<n-1>::infinity_norm(*this);
+    }
+
+    //! simplified infinity norm (uses Manhattan norm for complex values)
+    double infinity_norm_real () const
+    {
+      return fvmeta_infinity_norm_real<n-1>::infinity_norm_real(*this);
+    }
+
+
+    //===== sizes
+
+    //! number of blocks in the vector (are of size 1 here)
+    int N () const
+    {
+      return n;
+    }
+
+    //! dimension of the vector space
+    int dim () const
+    {
+      return n;
+    }
+
+  private:
+    // the data, very simply a built in array
+    K p[n];
+  };
+
+  /** @} end documentation */
+
+} // end namespace
+
+#endif

istl/vbvector.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_VBVECTOR_HH__
+#define __DUNE_VBVECTOR_HH__
+
+#include <math.h>
+#include <complex>
+
+#include "../common/exceptions.hh"
+
+#include "spmv.hh"
+#include "bvector.hh"
+
+/*! \file __FILE__
+
+   This file implements a vector space as a tensor product of
+   a given vector space. The number of components can be given at
+   run-time as well as the size of the blocks.
+ */
+
+namespace Dune {
+
+  /** @defgroup SPMV Dune Sparse Matrix Vector Templates
+              @addtogroup SPMV
+              @{
+   */
+
+  /**! Construct a vector space out of a tensor product of other
+           vector spaces. The number of components is given at run-time.
+
+   */
+  template<class B, class A=SPMV_Allocator>
+  class VariableBlockVector
+  {
+  public:
+
+    //===== type definitions and constants
+
+    //! export the type representing the field
+    typedef typename B::field_type field_type;
+
+    //! export the allocator type
+    typedef A allocator_type;
+
+    //! export the type representing the components
+    typedef BlockVector<B> block_type;
+
+    //! increment block level counter, yes, it is two levels !
+    enum {blocklevel = B::blocklevel+2};
+
+
+    //===== constructors and such
+
+    //! constructor without arguments makes empty vector, put subobject in window mode
+    VariableBlockVector () : window(true)
+    {
+      // nothing is known ...
+      nblocks = 0;
+      s = 0;
+      e = 0;
+      n = 0;
+      p = 0;
+      initialized = false;
+    }
+
+    //! make vector with _n components, put subobject in window mode
+    VariableBlockVector (int _nblocks) : window(true)
+    {
+      // now the s and e array can be allocated
+      nblocks = _nblocks;
+      if (nblocks>0)
+      {
+        s = A::template malloc<int>(nblocks);
+        e = A::template malloc<int>(nblocks);
+      }
+      else
+      {
+        nblocks = 0;
+        s = 0;
+        e = 0;
+      }
+
+      // the p array is still not allocated
+      n = 0;
+      p = 0;
+
+      // and the vector not fully usable
+      initialized = false;
+    }
+
+    //! copy constructor, put subobject in window mode
+    VariableBlockVector (const VariableBlockVector& a) : window(true)
+    {
+      // allocate s and e arrays
+      nblocks = a.nblocks;
+      if (nblocks>0)
+      {
+        s = A::template malloc<int>(nblocks);
+        e = A::template malloc<int>(nblocks);
+      }
+      else
+      {
+        nblocks = 0;
+        s = 0;
+        e = 0;
+      }
+
+      // allocate p array
+      n = a.n;
+      if (n>0)
+        p = A::template malloc<B>(n);
+      else
+      {
+        n = 0;
+        p = 0;
+      }
+
+      // and copy the data
+      for (int i=0; i<nblocks; i++)
+      {
+        s[i]=a.s[i];
+        e[i]=a.e[i];
+      }
+      for (int i=0; i<n; i++) p[i]=a.p[i];
+
+      // and we have a usable vector
+      initialized = true;
+    }
+
+    //! free dynamic memory in reverse order
+    ~VariableBlockVector ()
+    {
+      if (n>0) A::template free<B>(p);
+      if (nblocks>0)
+      {
+        A::template free<int>(e);
+        A::template free<int>(s);
+      }
+    }
+
+    //! same effect as constructor with same argument
+    void resize (int _nblocks)
+    {
+      // deallocate in reverse order
+      if (n>0) A::template free<B>(p);
+      if (nblocks>0)
+      {
+        A::template free<int>(e);
+        A::template free<int>(s);
+      }
+
+      // now the s and e array can be allocated
+      nblocks = _nblocks;
+      if (nblocks>0)
+      {
+        s = A::template malloc<int>(nblocks);
+        e = A::template malloc<int>(nblocks);
+      }
+      else
+      {
+        nblocks = 0;
+        s = 0;
+        e = 0;
+      }
+
+      // the p array is still not allocated
+      n = 0;
+      p = 0;
+
+      // and the vector not fully usable
+      initialized = false;
+    }
+
+    //! assignment
+    VariableBlockVector& operator= (const VariableBlockVector& a)
+    {
+      if (&a!=this)     // check if this and a are different objects
+      {
+        // reallocate arrays if necessary
+        if (n!=a.n || nblocks!=a.nblocks)
+        {
+          if (n>0) A::template free<B>(p);
+          if (nblocks>0)
+          {
+            A::template free<int>(e);
+            A::template free<int>(s);
+          }
+          nblocks = a.nblocks;
+          if (nblocks>0)
+          {
+            s = A::template malloc<int>(nblocks);
+            e = A::template malloc<int>(nblocks);
+          }
+          else
+          {
+            nblocks = 0;
+            s = 0;
+            e = 0;
+          }
+          n = a.n;
+          if (n>0)
+            p = A::template malloc<B>(n);
+          else
+          {
+            n = 0;
+            p = 0;
+          }
+        }
+
+        // and copy the data
+        for (int i=0; i<nblocks; i++)
+        {
+          s[i]=a.s[i];
+          e[i]=a.e[i];
+        }
+        for (int i=0; i<n; i++) p[i]=a.p[i];
+      }
+
+      // and we have a usable vector
+      initialized = true;
+
+      return *this;     // Gebe Referenz zurueck damit a=b=c; klappt
+    }
+
+
+    //===== the creation interface
+
+    //! Iterator class for sequential creation of blocks
+    class CreateIterator
+    {
+    public:
+      //! constructor
+      CreateIterator (VariableBlockVector& _v, int _i) : v(_v)
+      {
+        i = _i;
+      }
+
+      //! prefix increment
+      CreateIterator& operator++()
+      {
+        ++i;
+        return *this;
+      }
+
+      //! inequality
+      bool operator!= (const CreateIterator& it) const
+      {
+        return (i!=it.i) || (&v!=&it.v);
+      }
+
+      //! equality
+      bool operator== (const CreateIterator& it) const
+      {
+        return (i==it.i) && (&v==&it.v);
+      }
+
+      //! dereferencing
+      int blockindex ()
+      {
+        return i;
+      }
+
+      //! set size of current block
+      void setblocksize (int k)
+      {
+        // do nothing if vector is already built (requires clear first)
+#ifdef DUNE_ISTL_WITH_CHECKING
+        if (v.initialized) DUNE_THROW(ISTLError,"tried to set block size for initialized vector");
+#endif
+
+        // set the blocks size
+        v.s[i] = k;
+
+        // compute entry into big array
+        if (i==0)
+          v.e[i] = 0;
+        else
+          v.e[i] = v.e[i-1]+v.s[i-1];
+
+        // if last block, finish off
+        if (i==v.nblocks-1)
+        {
+          // allocate p array
+          v.n = v.e[i]+v.s[i];
+          if (v.n>0)
+            v.p = A::template malloc<B>(v.n);
+          else
+          {
+            v.n = 0;
+            v.p = 0;
+          }
+
+          // and the vector is ready
+          v.initialized = true;
+        }
+      }
+
+    private:
+      VariableBlockVector& v;     // my vector
+      int i;                      // current block to be defined
+    };
+
+    // CreateIterator wants to set all the arrays ...
+    friend class CreateIterator;
+
+    //! get initial create iterator
+    CreateIterator createbegin ()
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (initialized) DUNE_THROW(ISTLError,"no CreateIterator in initialized state");
+#endif
+      return CreateIterator(*this,0);
+    }
+
+    //! get create iterator pointing to one after the last block
+    CreateIterator createend ()
+    {
+      return CreateIterator(*this,nblocks);
+    }
+
+    //===== assignment from scalar
+    VariableBlockVector& operator= (const field_type& k)
+    {
+      for (int i=0; i<n; i++)
+        p[i] = k;
+      return *this;
+    }
+
+    //===== access to components
+
+    //! random access to blocks
+    block_type& operator[] (int i)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=nblocks) DUNE_THROW(ISTLError,"index out of range");
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+#endif
+      window.set(s[i],p+e[i]);
+      return window;
+    }
+
+    //! same for read only access
+    const block_type& operator[] (int i) const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=nblocks) DUNE_THROW(ISTLError,"index out of range");
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+#endif
+      window.set(s[i],p+e[i]);
+      return window;
+    }
+
+    //! random access via round brackets supplied to be consistent with matrix
+    block_type& operator() (int i)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=nblocks) DUNE_THROW(ISTLError,"index out of range");
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+#endif
+      window.set(s[i],p+e[i]);
+      return window;
+    }
+
+    //! same for read only access
+    const block_type& operator() (int i) const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (i<0 || i>=nblocks) DUNE_THROW(ISTLError,"index out of range");
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+#endif
+      window.set(s[i],p+e[i]);
+      return window;
+    }
+
+    //! Iterator class for sequential access
+    class Iterator
+    {
+    public:
+      //! constructor
+      Iterator (VariableBlockVector& _v, int _i, B* _p) : v(_v),window(true)
+      {
+        i = _i;
+        p = _p;
+      }
+
+      //! prefix increment
+      Iterator& operator++()
+      {
+        p += v.s[i];
+        ++i;
+        window.set(s[i],p);
+        return *this;
+      }
+
+      //! equality
+      bool operator== (const Iterator& i) const
+      {
+        return p==i.p;
+      }
+
+      //! inequality
+      bool operator!= (const Iterator& i) const
+      {
+        return p!=i.p;
+      }
+
+      //! dereferencing
+      block_type& operator* ()
+      {
+        return window;
+      }
+
+      //! arrow
+      block_type* operator-> ()
+      {
+        return &window;
+      }
+
+    private:
+      mutable block_type window;     // provides a window into the vector
+      int i;
+      B* p;
+      VariableBlockVector& v;
+    };
+
+    // Iterator wants to see all the arrays ...
+    friend class Iterator;
+
+    //! begin iterator
+    Iterator begin () const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+#endif
+      return Iterator(*this,0,p);
+    }
+
+    //! end iterator
+    Iterator end () const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+#endif
+      return Iterator(*this,nblocks,p+n);
+    }
+
+
+    //===== vector space arithmetic
+
+    //! vector space addition
+    VariableBlockVector& operator+= (const VariableBlockVector& y)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+      if (n!=y.n) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+      for (int i=0; i<n; ++i) p[i] += y.p[i];
+      return *this;
+    }
+
+    //! vector space subtraction
+    VariableBlockVector& operator-= (const VariableBlockVector& y)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+      if (n!=y.n) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+      for (int i=0; i<n; ++i) p[i] -= y.p[i];
+      return *this;
+    }
+
+    //! vector space multiplication with scalar
+    VariableBlockVector& operator*= (const field_type& k)
+    {
+      for (int i=0; i<n; ++i) p[i] *= k;
+      return *this;
+    }
+
+    //! vector space division by scalar
+    VariableBlockVector& operator/= (const field_type& k)
+    {
+      for (int i=0; i<n; ++i) p[i] /= k;
+      return *this;
+    }
+
+    //! vector space axpy operation
+    VariableBlockVector& axpy (const field_type& a, const VariableBlockVector& y)
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+      if (n!=y.n) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+      for (int i=0; i<n; ++i) p[i].axpy(a,y.p[i]);
+      return *this;
+    }
+
+
+    //===== Euclidean scalar product
+
+    //! scalar product
+    field_type operator* (const VariableBlockVector& y) const
+    {
+#ifdef DUNE_ISTL_WITH_CHECKING
+      if (!initialized) DUNE_THROW(ISTLError,"tried to access uninitialized vector");
+      if (n!=y.n) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+      field_type sum=0;
+      for (int i=0; i<n; ++i) sum += p[i]*y.p[i];
+      return sum;
+    }
+
+
+    //===== norms
+
+    //! one norm (sum over absolute values of entries)
+    double one_norm () const
+    {
+      double sum=0;
+      for (int i=0; i<n; ++i) sum += p[i].one_norm();
+      return sum;
+    }
+
+    //! simplified one norm (uses Manhattan norm for complex values)
+    double one_norm_real () const
+    {
+      double sum=0;
+      for (int i=0; i<n; ++i) sum += p[i].one_norm_real();
+      return sum;
+    }
+
+    //! two norm sqrt(sum over squared values of entries)
+    double two_norm () const
+    {
+      double sum=0;
+      for (int i=0; i<n; ++i) sum += p[i].two_norm2();
+      return sqrt(sum);
+    }
+
+    //! sqare of two norm (sum over squared values of entries), need for block recursion
+    double two_norm2 () const
+    {
+      double sum=0;
+      for (int i=0; i<n; ++i) sum += p[i].two_norm2();
+      return sum;
+    }
+
+    //! infinity norm (maximum of absolute values of entries)
+    double infinity_norm () const
+    {
+      double max=0;
+      for (int i=0; i<n; ++i) max = std::max(max,p[i].infinity_norm());
+      return max;
+    }
+
+    //! simplified infinity norm (uses Manhattan norm for complex values)
+    double infinity_norm_real () const
+    {
+      double max=0;
+      for (int i=0; i<n; ++i) max = std::max(max,p[i].infinity_norm_real());
+      return max;
+    }
+
+
+    //===== sizes
+
+    //! number of blocks in the vector (are of size 1 here)
+    int N () const
+    {
+      return n;
+    }
+
+    //! dimension of the vector space
+    int dim () const
+    {
+      int d=0;
+      for (int i=0; i<n; i++)
+        d += p[i].dim();
+      return d;
+    }
+
+
+  private:
+    int nblocks;    // number of blocks in vector
+    int* s;         // s[nblocks] gives number of elements in each block
+    int* e;         // e[nblocks] says that p[e[i]] is the first element of block i
+    int n;          // number of elements in p array
+    B *p;           // pointer to dynamically allocated built-in array
+    mutable block_type window;     // provides a window into the vector
+    bool initialized;              // true if vector has been initialized
+  };
+
+
+
+  /** @} end documentation */
+
+} // end namespace
+
+#endif