From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Mon, 2 Dec 2002 17:32:38 +0000 (+0000)
Subject: Should not have deleted this file.
X-Git-Tag: v8.0.0~17192
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=45348165d5d8a564b309f46a1edd496a8b573196;p=dealii.git

Should not have deleted this file.


git-svn-id: https://svn.dealii.org/trunk@6794 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h
new file mode 100644
index 0000000000..29c71926bc
--- /dev/null
+++ b/deal.II/lac/include/lac/sparse_matrix.h
@@ -0,0 +1,1491 @@
+//----------------------------  sparse_matrix.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_matrix.h  ---------------------------
+#ifndef __deal2__sparse_matrix_h
+#define __deal2__sparse_matrix_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+#include <lac/sparsity_pattern.h>
+
+template<typename number> class Vector;
+template<typename number> class FullMatrix;
+
+/**
+ * Sparse matrix.
+ *
+ *
+ * @sect2{On template instantiations}
+ *
+ * Member functions of this class are either implemented in this file
+ * or in a file of the same name with suffix ``.templates.h''. For the
+ * most common combinations of the template parameters, instantiations
+ * of this class are provided in a file with suffix ``.cc'' in the
+ * ``source'' directory. If you need an instantiation that is not
+ * listed there, you have to include this file along with the
+ * corresponding ``.templates.h'' file and instantiate the respective
+ * class yourself.
+ *
+ * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998
+ */
+template <typename number>
+class SparseMatrix : public Subscriptor
+{
+  public:
+				     /**
+				      * Accessor class for iterators
+				      */
+    class Accessor
+    {
+      public:
+					 /**
+					  * Constructor. Since we use
+					  * accessors only for read
+					  * access, a const matrix
+					  * pointer is sufficient.
+					  */
+	Accessor (const SparseMatrix<number>*,
+		  unsigned int row,
+		  unsigned short index);
+
+					 /**
+					  * Row number of the element
+					  * represented by this
+					  * object.
+					  */
+	unsigned int row() const;
+
+					 /**
+					  * Index in row of the element
+					  * represented by this
+					  * object.
+					  */
+	unsigned short index() const;
+
+					 /**
+					  * Column number of the
+					  * element represented by
+					  * this object.
+					  */
+	unsigned int column() const;
+
+					 /**
+					  * Value of this matrix entry.
+					  */
+	number value() const;
+	
+	protected:
+					 /**
+					  * The matrix accessed.
+					  */
+	const SparseMatrix<number>* matrix;
+
+					 /**
+					  * Current row number.
+					  */
+	unsigned int a_row;
+
+					 /**
+					  * Current index in row.
+					  */
+	unsigned short a_index;
+      };
+
+				     /**
+				      * STL conforming iterator.
+				      */
+    class const_iterator : private Accessor
+      {
+	public:
+					   /**
+					    * Constructor.
+					    */ 
+	const_iterator(const SparseMatrix<number>*,
+		       unsigned int row,
+		       unsigned short index);
+	  
+					   /**
+					    * Prefix increment.
+					    */
+	const_iterator& operator++ ();
+
+					   /**
+					    * Postfix increment.
+					    */
+	const_iterator& operator++ (int);
+
+					   /**
+					    * Dereferencing operator.
+					    */
+	const Accessor& operator* () const;
+
+					   /**
+					    * Dereferencing operator.
+					    */
+	const Accessor* operator-> () const;
+
+					   /**
+					    * Comparison. True, if
+					    * both iterators point to
+					    * the same matrix
+					    * position.
+					    */
+	bool operator == (const const_iterator&) const;
+					   /**
+					    * Inverse of @p{==}.
+					    */
+	bool operator != (const const_iterator&) const;
+
+					   /**
+					    * Comparison
+					    * operator. Result is true
+					    * if either the first row
+					    * number is smaller or if
+					    * the row numbers are
+					    * equal and the first
+					    * index is smaller.
+					    */
+	bool operator < (const const_iterator&) const;
+      };
+    
+				     /**
+				      * Type of matrix entries. In analogy to
+				      * the STL container classes.
+				      */
+    typedef number value_type;
+    
+				     /**
+				      * Constructor; initializes the matrix to
+				      * be empty, without any structure, i.e.
+				      * the matrix is not usable at all. This
+				      * constructor is therefore only useful
+				      * for matrices which are members of a
+				      * class. All other matrices should be
+				      * created at a point in the data flow
+				      * where all necessary information is
+				      * available.
+				      *
+				      * You have to initialize
+				      * the matrix before usage with
+				      * @p{reinit(SparsityPattern)}.
+				      */
+    SparseMatrix ();
+
+				     /**
+				      * Copy constructor. This constructor is
+				      * only allowed to be called if the matrix
+				      * to be copied is empty. This is for the
+				      * same reason as for the
+				      * @p{SparsityPattern}, see there for the
+				      * details.
+				      *
+				      * If you really want to copy a whole
+				      * matrix, you can do so by using the
+				      * @p{copy_from} function.
+				      */
+    SparseMatrix (const SparseMatrix &);
+
+				     /**
+				      * Constructor. Takes the given
+				      * matrix sparsity structure to
+				      * represent the sparsity pattern
+				      * of this matrix. You can change
+				      * the sparsity pattern later on
+				      * by calling the @p{reinit}
+				      * function.
+				      *
+				      * You have to make sure that the
+				      * lifetime of the sparsity
+				      * structure is at least as long
+				      * as that of this matrix or as
+				      * long as @p{reinit} is not
+				      * called with a new sparsity
+				      * structure.
+				      *
+				      * The constructor is marked
+				      * explicit so as to disallow
+				      * that someone passes a sparsity
+				      * pattern in place of a sparse
+				      * matrix to some function, where
+				      * an empty matrix would be
+				      * generated then.
+				      */
+    explicit SparseMatrix (const SparsityPattern &sparsity);
+    
+				     /**
+				      * Destructor. Free all memory, but do not
+				      * release the memory of the sparsity
+				      * structure.
+				      */
+    virtual ~SparseMatrix ();
+
+				     /** 
+				      * Pseudo operator only copying
+				      * empty objects.
+				      */
+    SparseMatrix<number>& operator = (const SparseMatrix<number> &);
+
+				     /**
+				      * Reinitialize the object but
+				      * keep to the sparsity pattern
+				      * previously used.  This may be
+				      * necessary if you @p{reinit}'d
+				      * the sparsity structure and
+				      * want to update the size of the
+				      * matrix.
+				      *
+				      * Note that memory is only
+				      * reallocated if the new size
+				      * exceeds the old size. If that
+				      * is not the case, the allocated
+				      * memory is not reduced. However,
+				      * if the sparsity structure is
+				      * empty (i.e. the dimensions are
+				      * zero), then all memory is
+				      * freed.
+				      *
+				      * If the sparsity pattern has
+				      * not changed, then the effect
+				      * of this function is simply to
+				      * reset all matrix entries to
+				      * zero.
+				      */
+    virtual void reinit ();
+
+				     /**
+				      * Reinitialize the sparse matrix
+				      * with the given sparsity
+				      * pattern. The latter tells the
+				      * matrix how many nonzero
+				      * elements there need to be
+				      * reserved.
+				      *
+				      * Regarding memory allocation,
+				      * the same applies as said
+				      * above.
+				      *
+				      * You have to make sure that the
+				      * lifetime of the sparsity
+				      * structure is at least as long
+				      * as that of this matrix or as
+				      * long as @p{reinit} is not called
+				      * with a new sparsity structure.
+				      *
+				      * The elements of the matrix are
+				      * set to zero by this function.
+				      */
+    virtual void reinit (const SparsityPattern &sparsity);
+
+				     /**
+				      * Release all memory and return
+				      * to a state just like after
+				      * having called the default
+				      * constructor. It also forgets
+				      * the sparsity pattern it was
+				      * previously tied to.
+				      */
+    virtual void clear ();
+    
+				     /**
+				      * Return whether the object is
+				      * empty. It is empty if either
+				      * both dimensions are zero or no
+				      * @p{SparsityPattern} is
+				      * associated.
+				      */
+    bool empty () const;
+
+				     /**
+				      * Return the dimension of the
+				      * image space.  To remember: the
+				      * matrix is of dimension
+				      * $m \times n$.
+				      */
+    unsigned int m () const;
+    
+				     /**
+				      * Return the dimension of the
+				      * range space.  To remember: the
+				      * matrix is of dimension
+				      * $m \times n$.
+				      */
+    unsigned int n () const;
+
+				     /**
+				      * Return the number of nonzero
+				      * elements of this
+				      * matrix. Actually, it returns
+				      * the number of entries in the
+				      * sparsity pattern; if any of
+				      * the entries should happen to
+				      * be zero, it is counted anyway.
+				      */
+    unsigned int n_nonzero_elements () const;
+
+				     /**
+				      * Return the number of actually
+				      * nonzero elements of this
+				      * matrix.
+				      *
+				      * Note, that this function does
+				      * (in contrary to the
+				      * @p{n_nonzero_elements}) NOT
+				      * count all entries of the
+				      * sparsity pattern but only the
+				      * ones that are nonzero.
+				      */
+    unsigned int n_actually_nonzero_elements () const;
+    
+				     /**
+				      * Set the element @p{(i,j)} to @p{value}.
+				      * Throws an error if the entry does
+				      * not exist. Still, it is allowed to store
+				      * zero values in non-existent fields.
+				      */
+    void set (const unsigned int i, const unsigned int j,
+	      const number value);
+    
+				     /**
+				      * Add @p{value} to the element
+				      * @p{(i,j)}.  Throws an error if
+				      * the entry does not
+				      * exist. Still, it is allowed to
+				      * store zero values in
+				      * non-existent fields.
+				      */
+    void add (const unsigned int i, const unsigned int j,
+	      const number value);
+
+				     /**
+				      * Symmetrize the matrix by
+				      * forming the mean value between
+				      * the existing matrix and its
+				      * transpose, $A = \frac 12(A+A^T)$.
+				      *
+				      * This operation assumes that
+				      * the underlying sparsity
+				      * pattern represents a symmetric
+				      * object. If this is not the
+				      * case, then the result of this
+				      * operation will not be a
+				      * symmetric matrix, since it
+				      * only explicitly symmetrizes
+				      * by looping over the lower left
+				      * triangular part for efficiency
+				      * reasons; if there are entries
+				      * in the upper right triangle,
+				      * then these elements are missed
+				      * in the
+				      * symmetrization. Symmetrization
+				      * of the sparsity pattern can be
+				      * obtain by the
+				      * @ref{SparsityPattern}@p{::symmetrize}
+				      * function.
+				      */
+    void symmetrize ();
+    
+				     /**
+				      * Copy the given matrix to this
+				      * one.  The operation throws an
+				      * error if the sparsity patterns
+				      * of the two involved matrices
+				      * do not point to the same
+				      * object, since in this case the
+				      * copy operation is
+				      * cheaper. Since this operation
+				      * is notheless not for free, we
+				      * do not make it available
+				      * through @p{operator =}, since
+				      * this may lead to unwanted
+				      * usage, e.g. in copy arguments
+				      * to functions, which should
+				      * really be arguments by
+				      * reference.
+				      *
+				      * The source matrix may be a matrix
+				      * of arbitrary type, as long as its
+				      * data type is convertible to the
+				      * data type of this matrix.
+				      *
+				      * The function returns a reference to
+				      * @p{this}.
+				      */
+    template <typename somenumber>
+    SparseMatrix<number> &
+    copy_from (const SparseMatrix<somenumber> &source);
+
+				     /**
+				      * This function is complete
+				      * analogous to the
+				      * @ref{SparsityPattern}@p{::copy_from}
+				      * function in that it allows to
+				      * initialize a whole matrix in
+				      * one step. See there for more
+				      * information on argument types
+				      * and their meaning. You can
+				      * also find a small example on
+				      * how to use this function
+				      * there.
+				      *
+				      * The only difference to the
+				      * cited function is that the
+				      * objects which the inner
+				      * iterator points to need to be
+				      * of type @p{std::pair<unsigned int, value},
+				      * where @p{value}
+				      * needs to be convertible to the
+				      * element type of this class, as
+				      * specified by the @p{number}
+				      * template argument.
+				      *
+				      * Previous content of the matrix
+				      * is overwritten. Note that the
+				      * entries specified by the input
+				      * parameters need not
+				      * necessarily cover all elements
+				      * of the matrix. Elements not
+				      * covered remain untouched.
+				      */
+    template <typename ForwardIterator>
+    void copy_from (const ForwardIterator begin,
+		    const ForwardIterator end);    
+
+				     /**
+				      * Copy the nonzero entries of a
+				      * full matrix into this
+				      * object. Previous content is
+				      * deleted. Note that the
+				      * underlying sparsity pattern
+				      * must be appropriate to hold
+				      * the nonzero entries of the
+				      * full matrix.
+				      */
+    template <typename somenumber>
+    void copy_from (const FullMatrix<somenumber> &matrix);
+    
+				     /**
+				      * Add @p{matrix} scaled by
+				      * @p{factor} to this matrix. The
+				      * function throws an error if
+				      * the sparsity patterns of the
+				      * two involved matrices do not
+				      * point to the same object,
+				      * since in this case the
+				      * operation is cheaper.
+				      *
+				      * The source matrix may be a matrix
+				      * of arbitrary type, as long as its
+				      * data type is convertible to the
+				      * data type of this matrix.
+				      */
+    template <typename somenumber>
+    void add_scaled (const number factor,
+		     const SparseMatrix<somenumber> &matrix);
+    
+				     /**
+				      * Return the value of the entry
+				      * (i,j).  This may be an
+				      * expensive operation and you
+				      * should always take care where
+				      * to call this function.  In
+				      * order to avoid abuse, this
+				      * function throws an exception
+				      * if the required element does
+				      * not exist in the matrix.
+				      *
+				      * In case you want a function
+				      * that returns zero instead (for
+				      * entries that are not in the
+				      * sparsity pattern of the
+				      * matrix), use the @p{el}
+				      * function.
+				      */
+    number operator () (const unsigned int i,
+			const unsigned int j) const;
+
+				     /**
+				      * This function is mostly like
+				      * @p{operator()} in that it
+				      * returns the value of the
+				      * matrix entry @p{(i,j)}. The only
+				      * difference is that if this
+				      * entry does not exist in the
+				      * sparsity pattern, then instead
+				      * of raising an exception, zero
+				      * is returned. While this may be
+				      * convenient in some cases, note
+				      * that it is simple to write
+				      * algorithms that are slow
+				      * compared to an optimal
+				      * solution, since the sparsity
+				      * of the matrix is not used.
+				      */
+    number el (const unsigned int i,
+	       const unsigned int j) const;
+
+				     /**
+				      * Return the main diagonal element in
+				      * the @p{i}th row. This function throws an
+				      * error if the matrix is not square.
+				      *
+				      * This function is considerably
+				      * faster than the @p{operator()},
+				      * since for square matrices, the
+				      * diagonal entry is always the
+				      * first to be stored in each row
+				      * and access therefore does not
+				      * involve searching for the
+				      * right column number.
+				      */
+    number diag_element (const unsigned int i) const;
+
+				     /**
+				      * Same as above, but return a
+				      * writeable reference. You're
+				      * sure you know what you do?
+				      */
+    number & diag_element (const unsigned int i);
+
+				     /**
+				      * Access to values in internal
+				      * mode.  Returns the value of
+				      * the @p{index}th entry in
+				      * @p{row}. Here, @p{index} refers to
+				      * the internal representation of
+				      * the matrix, not the column. Be
+				      * sure to understand what you are
+				      * doing here.
+				      */
+    number raw_entry (const unsigned int row,
+		      const unsigned int index) const;
+    
+    				     /**
+				      * This is for hackers. Get
+				      * access to the @p{i}th element of
+				      * this matrix. The elements are
+				      * stored in a consecutive way,
+				      * refer to the @p{SparsityPattern}
+				      * class for more details.
+				      *
+				      * You should use this interface
+				      * very carefully and only if you
+				      * are absolutely sure to know
+				      * what you do. You should also
+				      * note that the structure of
+				      * these arrays may change over
+				      * time.  If you change the
+				      * layout yourself, you should
+				      * also rename this function to
+				      * avoid programs relying on
+				      * outdated information!
+				      */
+    number global_entry (const unsigned int i) const;
+
+				     /**
+				      * Same as above, but with write
+				      * access.  You certainly know
+				      * what you do?
+				      */
+    number & global_entry (const unsigned int i);
+
+				     /**
+				      * Matrix-vector multiplication:
+				      * let $dst = M*src$ with $M$
+				      * being this matrix.
+				      */
+    template <typename somenumber>
+    void vmult (Vector<somenumber>       &dst,
+		const Vector<somenumber> &src) const;
+    
+				     /**
+				      * Matrix-vector multiplication:
+				      * let $dst = M^T*src$ with $M$
+				      * being this matrix. This
+				      * function does the same as
+				      * @p{vmult} but takes the
+				      * transposed matrix.
+				      */
+    template <typename somenumber>
+    void Tvmult (Vector<somenumber>       &dst,
+		 const Vector<somenumber> &src) const;
+  
+				     /**
+				      * Adding Matrix-vector
+				      * multiplication. Add $M*src$ on
+				      * $dst$ with $M$ being this
+				      * matrix.
+				      */
+    template <typename somenumber>
+    void vmult_add (Vector<somenumber>       &dst,
+		    const Vector<somenumber> &src) const;
+    
+				     /**
+				      * Adding Matrix-vector
+				      * multiplication. Add $M^T*src$
+				      * to $dst$ with $M$ being this
+				      * matrix. This function does the
+				      * same as @p{vmult_add} but takes
+				      * the transposed matrix.
+				      */
+    template <typename somenumber>
+    void Tvmult_add (Vector<somenumber>       &dst,
+		     const Vector<somenumber> &src) const;
+  
+				     /**
+				      * Return the square of the norm
+				      * of the vector $v$ with respect
+				      * to the norm induced by this
+				      * matrix,
+				      * i.e. $\left(v,Mv\right)$. This
+				      * is useful, e.g. in the finite
+				      * element context, where the
+				      * $L_2$ norm of a function
+				      * equals the matrix norm with
+				      * respect to the mass matrix of
+				      * the vector representing the
+				      * nodal values of the finite
+				      * element function.
+				      *
+				      * Obviously, the matrix needs to
+				      * be square for this operation.
+				      */
+    template <typename somenumber>
+    somenumber matrix_norm_square (const Vector<somenumber> &v) const;
+
+				     /**
+				      * Compute the matrix scalar
+				      * product $\left(u,Mv\right)$.
+				      */
+    template <typename somenumber>
+    somenumber matrix_scalar_product (const Vector<somenumber> &u,
+				      const Vector<somenumber> &v) const;
+    
+    				     /**
+				      * Return the l1-norm of the matrix, that is
+				      * $|M|_1=max_{all columns j}\sum_{all 
+				      * rows i} |M_ij|$,
+				      * (max. sum of columns).
+				      * This is the
+				      * natural matrix norm that is compatible
+				      * to the l1-norm for vectors, i.e.
+				      * $|Mv|_1\leq |M|_1 |v|_1$.
+				      * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
+				      */
+    number l1_norm () const;
+
+    				     /**
+				      * Return the linfty-norm of the
+				      * matrix, that is
+				      * $|M|_infty=max_{all rows i}\sum_{all 
+				      * columns j} |M_ij|$,
+				      * (max. sum of rows).
+				      * This is the
+				      * natural matrix norm that is compatible
+				      * to the linfty-norm of vectors, i.e.
+				      * $|Mv|_infty \leq |M|_infty |v|_infty$.
+				      * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
+				      */
+    number linfty_norm () const;
+
+				     /**
+				      * Compute the residual of an
+				      * equation @p{Mx=b}, where the
+				      * residual is defined to be
+				      * @p{r=b-Mx} with @p{x} typically
+				      * being an approximate of the
+				      * true solution of the
+				      * equation. Write the residual
+				      * into @p{dst}. The l2 norm of
+				      * the residual vector is
+				      * returned.
+				      */
+    template <typename somenumber>
+    somenumber residual (Vector<somenumber>       &dst,
+			 const Vector<somenumber> &x,
+			 const Vector<somenumber> &b) const;
+    
+				     /**
+				      * Apply the Jacobi
+				      * preconditioner, which
+				      * multiplies every element of
+				      * the @p{src} vector by the
+				      * inverse of the respective
+				      * diagonal element and
+				      * multiplies the result with the
+				      * damping factor @p{omega}.
+				      */
+    template <typename somenumber>
+    void precondition_Jacobi (Vector<somenumber>       &dst,
+			      const Vector<somenumber> &src,
+			      const number              omega = 1.) const;
+
+				     /**
+				      * Apply SSOR preconditioning to
+				      * @p{src}.
+				      */
+    template <typename somenumber>
+    void precondition_SSOR (Vector<somenumber>       &dst,
+			    const Vector<somenumber> &src,
+			    const number              om = 1.) const;
+
+				     /**
+				      * Apply SOR preconditioning matrix to @p{src}.
+				      * The result of this method is
+				      * $dst = (om D - L)^{-1} src$.
+				      */
+    template <typename somenumber>
+    void precondition_SOR (Vector<somenumber>       &dst,
+			   const Vector<somenumber> &src,
+ 			   const number              om = 1.) const;
+    
+				     /**
+				      * Apply transpose SOR preconditioning matrix to @p{src}.
+				      * The result of this method is
+				      * $dst = (om D - U)^{-1} src$.
+				      */
+    template <typename somenumber>
+    void precondition_TSOR (Vector<somenumber>       &dst,
+			    const Vector<somenumber> &src,
+			    const number              om = 1.) const;
+    
+				     /**
+				      * Perform SSOR preconditioning
+				      * in-place.  Apply the
+				      * preconditioner matrix without
+				      * copying to a second vector.
+				      * @p{omega} is the relaxation
+				      * parameter.
+				      */
+    template <typename somenumber>
+    void SSOR (Vector<somenumber> &v,
+	       const number        omega = 1.) const;
+
+				     /**
+				      * Perform an SOR preconditioning in-place.
+				      * The result is $v = (\omega D - L)^{-1} v$.
+				      * @p{omega} is the damping parameter.
+				      */
+    template <typename somenumber>
+    void SOR (Vector<somenumber> &v,
+	      const number        om = 1.) const;
+
+				     /**
+				      * Perform a transpose SOR preconditioning in-place.
+				      * The result is $v = (\omega D - L)^{-1} v$.
+				      * @p{omega} is the damping parameter.
+				      */
+    template <typename somenumber>
+    void TSOR (Vector<somenumber> &v,
+	      const number        om = 1.) const;
+
+				     /**
+				      * Do one SOR step on @p{v}.
+				      * Performs a direct SOR step
+				      * with right hand side @p{b}.
+				      */
+    template <typename somenumber>
+    void SOR_step (Vector<somenumber> &v,
+		   const Vector<somenumber> &b,
+		   const number        om = 1.) const;
+
+				     /**
+				      * Do one adjoint SOR step on
+				      * @p{v}.  Performs a direct TSOR
+				      * step with right hand side @p{b}.
+				      */
+    template <typename somenumber>
+    void TSOR_step (Vector<somenumber> &v,
+		    const Vector<somenumber> &b,
+		    const number        om = 1.) const;
+
+				     /**
+				      * Do one adjoint SSOR step on
+				      * @p{v}.  Performs a direct SSOR
+				      * step with right hand side @p{b}
+				      * by performing TSOR after SOR.
+				      */
+    template <typename somenumber>
+    void SSOR_step (Vector<somenumber> &v,
+		    const Vector<somenumber> &b,
+		    const number        om = 1.) const;
+
+				     /**
+				      * Return a (constant) reference
+				      * to the underlying sparsity
+				      * pattern of this matrix.
+				      *
+				      * Though the return value is
+				      * declared @p{const}, you should
+				      * be aware that it may change if
+				      * you call any nonconstant
+				      * function of objects which
+				      * operate on it.
+				      */
+    const SparsityPattern & get_sparsity_pattern () const;
+
+				     /**
+				      * STL-like iterator with the
+				      * first entry.
+				      */
+    const_iterator begin () const;
+
+				     /**
+				      * Final iterator.
+				      */
+    const_iterator end () const;
+    
+				     /**
+				      * STL-like iterator with the
+				      * first entry of row @p{r}.
+				      */
+    const_iterator begin (unsigned int r) const;
+
+				     /**
+				      * Final iterator of row @p{r}.
+				      */
+    const_iterator end (unsigned int r) const;
+    
+				     /**
+				      * Print the matrix to the given
+				      * stream, using the format
+				      * @p{(line,col) value}, i.e. one
+				      * nonzero entry of the matrix
+				      * per line.
+				      */
+    void print (std::ostream &out) const;
+
+				     /**
+				      * Print the matrix in the usual
+				      * format, i.e. as a matrix and
+				      * not as a list of nonzero
+				      * elements. For better
+				      * readability, elements not in
+				      * the matrix are displayed as
+				      * empty space, while matrix
+				      * elements which are explicitly
+				      * set to zero are displayed as
+				      * such.
+				      *
+				      * The parameters allow for a
+				      * flexible setting of the output
+				      * format: @p{precision} and
+				      * @p{scientific} are used to
+				      * determine the number format,
+				      * where @p{scientific} = @p{false}
+				      * means fixed point notation.  A
+				      * zero entry for @p{width} makes
+				      * the function compute a width,
+				      * but it may be changed to a
+				      * positive value, if output is
+				      * crude.
+				      *
+				      * Additionally, a character for
+				      * an empty value may be
+				      * specified.
+				      *
+				      * Finally, the whole matrix can
+				      * be multiplied with a common
+				      * denominator to produce more
+				      * readable output, even
+				      * integers.
+				      *
+				      * This function
+				      * may produce @em{large} amounts of
+				      * output if applied to a large matrix!
+				      */
+    void print_formatted (std::ostream       &out,
+			  const unsigned int  precision   = 3,
+			  const bool          scientific  = true,
+			  const unsigned int  width       = 0,
+			  const char         *zero_string = " ",
+			  const double        denominator = 1.) const;
+
+				     /**
+				      * Write the data of this object
+				      * en bloc to a file. This is
+				      * done in a binary mode, so the
+				      * output is neither readable by
+				      * humans nor (probably) by other
+				      * computers using a different
+				      * operating system of number
+				      * format.
+				      *
+				      * The purpose of this function
+				      * is that you can swap out
+				      * matrices and sparsity pattern
+				      * if you are short of memory,
+				      * want to communicate between
+				      * different programs, or allow
+				      * objects to be persistent
+				      * across different runs of the
+				      * program.
+				      */
+    void block_write (std::ostream &out) const;
+
+				     /**
+				      * Read data that has previously
+				      * been written by
+				      * @p{block_write} en block from
+				      * a file. This is done using the
+				      * inverse operations to the
+				      * above function, so it is
+				      * reasonably fast because the
+				      * bitstream is not interpreted
+				      * except for a few numbers up
+				      * front.
+				      *
+				      * The object is resized on this
+				      * operation, and all previous
+				      * contents are lost. Note,
+				      * however, that no checks are
+				      * performed whether new data and
+				      * the underlying
+				      * @ref{SparsityPattern} object
+				      * fit together. It is your
+				      * responsibility to make sure
+				      * that the sparsity pattern and
+				      * the data to be read match.
+				      *
+				      * A primitive form of error
+				      * checking is performed which
+				      * will recognize the bluntest
+				      * attempts to interpret some
+				      * data as a vector stored
+				      * bitwise to a file, but not
+				      * more.
+				      */
+    void block_read (std::istream &in);
+
+				     /**
+				      * Determine an estimate for the
+				      * memory consumption (in bytes)
+				      * of this object.
+				      */
+    unsigned int memory_consumption () const;
+    
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcNotCompressed);
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcMatrixNotInitialized);
+				     /**
+				      * Exception
+				      */
+    DeclException2 (ExcInvalidIndex,
+		    int, int,
+		    << "The entry with index <" << arg1 << ',' << arg2
+		    << "> does not exist.");
+				     /**
+				      * Exception
+				      */
+    DeclException1 (ExcInvalidIndex1,
+		    int,
+		    << "The index " << arg1 << " is not in the allowed range.");
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcMatrixNotSquare);
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcDifferentSparsityPatterns);
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcInvalidConstructorCall);
+				     /**
+				      * Exception
+				      */
+    DeclException2 (ExcIteratorRange,
+		    int, int,
+		    << "The iterators denote a range of " << arg1
+		    << " elements, but the given number of rows was " << arg2);
+    
+  private:
+				     /**
+				      * Pointer to the sparsity
+				      * pattern used for this
+				      * matrix. In order to guarantee
+				      * that it is not deleted while
+				      * still in use, we subscribe to
+				      * it using the @p{SmartPointer}
+				      * class.
+				      */
+    SmartPointer<const SparsityPattern> cols;
+
+				     /**
+				      * Array of values for all the
+				      * nonzero entries. The position
+				      * within the matrix, i.e.  the
+				      * row and column number for a
+				      * given entry can only be
+				      * deduced using the sparsity
+				      * pattern. The same holds for
+				      * the more common operation of
+				      * finding an entry by its
+				      * coordinates.
+				      */
+    number *val;
+
+				     /**
+				      * Allocated size of
+				      * @p{val}. This can be larger
+				      * than the actually used part if
+				      * the size of the matrix was
+				      * reduced somewhen in the past
+				      * by associating a sparsity
+				      * pattern with a smaller size to
+				      * this object, using the
+				      * @p{reinit} function.
+				      */
+    unsigned int max_len;
+
+				     /**
+				      * Version of @p{vmult} which only
+				      * performs its actions on the
+				      * region defined by
+				      * @p{[begin_row,end_row)}. This
+				      * function is called by @p{vmult}
+				      * in the case of enabled
+				      * multithreading.
+				      */
+    template <typename somenumber>
+    void threaded_vmult (Vector<somenumber>       &dst,
+			 const Vector<somenumber> &src,
+			 const unsigned int        begin_row,
+			 const unsigned int        end_row) const;
+
+				     /**
+				      * Version of
+				      * @p{matrix_norm_square} which
+				      * only performs its actions on
+				      * the region defined by
+				      * @p{[begin_row,end_row)}. This
+				      * function is called by
+				      * @p{matrix_norm_square} in the
+				      * case of enabled
+				      * multithreading.
+				      */
+    template <typename somenumber>
+    void threaded_matrix_norm_square (const Vector<somenumber> &v,
+				      const unsigned int        begin_row,
+				      const unsigned int        end_row,
+				      somenumber               *partial_sum) const;
+
+    				     /**
+				      * Version of
+				      * @p{matrix_scalar_product} which
+				      * only performs its actions on
+				      * the region defined by
+				      * @p{[begin_row,end_row)}. This
+				      * function is called by
+				      * @p{matrix_scalar_product} in the
+				      * case of enabled
+				      * multithreading.
+				      */
+    template <typename somenumber>
+    void threaded_matrix_scalar_product (const Vector<somenumber> &u,
+					 const Vector<somenumber> &v,
+					 const unsigned int        begin_row,
+					 const unsigned int        end_row,
+					 somenumber               *partial_sum) const;
+
+				     /**
+				      * Version of @p{residual} which
+				      * only performs its actions on
+				      * the region defined by
+				      * @p{[begin_row,end_row)} (these
+				      * numbers are the components of
+				      * @p{interval}). This function is
+				      * called by @p{residual} in the
+				      * case of enabled
+				      * multithreading.
+				      */
+    template <typename somenumber>
+    void threaded_residual (Vector<somenumber>       &dst,
+			    const Vector<somenumber> &u,
+			    const Vector<somenumber> &b,
+			    const std::pair<unsigned int,unsigned int> interval,
+			    somenumber               *partial_norm) const;
+
+				     // make all other sparse matrices
+				     // friends
+    template <typename somenumber> friend class SparseMatrix;
+};
+
+
+/*---------------------- Inline functions -----------------------------------*/
+
+
+
+template <typename number>
+inline
+unsigned int SparseMatrix<number>::m () const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  return cols->rows;
+};
+
+
+template <typename number>
+inline
+unsigned int SparseMatrix<number>::n () const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  return cols->cols;
+};
+
+
+template <typename number>
+inline
+void SparseMatrix<number>::set (const unsigned int i,
+				const unsigned int j,
+				const number value)
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+				   // it is allowed to set elements of
+				   // the matrix that are not part of
+				   // the sparsity pattern, if the
+				   // value to which we set it is zero
+  const unsigned int index = cols->operator()(i,j);
+  Assert ((index != SparsityPattern::invalid_entry) ||
+	  (value == 0.),
+	  ExcInvalidIndex(i,j));
+
+  if (index != SparsityPattern::invalid_entry)
+    val[index] = value;
+};
+
+
+
+template <typename number>
+inline
+void SparseMatrix<number>::add (const unsigned int i,
+				const unsigned int j,
+				const number value)
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+
+  const unsigned int index = cols->operator()(i,j);
+  Assert ((index != SparsityPattern::invalid_entry) ||
+	  (value == 0.),
+	  ExcInvalidIndex(i,j));
+
+  if (value != 0.)
+    val[index] += value;
+};
+
+
+
+template <typename number>
+inline
+number SparseMatrix<number>::operator () (const unsigned int i,
+					  const unsigned int j) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry,
+	  ExcInvalidIndex(i,j));
+  return val[cols->operator()(i,j)];
+};
+
+
+
+template <typename number>
+inline
+number SparseMatrix<number>::el (const unsigned int i,
+				 const unsigned int j) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  const unsigned int index = cols->operator()(i,j);
+
+  if (index != SparsityPattern::invalid_entry)
+    return val[index];
+  else
+    return 0;
+};
+
+
+
+template <typename number>
+inline
+number SparseMatrix<number>::diag_element (const unsigned int i) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (i<m(), ExcInvalidIndex1(i));
+  
+				   // Use that the first element in each
+				   // row of a square matrix is the main
+				   // diagonal
+  return val[cols->rowstart[i]];
+};
+
+
+
+template <typename number>
+inline
+number & SparseMatrix<number>::diag_element (const unsigned int i)
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (i<m(), ExcInvalidIndex1(i));
+  
+				   // Use that the first element in each
+				   // row of a square matrix is the main
+				   // diagonal
+  return val[cols->rowstart[i]];
+};
+
+
+
+template <typename number>
+inline
+number
+SparseMatrix<number>::raw_entry (const unsigned int row,
+				 const unsigned int index) const
+{
+  Assert(row<cols->rows, ExcIndexRange(row,0,cols->rows));
+  Assert(index<cols->row_length(row),
+	 ExcIndexRange(index,0,cols->row_length(row)));
+
+  return val[cols->rowstart[row]+index];
+};
+
+
+
+template <typename number>
+inline
+number SparseMatrix<number>::global_entry (const unsigned int j) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (j < cols->n_nonzero_elements(),
+	  ExcIndexRange (j, 0, cols->n_nonzero_elements()));
+  
+  return val[j];
+};
+
+
+
+template <typename number>
+inline
+number & SparseMatrix<number>::global_entry (const unsigned int j)
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (j < cols->n_nonzero_elements(),
+	  ExcIndexRange (j, 0, cols->n_nonzero_elements()));
+
+  return val[j];
+};
+
+
+
+template <typename number>
+template <typename ForwardIterator>
+void
+SparseMatrix<number>::copy_from (const ForwardIterator begin,
+				 const ForwardIterator end)
+{
+  Assert (static_cast<unsigned int>(std::distance (begin, end)) == m(),
+	  ExcIteratorRange (std::distance (begin, end), m()));
+
+				   // for use in the inner loop, we
+				   // define a typedef to the type of
+				   // the inner iterators
+  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
+  unsigned int row=0;
+  for (ForwardIterator i=begin; i!=end; ++i, ++row)
+    {
+      const inner_iterator end_of_row = i->end();
+      for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
+					 // write entries
+	set (row, j->first, j->second);
+    };
+};
+
+
+//----------------------------------------------------------------------//
+
+template <typename number>
+inline
+SparseMatrix<number>::Accessor::Accessor (
+  const SparseMatrix<number>* matrix,
+  unsigned int r,
+  unsigned short i)
+		:
+		matrix(matrix),
+		a_row(r),
+		a_index(i)
+{}
+
+
+template <typename number>
+inline
+unsigned int
+SparseMatrix<number>::Accessor::row() const
+{
+  return a_row;
+}
+
+
+template <typename number>
+inline
+unsigned int
+SparseMatrix<number>::Accessor::column() const
+{
+  const SparsityPattern& pat = matrix->get_sparsity_pattern();
+  return pat.get_column_numbers()[pat.get_rowstart_indices()[a_row]+a_index];
+}
+
+
+template <typename number>
+inline
+unsigned short
+SparseMatrix<number>::Accessor::index() const
+{
+  return a_index;
+}
+
+
+
+template <typename number>
+inline
+number
+SparseMatrix<number>::Accessor::value() const
+{
+  return matrix->raw_entry(a_row, a_index);
+}
+
+
+template <typename number>
+inline
+SparseMatrix<number>::const_iterator::const_iterator(
+  const SparseMatrix<number>* matrix,
+  unsigned int r,
+  unsigned short i)
+		:
+		Accessor(matrix, r, i)
+{}
+
+
+template <typename number>
+inline
+typename SparseMatrix<number>::const_iterator&
+SparseMatrix<number>::const_iterator::operator++ ()
+{
+  Assert (a_row < matrix->m(), ExcIteratorPastEnd());
+  
+  ++a_index;
+  if (a_index >= matrix->get_sparsity_pattern().row_length(a_row))
+    {
+      a_index = 0;
+      a_row++;
+    }
+  return *this;
+}
+
+
+template <typename number>
+inline
+const typename SparseMatrix<number>::Accessor&
+SparseMatrix<number>::const_iterator::operator* () const
+{
+  return *this;
+}
+
+
+template <typename number>
+inline
+const typename SparseMatrix<number>::Accessor*
+SparseMatrix<number>::const_iterator::operator-> () const
+{
+  return this;
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrix<number>::const_iterator::operator == (
+  const const_iterator& other) const
+{
+  return (row() == other->row() && index() == other->index());
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrix<number>::const_iterator::operator != (
+  const const_iterator& other) const
+{
+  return ! (*this == other);
+}
+
+
+template <typename number>
+inline
+bool
+SparseMatrix<number>::const_iterator::operator < (
+  const const_iterator& other) const
+{
+  return (row() < other->row() ||
+	  (row() == other->row() && index() < other->index()));
+}
+
+
+template <typename number>
+inline
+typename SparseMatrix<number>::const_iterator
+SparseMatrix<number>::begin () const
+{
+  return const_iterator(this, 0, 0);
+}
+
+template <typename number>
+inline
+typename SparseMatrix<number>::const_iterator
+SparseMatrix<number>::end () const
+{
+  return const_iterator(this, m(), 0);
+}
+
+template <typename number>
+inline
+typename SparseMatrix<number>::const_iterator
+SparseMatrix<number>::begin (unsigned int r) const
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return const_iterator(this, r, 0);
+}
+
+template <typename number>
+inline
+typename SparseMatrix<number>::const_iterator
+SparseMatrix<number>::end (unsigned int r) const
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return const_iterator(this, r+1, 0);
+}
+
+
+
+
+/*----------------------------   sparse_matrix.h     ---------------------------*/
+
+#endif
+/*----------------------------   sparse_matrix.h     ---------------------------*/