From 1f7afc93624dc4b6c0a8ea8f8c4b2732f8ff439c Mon Sep 17 00:00:00 2001
From: Peter Bastian <peter@dune-project.org>
Date: Thu, 14 Oct 2004 13:40:33 +0000
Subject: [PATCH] restructuring. source files now in istl/
[[Imported from SVN: r46]]
---
istl/allocator.hh | 62 ++
istl/basearray.hh | 889 +++++++++++++++++++++++++++
istl/bcrsmatrix.hh | 1242 ++++++++++++++++++++++++++++++++++++++
istl/bvector.hh | 850 ++++++++++++++++++++++++++
istl/doc/istl.tex | 29 +-
istl/fmatrix.hh | 1227 +++++++++++++++++++++++++++++++++++++
istl/fvector.hh | 885 +++++++++++++++++++++++++++
istl/gsetc.hh | 191 ++++++
istl/io.hh | 190 ++++++
istl/istl.hh | 38 ++
istl/istlexception.hh | 25 +
istl/precision.hh | 69 +++
istl/tutorial/example.cc | 45 +-
istl/vbcrsmatrix.hh | 31 +
istl/vbvector.hh | 701 +++++++++++++++++++++
15 files changed, 6446 insertions(+), 28 deletions(-)
create mode 100644 istl/allocator.hh
create mode 100644 istl/basearray.hh
create mode 100644 istl/bcrsmatrix.hh
create mode 100644 istl/bvector.hh
create mode 100644 istl/fmatrix.hh
create mode 100644 istl/fvector.hh
create mode 100644 istl/gsetc.hh
create mode 100644 istl/io.hh
create mode 100644 istl/istl.hh
create mode 100644 istl/istlexception.hh
create mode 100644 istl/precision.hh
create mode 100644 istl/vbcrsmatrix.hh
create mode 100644 istl/vbvector.hh
diff --git a/istl/allocator.hh b/istl/allocator.hh
new file mode 100644
index 000000000..20f498e67
--- /dev/null
+++ b/istl/allocator.hh
@@ -0,0 +1,62 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_ALLOCATOR_HH__
+#define __DUNE_ALLOCATOR_HH__
+
+#include <stdlib.h>
+
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ /** The default allocator for the sparse matrix vector classes:
+
+ - uses malloc and free
+ - member templates for type safety to the outside
+ - is a singleton
+ - illustrates state handling through counter
+ - throws std::bad_alloc just as new does
+
+ */
+ class ISTLAllocator { // uses standard malloc and free
+ public:
+ //! allocate array of nmemb objects of type T
+ template<class T>
+ static T* malloc (size_t nmemb)
+ {
+ T* p = static_cast<T*>(std::malloc(nmemb*sizeof(T)));
+ if (p==NULL) throw std::bad_alloc();
+ count++;
+ return p;
+ }
+
+ //! release memory previously allocated with malloc member
+ template<class T>
+ static void free (T* p)
+ {
+ std::free(p);
+ count--;
+ }
+
+ //! return number of allocated objects
+ unsigned int nobjects ()
+ {
+ return count;
+ }
+
+ private:
+ // just to demonstrate some state information
+ static unsigned int count;
+ };
+
+ unsigned int ISTLAllocator::count=0;
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/basearray.hh b/istl/basearray.hh
new file mode 100644
index 000000000..48beac12c
--- /dev/null
+++ b/istl/basearray.hh
@@ -0,0 +1,889 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_BASEARRAY_HH__
+#define __DUNE_BASEARRAY_HH__
+
+#include <math.h>
+#include <complex>
+
+#include "istlexception.hh"
+#include "allocator.hh"
+
+/*! \file __FILE__
+
+ This file implements several basic array containers.
+ */
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ /** A simple array container for objects of type B implementing
+
+ - iterator access
+ - const_iterator access
+ - random access
+
+ This container has *NO* memory management at all,
+ copy constuctor, assignment and destructor are all default.
+
+ The constructor is made protected to emphasize that objects
+ are only usably in derived classes.
+
+ 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=ISTLAllocator>
+ class base_array_unmanaged
+ {
+ public:
+
+ //===== type definitions and constants
+
+ //! export the type representing the components
+ typedef B member_type;
+
+ //! export the allocator type
+ typedef A allocator_type;
+
+
+ //===== 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];
+ }
+
+ // forward declaration
+ class const_iterator;
+
+ //! iterator class for sequential access
+ class iterator
+ {
+ public:
+ //! constructor
+ iterator ()
+ {
+ p = 0;
+ i = 0;
+ }
+
+ iterator (B* _p, int _i) : p(_p), i(_i)
+ { }
+
+ //! prefix increment
+ iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! equality
+ bool operator== (const const_iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const const_iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ B& operator* ()
+ {
+ return p[i];
+ }
+
+ //! arrow
+ B* operator-> ()
+ {
+ return p+i;
+ }
+
+ // return index corresponding to pointer
+ int index () const
+ {
+ return i;
+ }
+
+ friend class const_iterator;
+
+ private:
+ B* p;
+ int i;
+ };
+
+ //! begin iterator
+ iterator begin ()
+ {
+ return iterator(p,0);
+ }
+
+ //! end iterator
+ iterator end ()
+ {
+ return iterator(p,n);
+ }
+
+ //! begin iterator
+ iterator rbegin ()
+ {
+ return iterator(p,n-1);
+ }
+
+ //! end iterator
+ iterator rend ()
+ {
+ return iterator(p,-1);
+ }
+
+ //! random access returning iterator (end if not contained)
+ iterator find (int i)
+ {
+ if (i>=0 && i<n)
+ return iterator(p,i);
+ else
+ return iterator(p,n);
+ }
+
+ //! const_iterator class for sequential access
+ class const_iterator
+ {
+ public:
+ //! constructor
+ const_iterator ()
+ {
+ p = 0;
+ i = 0;
+ }
+
+ const_iterator (const B* _p, int _i) : p(_p), i(_i)
+ { }
+
+ const_iterator (const iterator& it) : p(it.p), i(it.i)
+ { }
+
+ //! prefix increment
+ const_iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ const_iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const const_iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const const_iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! equality
+ bool operator== (const iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ const B& operator* () const
+ {
+ return p[i];
+ }
+
+ //! arrow
+ const B* operator-> () const
+ {
+ return p+i;
+ }
+
+ // return index corresponding to pointer
+ int index () const
+ {
+ return i;
+ }
+
+ friend class iterator;
+
+ private:
+ const B* p;
+ int i;
+ };
+
+ //! begin const_iterator
+ const_iterator begin () const
+ {
+ return const_iterator(p,0);
+ }
+
+ //! end const_iterator
+ const_iterator end () const
+ {
+ return const_iterator(p,n);
+ }
+
+ //! begin const_iterator
+ const_iterator rbegin () const
+ {
+ return const_iterator(p,n-1);
+ }
+
+ //! end const_iterator
+ const_iterator rend () const
+ {
+ return const_iterator(p,-1);
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in the array (are of size 1 here)
+ int size () const
+ {
+ return n;
+ }
+
+ protected:
+ //! makes empty array
+ base_array_unmanaged ()
+ {
+ n = 0;
+ p = 0;
+ }
+
+ int n; // number of elements in array
+ B *p; // pointer to dynamically allocated built-in array
+ };
+
+
+
+ /** Extend base_array_unmanaged by functions to manipulate
+
+ This container has *NO* memory management at all,
+ copy constuctor, assignment and destructor are all default.
+ A container can be constructed as empty or from a given pointer
+ and size. This class is used to implement a view into a larger
+ array.
+
+ You can copy such an object to a base_array to make a real copy.
+
+ 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=ISTLAllocator>
+ class base_array_window : public base_array_unmanaged<B,A>
+ {
+ public:
+
+ //===== type definitions and constants
+
+ //! export the type representing the components
+ typedef B member_type;
+
+ //! export the allocator type
+ typedef A allocator_type;
+
+ //! make iterators available as types
+ typedef typename base_array_unmanaged<B,A>::iterator iterator;
+
+ //! make iterators available as types
+ typedef typename base_array_unmanaged<B,A>::const_iterator const_iterator;
+
+
+ //===== constructors and such
+
+ //! makes empty array
+ base_array_window ()
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ //! make array from given pointer and size
+ base_array_window (B* _p, int _n)
+ {
+ this->n = _n;
+ this->p = _p;
+ }
+
+ //===== window manipulation methods
+
+ //! set pointer and length
+ void set (int _n, B* _p)
+ {
+ this->n = _n;
+ this->p = _p;
+ }
+
+ //! advance pointer by newsize elements and then set size to new size
+ void advance (int newsize)
+ {
+ this->p += this->n;
+ this->n = newsize;
+ }
+
+ //! increment pointer by offset and set size
+ void move (int offset, int newsize)
+ {
+ this->p += offset;
+ this->n = newsize;
+ }
+
+ //! increment pointer by offset, leave size
+ void move (int offset)
+ {
+ this->p += offset;
+ }
+
+ //! return the pointer
+ B* getptr ()
+ {
+ return this->p;
+ }
+ };
+
+
+
+ /** This container extends base_array_unmanaged by memory management
+ with the usual copy semantics providing the full range of
+ copy constructor, destructor and assignement operators.
+
+ You can make
+
+ - empty array
+ - array with n components dynamically allocated
+ - resize an array with complete loss of data
+ - assign/construct from a base_array_window to make a copy of the data
+
+ 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=ISTLAllocator>
+ class base_array : public base_array_unmanaged<B,A>
+ {
+ public:
+
+ //===== type definitions and constants
+
+ //! export the type representing the components
+ typedef B member_type;
+
+ //! export the allocator type
+ typedef A allocator_type;
+
+ //! make iterators available as types
+ typedef typename base_array_unmanaged<B,A>::iterator iterator;
+
+ //! make iterators available as types
+ typedef typename base_array_unmanaged<B,A>::const_iterator const_iterator;
+
+ //===== constructors and such
+
+ //! makes empty array
+ base_array ()
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ //! make array with _n components
+ base_array (int _n)
+ {
+ this->n = _n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+ }
+
+ //! copy constructor
+ base_array (const base_array& a)
+ {
+ // allocate memory with same size as a
+ this->n = a.n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ // and copy elements
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+
+ //! construct from base class object
+ base_array (const base_array_unmanaged<B,A>& _a)
+ {
+ const base_array& a = static_cast<const base_array&>(_a);
+
+ // allocate memory with same size as a
+ this->n = a.n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ // and copy elements
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+
+
+ //! free dynamic memory
+ ~base_array ()
+ {
+ if (this->n>0) A::template free<B>(this->p);
+ }
+
+ //! reallocate array to given size, any data is lost
+ void resize (int _n)
+ {
+ if (this->n==_n) return;
+
+ if (this->n>0) A::template free<B>(this->p);
+ this->n = _n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+ }
+
+ //! assignment
+ base_array& operator= (const base_array& a)
+ {
+ if (&a!=this) // check if this and a are different objects
+ {
+ // adjust size of array
+ if (this->n!=a.n) // check if size is different
+ {
+ if (this->n>0) A::template free<B>(this->p); // delete old memory
+ this->n = a.n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+ }
+ // copy data
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+ return *this;
+ }
+
+ //! assign from base class object
+ base_array& operator= (const base_array_unmanaged<B,A>& a)
+ {
+ return this->operator=(static_cast<const base_array&>(a));
+ }
+
+ };
+
+
+
+
+ /** A simple array container with non-consecutive index set.
+ Elements in the array are of type B. This class provides
+
+ - iterator access
+ - const_iterator access
+ - random access in log(n) steps using binary search
+ - find returning iterator
+
+ This container has *NO* memory management at all,
+ copy constuctor, assignment and destructor are all default.
+
+ The constructor is made protected to emphasize that objects
+ are only usably in derived classes.
+
+ 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=ISTLAllocator>
+ class compressed_base_array_unmanaged
+ {
+ public:
+
+ //===== type definitions and constants
+
+ //! export the type representing the components
+ typedef B member_type;
+
+ //! export the allocator type
+ typedef A allocator_type;
+
+
+ //===== access to components
+
+ //! random access to blocks, assumes ascending ordering
+ B& operator[] (int i)
+ {
+ int l=0, r=n-1;
+ while (l<r)
+ {
+ int q = (l+r)/2;
+ if (i <= j[q]) r=q;
+ else l = q+1;
+ }
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (j[l]!=i) DUNE_THROW(ISTLError,"index not in compressed array");
+#endif
+ return p[l];
+ }
+
+ //! same for read only access, assumes ascending ordering
+ const B& operator[] (int i) const
+ {
+ int l=0, r=n-1;
+ while (l<r)
+ {
+ int q = (l+r)/2;
+ if (i <= j[q]) r=q;
+ else l = q+1;
+ }
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (j[l]!=i) DUNE_THROW(ISTLError,"index not in compressed array");
+#endif
+ return p[l];
+ }
+
+ // forward declaration
+ class const_iterator;
+
+ //! iterator class for sequential access
+ class iterator
+ {
+ public:
+ //! constructor
+ iterator ()
+ {
+ p = 0;
+ j = 0;
+ i = 0;
+ }
+
+ iterator (B* _p, int* _j, int _i) : p(_p), j(_j), i(_i)
+ { }
+
+ //! prefix increment
+ iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const iterator& it) const
+ {
+ // return (p+i)==(it.p+it.i);
+ return (i)==(it.i);
+ }
+
+ //! inequality
+ bool operator!= (const iterator& it) const
+ {
+ // return (p+i)!=(it.p+it.i);
+ return (i)!=(it.i);
+ }
+
+ //! equality
+ bool operator== (const const_iterator& it) const
+ {
+ // return (p+i)==(it.p+it.i);
+ return (i)==(it.i);
+ }
+
+ //! inequality
+ bool operator!= (const const_iterator& it) const
+ {
+ // return (p+i)!=(it.p+it.i);
+ return (i)!=(it.i);
+ }
+
+ //! dereferencing
+ B& operator* ()
+ {
+ return p[i];
+ }
+
+ //! arrow
+ B* operator-> ()
+ {
+ return p+i;
+ }
+
+ // return index corresponding to pointer
+ int index () const
+ {
+ return j[i];
+ }
+
+ // return index corresponding to pointer
+ void setindex (int k)
+ {
+ return j[i] = k;
+ }
+
+ friend class const_iterator;
+
+ private:
+ B* p;
+ int* j;
+ int i;
+ };
+
+ //! begin iterator
+ iterator begin ()
+ {
+ return iterator(p,j,0);
+ }
+
+ //! end iterator
+ iterator end ()
+ {
+ return iterator(p,j,n);
+ }
+
+ //! begin iterator
+ iterator rbegin ()
+ {
+ return iterator(p,j,n-1);
+ }
+
+ //! end iterator
+ iterator rend ()
+ {
+ return iterator(p,j,-1);
+ }
+
+ //! const_iterator class for sequential access
+ class const_iterator
+ {
+ public:
+ //! constructor
+ const_iterator ()
+ {
+ p = 0;
+ j = 0;
+ i = 0;
+ }
+
+ const_iterator (const B* _p, const int* _j, int _i) : p(_p), j(_j), i(_i)
+ { }
+
+ const_iterator (const iterator& it) : p(it.p), j(it.j), i(it.i)
+ { }
+
+ //! prefix increment
+ const_iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ const_iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const const_iterator& it) const
+ {
+ // return (p+i)==(it.p+it.i);
+ return (i)==(it.i);
+ }
+
+ //! inequality
+ bool operator!= (const const_iterator& it) const
+ {
+ // return (p+i)!=(it.p+it.i);
+ return (i)!=(it.i);
+ }
+
+ //! equality
+ bool operator== (const iterator& it) const
+ {
+ // return (p+i)==(it.p+it.i);
+ return (i)==(it.i);
+ }
+
+ //! inequality
+ bool operator!= (const iterator& it) const
+ {
+ // return (p+i)!=(it.p+it.i);
+ return (i)!=(it.i);
+ }
+
+ //! dereferencing
+ const B& operator* () const
+ {
+ return p[i];
+ }
+
+ //! arrow
+ const B* operator-> () const
+ {
+ return p+i;
+ }
+
+ // return index corresponding to pointer
+ int index () const
+ {
+ return j[i];
+ }
+
+ friend class iterator;
+
+ private:
+ const B* p;
+ const int* j;
+ int i;
+ };
+
+ //! begin const_iterator
+ const_iterator begin () const
+ {
+ return const_iterator(p,j,0);
+ }
+
+ //! end const_iterator
+ const_iterator end () const
+ {
+ return const_iterator(p,j,n);
+ }
+
+ //! begin const_iterator
+ const_iterator rbegin () const
+ {
+ return const_iterator(p,j,n-1);
+ }
+
+ //! end const_iterator
+ const_iterator rend () const
+ {
+ return const_iterator(p,j,-1);
+ }
+
+ //! random access returning iterator (end if not contained)
+ iterator find (int i)
+ {
+ int l=0, r=n-1;
+ while (l<r)
+ {
+ int q = (l+r)/2;
+ if (i <= j[q]) r=q;
+ else l = q+1;
+ }
+
+ if (i==j[l])
+ return iterator(p,j,l);
+ else
+ return iterator(p,j,n);
+ }
+
+ //===== sizes
+
+ //! number of blocks in the array (are of size 1 here)
+ int size () const
+ {
+ return n;
+ }
+
+ protected:
+ //! makes empty array
+ compressed_base_array_unmanaged ()
+ {
+ n = 0;
+ p = 0;
+ j = 0;
+ }
+
+ int n; // number of elements in array
+ B *p; // pointer to dynamically allocated built-in array
+ int* j; // the index set
+ };
+
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/bcrsmatrix.hh b/istl/bcrsmatrix.hh
new file mode 100644
index 000000000..2334e508c
--- /dev/null
+++ b/istl/bcrsmatrix.hh
@@ -0,0 +1,1242 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_BCRSMATRIX_HH__
+#define __DUNE_BCRSMATRIX_HH__
+
+#include <math.h>
+#include <complex>
+#include <set>
+
+#include "istlexception.hh"
+#include "allocator.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 ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+
+ /**! implements a block compressed row storage scheme. The block
+ type B can be any type implementing the matrix interface.
+
+ Different ways to build up a compressed row
+ storage matrix are supported:
+
+ 1. Row-wise scheme
+
+ Rows are built up in sequential order. Size of the row and
+ the column indices are defined. A row can be used as soon as it
+ is initialized. With respect to memory there are two variants of
+ this scheme: (a) number of non-zeroes known in advance (application
+ finite difference schemes), (b) number of non-zeroes not known
+ in advance (application: Sparse LU, ILU(n)).
+
+ 2. Random scheme
+
+ For general finite element implementations the number of rows n
+ is known, the number of non-zeroes might also be known (e.g.
+ #edges + #nodes for P1) but the size of a row and the indices of a row
+ can not be defined in sequential order.
+
+ 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=ISTLAllocator>
+ class BCRSMatrix
+ {
+ 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;
+
+ //! implement row_type with compressed vector
+ typedef CompressedBlockVectorWindow<B,A> row_type;
+
+ //! increment block level counter
+ enum {blocklevel = B::blocklevel+1};
+
+ //! we support two modes
+ enum BuildMode {row_wise, random, unknown};
+
+
+ //===== random access interface to rows of the matrix
+
+ //! random access to the rows
+ row_type& operator[] (int i)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (!ready) DUNE_THROW(ISTLError,"row not initialized yet");
+ if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ return r[i];
+ }
+
+ //! same for read only access
+ const row_type& operator[] (int i) const
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (!ready) DUNE_THROW(ISTLError,"row not initialized yet");
+ if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ return r[i];
+ }
+
+
+ //===== iterator interface to rows of the matrix
+
+ // forward declaration
+ class ConstIterator;
+
+ //! Iterator access to rows
+ class Iterator
+ {
+ public:
+ //! constructor
+ Iterator (row_type* _p, int _i)
+ {
+ p = _p;
+ i = _i;
+ }
+
+ //! empty constructor, use with care!
+ Iterator ()
+ { }
+
+ //! prefix increment
+ Iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ Iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ // return (p+i)==(it.p+it.i);
+ return (i)==(it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ // return (p+i)!=(it.p+it.i);
+ return (i)!=(it.i);
+ }
+
+ //! less than
+ bool operator< (const Iterator& it) const
+ {
+ return (i)<(it.i);
+ }
+
+ //! dereferencing
+ row_type& operator* ()
+ {
+ return p[i];
+ }
+
+ //! arrow
+ row_type* operator-> ()
+ {
+ return p+i;
+ }
+
+ //! return index
+ int index ()
+ {
+ return i;
+ }
+
+ friend class ConstIterator;
+
+ private:
+ row_type* p;
+ int i;
+ };
+
+ //! begin iterator
+ Iterator begin ()
+ {
+ return Iterator(r,0);
+ }
+
+ //! end iterator
+ Iterator end ()
+ {
+ return Iterator(r,n);
+ }
+
+ //! begin iterator
+ Iterator rbegin ()
+ {
+ return Iterator(r,n-1);
+ }
+
+ //! end iterator
+ Iterator rend ()
+ {
+ return Iterator(r,-1);
+ }
+
+ //! rename the iterators for easier access (no const iterators?)
+ typedef Iterator RowIterator;
+ typedef typename row_type::Iterator ColIterator;
+
+
+ //! Iterator access to rows
+ class ConstIterator
+ {
+ public:
+ //! constructor
+ ConstIterator (const row_type* _p, int _i) : p(_p), i(_i)
+ { }
+
+ //! empty constructor, use with care!
+ ConstIterator ()
+ {
+ p = 0;
+ i = 0;
+ }
+
+ //! prefix increment
+ ConstIterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ ConstIterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const ConstIterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const ConstIterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ const row_type& operator* () const
+ {
+ return p[i];
+ }
+
+ //! arrow
+ const row_type* operator-> () const
+ {
+ return p+i;
+ }
+
+ //! return index
+ int index () const
+ {
+ return i;
+ }
+
+ friend class Iterator;
+
+ private:
+ const row_type* p;
+ int i;
+ };
+
+ //! begin iterator
+ ConstIterator begin () const
+ {
+ return ConstIterator(r,0);
+ }
+
+ //! end iterator
+ ConstIterator end () const
+ {
+ return ConstIterator(r,n);
+ }
+
+ //! begin iterator
+ ConstIterator rbegin () const
+ {
+ return ConstIterator(r,n-1);
+ }
+
+ //! end iterator
+ ConstIterator rend () const
+ {
+ return ConstIterator(r,-1);
+ }
+
+ //! rename the iterators for easier access
+ typedef ConstIterator ConstRowIterator;
+ typedef typename row_type::ConstIterator ConstColIterator;
+
+ //===== constructors & resizers
+
+ //! an empty matrix
+ BCRSMatrix ()
+ {
+ // the state
+ build_mode = unknown;
+ ready = false;
+
+ // store parameters
+ n = 0;
+ m = 0;
+ nnz = 0;
+
+ // allocate rows
+ r = 0;
+
+ // allocate a and j array
+ a = 0;
+ j = 0;
+ }
+
+ //! matrix with known number of nonzeroes
+ BCRSMatrix (int _n, int _m, int _nnz, BuildMode bm)
+ {
+ // the state
+ build_mode = bm;
+ ready = false;
+
+ // store parameters
+ n = _n;
+ m = _m;
+ nnz = _nnz;
+
+ // allocate rows
+ if (n>0)
+ {
+ r = A::template malloc<row_type>(n);
+ }
+ else
+ {
+ r = 0;
+ }
+
+ // allocate a and j array
+ if (nnz>0)
+ {
+ a = A::template malloc<B>(nnz);
+ j = A::template malloc<int>(nnz);
+ }
+ else
+ {
+ a = 0;
+ j = 0;
+ }
+ }
+
+ //! matrix with unknown number of nonzeroes
+ BCRSMatrix (int _n, int _m, BuildMode bm)
+ {
+ // the state
+ build_mode = bm;
+ ready = false;
+
+ // store parameters
+ n = _n;
+ m = _m;
+ nnz = 0;
+
+ // allocate rows
+ if (n>0)
+ {
+ r = A::template malloc<row_type>(n);
+ }
+ else
+ {
+ r = 0;
+ }
+
+ // allocate a and j array
+ a = 0;
+ j = 0;
+ }
+
+ //! copy constructor
+ BCRSMatrix (const BCRSMatrix& Mat)
+ {
+ // deep copy in global array
+
+ // copy sizes
+ n = Mat.n;
+ m = Mat.m;
+ nnz = Mat.nnz;
+
+ // in case of row-wise allocation
+ if (nnz<=0)
+ {
+ nnz = 0;
+ for (int i=0; i<n; i++)
+ nnz += Mat.r[i].getsize();
+ }
+
+ // allocate rows
+ if (n>0)
+ {
+ r = A::template malloc<row_type>(n);
+ }
+ else
+ {
+ r = 0;
+ }
+
+ // allocate a and j array
+ if (nnz>0)
+ {
+ a = A::template malloc<B>(nnz);
+ j = A::template malloc<int>(nnz);
+ }
+ else
+ {
+ a = 0;
+ j = 0;
+ }
+
+ // build window structure
+ for (int i=0; i<n; i++)
+ {
+ // set row i
+ r[i].setsize(Mat.r[i].getsize());
+ if (i==0)
+ {
+ r[i].setptr(a);
+ r[i].setindexptr(j);
+ }
+ else
+ {
+ r[i].setptr( r[i-1].getptr()+r[i-1].getsize() );
+ r[i].setindexptr( r[i-1].getindexptr()+r[i-1].getsize() );
+ }
+ }
+
+ // copy data
+ for (int i=0; i<n; i++) r[i] = Mat.r[i];
+
+ // finish off
+ build_mode = row_wise; // dummy
+ ready = true;
+ }
+
+ //! destructor
+ ~BCRSMatrix ()
+ {
+ if (nnz>0)
+ {
+ // a,j have been allocated as one long vector
+ A::template free<int>(j);
+ A::template free<B>(a);
+ }
+ else
+ {
+ // check if memory for rows have been allocated individually
+ for (int i=0; i<n; i++)
+ if (r[i].getsize()>0)
+ {
+ A::template free<int>(r[i].getindexptr());
+ A::template free<B>(r[i].getptr());
+ }
+ }
+
+ // deallocate the rows
+ if (n>0) A::template free<row_type>(r);
+ }
+
+ //! assignment
+ BCRSMatrix& operator= (const BCRSMatrix& Mat)
+ {
+ // return immediately when self-assignment
+ if (&Mat==this) return *this;
+
+ // make it simple: ALWAYS throw away memory for a and j
+ if (nnz>0)
+ {
+ // a,j have been allocated as one long vector
+ A::template free<int>(j);
+ A::template free<B>(a);
+ }
+ else
+ {
+ // check if memory for rows have been allocated individually
+ for (int i=0; i<n; i++)
+ if (r[i].getsize()>0)
+ {
+ A::template free<int>(r[i].getindexptr());
+ A::template free<B>(r[i].getptr());
+ }
+ }
+
+ // reallocate the rows if required
+ if (n!=Mat.n)
+ {
+ // free rows
+ A::template free<row_type>(r);
+
+ // allocate rows
+ n = Mat.n;
+ if (n>0)
+ {
+ r = A::template malloc<row_type>(n);
+ }
+ else
+ {
+ r = 0;
+ }
+ }
+
+ // allocate a,j
+ m = Mat.m;
+ nnz = Mat.nnz;
+
+ // in case of row-wise allocation
+ if (nnz<=0)
+ {
+ nnz = 0;
+ for (int i=0; i<n; i++)
+ nnz += Mat.r[i].getsize();
+ }
+
+ // allocate a and j array
+ if (nnz>0)
+ {
+ a = A::template malloc<B>(nnz);
+ j = A::template malloc<int>(nnz);
+ }
+ else
+ {
+ a = 0;
+ j = 0;
+ }
+
+ // build window structure
+ for (int i=0; i<n; i++)
+ {
+ // set row i
+ r[i].setsize(Mat.r[i].getsize());
+ if (i==0)
+ {
+ r[i].setptr(a);
+ r[i].setindexptr(j);
+ }
+ else
+ {
+ r[i].setptr( r[i-1].getptr()+r[i-1].getsize() );
+ r[i].setindexptr( r[i-1].getindexptr()+r[i-1].getsize() );
+ }
+ }
+
+ // copy data
+ for (int i=0; i<n; i++) r[i] = Mat.r[i];
+
+ // finish off
+ build_mode = row_wise; // dummy
+ ready = true;
+ }
+
+ BCRSMatrix& operator= (const field_type& k)
+ {
+ for (int i=0; i<n; i++) r[i] = k;
+ }
+
+ //===== row-wise creation interface
+
+ //! Iterator class for sequential creation of blocks
+ class CreateIterator
+ {
+ public:
+ //! constructor
+ CreateIterator (BCRSMatrix& _Mat, int _i) : Mat(_Mat)
+ {
+ if (Mat.build_mode!=row_wise)
+ DUNE_THROW(ISTLError,"creation only allowed for uninitialized matrix");
+
+ i = _i;
+ nnz = 0;
+ }
+
+ //! prefix increment
+ CreateIterator& operator++()
+ {
+ // this should only be called if matrix is in creation
+ if (Mat.ready)
+ DUNE_THROW(ISTLError,"matrix already built up");
+
+ // row i is defined through the pattern
+ // get memory for the row and initialize the j array
+ // this depends on the allocation mode
+
+ // compute size of the row
+ int s = pattern.size();
+
+ // update number of nonzeroes including this row
+ nnz += s;
+
+ // alloc memory / set window
+ if (Mat.nnz>0)
+ {
+ // memory is allocated in one long array
+
+ // check if that memory is sufficient
+ if (nnz>Mat.nnz)
+ DUNE_THROW(ISTLError,"allocated nnz too small");
+
+ // set row i
+ if (i==0)
+ Mat.r[i].set(s,Mat.a,Mat.j);
+ else
+ Mat.r[i].set(s,Mat.r[i-1].getptr()+Mat.r[i-1].getsize(),
+ Mat.r[i-1].getindexptr()+Mat.r[i-1].getsize());
+ }
+ else
+ {
+ // memory is allocated individually per row
+ // allocate and set row i
+ if (s>0)
+ {
+ B* a = A::template malloc<B>(s);
+ int* j = A::template malloc<int>(s);
+ Mat.r[i].set(s,a,j);
+ }
+ else
+ Mat.r[i].set(0,0,0);
+ }
+
+ // initialize the j array for row i from pattern
+ int k=0;
+ int *j = Mat.r[i].getindexptr();
+ for (std::set<int>::const_iterator it=pattern.begin(); it!=pattern.end(); ++it)
+ j[k++] = *it;
+
+ // now go to next row
+ i++;
+ pattern.clear();
+
+ // check if this was last row
+ if (i==Mat.n)
+ {
+ Mat.ready = true;
+ }
+ // done
+ return *this;
+ }
+
+ //! inequality
+ bool operator!= (const CreateIterator& it) const
+ {
+ return (i!=it.i) || (&Mat!=&it.Mat);
+ }
+
+ //! equality
+ bool operator== (const CreateIterator& it) const
+ {
+ return (i==it.i) && (&Mat==&it.Mat);
+ }
+
+ //! dereferencing
+ int index ()
+ {
+ return i;
+ }
+
+ //! put column index in row
+ void insert (int j)
+ {
+ pattern.insert(j);
+ }
+
+ //! return true if column index is in row
+ bool contains (int j)
+ {
+ if (pattern.find(j)!=pattern.end())
+ return true;
+ else
+ return false;
+ }
+
+ private:
+ BCRSMatrix& Mat; // the matrix we are defining
+ int i; // current row to be defined
+ int nnz; // count total number of nonzeros
+ std::set<int> pattern; // used to compile entries in a row
+ };
+
+
+ //! allow CreateIterator to access internal data
+ friend class CreateIterator;
+
+ //! get initial create iterator
+ CreateIterator createbegin ()
+ {
+ return CreateIterator(*this,0);
+ }
+
+ //! get create iterator pointing to one after the last block
+ CreateIterator createend ()
+ {
+ return CreateIterator(*this,n);
+ }
+
+
+ //===== random creation interface
+
+ //! set number of indices in row i to s
+ void setrowsize (int i, int s)
+ {
+ if (build_mode!=random)
+ DUNE_THROW(ISTLError,"requires random build mode");
+ if (ready)
+ DUNE_THROW(ISTLError,"matrix already built up");
+
+ r[i].setsize(s);
+ }
+
+ //! increment size of row i by 1
+ void incrementrowsize (int i)
+ {
+ if (build_mode!=random)
+ DUNE_THROW(ISTLError,"requires random build mode");
+ if (ready)
+ DUNE_THROW(ISTLError,"matrix already built up");
+
+ r[i].setsize(r[i].getsize()+1);
+ }
+
+ //! indicate that size of all rows is defined
+ void endrowsizes ()
+ {
+ if (build_mode!=random)
+ DUNE_THROW(ISTLError,"requires random build mode");
+ if (ready)
+ DUNE_THROW(ISTLError,"matrix already built up");
+
+ // compute total size, check positivity
+ int total=0;
+ for (int i=0; i<n; i++)
+ {
+ if (r[i].getsize()<=0)
+ DUNE_THROW(ISTLError,"rowsize must be positive");
+ total += r[i].getsize();
+ }
+
+ // allocate/check memory
+ if (nnz>0)
+ {
+ // there is already memory
+ if (total>nnz)
+ DUNE_THROW(ISTLError,"nnz too small");
+ }
+ else
+ {
+ // we allocate the memory here
+ nnz = total;
+ if (nnz>0)
+ {
+ a = A::template malloc<B>(nnz);
+ j = A::template malloc<int>(nnz);
+ }
+ else
+ {
+ a = 0;
+ j = 0;
+ }
+ }
+
+ // set the window pointers correctly
+ for (int i=0; i<n; i++)
+ {
+ // set row i
+ if (i==0)
+ {
+ r[i].setptr(a);
+ r[i].setindexptr(j);
+ }
+ else
+ {
+ r[i].setptr( r[i-1].getptr()+r[i-1].getsize() );
+ r[i].setindexptr( r[i-1].getindexptr()+r[i-1].getsize() );
+ }
+ }
+
+ // initialize j array with m (an invalid column index)
+ // this indicates an unused entry
+ for (int k=0; k<nnz; k++)
+ j[k] = m;
+ }
+
+ //! add index (row,col) to the matrix
+ void addindex (int row, int col)
+ {
+ if (build_mode!=random)
+ DUNE_THROW(ISTLError,"requires random build mode");
+ if (ready)
+ DUNE_THROW(ISTLError,"matrix already built up");
+
+ // get row
+ int* p = r[row].getindexptr();
+ int s = r[row].getsize();
+
+ // binary search for col
+ int l=0, r=s-1;
+ while (l<r)
+ {
+ int q = (l+r)/2;
+ if (col <= p[q]) r=q;
+ else l = q+1;
+ }
+
+ // check if index is already in row
+ if (p[l]==col) return;
+
+ // no, find first free entry by binary search
+ l=0; r=s-1;
+ while (l<r)
+ {
+ int q = (l+r)/2;
+ if (m <= p[q]) r=q;
+ else l = q+1;
+ }
+ if (p[l]!=m)
+ DUNE_THROW(ISTLError,"row is too small");
+
+ // now p[l] is the leftmost free entry
+ p[l] = col;
+
+ // use insertion sort to move index to correct position
+ for (int i=l-1; i>=0; i--)
+ if (p[i]>p[i+1])
+ std::swap(p[i],p[i+1]);
+ else
+ break;
+ }
+
+ //! indicate that all indices are defined, check consistency
+ void endindices ()
+ {
+ if (build_mode!=random)
+ DUNE_THROW(ISTLError,"requires random build mode");
+ if (ready)
+ DUNE_THROW(ISTLError,"matrix already built up");
+
+ // check if there are undefined indices
+ for (int k=0; k<nnz; k++)
+ if (j[k]<0 || j[k]>=m)
+ DUNE_THROW(ISTLError,"undefined index detected");
+
+ // if not, set matrix to ready
+ ready = true;
+ }
+
+ //===== vector space arithmetic
+
+ //! vector space multiplication with scalar
+ BCRSMatrix& operator*= (const field_type& k)
+ {
+ if (nnz>0)
+ {
+ // process 1D array
+ for (int i=0; i<nnz; i++)
+ a[i] *= k;
+ }
+ else
+ {
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j) *= k;
+ }
+ }
+
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ BCRSMatrix& operator/= (const field_type& k)
+ {
+ if (nnz>0)
+ {
+ // process 1D array
+ for (int i=0; i<nnz; i++)
+ a[i] /= k;
+ }
+ else
+ {
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j) /= k;
+ }
+ }
+
+ return *this;
+ }
+
+
+ //===== linear maps
+
+ //! y += A x
+ template<class X, class Y>
+ void umv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).umv(x[j.index()],y[i.index()]);
+ }
+ }
+
+ //! y -= A x
+ template<class X, class Y>
+ void mmv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).mmv(x[j.index()],y[i.index()]);
+ }
+ }
+
+ //! y += alpha A x
+ template<class X, class Y>
+ void usmv (const field_type& alpha, X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).usmv(alpha,x[j.index()],y[i.index()]);
+ }
+ }
+
+ //! y += A^T x
+ template<class X, class Y>
+ void umtv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).umtv(x[i.index()],y[j.index()]);
+ }
+ }
+
+ //! y -= A^T x
+ template<class X, class Y>
+ void mmtv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).mmtv(x[i.index()],y[j.index()]);
+ }
+ }
+
+ //! y += alpha A^T x
+ template<class X, class Y>
+ void usmtv (const field_type& alpha, X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).usmtv(alpha,x[i.index()],y[j.index()]);
+ }
+ }
+
+ //! y += A^H x
+ template<class X, class Y>
+ void umhv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).umhv(x[i.index()],y[j.index()]);
+ }
+ }
+
+ //! y -= A^H x
+ template<class X, class Y>
+ void mmhv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).mmhv(x[i.index()],y[j.index()]);
+ }
+ }
+
+ //! y += alpha A^H x
+ template<class X, class Y>
+ void usmhv (const field_type& alpha, X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ (*j).usmhv(alpha,x[i.index()],y[j.index()]);
+ }
+ }
+
+
+ //===== norms
+
+ //! square of frobenius norm, need for block recursion
+ double frobenius_norm2 () const
+ {
+ double sum=0;
+
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ sum += (*j).frobenius_norm2();
+ }
+
+ return sum;
+ }
+
+ //! frobenius norm: sqrt(sum over squared values of entries)
+ double frobenius_norm () const
+ {
+ return sqrt(frobenius_norm2());
+ }
+
+ //! infinity norm (row sum norm, how to generalize for blocks?)
+ double infinity_norm () const
+ {
+ double max=0;
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ double sum=0;
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ sum += (*j).infinity_norm();
+ max = std::max(max,sum);
+ }
+ return max;
+ }
+
+ //! simplified infinity norm (uses Manhattan norm for complex values)
+ double infinity_norm_real () const
+ {
+ double max=0;
+ RowIterator endi=end();
+ for (RowIterator i=begin(); i!=endi; ++i)
+ {
+ double sum=0;
+ ColIterator endj = (*i).end();
+ for (ColIterator j=(*i).begin(); j!=endj; ++j)
+ sum += (*j).infinity_norm_real();
+ max = std::max(max,sum);
+ }
+ return max;
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in row direction
+ int N () const
+ {
+ return n;
+ }
+
+ //! number of blocks in column direction
+ int M () const
+ {
+ return m;
+ }
+
+ //! row dimension of block r
+ int rowdim (int i) const
+ {
+ return r[i].getptr()->rowdim();
+ }
+
+ //! col dimension of block c
+ int coldim (int c) const
+ {
+ // find an entry in column j
+ B* entry=0;
+ if (nnz>0)
+ {
+ for (int k=0; k<nnz; k++)
+ if (j[k]==c) {
+ return a[k].coldim();
+ }
+ }
+ else
+ {
+ for (int i=0; i<n; i++)
+ {
+ int* j = r[i].getindexptr();
+ B* a = r[i].getptr();
+ for (int k=0; k<r[i].getsize(); k++)
+ if (j[k]==c) {
+ return a[k].coldim();
+ }
+ }
+ }
+
+ // not found
+ return 0;
+ }
+
+ //! dimension of the destination vector space
+ int rowdim () const
+ {
+ int nn=0;
+ for (int i=0; i<n; i++)
+ nn += rowdim(i);
+ return nn;
+ }
+
+ //! dimension of the source vector space
+ int coldim () const
+ {
+ int mm=0;
+ for (int i=0; i<m; i++)
+ mm += coldim(i);
+ return mm;
+ }
+
+ //===== query
+
+ // return true when (i,j) is in pattern
+ bool exists (int i, int j)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+ if (j<0 || i>=m) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ if (r[i].find(j)!=r[i].end())
+ return true;
+ else
+ return false;
+ }
+
+
+ private:
+ // state information
+ BuildMode build_mode; // row wise or whole matrix
+ bool ready; // true if matrix is ready to use
+
+ // size of the matrix
+ int n; // number of rows
+ int m; // number of columns
+ int nnz; // number of nonzeros allocated in the a and j array below
+ // zero means that memory is allocated seperately for each row.
+
+ // the rows are dynamically allocated
+ row_type* r; // [n] the individual rows having pointers into a,j arrays
+
+ // dynamically allocated memory
+ B* a; // [nnz] non-zero entries of the matrix in row-wise ordering
+ int* j; // [nnz] column indices of entries
+ };
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/bvector.hh b/istl/bvector.hh
new file mode 100644
index 000000000..7ec4cb698
--- /dev/null
+++ b/istl/bvector.hh
@@ -0,0 +1,850 @@
+// -*- 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 "istlexception.hh"
+#include "allocator.hh"
+#include "basearray.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 ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+
+ /** block_vector_unmanaged extends the base_array_unmanaged by
+ vector operations such as addition and scalar multiplication.
+ No memory management is added.
+
+ 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=ISTLAllocator>
+ class block_vector_unmanaged : public base_array_unmanaged<B,A>
+ {
+ 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;
+
+ //! make iterators available as types
+ typedef typename base_array_unmanaged<B,A>::iterator Iterator;
+
+ //! make iterators available as types
+ typedef typename base_array_unmanaged<B,A>::const_iterator ConstIterator;
+
+
+ //===== assignment from scalar
+
+ block_vector_unmanaged& operator= (const field_type& k)
+ {
+ for (int i=0; i<this->n; i++)
+ (*this)[i] = k;
+ return *this;
+ }
+
+
+ //===== vector space arithmetic
+
+ //! vector space addition
+ block_vector_unmanaged& operator+= (const block_vector_unmanaged& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+ for (int i=0; i<this->n; ++i) (*this)[i] += y[i];
+ return *this;
+ }
+
+ //! vector space subtraction
+ block_vector_unmanaged& operator-= (const block_vector_unmanaged& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+ for (int i=0; i<this->n; ++i) (*this)[i] -= y[i];
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ block_vector_unmanaged& operator*= (const field_type& k)
+ {
+ for (int i=0; i<this->n; ++i) (*this)[i] *= k;
+ return *this;
+ }
+
+ //! vector space division by scalar
+ block_vector_unmanaged& operator/= (const field_type& k)
+ {
+ for (int i=0; i<this->n; ++i) (*this)[i] /= k;
+ return *this;
+ }
+
+ //! vector space axpy operation
+ block_vector_unmanaged& axpy (const field_type& a, const block_vector_unmanaged& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+ for (int i=0; i<this->n; ++i) (*this)[i].axpy(a,y[i]);
+ return *this;
+ }
+
+
+ //===== Euclidean scalar product
+
+ //! scalar product
+ field_type operator* (const block_vector_unmanaged& y) const
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (this->n!=y.N()) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+ field_type sum=0;
+ for (int i=0; i<this->n; ++i) sum += (*this)[i]*y[i];
+ return sum;
+ }
+
+
+ //===== norms
+
+ //! one norm (sum over absolute values of entries)
+ double one_norm () const
+ {
+ double sum=0;
+ for (int i=0; i<this->n; ++i) sum += (*this)[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<this->n; ++i) sum += (*this)[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<this->n; ++i) sum += (*this)[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<this->n; ++i) sum += (*this)[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<this->n; ++i) max = std::max(max,(*this)[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<this->n; ++i) max = std::max(max,(*this)[i].infinity_norm_real());
+ return max;
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in the vector (are of size 1 here)
+ int N () const
+ {
+ return this->n;
+ }
+
+ //! dimension of the vector space
+ int dim () const
+ {
+ int d=0;
+ for (int i=0; i<this->n; i++)
+ d += (*this)[i].dim();
+ return d;
+ }
+
+ protected:
+ //! make constructor protected, so only derived classes can be instantiated
+ block_vector_unmanaged () : base_array_unmanaged<B,A>()
+ { }
+ };
+
+
+ /** BlockVector adds memory management with ordinary
+ copy semantics to the block_vector_unmanaged template.
+
+ 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=ISTLAllocator>
+ class BlockVector : public block_vector_unmanaged<B,A>
+ {
+ 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};
+
+ //! make iterators available as types
+ typedef typename block_vector_unmanaged<B,A>::Iterator Iterator;
+
+ //! make iterators available as types
+ typedef typename block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
+
+ //===== constructors and such
+
+ //! makes empty vector
+ BlockVector () : block_vector_unmanaged<B,A>()
+ { }
+
+ //! make vector with _n components
+ BlockVector (int _n)
+ {
+ this->n = _n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+ }
+
+ //! copy constructor
+ BlockVector (const BlockVector& a)
+ {
+ // allocate memory with same size as a
+ this->n = a.n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ // and copy elements
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+
+ //! construct from base class object
+ BlockVector (const block_vector_unmanaged<B,A>& _a)
+ {
+ // upcast, because protected data inaccessible
+ const BlockVector& a = static_cast<const BlockVector&>(_a);
+
+ // allocate memory with same size as a
+ this->n = a.n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ // and copy elements
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+
+
+ //! free dynamic memory
+ ~BlockVector ()
+ {
+ if (this->n>0) A::template free<B>(this->p);
+ }
+
+ //! reallocate vector to given size, any data is lost
+ void resize (int _n)
+ {
+ if (this->n==_n) return;
+
+ if (this->n>0) A::template free<B>(this->p);
+ this->n = _n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+ }
+
+ //! assignment
+ BlockVector& operator= (const BlockVector& a)
+ {
+ if (&a!=this) // check if this and a are different objects
+ {
+ // adjust size of vector
+ if (this->n!=a.n) // check if size is different
+ {
+ if (this->n>0) A::template free<B>(this->p); // delete old memory
+ this->n = a.n;
+ if (this->n>0)
+ this->p = A::template malloc<B>(this->n);
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+ }
+ // copy data
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+ return *this;
+ }
+
+ //! assign from base class object
+ BlockVector& operator= (const block_vector_unmanaged<B,A>& a)
+ {
+ // forward to regular assignement operator
+ return this->operator=(static_cast<const BlockVector&>(a));
+ }
+
+ //! assign from scalar
+ BlockVector& operator= (const field_type& k)
+ {
+ // forward to operator= in base class
+ (static_cast<block_vector_unmanaged<B,A>&>(*this)) = k;
+ return *this;
+ }
+ };
+
+
+ /** BlockVectorWindow adds window manipulation functions
+ to the block_vector_unmanaged template.
+
+ This class has no memory management. It assumes that the storage
+ for the entries of the vector is maintained outside of this class.
+
+ But you can copy objects of this class and of the base class
+ with reference semantics.
+
+ Assignment copies the data, if the format is incopmpatible with
+ the argument an exception is thrown in debug mode.
+
+ 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=ISTLAllocator>
+ class BlockVectorWindow : public block_vector_unmanaged<B,A>
+ {
+ 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};
+
+ //! make iterators available as types
+ typedef typename block_vector_unmanaged<B,A>::Iterator Iterator;
+
+ //! make iterators available as types
+ typedef typename block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
+
+
+ //===== constructors and such
+ //! makes empty array
+ BlockVectorWindow () : block_vector_unmanaged<B,A>()
+ { }
+
+ //! make array from given pointer and size
+ BlockVectorWindow (B* _p, int _n)
+ {
+ this->n = _n;
+ this->p = _p;
+ }
+
+ //! copy constructor, this has reference semantics!
+ BlockVectorWindow (const BlockVectorWindow& a)
+ {
+ this->n = a.n;
+ this->p = a.p;
+ }
+
+ //! construct from base class object with reference semantics!
+ BlockVectorWindow (const block_vector_unmanaged<B,A>& _a)
+ {
+ // cast needed to access protected data
+ const BlockVectorWindow& a = static_cast<const BlockVectorWindow&>(_a);
+
+ // make me point to the other's data
+ this->n = a.n;
+ this->p = a.p;
+ }
+
+
+ //! assignment
+ BlockVectorWindow& operator= (const BlockVectorWindow& a)
+ {
+ // check correct size
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (this->n!=a.N()) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+
+ if (&a!=this) // check if this and a are different objects
+ {
+ // copy data
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+ return *this;
+ }
+
+ //! assign from base class object
+ BlockVectorWindow& operator= (const block_vector_unmanaged<B,A>& a)
+ {
+ // forward to regular assignment operator
+ return this->operator=(static_cast<const BlockVectorWindow&>(a));
+ }
+
+ //! assign from scalar
+ BlockVectorWindow& operator= (const field_type& k)
+ {
+ (static_cast<block_vector_unmanaged<B,A>&>(*this)) = k;
+ return *this;
+ }
+
+
+ //===== window manipulation methods
+
+ //! set size and pointer
+ void set (int _n, B* _p)
+ {
+ this->n = _n;
+ this->p = _p;
+ }
+
+ //! set size only
+ void setsize (int _n)
+ {
+ this->n = _n;
+ }
+
+ //! set pointer only
+ void setptr (B* _p)
+ {
+ this->p = _p;
+ }
+
+ //! get pointer
+ B* getptr ()
+ {
+ return this->p;
+ }
+
+ //! get size
+ int getsize ()
+ {
+ return this->n;
+ }
+ };
+
+
+
+ /** compressed_block_vector_unmanaged extends the compressed base_array_unmanaged by
+ vector operations such as addition and scalar multiplication.
+ No memory management is added.
+
+ 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=ISTLAllocator>
+ class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,A>
+ {
+ 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;
+
+ //! make iterators available as types
+ typedef typename compressed_base_array_unmanaged<B,A>::iterator Iterator;
+
+ //! make iterators available as types
+ typedef typename compressed_base_array_unmanaged<B,A>::const_iterator ConstIterator;
+
+
+ //===== assignment from scalar
+
+ compressed_block_vector_unmanaged& operator= (const field_type& k)
+ {
+ for (int i=0; i<this->n; i++)
+ (this->p)[i] = k;
+ return *this;
+ }
+
+
+ //===== vector space arithmetic
+
+ //! vector space addition
+ template<class V>
+ compressed_block_vector_unmanaged& operator+= (const V& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
+#endif
+ for (int i=0; i<this->n; ++i) (this->p)[i] += y[(this->j)[i]];
+ return *this;
+ }
+
+ //! vector space subtraction
+ template<class V>
+ compressed_block_vector_unmanaged& operator-= (const V& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
+#endif
+ for (int i=0; i<this->n; ++i) (this->p)[i] -= y[(this->j)[i]];
+ return *this;
+ }
+
+ //! vector space axpy operation
+ template<class V>
+ compressed_block_vector_unmanaged& axpy (const field_type& a, const V& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
+#endif
+ for (int i=0; i<this->n; ++i) (this->p)[i].axpy(a,y[(this->j)[i]]);
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ compressed_block_vector_unmanaged& operator*= (const field_type& k)
+ {
+ for (int i=0; i<this->n; ++i) (this->p)[i] *= k;
+ return *this;
+ }
+
+ //! vector space division by scalar
+ compressed_block_vector_unmanaged& operator/= (const field_type& k)
+ {
+ for (int i=0; i<this->n; ++i) (this->p)[i] /= k;
+ return *this;
+ }
+
+
+ //===== Euclidean scalar product
+
+ //! scalar product
+ field_type operator* (const compressed_block_vector_unmanaged& y) const
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (!includesindexset(y)) DUNE_THROW(ISTLError,"index set mismatch");
+#endif
+ field_type sum=0;
+ for (int i=0; i<this->n; ++i)
+ sum += (this->p)[i] * y[(this->j)[i]];
+ return sum;
+ }
+
+
+ //===== norms
+
+ //! one norm (sum over absolute values of entries)
+ double one_norm () const
+ {
+ double sum=0;
+ for (int i=0; i<this->n; ++i) sum += (this->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<this->n; ++i) sum += (this->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<this->n; ++i) sum += (this->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<this->n; ++i) sum += (this->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<this->n; ++i) max = std::max(max,(this->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<this->n; ++i) max = std::max(max,(this->p)[i].infinity_norm_real());
+ return max;
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in the vector (are of size 1 here)
+ int N () const
+ {
+ return this->n;
+ }
+
+ //! dimension of the vector space
+ int dim () const
+ {
+ int d=0;
+ for (int i=0; i<this->n; i++)
+ d += (this->p)[i].dim();
+ return d;
+ }
+
+ protected:
+ //! make constructor protected, so only derived classes can be instantiated
+ compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,A>()
+ { }
+
+ //! return true if index sets coincide
+ template<class V>
+ bool includesindexset (const V& y)
+ {
+ typename V::Iterator e=y.end();
+ for (int i=0; i<this->n; i++)
+ if (y.find((this->j)[i])==e)
+ return false;
+ return true;
+ }
+ };
+
+
+ /** CompressedBlockVectorWindow adds window manipulation functions
+ to the compressed_block_vector_unmanaged template.
+
+ This class has no memory management. It assumes that the storage
+ for the entries of the vector and its index set is maintained outside of this class.
+
+ But you can copy objects of this class and of the base class
+ with reference semantics.
+
+ Assignment copies the data, if the format is incopmpatible with
+ the argument an exception is thrown in debug mode.
+
+ 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=ISTLAllocator>
+ class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,A>
+ {
+ 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};
+
+ //! make iterators available as types
+ typedef typename compressed_block_vector_unmanaged<B,A>::Iterator Iterator;
+
+ //! make iterators available as types
+ typedef typename compressed_block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
+
+
+ //===== constructors and such
+ //! makes empty array
+ CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,A>()
+ { }
+
+ //! make array from given pointers and size
+ CompressedBlockVectorWindow (B* _p, int* _j, int _n)
+ {
+ this->n = _n;
+ this->p = _p;
+ this->j = _j;
+ }
+
+ //! copy constructor, this has reference semantics!
+ CompressedBlockVectorWindow (const CompressedBlockVectorWindow& a)
+ {
+ this->n = a.n;
+ this->p = a.p;
+ this->j = a.j;
+ }
+
+ //! construct from base class object with reference semantics!
+ CompressedBlockVectorWindow (const compressed_block_vector_unmanaged<B,A>& _a)
+ {
+ // cast needed to access protected data (upcast)
+ const CompressedBlockVectorWindow& a = static_cast<const CompressedBlockVectorWindow&>(_a);
+
+ // make me point to the other's data
+ this->n = a.n;
+ this->p = a.p;
+ this->j = a.j;
+ }
+
+
+ //! assignment
+ CompressedBlockVectorWindow& operator= (const CompressedBlockVectorWindow& a)
+ {
+ // check correct size
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (this->n!=a.N()) DUNE_THROW(ISTLError,"vector size mismatch");
+#endif
+
+ if (&a!=this) // check if this and a are different objects
+ {
+ // copy data
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ for (int i=0; i<this->n; i++) this->j[i]=a.j[i];
+ }
+ return *this;
+ }
+
+ //! assign from base class object
+ CompressedBlockVectorWindow& operator= (const compressed_block_vector_unmanaged<B,A>& a)
+ {
+ // forward to regular assignment operator
+ return this->operator=(static_cast<const CompressedBlockVectorWindow&>(a));
+ }
+
+ //! assign from scalar
+ CompressedBlockVectorWindow& operator= (const field_type& k)
+ {
+ (static_cast<compressed_block_vector_unmanaged<B,A>&>(*this)) = k;
+ return *this;
+ }
+
+
+ //===== window manipulation methods
+
+ //! set size and pointer
+ void set (int _n, B* _p, int* _j)
+ {
+ this->n = _n;
+ this->p = _p;
+ this->j = _j;
+ }
+
+ //! set size only
+ void setsize (int _n)
+ {
+ this->n = _n;
+ }
+
+ //! set pointer only
+ void setptr (B* _p)
+ {
+ this->p = _p;
+ }
+
+ //! set pointer only
+ void setindexptr (int* _j)
+ {
+ this->j = _j;
+ }
+
+ //! get pointer
+ B* getptr ()
+ {
+ return this->p;
+ }
+
+ //! get pointer
+ int* getindexptr ()
+ {
+ return this->j;
+ }
+
+ //! get size
+ int getsize ()
+ {
+ return this->n;
+ }
+ };
+
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/doc/istl.tex b/istl/doc/istl.tex
index acdd177b2..3d260d1ee 100644
--- a/istl/doc/istl.tex
+++ b/istl/doc/istl.tex
@@ -13,7 +13,7 @@
\usepackage[a4paper,body={148mm,240mm,nohead}]{geometry}
\usepackage[ansinew]{inputenc}
\usepackage{listings}
-\lstset{language=C++, indent=20pt, basicstyle=\ttfamily,
+\lstset{language=C++, basicstyle=\ttfamily,
stringstyle=\ttfamily, commentstyle=\it, extendedchars=true}
\newif\ifpdf
@@ -235,7 +235,8 @@ access is provided via \lstinline!operator[](int i)! where the indices
are in the range $0,\ldots,n-1$ with the number of blocks $n$ given by
the \lstinline!N! method. Here is a code fragment for illustration:
\begin{lstlisting}{}
-Dune::BlockVector<Dune::FieldVector<std::complex<double>,2> > v(20);
+typedef Dune::FieldVector<std::complex<double>,2> BType;
+Dune::BlockVector<BType> v(20);
v[1] = 3.14;
v[3][0] = 2.56;
v[3][1] = std::complex<double>(1,-1);
@@ -299,8 +300,8 @@ in the section on memory management.
\texttt{X::operator[](int)} & \texttt{const field\_type\&} & \\
\texttt{X::begin()} & \texttt{Iterator} & \\
\texttt{X::end()} & \texttt{Iterator} & \\
-\texttt{X::begin()} & \texttt{ConstIterator} & \\
-\texttt{X::end()} & \texttt{ConstIterator} & \\
+\texttt{X::rbegin()} & \texttt{Iterator} & for reverse iteration \\
+\texttt{X::rend()} & \texttt{Iterator} & \\
\texttt{X::find(int)} & \texttt{Iterator} & \\
\hline
\texttt{X::operator+=(X\&)} & \texttt{X\&} & $x = x+y$\\
@@ -327,7 +328,7 @@ in the section on memory management.
\label{Fig:VectorMembers}
\end{figure}
-\subsection{Object memory model and memory management}
+\subsection[Memory model]{Object memory model and memory management}
The memory model for all ISTL objects is deep copy as in the Standard
Template Library and in contrast to the Matrix Template
@@ -346,7 +347,9 @@ structure can be generated simply with the copy constructor.
\subsection{Matrix classes}
-\subsection{Matrices are containers of containers}
+\subsection[Matrix containers]{Matrices are containers of containers}
+
+\subsection{Precision control}
\subsection{Operations}
@@ -364,6 +367,8 @@ structure can be generated simply with the copy constructor.
\texttt{M::blocklevel} & \texttt{int} & block levels inside\\
\texttt{M::RowIterator} & T & over rows\\
\texttt{M::ColIterator} & T & over columns\\
+\texttt{M::ConstRowIterator} & T & over rows\\
+\texttt{M::ConstColIterator} & T & over columns\\
\hline
\texttt{M::M()} & & empty matrix \\
\texttt{M::M(M\&)} & & deep copy\\
@@ -375,6 +380,8 @@ structure can be generated simply with the copy constructor.
\texttt{M::operator[](int)} & \texttt{const row\_type\&} & \\
\texttt{M::begin()} & \texttt{RowIterator} & \\
\texttt{M::end()} & \texttt{RowIterator} & \\
+\texttt{M::rbegin()} & \texttt{RowIterator} & reverse iteration\\
+\texttt{M::rend()} & \texttt{RowIterator} & \\
\hline
\texttt{M::operator*=(field\_type\&)} & \texttt{M\&} & $A = \alpha A$\\
\texttt{M::operator/=(field\_type\&)} & \texttt{M\&} & $A = \alpha^{-1} A$\\
@@ -389,8 +396,14 @@ structure can be generated simply with the copy constructor.
\texttt{M::mmhv(X\& x,Y\& y)} & & $y = y - A^Hx$\\
\texttt{M::usmhv(field\_type\&,X\& x,Y\& y)} & & $y = y + \alpha A^Hx$\\
\hline
-\texttt{M::frobenius\_norm()} & \texttt{double} &\small$\sqrt{\sum\limits_{i,j} (Re(a_{ij})^2+Im(a_{ij})^2)}$\\
-\texttt{M::frobenius\_norm2()} & \texttt{double} &\small$\sum\limits_{i,j} (Re(a_{ij})^2+Im(a_{ij})^2)$\\
+\texttt{M::solve(X\& x,Y\& b)} & & $x = A^{-1}b$\\
+\texttt{M::inverse(M\& B)} & & $B=A^{-1}$\\
+\texttt{M::leftmultiply(M\& B)} & \texttt{M\&} & $A = BA$\\
+\hline
+\texttt{M::frobenius\_norm()} & \texttt{double} & see text\\
+\texttt{M::frobenius\_norm2()} & \texttt{double} &see text\\
+\texttt{X::infinity\_norm()} & \texttt{double} & see text\\
+\texttt{X::infinity\_norm\_real()} & \texttt{double} &see text\\
\hline
\texttt{M::N()} & \texttt{int} & row blocks\\
\texttt{M::M()} & \texttt{int} & col blocks\\
diff --git a/istl/fmatrix.hh b/istl/fmatrix.hh
new file mode 100644
index 000000000..eb5ec5759
--- /dev/null
+++ b/istl/fmatrix.hh
@@ -0,0 +1,1227 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_FMATRIX_HH__
+#define __DUNE_FMATRIX_HH__
+
+#include <math.h>
+#include <complex>
+
+#include "istlexception.hh"
+#include "allocator.hh"
+#include "fvector.hh"
+#include "precision.hh"
+
+/*! \file __FILE__
+
+ This file implements a matrix constructed from a given type
+ representing a field and compile-time given number of rows and columns.
+ */
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ template<class K, int n, int m> class FieldMatrix;
+
+ // template meta program for assignment from scalar
+ template<int I>
+ struct fmmeta_assignscalar {
+ template<class T, class K>
+ static void assignscalar (T* x, const K& k)
+ {
+ fmmeta_assignscalar<I-1>::assignscalar(x,k);
+ x[I] = k;
+ }
+ };
+ template<>
+ struct fmmeta_assignscalar<0> {
+ template<class T, class K>
+ static void assignscalar (T* x, const K& k)
+ {
+ x[0] = k;
+ }
+ };
+
+ // template meta program for operator+=
+ template<int I>
+ struct fmmeta_plusequal {
+ template<class T>
+ static void plusequal (T& x, const T& y)
+ {
+ x[I] += y[I];
+ fmmeta_plusequal<I-1>::plusequal(x,y);
+ }
+ };
+ template<>
+ struct fmmeta_plusequal<0> {
+ template<class T>
+ static void plusequal (T& x, const T& y)
+ {
+ x[0] += y[0];
+ }
+ };
+
+ // template meta program for operator-=
+ template<int I>
+ struct fmmeta_minusequal {
+ template<class T>
+ static void minusequal (T& x, const T& y)
+ {
+ x[I] -= y[I];
+ fmmeta_minusequal<I-1>::minusequal(x,y);
+ }
+ };
+ template<>
+ struct fmmeta_minusequal<0> {
+ template<class T>
+ static void minusequal (T& x, const T& y)
+ {
+ x[0] -= y[0];
+ }
+ };
+
+ // template meta program for operator*=
+ template<int I>
+ struct fmmeta_multequal {
+ template<class T, class K>
+ static void multequal (T& x, const K& k)
+ {
+ x[I] *= k;
+ fmmeta_multequal<I-1>::multequal(x,k);
+ }
+ };
+ template<>
+ struct fmmeta_multequal<0> {
+ template<class T, class K>
+ static void multequal (T& x, const K& k)
+ {
+ x[0] *= k;
+ }
+ };
+
+ // template meta program for operator/=
+ template<int I>
+ struct fmmeta_divequal {
+ template<class T, class K>
+ static void divequal (T& x, const K& k)
+ {
+ x[I] /= k;
+ fmmeta_divequal<I-1>::divequal(x,k);
+ }
+ };
+ template<>
+ struct fmmeta_divequal<0> {
+ template<class T, class K>
+ static void divequal (T& x, const K& k)
+ {
+ x[0] /= k;
+ }
+ };
+
+ // template meta program for dot
+ template<int I>
+ struct fmmeta_dot {
+ template<class X, class Y, class K>
+ static K dot (const X& x, const Y& y)
+ {
+ return x[I]*y[I] + fmmeta_dot<I-1>::template dot<X,Y,K>(x,y);
+ }
+ };
+ template<>
+ struct fmmeta_dot<0> {
+ template<class X, class Y, class K>
+ static K dot (const X& x, const Y& y)
+ {
+ return x[0]*y[0];
+ }
+ };
+
+ // template meta program for umv(x,y)
+ template<int I>
+ struct fmmeta_umv {
+ template<class Mat, class X, class Y, int c>
+ static void umv (const Mat& A, const X& x, Y& y)
+ {
+ typedef typename Mat::row_type R;
+ typedef typename Mat::field_type K;
+ y[I] += fmmeta_dot<c>::template dot<R,X,K>(A[I],x);
+ fmmeta_umv<I-1>::template umv<Mat,X,Y,c>(A,x,y);
+ }
+ };
+ template<>
+ struct fmmeta_umv<0> {
+ template<class Mat, class X, class Y, int c>
+ static void umv (const Mat& A, const X& x, Y& y)
+ {
+ typedef typename Mat::row_type R;
+ typedef typename Mat::field_type K;
+ y[0] += fmmeta_dot<c>::template dot<R,X,K>(A[0],x);
+ }
+ };
+
+ // template meta program for mmv(x,y)
+ template<int I>
+ struct fmmeta_mmv {
+ template<class Mat, class X, class Y, int c>
+ static void mmv (const Mat& A, const X& x, Y& y)
+ {
+ typedef typename Mat::row_type R;
+ typedef typename Mat::field_type K;
+ y[I] -= fmmeta_dot<c>::template dot<R,X,K>(A[I],x);
+ fmmeta_mmv<I-1>::template mmv<Mat,X,Y,c>(A,x,y);
+ }
+ };
+ template<>
+ struct fmmeta_mmv<0> {
+ template<class Mat, class X, class Y, int c>
+ static void mmv (const Mat& A, const X& x, Y& y)
+ {
+ typedef typename Mat::row_type R;
+ typedef typename Mat::field_type K;
+ y[0] -= fmmeta_dot<c>::template dot<R,X,K>(A[0],x);
+ }
+ };
+
+ template<class K, int n, int m, class X, class Y>
+ inline void fm_mmv (const FieldMatrix<K,n,m>& A, const X& x, Y& y)
+ {
+ for (int i=0; i<n; i++)
+ for (int j=0; j<m; j++)
+ y[i] -= A[i][j]*x[j];
+ }
+
+ template<class K>
+ inline void fm_mmv (const FieldMatrix<K,1,1>& A, const FieldVector<K,1>& x, FieldVector<K,1>& y)
+ {
+ y[0] -= A[0][0]*x[0];
+ }
+
+ // template meta program for usmv(x,y)
+ template<int I>
+ struct fmmeta_usmv {
+ template<class Mat, class K, class X, class Y, int c>
+ static void usmv (const Mat& A, const K& alpha, const X& x, Y& y)
+ {
+ typedef typename Mat::row_type R;
+ y[I] += alpha*fmmeta_dot<c>::template dot<R,X,K>(A[I],x);
+ fmmeta_usmv<I-1>::template usmv<Mat,K,X,Y,c>(A,alpha,x,y);
+ }
+ };
+ template<>
+ struct fmmeta_usmv<0> {
+ template<class Mat, class K, class X, class Y, int c>
+ static void usmv (const Mat& A, const K& alpha, const X& x, Y& y)
+ {
+ typedef typename Mat::row_type R;
+ y[0] += alpha*fmmeta_dot<c>::template dot<R,X,K>(A[0],x);
+ }
+ };
+
+ // conjugate komplex does nothing for non-complex types
+ template<class K>
+ inline K fm_ck (const K& k)
+ {
+ return k;
+ }
+
+ // conjugate komplex
+ template<class K>
+ inline std::complex<K> fm_ck (const std::complex<K>& c)
+ {
+ return std::complex<K>(c.real(),-c.imag());
+ }
+
+
+ //! solve small system
+ template<class K, int n, class V>
+ void fm_solve (const FieldMatrix<K,n,n>& A, V& x, const V& b)
+ {
+ // Gaussian elimination with maximum column pivot
+ double norm=A.one_norm_real(); // for relative thresholds
+ double pivthres = std::max(ISTLPrecision::absolute_limit(),norm*ISTLPrecision::pivoting_limit());
+ double singthres = std::max(ISTLPrecision::absolute_limit(),norm*ISTLPrecision::singular_limit());
+ V& rhs = x; // use x to store rhs
+ rhs = b; // copy data
+
+ // elimination phase
+ for (int i=0; i<n; i++) // loop over all rows
+ {
+ double pivmax=fvmeta_absreal(A[i][i]);
+
+ // pivoting ?
+ if (pivmax<pivthres)
+ {
+ // compute maximum of row
+ int imax=i; double abs;
+ for (int k=i+1; k<n; k++)
+ if ((abs=fvmeta_absreal(A[k][i]))>pivmax)
+ {
+ pivmax = abs; imax = k;
+ }
+ // swap rows
+ if (imax!=i)
+ for (int j=i; j<n; j++)
+ std::swap(A[i][j],A[imax][j]);
+ }
+
+ // singular ?
+ if (pivmax<singthres)
+ DUNE_THROW(ISTLError,"matrix is singular");
+
+ // eliminate
+ for (int k=i+1; k<n; k++)
+ {
+ K factor = -A[k][i]/A[i][i];
+ for (int j=i+1; j<n; j++)
+ A[k][j] += factor*A[i][j];
+ rhs[k] += factor*rhs[i];
+ }
+ }
+
+ // backsolve
+ for (int i=n-1; i>=0; i--)
+ {
+ for (int j=i+1; j<n; j++)
+ rhs[i] -= A[i][j]*x[j];
+ x[i] = rhs[i]/A[i][i];
+ }
+ }
+
+ //! special case for 1x1 matrix, x and b may be identical
+ template<class K, class V>
+ inline void fm_solve (const FieldMatrix<K,1,1>& A, V& x, const V& b)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (fvmeta_absreal(A[0][0])<ISTLPrecision::absolute_limit())
+ DUNE_THROW(ISTLError,"matrix is singular");
+#endif
+ x[0] = b[0]/A[0][0];
+ }
+
+ //! special case for 2x2 matrix, x and b may be identical
+ template<class K, class V>
+ inline void fm_solve (const FieldMatrix<K,2,2>& A, V& x, const V& b)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ K detinv = A[0][0]*A[1][1]-A[0][1]*A[1][0];
+ if (fvmeta_absreal(detinv)<ISTLPrecision::absolute_limit())
+ DUNE_THROW(ISTLError,"matrix is singular");
+ detinv = 1/detinv;
+#else
+ K detinv = 1.0/(A[0][0]*A[1][1]-A[0][1]*A[1][0]);
+#endif
+
+ K temp = b[0];
+ x[0] = detinv*(A[1][1]*b[0]-A[0][1]*b[1]);
+ x[1] = detinv*(A[0][0]*b[1]-A[1][0]*temp);
+ }
+
+
+
+ //! compute inverse
+ template<class K, int n>
+ void fm_inverse (const FieldMatrix<K,n,n>& Ain, FieldMatrix<K,n,n>& Ainv)
+ {
+ FieldMatrix<K,n,n> A(Ain);
+ FieldMatrix<K,n,n>& L=A;
+ FieldMatrix<K,n,n>& U=A;
+
+ double norm=A.one_norm_real(); // for relative thresholds
+ double pivthres = std::max(ISTLPrecision::absolute_limit(),norm*ISTLPrecision::pivoting_limit());
+ double singthres = std::max(ISTLPrecision::absolute_limit(),norm*ISTLPrecision::singular_limit());
+
+ // LU decomposition of A
+ for (int i=0; i<n; i++) // loop over all rows
+ {
+ double pivmax=fvmeta_absreal(A[i][i]);
+
+ // pivoting ?
+ if (pivmax<pivthres)
+ {
+ // compute maximum of row
+ int imax=i; double abs;
+ for (int k=i+1; k<n; k++)
+ if ((abs=fvmeta_absreal(A[k][i]))>pivmax)
+ {
+ pivmax = abs; imax = k;
+ }
+ // swap rows
+ if (imax!=i)
+ for (int j=i; j<n; j++)
+ std::swap(A[i][j],A[imax][j]);
+ }
+
+ // singular ?
+ if (pivmax<singthres)
+ DUNE_THROW(ISTLError,"matrix is singular");
+
+ // eliminate
+ for (int k=i+1; k<n; k++)
+ {
+ K factor = -A[k][i]/A[i][i];
+ L[k][i] = factor;
+ for (int j=i+1; j<n; j++)
+ A[k][j] += factor*A[i][j];
+ }
+ }
+
+ // initialize inverse
+ Ainv = 0;
+ for (int i=0; i<n; i++) Ainv[i][i] = 1;
+
+ // L Y = I
+ for (int i=0; i<n; i++)
+ for (int k=0; k<n; k++)
+ {
+ for (int j=0; j<n; j++)
+ Ainv[i][k] -= L[i][j]*Ainv[j][k];
+ }
+
+ // U A^{-1} = Y
+ for (int i=n-1; i>=0; i--)
+ for (int k=0; k<n; k++)
+ {
+ for (int j=i+1; j<n; j++)
+ Ainv[i][k] -= U[i][j]*Ainv[j][k];
+ Ainv[i][k] = Ainv[i][k]/U[i][i];
+ }
+ }
+
+ //! compute inverse n=1
+ template<class K>
+ void fm_inverse (const FieldMatrix<K,1,1>& Ain, FieldMatrix<K,1,1>& Ainv)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (fvmeta_absreal(Ain[0][0])<ISTLPrecision::absolute_limit())
+ DUNE_THROW(ISTLError,"matrix is singular");
+#endif
+ Ainv[0][0] = 1/Ain[0][0];
+ }
+
+ //! compute inverse n=2
+ template<class K>
+ void fm_inverse (const FieldMatrix<K,2,2>& Ain, FieldMatrix<K,2,2>& Ainv)
+ {
+ K detinv = Ain[0][0]*Ain[1][1]-Ain[0][1]*Ain[1][0];
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (fvmeta_absreal(detinv)<ISTLPrecision::absolute_limit())
+ DUNE_THROW(ISTLError,"matrix is singular");
+#endif
+ detinv = 1/detinv;
+
+ Ainv[0][0] = Ain[1][1]*detinv;
+ Ainv[0][1] = -Ain[0][1]*detinv;
+ Ainv[1][0] = -Ain[1][0]*detinv;
+ Ainv[1][1] = Ain[0][0]*detinv;
+ }
+
+ //! left multiplication with matrix
+ template<class K, int n, int m>
+ void fm_leftmultiply (const FieldMatrix<K,n,n>& M, FieldMatrix<K,n,m>& A)
+ {
+ FieldMatrix<K,n,m> C(A);
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<m; j++)
+ {
+ A[i][j] = 0;
+ for (int k=0; k<n; k++)
+ A[i][j] += M[i][k]*C[k][j];
+ }
+ }
+
+ //! left multiplication with matrix, n=1
+ template<class K>
+ void fm_leftmultiply (const FieldMatrix<K,1,1>& M, FieldMatrix<K,1,1>& A)
+ {
+ A[0][0] *= M[0][0];
+ }
+
+ //! left multiplication with matrix, n=2
+ template<class K>
+ void fm_leftmultiply (const FieldMatrix<K,2,2>& M, FieldMatrix<K,2,2>& A)
+ {
+ FieldMatrix<K,2,2> C(A);
+
+ A[0][0] = M[0][0]*C[0][0] + M[0][1]*C[1][0];
+ A[0][1] = M[0][0]*C[0][1] + M[0][1]*C[1][1];
+ A[1][0] = M[1][0]*C[0][0] + M[1][1]*C[1][0];
+ A[1][1] = M[1][0]*C[0][1] + M[1][1]*C[1][1];
+ }
+
+ /**! Matrices represent linear maps from a vector space V to a vector space W.
+ This class represents such a linear map by storing a two-dimensional
+ array of numbers of a given field type K. The number of rows and
+ columns is given at compile time.
+
+ Implementation of all members uses template meta programs where appropriate
+ */
+ template<class K, int n, int m>
+ class FieldMatrix
+ {
+ 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 = 1};
+
+ //! Each row is implemented by a field vector
+ typedef FieldVector<K,m> row_type;
+
+ //! export size
+ enum {rows = n, cols = m};
+
+ //===== random access interface to rows of the matrix
+
+ //! random access to the rows
+ row_type& 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 row_type& 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 interface to rows of the matrix
+
+ // forward declaration
+ class ConstIterator;
+
+ //! Iterator access to rows
+ class Iterator
+ {
+ public:
+ //! constructor
+ Iterator (row_type* _p, int _i)
+ {
+ p = _p;
+ i = _i;
+ }
+
+ //! empty constructor, use with care!
+ Iterator ()
+ { }
+
+ //! prefix increment
+ Iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ Iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ row_type& operator* ()
+ {
+ return p[i];
+ }
+
+ //! arrow
+ row_type* operator-> ()
+ {
+ return p+i;
+ }
+
+ //! return index
+ int index ()
+ {
+ return i;
+ }
+
+ friend class ConstIterator;
+
+ private:
+ row_type* p;
+ int i;
+ };
+
+ //! begin iterator
+ Iterator begin ()
+ {
+ return Iterator(p,0);
+ }
+
+ //! end iterator
+ Iterator end ()
+ {
+ return Iterator(p,n);
+ }
+
+ //! begin iterator
+ Iterator rbegin ()
+ {
+ return Iterator(p,n-1);
+ }
+
+ //! end iterator
+ Iterator rend ()
+ {
+ return Iterator(p,-1);
+ }
+
+ //! rename the iterators for easier access
+ typedef Iterator RowIterator;
+ typedef typename row_type::Iterator ColIterator;
+
+
+ //! Iterator access to rows
+ class ConstIterator
+ {
+ public:
+ //! constructor
+ ConstIterator (const row_type* _p, int _i) : p(_p), i(_i)
+ { }
+
+ //! empty constructor, use with care!
+ ConstIterator ()
+ {
+ p = 0;
+ i = 0;
+ }
+
+ //! prefix increment
+ ConstIterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ ConstIterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const ConstIterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const ConstIterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ const row_type& operator* () const
+ {
+ return p[i];
+ }
+
+ //! arrow
+ const row_type* operator-> () const
+ {
+ return p+i;
+ }
+
+ //! return index
+ int index () const
+ {
+ return i;
+ }
+
+ friend class Iterator;
+
+ private:
+ const row_type* p;
+ int i;
+ };
+
+ //! begin iterator
+ ConstIterator begin () const
+ {
+ return ConstIterator(p,0);
+ }
+
+ //! end iterator
+ ConstIterator end () const
+ {
+ return ConstIterator(p,n);
+ }
+
+ //! begin iterator
+ ConstIterator rbegin () const
+ {
+ return ConstIterator(p,n-1);
+ }
+
+ //! end iterator
+ ConstIterator rend () const
+ {
+ return ConstIterator(p,-1);
+ }
+
+ //! rename the iterators for easier access
+ typedef ConstIterator ConstRowIterator;
+ typedef typename row_type::ConstIterator ConstColIterator;
+
+
+ //===== assignment from scalar
+ FieldMatrix& operator= (const K& k)
+ {
+ fmmeta_assignscalar<n-1>::assignscalar(p,k);
+ return *this;
+ }
+
+ //===== vector space arithmetic
+
+ //! vector space addition
+ FieldMatrix& operator+= (const FieldMatrix& y)
+ {
+ fmmeta_plusequal<n-1>::plusequal(*this,y);
+ return *this;
+ }
+
+ //! vector space subtraction
+ FieldMatrix& operator-= (const FieldMatrix& y)
+ {
+ fmmeta_minusequal<n-1>::minusequal(*this,y);
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ FieldMatrix& operator*= (const K& k)
+ {
+ fmmeta_multequal<n-1>::multequal(*this,k);
+ return *this;
+ }
+
+ //! vector space division by scalar
+ FieldMatrix& operator/= (const K& k)
+ {
+ fmmeta_divequal<n-1>::divequal(*this,k);
+ return *this;
+ }
+
+ //===== linear maps
+
+ //! y += A x
+ template<class X, class Y>
+ void umv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ fmmeta_umv<n-1>::template umv<FieldMatrix,X,Y,m-1>(*this,x,y);
+ }
+
+ //! y += A^T x
+ template<class X, class Y>
+ void umtv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<m; j++)
+ y[j] += p[i][j]*x[i];
+ }
+
+ //! y += A^H x
+ template<class X, class Y>
+ void umhv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<m; j++)
+ y[j] += fm_ck(p[i][j])*x[i];
+ }
+
+ //! y -= A x
+ template<class X, class Y>
+ void mmv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ fmmeta_mmv<n-1>::template mmv<FieldMatrix,X,Y,m-1>(*this,x,y);
+ //fm_mmv(*this,x,y);
+ }
+
+ //! y -= A^T x
+ template<class X, class Y>
+ void mmtv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<m; j++)
+ y[j] -= p[i][j]*x[i];
+ }
+
+ //! y -= A^H x
+ template<class X, class Y>
+ void mmhv (X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<m; j++)
+ y[j] -= fm_ck(p[i][j])*x[i];
+ }
+
+ //! y += alpha A x
+ template<class X, class Y>
+ void usmv (const K& alpha, X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ fmmeta_usmv<n-1>::template usmv<FieldMatrix,K,X,Y,m-1>(*this,alpha,x,y);
+ }
+
+ //! y += alpha A^T x
+ template<class X, class Y>
+ void usmtv (const K& alpha, X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<m; j++)
+ y[j] += alpha*p[i][j]*x[i];
+ }
+
+ //! y += alpha A^H x
+ template<class X, class Y>
+ void usmhv (const K& alpha, X& x, Y& y)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (x.N()!=N()) DUNE_THROW(ISTLError,"index out of range");
+ if (y.N()!=M()) DUNE_THROW(ISTLError,"index out of range");
+#endif
+
+ for (int i=0; i<n; i++)
+ for (int j=0; j<m; j++)
+ y[j] += alpha*fm_ck(p[i][j])*x[i];
+ }
+
+ //===== norms
+
+ //! frobenius norm: sqrt(sum over squared values of entries)
+ double frobenius_norm () const
+ {
+ double sum=0;
+ for (int i=0; i<n; ++i) sum += p[i].two_norm2();
+ return sqrt(sum);
+ }
+
+ //! square of frobenius norm, need for block recursion
+ double frobenius_norm2 () const
+ {
+ double sum=0;
+ for (int i=0; i<n; ++i) sum += p[i].two_norm2();
+ return sum;
+ }
+
+ //! infinity norm (row sum norm, how to generalize for blocks?)
+ double infinity_norm () const
+ {
+ double max=0;
+ for (int i=0; i<n; ++i) max = std::max(max,p[i].one_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].one_norm_real());
+ return max;
+ }
+
+ //===== solve
+
+ //! Solve system A x = b
+ template<class V>
+ void solve (V& x, const V& b) const
+ {
+ fm_solve(*this,x,b);
+ }
+
+ //! compute inverse
+ void inverse (FieldMatrix& Ainv) const
+ {
+ fm_inverse(*this,Ainv);
+ }
+
+ //! left multiplication
+ FieldMatrix& leftmultiply (const FieldMatrix<K,n,n>& M)
+ {
+ fm_leftmultiply(M,*this);
+ return *this;
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in row direction
+ int N () const
+ {
+ return n;
+ }
+
+ //! number of blocks in column direction
+ int M () const
+ {
+ return m;
+ }
+
+ //! row dimension of block r
+ int rowdim (int r) const
+ {
+ return 1;
+ }
+
+ //! col dimension of block c
+ int coldim (int c) const
+ {
+ return 1;
+ }
+
+ //! dimension of the destination vector space
+ int rowdim () const
+ {
+ return n;
+ }
+
+ //! dimension of the source vector space
+ int coldim () const
+ {
+ return m;
+ }
+
+ //===== query
+
+ // return true when (i,j) is in pattern
+ bool exists (int i, int j)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (i<0 || i>=n) DUNE_THROW(ISTLError,"index out of range");
+ if (j<0 || i>=m) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ return true;
+ }
+
+ //===== conversion operator
+
+ operator K () {return p[0][0];}
+
+
+ private:
+ // the data, very simply a built in array with rowwise ordering
+ row_type p[n];
+ };
+
+
+ /**! Matrices of size 1x1 are treated in a special way
+ */
+ template<class K>
+ class K11Matrix
+ {
+ 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 = 1};
+
+ //! export size
+ enum {rows = 1, cols = 1};
+
+ //===== random access interface to rows of the matrix
+
+ //! random access to data
+ K& operator() ()
+ {
+ return a;
+ }
+
+ //! same for read only access
+ const K& operator() () const
+ {
+ return a;
+ }
+
+ //===== assignment from scalar
+
+ K11Matrix& operator= (const K& k)
+ {
+ a = k;
+ return *this;
+ }
+
+ //===== vector space arithmetic
+
+ //! vector space addition
+ K11Matrix& operator+= (const K11Matrix& y)
+ {
+ a += y.a;
+ return *this;
+ }
+
+ //! vector space subtraction
+ K11Matrix& operator-= (const K11Matrix& y)
+ {
+ a -= y.a;
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ K11Matrix& operator*= (const K& k)
+ {
+ a *= k;
+ return *this;
+ }
+
+ //! vector space division by scalar
+ K11Matrix& operator/= (const K& k)
+ {
+ a /= k;
+ return *this;
+ }
+
+ //===== linear maps
+
+ //! y += A x
+ void umv (const K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p += a * x.p;
+ }
+
+ //! y += A^T x
+ void umtv (K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p += a * x.p;
+ }
+
+ //! y += A^H x
+ void umhv (K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p += fm_ck(a) * x.p;
+ }
+
+ //! y -= A x
+ void mmv (K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p -= a * x.p;
+ }
+
+ //! y -= A^T x
+ void mmtv (K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p -= a * x.p;
+ }
+
+ //! y -= A^H x
+ void mmhv (K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p -= fm_ck(a) * x.p;
+ }
+
+ //! y += alpha A x
+ void usmv (const K& alpha, K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p += alpha * a * x.p;
+ }
+
+ //! y += alpha A^T x
+ void usmtv (const K& alpha, K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p += alpha * a * x.p;
+ }
+
+ //! y += alpha A^H x
+ void usmhv (const K& alpha, K1Vector<K>& x, K1Vector<K>& y)
+ {
+ y.p += alpha * fm_ck(a) * x.p;
+ }
+
+ //===== norms
+
+ //! frobenius norm: sqrt(sum over squared values of entries)
+ double frobenius_norm () const
+ {
+ return sqrt(fvmeta_abs2(a));
+ }
+
+ //! square of frobenius norm, need for block recursion
+ double frobenius_norm2 () const
+ {
+ return fvmeta_abs2(a);
+ }
+
+ //! infinity norm (row sum norm, how to generalize for blocks?)
+ double infinity_norm () const
+ {
+ return fvmeta_abs(a);
+ }
+
+ //! simplified infinity norm (uses Manhattan norm for complex values)
+ double infinity_norm_real () const
+ {
+ return fvmeta_abs_real(a);
+ }
+
+ //===== solve
+
+ //! Solve system A x = b
+ void solve (K1Vector<K>& x, const K1Vector<K>& b) const
+ {
+ x.p = b.p/a;
+ }
+
+ //! compute inverse
+ void inverse (K11Matrix& Ainv) const
+ {
+ Ainv.a = 1/a;
+ }
+
+ //! left multiplication
+ K11Matrix& leftmultiply (const K11Matrix<K>& M)
+ {
+ a *= M.a;
+ return *this;
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in row direction
+ int N () const
+ {
+ return 1;
+ }
+
+ //! number of blocks in column direction
+ int M () const
+ {
+ return 1;
+ }
+
+ //! row dimension of block r
+ int rowdim (int r) const
+ {
+ return 1;
+ }
+
+ //! col dimension of block c
+ int coldim (int c) const
+ {
+ return 1;
+ }
+
+ //! dimension of the destination vector space
+ int rowdim () const
+ {
+ return 1;
+ }
+
+ //! dimension of the source vector space
+ int coldim () const
+ {
+ return 1;
+ }
+
+ //===== query
+
+ // return true when (i,j) is in pattern
+ bool exists (int i, int j)
+ {
+ return true;
+ }
+
+ //===== conversion operator
+
+ operator K () {return a;}
+
+
+ private:
+ // the data, very simply a built in array with rowwise ordering
+ K a;
+ };
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/fvector.hh b/istl/fvector.hh
new file mode 100644
index 000000000..17455d2d7
--- /dev/null
+++ b/istl/fvector.hh
@@ -0,0 +1,885 @@
+// -*- 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 "istlexception.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 ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ // 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<class K, int n>
+ static void plusequal (FieldVector<K,n>& x, const K& y)
+ {
+ fvmeta_plusequal<I-1>::plusequal(x,y);
+ x[I] += y;
+ }
+ };
+ 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<class K, int n>
+ static void plusequal (FieldVector<K,n>& x, const K& y)
+ {
+ x[0] += y;
+ }
+ };
+
+ // 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<class K, int n>
+ static void minusequal (FieldVector<K,n>& x, const K& y)
+ {
+ fvmeta_minusequal<I-1>::minusequal(x,y);
+ x[I] -= y;
+ }
+ };
+ 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<class K, int n>
+ static void minusequal (FieldVector<K,n>& x, const K& y)
+ {
+ x[0] -= y;
+ }
+ };
+
+ // 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>
+ inline double fvmeta_abs (const K& k)
+ {
+ return (k >= 0) ? k : -k;
+ }
+
+ template<class K>
+ inline double fvmeta_abs (const std::complex<K>& c)
+ {
+ return sqrt(c.real()*c.real() + c.imag()*c.imag());
+ }
+
+ template<class K>
+ inline double fvmeta_absreal (const K& k)
+ {
+ return fvmeta_abs(k);
+ }
+
+ template<class K>
+ inline double fvmeta_absreal (const std::complex<K>& c)
+ {
+ return fvmeta_abs(c.real()) + fvmeta_abs(c.imag());
+ }
+
+ template<class K>
+ inline double fvmeta_abs2 (const K& k)
+ {
+ return k*k;
+ }
+
+ template<class K>
+ inline 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 = 1};
+
+ //! export size
+ enum {size = n};
+
+
+ //===== 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];
+ }
+
+ // forward declaration
+ class ConstIterator;
+
+ //! Iterator class for sequential access
+ class Iterator
+ {
+ public:
+ //! constructor
+ Iterator (K* _p, int _i)
+ {
+ p = _p;
+ i = _i;
+ }
+
+ //! empty constructor, use with care
+ Iterator ()
+ { }
+
+ //! prefix increment
+ Iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ Iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ K& operator* ()
+ {
+ return p[i];
+ }
+
+ //! arrow
+ K* operator-> ()
+ {
+ return p+i;
+ }
+
+ //! return index
+ int index ()
+ {
+ return i;
+ }
+
+ friend class ConstIterator;
+
+ private:
+ K* p;
+ int i;
+ };
+
+ //! begin iterator
+ Iterator begin ()
+ {
+ return Iterator(p,0);
+ }
+
+ //! end iterator
+ Iterator end ()
+ {
+ return Iterator(p,n);
+ }
+
+ //! begin iterator
+ Iterator rbegin ()
+ {
+ return Iterator(p,n-1);
+ }
+
+ //! end iterator
+ Iterator rend ()
+ {
+ return Iterator(p,-1);
+ }
+
+ //! return iterator to given element or end()
+ Iterator find (int i)
+ {
+ if (i>=0 && i<n)
+ return Iterator(p,i);
+ else
+ return Iterator(p,n);
+ }
+
+ //! ConstIterator class for sequential access
+ class ConstIterator
+ {
+ public:
+ //! constructor
+ ConstIterator ()
+ {
+ p = 0;
+ i = 0;
+ }
+
+ //! constructor from pointer
+ ConstIterator (const K* _p, int _i) : p(_p), i(_i)
+ { }
+
+ //! constructor from non_const iterator
+ ConstIterator (const Iterator& it) : p(it.p), i(it.i)
+ { }
+
+ //! prefix increment
+ ConstIterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ ConstIterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const ConstIterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const ConstIterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ const K& operator* () const
+ {
+ return p[i];
+ }
+
+ //! arrow
+ const K* operator-> () const
+ {
+ return p+i;
+ }
+
+ // return index corresponding to pointer
+ int index () const
+ {
+ return i;
+ }
+
+ friend class Iterator;
+
+ private:
+ const K* p;
+ int i;
+ };
+
+ //! begin ConstIterator
+ ConstIterator begin () const
+ {
+ return ConstIterator(p,0);
+ }
+
+ //! end ConstIterator
+ ConstIterator end () const
+ {
+ return ConstIterator(p,n);
+ }
+
+ //! begin ConstIterator
+ ConstIterator rbegin () const
+ {
+ return ConstIterator(p,n-1);
+ }
+
+ //! end ConstIterator
+ ConstIterator rend () const
+ {
+ return ConstIterator(p,-1);
+ }
+
+
+ //===== 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 add scalar to all comps
+ FieldVector& operator+= (const K& k)
+ {
+ fvmeta_plusequal<n-1>::plusequal(*this,k);
+ return *this;
+ }
+
+ //! vector space subtract scalar from all comps
+ FieldVector& operator-= (const K& k)
+ {
+ fvmeta_minusequal<n-1>::minusequal(*this,k);
+ 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;
+ }
+
+ //===== conversion operator
+
+ operator K () {return *p;}
+
+ private:
+ // the data, very simply a built in array
+ K p[n];
+ };
+
+ // forward declarations
+ template<class K> class K1Vector;
+ template<class K> class K11Matrix;
+
+ /**! Vectors containing only one component
+ */
+ template<class K>
+ class K1Vector
+ {
+ public:
+ friend class K11Matrix<K>;
+
+ //===== 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 = 1};
+
+ //! export size
+ enum {size = 1};
+
+ //===== construction
+
+ K1Vector ()
+ { }
+
+ K1Vector (const K& k)
+ {
+ p = k;
+ }
+
+ //===== assignment from basetype
+ K1Vector& operator= (const K& k)
+ {
+ p = k;
+ return *this;
+ }
+
+ //===== access to data
+
+ //! random access
+ K& operator() (void)
+ {
+ return p;
+ }
+
+ //! same for read only access
+ const K& operator() (void) const
+ {
+ return p;
+ }
+
+ //===== vector space arithmetic
+
+ //! vector space addition
+ K1Vector& operator+= (const K1Vector& y)
+ {
+ p += y.p;
+ return *this;
+ }
+
+ //! vector space subtraction
+ K1Vector& operator-= (const K1Vector& y)
+ {
+ p -= y.p;
+ return *this;
+ }
+
+ //! vector space add scalar to each comp
+ K1Vector& operator+= (const K& k)
+ {
+ p += k;
+ return *this;
+ }
+
+ //! vector space subtract scalar from each comp
+ K1Vector& operator-= (const K& k)
+ {
+ p -= k;
+ return *this;
+ }
+
+ //! vector space multiplication with scalar
+ K1Vector& operator*= (const K& k)
+ {
+ p *= k;
+ return *this;
+ }
+
+ //! vector space division by scalar
+ K1Vector& operator/= (const K& k)
+ {
+ p /= k;
+ return *this;
+ }
+
+ //! vector space axpy operation
+ K1Vector& axpy (const K& a, const K1Vector& y)
+ {
+ p += a*y.p;
+ return *this;
+ }
+
+
+ //===== Euclidean scalar product
+
+ //! scalar product
+ const K operator* (const K1Vector& y) const
+ {
+ return p*y.p;
+ }
+
+
+ //===== norms
+
+ //! one norm (sum over absolute values of entries)
+ double one_norm () const
+ {
+ return fvmeta_abs(p);
+ }
+
+ //! simplified one norm (uses Manhattan norm for complex values)
+ double one_norm_real () const
+ {
+ return fvmeta_abs_real(p);
+ }
+
+ //! two norm sqrt(sum over squared values of entries)
+ double two_norm () const
+ {
+ return sqrt(fvmeta_abs2(p));
+ }
+
+ //! sqare of two norm (sum over squared values of entries), need for block recursion
+ double two_norm2 () const
+ {
+ return fvmeta_abs2(p);
+ }
+
+ //! infinity norm (maximum of absolute values of entries)
+ double infinity_norm () const
+ {
+ return fvmeta_abs(p);
+ }
+
+ //! simplified infinity norm (uses Manhattan norm for complex values)
+ double infinity_norm_real () const
+ {
+ return fvmeta_abs_real(p);
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in the vector (are of size 1 here)
+ int N () const
+ {
+ return 1;
+ }
+
+ //! dimension of the vector space
+ int dim () const
+ {
+ return 1;
+ }
+
+ //===== conversion operator
+
+ operator K () {return p;}
+
+ private:
+ // the data
+ K p;
+ };
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/gsetc.hh b/istl/gsetc.hh
new file mode 100644
index 000000000..700b51954
--- /dev/null
+++ b/istl/gsetc.hh
@@ -0,0 +1,191 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_GSETC_HH__
+#define __DUNE_GSETC_HH__
+
+#include <math.h>
+#include <complex>
+#include <iostream>
+#include <iomanip>
+#include <string>
+
+#include "istlexception.hh"
+#include "fvector.hh"
+#include "fmatrix.hh"
+
+/*! \file __FILE__
+
+ Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc.
+ in a generic way
+ */
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ //! return x = x + w D(A)^{-1} (b-Ax)
+ template<class M, class X, class Y, class K>
+ void bjac (M& A, X& x, Y& b, K w)
+ {
+ // iterator types
+ typedef typename M::RowIterator rowiterator;
+ typedef typename M::ColIterator coliterator;
+
+ // cast to field type
+ typename X::field_type omega=w;
+
+ // compute defect
+ Y d(b); // d=b
+ A.mmv(x,d); // d=d-Ax
+
+ // invert block diagonal, overwrite d with update
+ for (rowiterator i=A.begin(); i!=A.end(); ++i)
+ {
+ coliterator j = (*i).find(i.index());
+ (*j).solve(d[i.index()],d[i.index()]);
+ }
+
+ // update
+ x.axpy(omega,d);
+ }
+
+ //! return v = D(A)^{-1} d
+ template<class M, class X, class Y>
+ void bjac_update (M& A, X& v, Y& d)
+ {
+ // iterator types
+ typedef typename M::RowIterator rowiterator;
+ typedef typename M::ColIterator coliterator;
+
+ // invert block diagonal, overwrite d with update
+ rowiterator endi=A.end();
+ for (rowiterator i=A.begin(); i!=endi; ++i)
+ {
+ coliterator j = (*i).find(i.index());
+ (*j).solve(v[i.index()],d[i.index()]);
+ }
+ }
+
+ //! return x = x + L(A)^{-1} (b-Ax), implemented without temporary vector
+ template<class M, class X, class Y>
+ void bgs (M& A, X& x, Y& b)
+ {
+ // iterator types
+ typedef typename M::RowIterator rowiterator;
+ typedef typename M::ColIterator coliterator;
+ typedef typename Y::block_type bblock;
+
+ // local solve at each block and immediate update
+ coliterator diag;
+ for (rowiterator i=A.begin(); i!=A.end(); ++i)
+ {
+ // compute right hand side
+ bblock rhs(b[i.index()]);
+ for (typename M::ColIterator j=(*i).begin(); j!=(*i).end(); ++j)
+ if (j.index()!=i.index())
+ (*j).mmv(x[j.index()],rhs);
+ else
+ diag = j;
+
+ // solve locally
+ (*diag).solve(x[i.index()],rhs);
+ }
+ }
+
+ //! return v = L(A)^{-1} d
+ template<class M, class X, class Y>
+ void bgs_update (M& A, X& v, Y& d)
+ {
+ // iterator types
+ typedef typename M::RowIterator rowiterator;
+ typedef typename M::ColIterator coliterator;
+ typedef typename Y::block_type bblock;
+
+ // local solve at each block and immediate update
+ rowiterator endi=A.end();
+ for (rowiterator i=A.begin(); i!=endi; ++i)
+ {
+ bblock rhs(d[i.index()]);
+ coliterator j;
+ for (j=(*i).begin(); j.index()<i.index(); ++j)
+ (*j).mmv(v[j.index()],rhs);
+ (*j).solve(v[i.index()],rhs);
+ }
+ }
+
+ //! return x = x + w L(A)^{-1} (b-Ax), implemented without temporary vector
+ template<class M, class X, class Y, class K>
+ void bsor (M& A, X& x, Y& b, K w)
+ {
+ // iterator types
+ typedef typename M::RowIterator rowiterator;
+ typedef typename M::ColIterator coliterator;
+ typedef typename X::block_type xblock;
+ typedef typename Y::block_type bblock;
+
+ // cast to field type
+ typename X::field_type omega=w;
+
+ // local solve at each block and immediate update
+ coliterator diag;
+ for (rowiterator i=A.begin(); i!=A.end(); ++i)
+ {
+ // compute right hand side
+ bblock rhs(b[i.index()]);
+ for (typename M::ColIterator j=(*i).begin(); j!=(*i).end(); ++j)
+ if (j.index()!=i.index())
+ (*j).mmv(x[j.index()],rhs);
+ else
+ diag = j;
+
+ // solve locally
+ xblock xnew;
+ (*diag).solve(xnew,rhs);
+
+ // damped update
+ x[i.index()] *= 1-omega;
+ x[i.index()].axpy(omega,xnew);
+ }
+ }
+
+ //! return v = w L(A)^{-1} d
+ template<class M, class X, class Y, class K>
+ void bsor_update (M& A, X& v, Y& d, K w)
+ {
+ // iterator types
+ typedef typename M::RowIterator rowiterator;
+ typedef typename M::ColIterator coliterator;
+ typedef typename Y::block_type bblock;
+
+ // cast to field type
+ typename X::field_type omega=w;
+
+ // local solve at each block and immediate update
+ coliterator diag;
+ for (rowiterator i=A.begin(); i!=A.end(); ++i)
+ {
+ // compute right hand side
+ bblock rhs(d[i.index()]);
+ for (typename M::ColIterator j=(*i).begin(); j.index()<=i.index(); ++j)
+ if (j.index()<i.index())
+ (*j).mmv(v[j.index()],rhs);
+ else
+ diag = j;
+
+ // solve locally
+ (*diag).solve(v[i.index()],rhs);
+
+ // damp immediately
+ v[i.index] *= omega;
+ }
+ }
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/io.hh b/istl/io.hh
new file mode 100644
index 000000000..1959aae88
--- /dev/null
+++ b/istl/io.hh
@@ -0,0 +1,190 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_ISTLIO_HH__
+#define __DUNE_ISTLIO_HH__
+
+#include <math.h>
+#include <complex>
+#include <iostream>
+#include <iomanip>
+#include <string>
+
+#include "istlexception.hh"
+#include "fvector.hh"
+#include "fmatrix.hh"
+
+/*! \file __FILE__
+
+ Some generic functions for pretty printing vectors and matrices
+ */
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+
+ //==== pretty printing of vectors
+
+ // recursively print all the blocks
+ template<class V>
+ void recursive_printvector (std::ostream& s, V& v, std::string rowtext, int& counter,
+ int columns, int width, int precision)
+ {
+ for (typename V::Iterator i=v.begin(); i!=v.end(); ++i)
+ recursive_printvector(s,*i,rowtext,counter,columns,width,precision);
+ }
+
+ // specialization for FieldVector
+ template<class K, int n>
+ void recursive_printvector (std::ostream& s, FieldVector<K,n>& v, std::string rowtext, int& counter,
+ int columns, int width, int precision)
+ {
+ // we now can print n numbers
+ for (int i=0; i<n; i++)
+ {
+ if (counter%columns==0)
+ {
+ s << rowtext; // start a new row
+ s << " "; // space in front of each entry
+ s.width(4); // set width for counter
+ s << counter; // number of first entry in a line
+ }
+ s << " "; // space in front of each entry
+ s.width(width); // set width for each entry anew
+ s << v[i]; // yeah, the number !
+ counter++; // increment the counter
+ if (counter%columns==0)
+ s << std::endl; // start a new line
+ }
+ }
+
+
+ template<class V>
+ void printvector (std::ostream& s, V& v, std::string title, std::string rowtext,
+ int columns=1, int width=10, int precision=2)
+ {
+ // count the numbers printed to make columns
+ int counter=0;
+
+ // set the output format
+ s.setf(std::ios_base::scientific, std::ios_base::floatfield);
+ int oldprec = s.precision();
+ s.precision(precision);
+
+ // print title
+ s << title << " [blocks=" << v.N() << ",dimension=" << v.dim() << "]" << std::endl;
+
+ // print data from all blocks
+ recursive_printvector(s,v,rowtext,counter,columns,width,precision);
+
+ // check if new line is required
+ if (counter%columns!=0)
+ s << std::endl;
+
+ // reset the output format
+ s.precision(oldprec);
+ s.setf(std::ios_base::fixed, std::ios_base::floatfield);
+ }
+
+
+ //==== pretty printing of matrices
+
+
+ //! print a row of zeros for a non-existing block
+ void fill_row (std::ostream& s, int m, int width, int precision)
+ {
+ for (int j=0; j<m; j++)
+ {
+ s << " "; // space in front of each entry
+ s.width(width); // set width for each entry anew
+ s << "0"; // yeah, the number !
+ }
+ }
+
+ //! print one row of a matrix
+ template<class K, int n, int m>
+ void print_row (std::ostream& s, FieldMatrix<K,n,m>& A, int I, int J, int therow, int width, int precision)
+ {
+ for (int i=0; i<n; i++)
+ if (I+i==therow)
+ for (int j=0; j<m; j++)
+ {
+ s << " "; // space in front of each entry
+ s.width(width); // set width for each entry anew
+ s << A[i][j]; // yeah, the number !
+ }
+ }
+
+ //! print one row of a matrix
+ template<class M>
+ void print_row (std::ostream& s, M& A, int I, int J, int therow, int width, int precision)
+ {
+ int i0=I;
+ for (int i=0; i<A.N(); i++)
+ {
+ if (therow>=i0 && therow<i0+A.rowdim(i))
+ {
+ // the row is in this block row !
+ int j0=J;
+ for (int j=0; j<A.M(); j++)
+ {
+ // find this block
+ typename M::ColIterator it = A[i].find(j);
+
+ // print row or filler
+ if (it!=A[i].end())
+ print_row(s,*it,i0,j0,therow,width,precision);
+ else
+ fill_row(s,A.coldim(j),width,precision);
+
+ // advance columns
+ j0 += A.coldim(j);
+ }
+ }
+ // advance rows
+ i0 += A.rowdim(i);
+ }
+ }
+
+ template<class M>
+ void printmatrix (std::ostream& s, M& A, std::string title, std::string rowtext,
+ int width=10, int precision=2)
+ {
+ // set the output format
+ s.setf(std::ios_base::scientific, std::ios_base::floatfield);
+ int oldprec = s.precision();
+ s.precision(precision);
+
+ // print title
+ s << title
+ << " [n=" << A.N()
+ << ",m=" << A.M()
+ << ",rowdim=" << A.rowdim()
+ << ",coldim=" << A.coldim()
+ << "]" << std::endl;
+
+ // print all rows
+ for (int i=0; i<A.rowdim(); i++)
+ {
+ s << rowtext; // start a new row
+ s << " "; // space in front of each entry
+ s.width(4); // set width for counter
+ s << i; // number of first entry in a line
+ print_row(s,A,0,0,i,width,precision); // generic print
+ s << std::endl; // start a new line
+ }
+
+ // reset the output format
+ s.precision(oldprec);
+ s.setf(std::ios_base::fixed, std::ios_base::floatfield);
+ }
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/istl.hh b/istl/istl.hh
new file mode 100644
index 000000000..00e0db6ba
--- /dev/null
+++ b/istl/istl.hh
@@ -0,0 +1,38 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_ISTL_HH__
+#define __DUNE_ISTL_HH__
+
+/** @defgroup ISTL Dune Sparse Matrix Vector Templates
+
+ The SPMV module defines a set of classes for sparse
+ matrices and vectors with the following features:
+
+ - They are efficient (including access to single elements)
+
+ - They provide a recursive block structure for a flexible
+ construction of iterative linear solvers
+
+ - They are very flexible
+ */
+
+// ISTL includes
+#include "fvector.hh"
+#include "bvector.hh"
+#include "vbvector.hh"
+#include "fmatrix.hh"
+
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/istlexception.hh b/istl/istlexception.hh
new file mode 100644
index 000000000..b36dea778
--- /dev/null
+++ b/istl/istlexception.hh
@@ -0,0 +1,25 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_ISTLEXC_HH__
+#define __DUNE_ISTLEXC_HH__
+
+#include <stdlib.h>
+
+#include "../common/exceptions.hh"
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ //! derive error class from the base class in common
+ class ISTLError : public Exception {};
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/precision.hh b/istl/precision.hh
new file mode 100644
index 000000000..75da98365
--- /dev/null
+++ b/istl/precision.hh
@@ -0,0 +1,69 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_PRECISION_HH__
+#define __DUNE_PRECISION_HH__
+
+#include <stdlib.h>
+
+namespace Dune {
+
+ /** @defgroup ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ class ISTLPrecision { // uses standard malloc and free
+ public:
+ //! return threshold to do pivoting
+ static double pivoting_limit ()
+ {
+ return _pivoting;
+ }
+
+ //! set pivoting threshold
+ static void set_pivoting_limit (double pivthres)
+ {
+ _pivoting = pivthres;
+ }
+
+ //! return threshold to declare matrix singular
+ static double singular_limit ()
+ {
+ return _singular;
+ }
+
+ //! set singular threshold
+ static void set_singular_limit (double singthres)
+ {
+ _singular = singthres;
+ }
+
+ //! return threshold to declare matrix singular
+ static double absolute_limit ()
+ {
+ return _absolute;
+ }
+
+ //! set singular threshold
+ static void set_absolute_limit (double absthres)
+ {
+ _absolute = absthres;
+ }
+
+ private:
+ // just to demonstrate some state information
+ static double _pivoting;
+ static double _singular;
+ static double _absolute;
+ };
+
+ double ISTLPrecision::_pivoting = 1E-8;
+ double ISTLPrecision::_singular = 1E-14;
+ double ISTLPrecision::_absolute = 1E-80;
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/tutorial/example.cc b/istl/tutorial/example.cc
index 60fdeb47d..0bd4059d6 100644
--- a/istl/tutorial/example.cc
+++ b/istl/tutorial/example.cc
@@ -92,8 +92,28 @@ void test_basearray ()
c.move(2,5);
}
+template<class V>
+void f (V& v)
+{
+ typedef typename V::Iterator iterator;
+ for (iterator i=v.begin(); i!=v.end(); ++i)
+ *i = i.index();
+
+ typedef typename V::ConstIterator const_iterator;
+ for (const_iterator i=v.begin(); i!=v.end(); ++i)
+ std::cout << (*i).two_norm() << std::endl;
+}
+
void test_BlockVector ()
{
+ Dune::BlockVector<Dune::FieldVector<std::complex<double>,2> > v(20);
+
+ v[1] = 3.14;
+ v[3][0] = 2.56;
+ v[3][1] = std::complex<double>(1,-1);
+
+ f(v);
+
typedef Dune::FieldVector<double,1> R1;
const int n=480;
@@ -323,9 +343,8 @@ void test_Iter ()
{
// block types
Timer t;
- const int BlockSize=1;
- typedef Dune::FieldVector<double,BlockSize> VB;
- typedef Dune::FieldMatrix<double,BlockSize,BlockSize> MB;
+ typedef Dune::K1Vector<double> VB;
+ typedef Dune::K11Matrix<double> MB;
// a fake discretization
t.start();
@@ -364,9 +383,9 @@ void test_Iter ()
Dune::BlockVector<VB> v(x); // memory for update
v = 1.23E-4;
double w=1.0; // damping factor
- for (int k=1; k<=50; k++)
+ for (int k=1; k<=10; k++)
{
- bgs_update(A,v,d); // compute update
+ bgs_update_lowlevel(A,v,d); // compute update
x.axpy(w,v); // update solution
A.usmv(-w,v,d); // update defect
std::cout << k << " " << d.two_norm() << std::endl;
@@ -376,20 +395,6 @@ void test_Iter ()
return;
// printvector(std::cout,x,"solution","entry",1,10,2);
-
- // a simple iteration
- x = 0; // initial guess
- d=b; A.mmv(x,d);
- std::cout.setf(std::ios_base::scientific, std::ios_base::floatfield);
- std::cout.precision(8);
- std::cout << 0 << " " << d.two_norm() << std::endl;
- for (int k=1; k<=50; k++)
- {
- bsor(A,x,b,1.0);
- d=b; A.mmv(x,d);
- std::cout << k << " " << d.two_norm() << std::endl;
- }
- // printvector(std::cout,x,"solution","entry",1,10,2);
}
int main (int argc , char ** argv)
@@ -400,7 +405,7 @@ int main (int argc , char ** argv)
// test_VariableBlockVector();
// test_FieldMatrix();
// test_BCRSMatrix();
- // test_IO();
+ // test_IO();
test_Iter();
}
catch (Dune::ISTLError& error)
diff --git a/istl/vbcrsmatrix.hh b/istl/vbcrsmatrix.hh
new file mode 100644
index 000000000..02e53c5e9
--- /dev/null
+++ b/istl/vbcrsmatrix.hh
@@ -0,0 +1,31 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_VBCRSMATRIX_HH__
+#define __DUNE_VBCRSMATRIX_HH__
+
+#include <math.h>
+#include <complex>
+
+#include "istlexection.hh"
+#include "allocator.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 ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
diff --git a/istl/vbvector.hh b/istl/vbvector.hh
new file mode 100644
index 000000000..6ecd45f57
--- /dev/null
+++ b/istl/vbvector.hh
@@ -0,0 +1,701 @@
+// -*- 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 <iostream>
+
+#include "istlexception.hh"
+#include "allocator.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 ISTL Iterative Solvers Template Library
+ @addtogroup ISTL
+ @{
+ */
+
+ /** implements a vector consisting of a number of blocks (to
+ be given at run-time) which themselves consist of a number
+ of blocks (also given at run-time) of the given type B.
+
+ VariableBlockVector is a container of containers!
+
+ */
+ template<class B, class A=ISTLAllocator>
+ class VariableBlockVector : public block_vector_unmanaged<B,A>
+ // this derivation gives us all the blas level 1 and norms
+ // on the large array. However, access operators have to be
+ // overwritten.
+ {
+ 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, note that this
+ is *not* the type refered to by the iterators and random access.
+ However, it can be used to copy blocks (which is its only purpose).
+ */
+ typedef BlockVector<B,A> block_type;
+
+ /** increment block level counter, yes, it is two levels because
+ VariableBlockVector is a container of containers
+ */
+ enum {blocklevel = B::blocklevel+2};
+
+ // just a shorthand
+ typedef BlockVectorWindow<B,A> window_type;
+
+
+ //===== constructors and such
+
+ /** constructor without arguments makes empty vector,
+ object cannot be used yet
+ */
+ VariableBlockVector () : block_vector_unmanaged<B,A>()
+ {
+ // nothing is known ...
+ nblocks = 0;
+ block = 0;
+ initialized = false;
+ }
+
+ /** make vector with given number of blocks, but size of each block is not yet known,
+ object cannot be used yet
+ */
+ VariableBlockVector (int _nblocks) : block_vector_unmanaged<B,A>()
+ {
+ // we can allocate the windows now
+ nblocks = _nblocks;
+ if (nblocks>0)
+ {
+ block = A::template malloc<window_type>(nblocks);
+ }
+ else
+ {
+ nblocks = 0;
+ block = 0;;
+ }
+
+ // Note: memory in base class still not allocated
+ // the vector not usable
+ initialized = false;
+ }
+
+ /** make vector with given number of blocks each having a constant size,
+ object is fully usable then
+ */
+ VariableBlockVector (int _nblocks, int m) : block_vector_unmanaged<B,A>()
+ {
+ // and we can allocate the big array in the base class
+ this->n = _nblocks*m;
+ if (this->n>0)
+ {
+ this->p = A::template malloc<B>(this->n);
+ }
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ // we can allocate the windows now
+ nblocks = _nblocks;
+ if (nblocks>0)
+ {
+ // alloc
+ block = A::template malloc<window_type>(nblocks);
+
+ // set the windows into the big array
+ for (int i=0; i<nblocks; ++i)
+ block[i].set(m,this->p+(i*m));
+ }
+ else
+ {
+ nblocks = 0;
+ block = 0;;
+ }
+
+ // and the vector is usable
+ initialized = true;
+ }
+
+ //! copy constructor, has copy semantics
+ VariableBlockVector (const VariableBlockVector& a)
+ {
+ // allocate the big array in the base class
+ this->n = a.n;
+ if (this->n>0)
+ {
+ // alloc
+ this->p = A::template malloc<B>(this->n);
+
+ // copy data
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ // we can allocate the windows now
+ nblocks = a.nblocks;
+ if (nblocks>0)
+ {
+ // alloc
+ block = A::template malloc<window_type>(nblocks);
+
+ // and we must set the windows
+ block[0].set(a.block[0].getsize(),this->p); // first block
+ for (int i=1; i<nblocks; ++i) // and the rest
+ block[i].set(a.block[i].getsize(),block[i-1].getptr()+block[i-1].getsize());
+ }
+ else
+ {
+ nblocks = 0;
+ block = 0;;
+ }
+
+ // and we have a usable vector
+ initialized = true;
+ }
+
+ //! free dynamic memory
+ ~VariableBlockVector ()
+ {
+ if (this->n>0) A::template free<B>(this->p);
+ if (nblocks>0) A::template free<window_type>(block);
+ }
+
+
+ //! same effect as constructor with same argument
+ void resize (int _nblocks)
+ {
+ // deallocate memory if necessary
+ if (this->n>0) A::template free<B>(this->p);
+ if (nblocks>0) A::template free<window_type>(block);
+ this->n = 0;
+ this->p = 0;
+
+ // we can allocate the windows now
+ nblocks = _nblocks;
+ if (nblocks>0)
+ {
+ block = A::template malloc<window_type>(nblocks);
+ }
+ else
+ {
+ nblocks = 0;
+ block = 0;;
+ }
+
+ // and the vector not fully usable
+ initialized = false;
+ }
+
+ //! same effect as constructor with same argument
+ void resize (int _nblocks, int m)
+ {
+ // deallocate memory if necessary
+ if (this->n>0) A::template free<B>(this->p);
+ if (nblocks>0) A::template free<window_type>(block);
+
+ // and we can allocate the big array in the base class
+ this->n = _nblocks*m;
+ if (this->n>0)
+ {
+ this->p = A::template malloc<B>(this->n);
+ }
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ // we can allocate the windows now
+ nblocks = _nblocks;
+ if (nblocks>0)
+ {
+ // alloc
+ block = A::template malloc<window_type>(nblocks);
+
+ // set the windows into the big array
+ for (int i=0; i<nblocks; ++i)
+ block[i].set(m,this->p+(i*m));
+ }
+ else
+ {
+ nblocks = 0;
+ block = 0;;
+ }
+
+ // and the vector is usable
+ initialized = true;
+ }
+
+ //! assignment
+ VariableBlockVector& operator= (const VariableBlockVector& a)
+ {
+ if (&a!=this) // check if this and a are different objects
+ {
+ // reallocate arrays if necessary
+ // Note: still the block sizes may vary !
+ if (this->n!=a.n || nblocks!=a.nblocks)
+ {
+ // deallocate memory if necessary
+ if (this->n>0) A::template free<B>(this->p);
+ if (nblocks>0) A::template free<window_type>(block);
+
+ // allocate the big array in the base class
+ this->n = a.n;
+ if (this->n>0)
+ {
+ // alloc
+ this->p = A::template malloc<B>(this->n);
+ }
+ else
+ {
+ this->n = 0;
+ this->p = 0;
+ }
+
+ // we can allocate the windows now
+ nblocks = a.nblocks;
+ if (nblocks>0)
+ {
+ // alloc
+ block = A::template malloc<window_type>(nblocks);
+ }
+ else
+ {
+ nblocks = 0;
+ block = 0;;
+ }
+ }
+
+ // copy block structure, might be different although
+ // sizes are the same !
+ if (nblocks>0)
+ {
+ block[0].set(a.block[0].getsize(),this->p); // first block
+ for (int i=1; i<nblocks; ++i) // and the rest
+ block[i].set(a.block[i].getsize(),block[i-1].getptr()+block[i-1].getsize());
+ }
+
+ // and copy the data
+ for (int i=0; i<this->n; i++) this->p[i]=a.p[i];
+ }
+
+ // and we have a usable vector
+ initialized = true;
+
+ return *this; // Gebe Referenz zurueck damit a=b=c; klappt
+ }
+
+
+ //===== assignment from scalar
+
+ //! assign from scalar
+ VariableBlockVector& operator= (const field_type& k)
+ {
+ (static_cast<block_vector_unmanaged<B,A>&>(*this)) = k;
+ return *this;
+ }
+
+
+ //===== the creation interface
+
+ //! Iterator class for sequential creation of blocks
+ class CreateIterator
+ {
+ public:
+ //! constructor
+ CreateIterator (VariableBlockVector& _v, int _i) : v(_v)
+ {
+ i = _i;
+ k = 0;
+ n = 0;
+ }
+
+ //! prefix increment
+ CreateIterator& operator++()
+ {
+ // we are at block i and the blocks size is known
+
+ // set the blocks size to current k
+ v.block[i].setsize(k);
+
+ // accumulate total size
+ n += k;
+
+ // go to next block
+ ++i;
+
+ // reset block size
+ k = 0;
+
+ // if we are past the last block, finish off
+ if (i==v.nblocks)
+ {
+ // now we can allocate the big array in the base class of v
+ v.n = n;
+ if (n>0)
+ {
+ // alloc
+ v.p = A::template malloc<B>(n);
+ }
+ else
+ {
+ v.n = 0;
+ v.p = 0;
+ }
+
+ // and we set the window pointer
+ if (v.nblocks>0)
+ {
+ v.block[0].setptr(v.p); // pointer tofirst block
+ for (int j=1; j<v.nblocks; ++j) // and the rest
+ v.block[j].setptr(v.block[j-1].getptr()+v.block[j-1].getsize());
+ }
+
+ // and the vector is ready
+ v.initialized = true;
+
+ //std::cout << "made vbvector with " << v.n << " components" << std::endl;
+ }
+
+ 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 index ()
+ {
+ return i;
+ }
+
+ //! set size of current block
+ void setblocksize (int _k)
+ {
+ k = _k;
+ }
+
+ private:
+ VariableBlockVector& v; // my vector
+ int i; // current block to be defined
+ int k; // size of current block to be defined
+ int n; // total number of elements to be allocated
+ };
+
+ // 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);
+ }
+
+
+ //===== access to components
+ // has to be overwritten from base class because it must
+ // return access to the windows
+
+ //! random access to blocks
+ window_type& operator[] (int i)
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (i<0 || i>=nblocks) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ return block[i];
+ }
+
+ //! same for read only access
+ const window_type& operator[] (int i) const
+ {
+#ifdef DUNE_ISTL_WITH_CHECKING
+ if (i<0 || i>=nblocks) DUNE_THROW(ISTLError,"index out of range");
+#endif
+ return block[i];
+ }
+
+ // forward declaration
+ class ConstIterator;
+
+ //! Iterator class for sequential access
+ class Iterator
+ {
+ public:
+ //! constructor, no arguments
+ Iterator ()
+ {
+ p = 0;
+ i = 0;
+ }
+
+ //! constructor
+ Iterator (window_type* _p, int _i) : p(_p), i(_i)
+ { }
+
+ //! prefix increment
+ Iterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ Iterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! equality
+ bool operator== (const ConstIterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const ConstIterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ window_type& operator* ()
+ {
+ return p[i];
+ }
+
+ //! arrow
+ window_type* operator-> ()
+ {
+ return p+i;
+ }
+
+ // return index corresponding to pointer
+ int index () const
+ {
+ return i;
+ }
+
+ friend class ConstIterator;
+
+ private:
+ window_type* p;
+ int i;
+ };
+
+ //! begin Iterator
+ Iterator begin ()
+ {
+ return Iterator(block,0);
+ }
+
+ //! end Iterator
+ Iterator end ()
+ {
+ return Iterator(block,nblocks);
+ }
+
+ //! begin Iterator
+ Iterator rbegin ()
+ {
+ return Iterator(block,nblocks-1);
+ }
+
+ //! end Iterator
+ Iterator rend ()
+ {
+ return Iterator(block,-1);
+ }
+
+ //! random access returning iterator (end if not contained)
+ Iterator find (int i)
+ {
+ if (i>=0 && i<nblocks)
+ return Iterator(block,i);
+ else
+ return Iterator(block,nblocks);
+ }
+
+ //! ConstIterator class for sequential access
+ class ConstIterator
+ {
+ public:
+ //! constructor
+ ConstIterator ()
+ {
+ p = 0;
+ i = 0;
+ }
+
+ //! constructor from pointer
+ ConstIterator (const window_type* _p, int _i) : p(_p), i(_i)
+ { }
+
+ //! constructor from non_const iterator
+ ConstIterator (const Iterator& it) : p(it.p), i(it.i)
+ { }
+
+ //! prefix increment
+ ConstIterator& operator++()
+ {
+ ++i;
+ return *this;
+ }
+
+ //! prefix decrement
+ ConstIterator& operator--()
+ {
+ --i;
+ return *this;
+ }
+
+ //! equality
+ bool operator== (const ConstIterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const ConstIterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! equality
+ bool operator== (const Iterator& it) const
+ {
+ return (p+i)==(it.p+it.i);
+ }
+
+ //! inequality
+ bool operator!= (const Iterator& it) const
+ {
+ return (p+i)!=(it.p+it.i);
+ }
+
+ //! dereferencing
+ const window_type& operator* () const
+ {
+ return p[i];
+ }
+
+ //! arrow
+ const window_type* operator-> () const
+ {
+ return p+i;
+ }
+
+ // return index corresponding to pointer
+ int index () const
+ {
+ return i;
+ }
+
+ friend class Iterator;
+
+ private:
+ const window_type* p;
+ int i;
+ };
+
+ //! begin ConstIterator
+ ConstIterator begin () const
+ {
+ return ConstIterator(block,0);
+ }
+
+ //! end ConstIterator
+ ConstIterator end () const
+ {
+ return ConstIterator(block,nblocks);
+ }
+
+ //! begin ConstIterator
+ ConstIterator rbegin () const
+ {
+ return ConstIterator(block,nblocks-1);
+ }
+
+ //! end ConstIterator
+ ConstIterator rend () const
+ {
+ return ConstIterator(block,-1);
+ }
+
+
+ //===== sizes
+
+ //! number of blocks in the vector (are of variable size here)
+ int N () const
+ {
+ return nblocks;
+ }
+
+
+ private:
+ int nblocks; // number of blocks in vector
+ window_type* block; // array of blocks pointing to the array in the base class
+ bool initialized; // true if vector has been initialized
+ };
+
+
+
+ /** @} end documentation */
+
+} // end namespace
+
+#endif
--
GitLab