]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
These files should not have been picked up.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 Dec 2002 17:29:16 +0000 (17:29 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 Dec 2002 17:29:16 +0000 (17:29 +0000)
git-svn-id: https://svn.dealii.org/trunk@6793 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_ilu.templates.h.x [deleted file]
deal.II/lac/include/lac/sparse_matrix.h [deleted file]
deal.II/lac/source/sparse_ilu.cc.x [deleted file]

diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h.x b/deal.II/lac/include/lac/sparse_ilu.templates.h.x
deleted file mode 100644 (file)
index 608d889..0000000
+++ /dev/null
@@ -1,353 +0,0 @@
-//----------------------------  sparse_ilu.templates.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_ilu.templates.h  ---------------------------
-#ifndef __deal2__sparse_ilu_templates_h
-#define __deal2__sparse_ilu_templates_h
-
-
-
-#include <base/config.h>
-#include <lac/vector.h>
-#include <lac/sparse_ilu.h>
-
-#include <algorithm>
-#include <cmath>
-
-
-template <typename number>
-SparseILU<number>::SparseILU () 
-{};
-
-
-
-template <typename number>
-SparseILU<number>::SparseILU (const SparsityPattern &sparsity) :
-               SparseMatrix<number> (sparsity)
-{};
-
-
-
-template <typename number>
-void SparseILU<number>::reinit ()
-{
-  SparseMatrix<number>::reinit ();
-};
-
-
-
-template <typename number>
-void SparseILU<number>::reinit (const SparsityPattern &sparsity)
-{
-  SparseMatrix<number>::reinit (sparsity);
-};
-
-
-
-template <typename number>
-template <typename somenumber>
-void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
-                                  const double                    strengthen_diagonal)
-{
-  Assert (matrix.m()==matrix.n(), ExcMatrixNotSquare ());
-  Assert (this->m()==this->n(),   ExcMatrixNotSquare ());
-  Assert (matrix.m()==this->m(),  ExcSizeMismatch(matrix.m(), this->m()));
-  
-  Assert (strengthen_diagonal>=0, ExcInvalidStrengthening (strengthen_diagonal));
-
-
-                                  // first thing: copy over all elements
-                                  // of @p{matrix} to the present object
-                                  //
-                                  // note that some elements in this
-                                  // matrix may not be in @p{matrix},
-                                  // so we need to preset our matrix
-                                  // by zeroes.
-  if (true)
-    {
-                                      // preset the elements
-      std::fill_n (&this->global_entry(0),
-                  this->n_nonzero_elements(),
-                  0);
-
-                                      // note: pointers to the sparsity
-                                      // pattern of the old matrix!
-      const unsigned int * const rowstart_indices
-       = matrix.get_sparsity_pattern().get_rowstart_indices();
-      const unsigned int * const column_numbers
-       = matrix.get_sparsity_pattern().get_column_numbers();
-      
-      for (unsigned int row=0; row<this->m(); ++row)
-       for (const unsigned int * col = &column_numbers[rowstart_indices[row]];
-            col != &column_numbers[rowstart_indices[row+1]]; ++col)
-         set (row, *col, matrix.global_entry(col-column_numbers));
-    };
-
-  if (strengthen_diagonal > 0)
-    for (unsigned int row=0; row<this->m(); ++row)
-      {
-                                        // get the length of the row
-                                        // (without the diagonal element)
-       const unsigned int
-         rowlength = (this->get_sparsity_pattern().get_rowstart_indices()[row+1]
-                      -
-                      this->get_sparsity_pattern().get_rowstart_indices()[row]
-                      -
-                      1);
-       
-                                        // get the global index of the first
-                                        // non-diagonal element in this row
-       const unsigned int rowstart
-         = this->get_sparsity_pattern().get_rowstart_indices()[row] + 1;
-       number * const diagonal_element = &this->global_entry(rowstart-1);
-
-       number rowsum = 0;
-       for (unsigned int global_index=rowstart;
-            global_index<rowstart+rowlength; ++global_index)
-         rowsum += std::fabs(this->global_entry(global_index));
-
-       *diagonal_element += strengthen_diagonal * rowsum;
-      };
-
-
-                                  // now work only on this
-                                  // matrix
-  const SparsityPattern             &sparsity = this->get_sparsity_pattern();
-  const unsigned int * const rowstart_indices = sparsity.get_rowstart_indices();
-  const unsigned int * const column_numbers   = sparsity.get_column_numbers();
-  
-/*
-  PSEUDO-ALGORITHM
-  (indices=0..N-1)
-  
-  for i=1..N-1
-    a[i-1,i-1] = a[i-1,i-1]^{-1}
-
-    for k=0..i-1
-      a[i,k] = a[i,k] * a[k,k]
-
-      for j=k+1..N-1
-        if (a[i,j] exists & a[k,j] exists)
-          a[i,j] -= a[i,k] * a[k,j]
-*/
-
-
-                                  // i := row
-  for (unsigned int row=1; row<this->m(); ++row)
-    {
-                                      // invert diagonal element of the
-                                      // previous row. this is a hack,
-                                      // which is possible since this
-                                      // element is not needed any more
-                                      // in the process of decomposition
-                                      // and since it makes the backward
-                                      // step when applying the decomposition
-                                      // significantly faster
-      AssertThrow((this->global_entry(rowstart_indices[row-1]) !=0),
-                 ExcDivideByZero());
-      
-      this->global_entry (rowstart_indices[row-1])
-       = 1./this->global_entry (rowstart_indices[row-1]);
-
-                                      // let k run over all lower-left
-                                      // elements of row i; skip
-                                      // diagonal element at start
-      const unsigned int * first_of_row
-       = &column_numbers[rowstart_indices[row]+1];
-      const unsigned int * first_after_diagonal
-       = std::lower_bound (&column_numbers[rowstart_indices[row]+1],
-                           &column_numbers[rowstart_indices[row+1]],
-                           row);
-
-                                      // k := *col_ptr
-      for (const unsigned int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++col_ptr)
-       {
-         const unsigned int global_index_ik = col_ptr-column_numbers;
-         this->global_entry(global_index_ik) *= this->diag_element(*col_ptr);
-
-                                          // now do the inner loop over
-                                          // j. note that we need to do
-                                          // it in the right order, i.e.
-                                          // taking into account that the
-                                          // columns are sorted within each
-                                          // row correctly, but excluding
-                                          // the main diagonal entry
-         const int global_index_ki = sparsity(*col_ptr,row);
-
-         if (global_index_ki != -1)
-           this->diag_element(row) -= this->global_entry(global_index_ik) *
-                                      this->global_entry(global_index_ki);
-
-         for (const unsigned int * j = col_ptr+1;
-              j<&column_numbers[rowstart_indices[row+1]];
-              ++j)
-           {
-//TODO:[WB] make code faster by using the following comment          
-                                              // note: this inner loop could
-                                              // be made considerable faster
-                                              // if we consulted the row
-                                              // with number *col_ptr,
-                                              // instead of always asking
-                                              // sparsity(*col_ptr,*j),
-                                              // since we traverse this
-                                              // row linearly. I just didn't
-                                              // have the time to figure out
-                                              // the details.
-                     const int global_index_ij = j - &column_numbers[0],
-                       global_index_kj = sparsity(*col_ptr,*j);
-             if ((global_index_ij != -1) &&
-                 (global_index_kj != -1))
-               this->global_entry(global_index_ij) -= this->global_entry(global_index_ik) *
-                                                      this->global_entry(global_index_kj);
-           };
-       };
-    };
-
-                                  // Here the very last diagonal
-                                  // element still has to be inverted
-                                  // because the for-loop doesn't do
-                                  // it...
- this->diag_element(this->m()-1) = 1./this->diag_element(this->m()-1);
-
-/*
-  OLD CODE, rather crude first implementation with an algorithm taken
-  from 'W. Hackbusch, G. Wittum: Incomplete Decompositions (ILU)-
-  Algorithms, Theory, and Applications', page 6.
-  
-  for (unsigned int k=0; k<m()-1; ++k)
-    for (unsigned int i=k+1; i<m(); ++i)
-      {
-                                        // get the global index
-                                        // of the element (i,k)
-       const int global_index_ik = get_sparsity_pattern()(i,k);
-
-                                        // if this element is zero,
-                                        // then we continue with the
-                                        // next i, since e would be
-                                        // zero and nothing would happen
-                                        // in this loop
-       if (global_index_ik == -1)
-         continue;
-       
-       const number e = global_entry(global_index_ik) / diag_element(k);
-       global_entry(global_index_ik) = e;
-
-       for (unsigned int j=k+1; j<m(); ++j)
-         {
-                                            // find out about a_kj
-                                            // if this does not exist,
-                                            // then the updates within
-                                            // this innermost loop would
-                                            // be zero, invariable of the
-                                            // fact of whether a_ij is a
-                                            // nonzero or a zero element
-           const int global_index_kj = get_sparsity_pattern()(k,j);
-           if (global_index_kj == -1)
-             continue;
-
-           const int global_index_ij = get_sparsity_pattern()(i,j);
-           if (global_index_ij != -1)
-             global_entry(global_index_ij) -= e*global_entry(global_index_kj);
-           else
-             diag_element(i) -= e*global_entry(global_index_kj);
-         };
-      };
-*/      
-};
-
-
-
-template <typename number>
-template <typename somenumber>
-void SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
-                                            const Vector<somenumber> &src) const 
-{
-  Assert (dst.size() == src.size(), ExcSizeMismatch(dst.size(), src.size()));
-  Assert (dst.size() == this->m(), ExcSizeMismatch(dst.size(), this->m()));
-  
-  const unsigned int N=dst.size();
-  const unsigned int * const rowstart_indices
-    = this->get_sparsity_pattern().get_rowstart_indices();
-  const unsigned int * const column_numbers
-    = this->get_sparsity_pattern().get_column_numbers();
-                                  // solve LUx=b in two steps:
-                                  // first Ly = b, then
-                                  //       Ux = y
-                                  //
-                                  // first a forward solve. since
-                                  // the diagonal values of L are
-                                  // one, there holds
-                                  // y_i = b_i
-                                  //       - sum_{j=0}^{i-1} L_{ij}y_j
-                                  // we split the y_i = b_i off and
-                                  // perform it at the outset of the
-                                  // loop
-  dst = src;
-  for (unsigned int row=0; row<N; ++row)
-    {
-                                      // get start of this row. skip the
-                                      // diagonal element
-      const unsigned int * const rowstart = &column_numbers[rowstart_indices[row]+1];
-                                      // find the position where the part
-                                      // right of the diagonal starts
-      const unsigned int * const first_after_diagonal
-       = std::lower_bound (rowstart,
-                           &column_numbers[rowstart_indices[row+1]],
-                           row);
-      
-      for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col)
-       dst(row) -= this->global_entry (col-column_numbers) * dst(*col);
-    };
-
-                                  // now the backward solve. same
-                                  // procedure, but we need not set
-                                  // dst before, since this is already
-                                  // done.
-                                  //
-                                  // note that we need to scale now,
-                                  // since the diagonal is not zero
-                                  // now
-  for (int row=N-1; row>=0; --row)
-    {
-                                      // get end of this row
-      const unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]];
-                                      // find the position where the part
-                                      // right of the diagonal starts
-      const unsigned int * const first_after_diagonal
-       = std::lower_bound (&column_numbers[rowstart_indices[row]+1],
-                           &column_numbers[rowstart_indices[row+1]],
-                           static_cast<unsigned int>(row));
-      
-      for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col)
-       dst(row) -= this->global_entry (col-column_numbers) * dst(*col);
-
-                                      // scale by the diagonal element.
-                                      // note that the diagonal element
-                                      // was stored inverted
-      dst(row) *= this->diag_element(row);
-    };
-};
-
-
-
-template <typename number>
-unsigned int
-SparseILU<number>::memory_consumption () const
-{
-  return SparseMatrix<number>::memory_consumption ();
-};
-
-
-
-/*----------------------------   sparse_ilu.templates.h     ---------------------------*/
-
-#endif
-/*----------------------------   sparse_ilu.templates.h     ---------------------------*/
diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h
deleted file mode 100644 (file)
index 29c7192..0000000
+++ /dev/null
@@ -1,1491 +0,0 @@
-//----------------------------  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     ---------------------------*/
diff --git a/deal.II/lac/source/sparse_ilu.cc.x b/deal.II/lac/source/sparse_ilu.cc.x
deleted file mode 100644 (file)
index 771c729..0000000
+++ /dev/null
@@ -1,38 +0,0 @@
-//----------------------------  sparse_ilu.cc  ---------------------------
-//    $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_ilu.cc  ---------------------------
-
-
-#include <lac/sparse_ilu.templates.h>
-
-
-// explicit instantiations
-template class SparseILU<double>;
-template void SparseILU<double>::decompose<double> (const SparseMatrix<double> &,
-                                                   const double);
-template void SparseILU<double>::apply_decomposition<double> (Vector<double> &,
-                                                             const Vector<double> &) const;
-template void SparseILU<double>::decompose<float> (const SparseMatrix<float> &,
-                                                  const double);
-template void SparseILU<double>::apply_decomposition<float> (Vector<float> &,
-                                                            const Vector<float> &) const;
-
-
-template class SparseILU<float>;
-template void SparseILU<float>::decompose<double> (const SparseMatrix<double> &,
-                                                  const double);
-template void SparseILU<float>::apply_decomposition<double> (Vector<double> &,
-                                                            const Vector<double> &) const;
-template void SparseILU<float>::decompose<float> (const SparseMatrix<float> &,
-                                                 const double);
-template void SparseILU<float>::apply_decomposition<float> (Vector<float> &,
-                                                           const Vector<float> &) const;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.