]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Remove the old-style matrices and vectors which were not templated.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Feb 1999 16:03:55 +0000 (16:03 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 26 Feb 1999 16:03:55 +0000 (16:03 +0000)
git-svn-id: https://svn.dealii.org/trunk@919 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/dblocksmatrix.h [deleted file]
deal.II/lac/include/lac/dfmatrix.h [deleted file]
deal.II/lac/include/lac/dsmatrix.h [deleted file]
deal.II/lac/include/lac/dvector.h [deleted file]
deal.II/lac/source/dblocksmatrix.cc [deleted file]
deal.II/lac/source/dfmatrix.cc [deleted file]
deal.II/lac/source/dsmatrix.cc [deleted file]
deal.II/lac/source/dvector.cc [deleted file]

diff --git a/deal.II/lac/include/lac/dblocksmatrix.h b/deal.II/lac/include/lac/dblocksmatrix.h
deleted file mode 100644 (file)
index 5501acc..0000000
+++ /dev/null
@@ -1,154 +0,0 @@
-/*----------------------------   dblocksmatrix.h     ---------------------------*/
-/*      $Id$                 */
-#ifndef __dblocksmatrix_H
-#define __dblocksmatrix_H
-/*----------------------------   dblocksmatrix.h     ---------------------------*/
-
-
-#include <lac/dsmatrix.h>
-#include <lac/dfmatrix.h>
-#include <vector.h>
-
-/**
- * Double precision block sparse matrix.
- * The block matrix assumes the matrix consisting of blocks on
- * the diagonal. These diagonal blocks and the elements below the
- * diagonal blocks are used in the #precondition_BlockSOR#.
- *
- * This block matrix structure is given e.g. for the DG method
- * for the transport equation and a downstream numbering.
- * If (as for this DG method) the matrix is empty above the
- * diagonal blocks BlockSOR is a direct solver.
- *
- * This first implementation of the BlockMatrix assumes the
- * matrix having blocks each of the same block size. Varying
- * block sizes within the matrix must still be implemented if needed.
- * @author Ralf Hartmann, 1999
- */
-class dBlockSMatrix: public dSMatrix 
-{
-  public:
-                                    /**
-                                     * Constructor
-                                     */
-    dBlockSMatrix();
-
-                                    /**
-                                     * Destructor
-                                     */
-    virtual ~dBlockSMatrix();
-
-                                    /**
-                                     * Call #dSMatrix::reinit()# and
-                                     * delete the inverse matrices if existent.
-                                     */
-
-    virtual void reinit();
-
-                                    /**
-                                     * Call #dSMatrix::reinit
-                                     * (const dSMatrixStruct &sparsity)# and
-                                     * delete the inverse matrices if existent.
-                                     */
-    virtual void reinit (const dSMatrixStruct &sparsity);
-
-                                    /**
-                                     * Call #dSMatrix::clear# and 
-                                     * delete the inverse matrices if existent.
-                                     */
-    virtual void clear ();
-
-                                    /**
-                                     * Stores the inverse matrices of
-                                     * the diagonal blocks matrices
-                                     * in #inverse#. This costs some 
-                                     * additional memory (for DG
-                                     * methods about 1/3 of that used for
-                                     * the matrix) but it
-                                     * makes the preconditioning much faster.
-                                     */
-    void invert_diagblocks();
-
-                                    /**
-                                     * Block SOR. Make sure that the right block size
-                                     * of the matrix is set by #set_block_size#
-                                     * before calling this function.
-                                     *
-                                     * BlockSOR will automatically use the
-                                     * inverse matrices if they exist, if not
-                                     * then BlockSOR will waste much time
-                                     * inverting the diagonal block
-                                     * matrices in each preconditioning step.
-                                     *
-                                     * For matrices which are
-                                     * empty above the diagonal blocks
-                                     * BlockSOR is a direct solver.
-                                     */
-    void precondition_BlockSOR (dVector &dst, const dVector &src) const;
-    
-                                    /**
-                                     * Set the right block size before calling
-                                     * #precondition_BlockSOR#.
-                                     * If block_size==1 BlockSOR is the same as SOR.
-                                     */
-    void set_block_size (const unsigned int bsize);
-
-                                    /**
-                                     * Gives back the size of the blocks.
-                                     */
-    unsigned int block_size() const;
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcWrongBlockSize,
-                   int, int,
-                   << "The blocksize " << arg1
-                   << " and the size of the matrix " << arg2
-                   << " do not match.");
-
-    DeclException2 (ExcWrongInverses,
-                   int, int,
-                   << "There are " << arg1
-                   << " inverse matrices but " << arg2
-                   << " cells.");
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInverseMatricesDoNotExist);
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInverseMatricesAlreadyExist);
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcBlockSizeNotSet);
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInternalError);
-    
-  private:
-                                    /**
-                                     * size of the blocks.
-                                     */
-    unsigned int blocksize;
-
-                                    /**
-                                     * stores the inverse matrices of
-                                     * the diagonal blocks matrices
-                                     */
-    vector<dFMatrix> inverse;
-};
-
-
-
-/*----------------------------   dblocksmatrix.h     ---------------------------*/
-/* end of #ifndef __dblocksmatrix_H */
-#endif
-/*----------------------------   dblocksmatrix.h     ---------------------------*/
diff --git a/deal.II/lac/include/lac/dfmatrix.h b/deal.II/lac/include/lac/dfmatrix.h
deleted file mode 100644 (file)
index a5fda7a..0000000
+++ /dev/null
@@ -1,568 +0,0 @@
-/*----------------------------   dfmatrix.h     ---------------------------*/
-/*      $Id$                 */
-#ifndef __dfmatrix_H
-#define __dfmatrix_H
-/*----------------------------   dfmatrix.h     ---------------------------*/
-
-// This file is part of the DEAL Library
-// DEAL is Copyright(1995) by
-// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
-// Revised by Wolfgang Bangerth
-
-
-#include <base/exceptions.h>
-
-
-// forward declarations
-class dVector;
-class iVector;
-
-
-
-/**
- *  Double precision full Matrix.
- *  Memory for Components is supplied explicitly <p>
- *  ( ! Amount of memory needs not to comply with actual dimension due to reinitializations ! ) <p>
- *  - all necessary methods for matrices are supplied <p>
- *  - operators available are '=' and '( )' <p>
- *  CONVENTIONS for used 'equations' : <p>
- *  - THIS matrix is always named 'A' <p>
- *  - matrices are always uppercase , vectors and scalars are lowercase <p>
- *  - Transp(A) used for transpose of matrix A
- *
- */
-class dFMatrix
-{
-                                    /**
-                                     * Component-array.
-                                     */
-    double* val;
-                                    /** 
-                                     * Dimension. Actual number of Columns
-                                     */
-    unsigned int dim_range;
-                                    /**
-                                     * Dimension. Actual number of Rows
-                                     */
-    unsigned int dim_image;
-                                    /**
-                                     * Dimension. Determines amount of reserved memory
-                                     */
-    unsigned int val_size;
-
-                                    /**
-                                     * Initialization   . initialize memory for Matrix <p>
-                                     * ( m rows , n columns )
-                                     */
-    void init (const unsigned int m, const unsigned int n);
-    
-                                    /**
-                                     * Return a read-write reference to the
-                                     * element #(i,j)#.
-                                     *
-                                     * This function does no bounds checking.
-                                     */
-    double& el (const unsigned int i, const unsigned int j);
-    
-                                    /**
-                                     * Return the value of the element #(i,j)#.
-                                     *
-                                     * This function does no bounds checking.
-                                     */
-    double el (const unsigned int i, const unsigned int j) const;
-    
-    
-  public:
-                                    /**@name 1: Basic Object-handling */
-                                    //@{
-                                    /**
-                                     * Constructor. Initialize the matrix as
-                                     * a square matrix with dimension #n#.
-                                     */
-    explicit dFMatrix (const unsigned int n = 1);
-    
-                                    /**
-                                     * Constructor. Initialize the matrix as
-                                     * a rectangular #m# times #n# matrix.
-                                     */
-    dFMatrix (const unsigned int m, const unsigned int n);
-    
-                                    /** 
-                                     * Copy constructor. Be very careful with
-                                     * this constructor, since it may take a
-                                     * huge amount of computing time for large
-                                     * matrices!!
-                                     */
-    dFMatrix (const dFMatrix&);
-
-                                    /**
-                                     * Destructor. Release all memory.
-                                     */
-    ~dFMatrix();
-    
-                                    /**
-                                     * Comparison operator. Be careful with
-                                     * this thing, it may eat up huge amounts
-                                     * of computing time! It is most commonly
-                                     * used for internal consistency checks
-                                     * of programs.
-                                     */
-    bool operator == (const dFMatrix &) const;
-
-                                    /**
-                                     *  A = B           . Copy all elements
-                                     */
-    dFMatrix& operator = (const dFMatrix& B);
-    
-    
-                                    /**
-                                     *  U(0-m,0-n) = s  . Fill all elements
-                                     */
-    void fill (const dFMatrix& src,
-              const unsigned int i=0, const unsigned int j=0);
-    
-                                    /**
-                                     * Change  Dimension.
-                                     * Set dimension to (m,n) <p>
-                                     * ( reinit rectangular matrix )
-                                     */
-    void reinit (const unsigned int m, const unsigned int n);
-    
-                                    /**
-                                     * Change  Dimension.
-                                     * Set dimension to (n,n) <p>
-                                     * ( reinit quadratic matrix )
-                                     */
-    void reinit (const unsigned int n);
-    
-                                    /**
-                                     * Adjust  Dimension.
-                                     * Set dimension to ( m(B),n(B) ) <p>
-                                     * ( adjust to dimensions of another matrix B )
-                                     */
-    void reinit (const dFMatrix &B);
-    
-                                    /**
-                                     * Return number of rows of this matrix.
-                                     * To remember: this matrix is an
-                                     * $m \times n$-matrix.
-                                     */
-    unsigned int m () const;
-    
-                                    /**
-                                     * Return number of columns of this matrix.
-                                     * To remember: this matrix is an
-                                     * $m \times n$-matrix.
-                                     */
-    unsigned int n () const;
-
-                                    /**
-                                     * Return whether the matrix contains only
-                                     * elements with value zero. This function
-                                     * is mainly for internal consistency
-                                     * check and should seldomly be used when
-                                     * not in debug mode since it uses quite
-                                     * some time.
-                                     */
-    bool all_zero () const;
-
-                                    //@}
-    
-    
-                                    /**@name 2: Data-Access
-                                     */
-                                    //@{
-                                    /**
-                                     *   Access Elements. returns element at relative 'address' i <p>
-                                     *   ( -> access to A(i/n , i mod n) )
-                                     */
-    double el (const unsigned int i) const;
-    
-                                    /**
-                                     * Return the value of the element #(i,j)#.
-                                     * Does the same as the #el(i,j)# function
-                                     * but does bounds checking.
-                                     */
-    double operator() (const unsigned int i, const unsigned int j) const;
-    
-                                    /**
-                                     * Return a read-write reference to
-                                     * the element #(i,j)#.
-                                     * Does the same as the #el(i,j)# function
-                                     * but does bounds checking.
-                                     */
-    double& operator() (const unsigned int i, const unsigned int j);
-    
-                                    /**
-                                     * Set all entries in the matrix to
-                                     * zero.
-                                     */
-    void clear ();
-                                    //@}
-    
-    
-                                    /**@name 3: Basic applications on matrices
-                                     */
-                                    //@{
-                                    /**
-                                     *  A+=B            . Simple addition
-                                     */
-    void add (const double s, const dFMatrix& B);
-
-                                    /**
-                                     * A+=Transp(B).
-                                     * Simple addition of the transpose of B to this
-                                     */
-    void Tadd (const double s, const dFMatrix& B);
-    
-                                    /**
-                                     * C=A*B.
-                                     * Matrix-matrix-multiplication 
-                                     */
-    void mmult (dFMatrix& C, const dFMatrix& B) const;
-    
-                                    /**
-                                     * C=Transp(A)*B.
-                                     * Matrix-matrix-multiplication using
-                                     * transpose of this
-                                     */
-    void Tmmult (dFMatrix& C, const dFMatrix& B) const;
-    
-                                    /**
-                                     *  w (+)= A*v.
-                                     *  Matrix-vector-multiplication ; <p>
-                                     *  ( application of this to a vector v )
-                                     *  flag adding=true : w+=A*v
-                                     */
-    void vmult (dVector& w, const dVector& v, const bool adding=false) const;
-    
-                                    /**
-                                     *  w (+)= Transp(A)*v.
-                                     *  Matrix-vector-multiplication ; <p>
-                                     *  (application of transpose of this to a vector v)
-                                     *  flag adding=true : w+=A*v
-                                     */
-    void Tvmult (dVector& w, const dVector& v, const bool adding=false) const;
-
-                                    /**
-                                     * Return 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.
-                                     *
-                                     * Note the order in which the matrix
-                                     * appears. For non-symmetric matrices
-                                     * there is a difference whether the
-                                     * matrix operates on the first
-                                     * or on the second operand of the
-                                     * scalar product.
-                                     *
-                                     * Obviously, the matrix needs to be square
-                                     * for this operation.
-                                     */
-    double matrix_norm (const dVector &v) const;
-
-                                    /**
-                                     * Build the matrix scalar product
-                                     * #u^T M v#. This function is mostly
-                                     * useful when building the cellwise
-                                     * scalar product of two functions in
-                                     * the finite element context.
-                                     */
-    double matrix_scalar_product (const dVector &u, const dVector &v) const;
-    
-                                    /**
-                                     * A=Inverse(A). Inversion of (*this) by
-                                     * Gauss-Jordan-algorithm
-                                     */
-    void gauss_jordan ();
-
-                                    /**
-                                      * Computes the determinant of a matrix.
-                                      * This is only implemented for one two and
-                                      * three dimensions, since for higher
-                                      * dimensions the numerical work explodes.
-                                      * Obviously, the matrix needs to be square
-                                      * for this function.
-                                      */
-    double determinant () const;
-
-                                    /**
-                                     * Compute the quadratic matrix norm.
-                                     * Return value is the root of the square
-                                     * sum of all matrix entries.
-                                     */
-    double norm2 () const;
-                                    /**
-                                     * Assign the inverse of the given
-                                     * matrix to #*this#. This function is
-                                     * only implemented (hardcoded) for
-                                     * square matrices of dimension one,
-                                     * two and three.
-                                     */
-    void invert (const dFMatrix &M);
-                                    //@}
-
-
-                                    /**@name 4: Basic applications on Rows or Columns
-                                     */
-                                    //@{
-                                    /**
-                                     *  A(i,1-n)+=s*A(j,1-n).
-                                     * Simple addition of rows of this
-                                     */
-    void add_row (const unsigned int i, const double s, const unsigned int j);
-
-                                    /**
-                                     *  A(i,1-n)+=s*A(j,1-n)+t*A(k,1-n).
-                                     *  Multiple addition of rows of this
-                                     */
-    void add_row (const unsigned int i,
-                 const double s, const unsigned int j,
-                 const double t, const unsigned int k);
-
-                                    /**
-                                     *  A(1-n,i)+=s*A(1-n,j).
-                                     *  Simple addition of columns of this
-                                     */
-    void add_col (const unsigned int i, const double s, const unsigned int j);
-
-                                    /**
-                                     *  A(1-n,i)+=s*A(1-n,j)+t*A(1-n,k).
-                                     *  Multiple addition of columns of this
-                                     */
-    void add_col (const unsigned int i,
-                 const double s, const unsigned int j,
-                 const double t, const unsigned int k);
-
-                                    /**
-                                     * Swap  A(i,1-n) <-> A(j,1-n).
-                                     * Swap rows i and j of this
-                                     */
-    void swap_row (const unsigned int i, const unsigned int j);
-
-                                    /**
-                                     *  Swap  A(1-n,i) <-> A(1-n,j).
-                                     *  Swap columns i and j of this
-                                     */
-    void swap_col (const unsigned int i, const unsigned int j);
-                                    //@}
-
-
-                                    /**@name 5: Mixed stuff. Including more
-                                     *  applications on matrices
-                                     */
-                                    //@{
-                                    /**
-                                     *  w=b-A*v.
-                                     *  Residual calculation , returns |w|
-                                     */
-    double residual (dVector& w, const dVector& v, const dVector& b) const;
-
-                                    /**
-                                     *  Inversion of lower triangle .
-                                     */
-    void forward (dVector& dst, const dVector& src) const;
-
-                                    /**
-                                     *  Inversion of upper triangle .
-                                     */
-    void backward (dVector& dst, const dVector& src) const;
-
-                                    /**
-                                     * QR - factorization of a matrix.
-                                     * The orthogonal transformation Q is
-                                     * applied to the vector y and this matrix. <p>
-                                     * After execution of householder, the upper
-                                     *  triangle contains the resulting matrix R, <p>
-                                     * the lower the incomplete factorization matrices.
-                                     * 
-                                     * #householder(src); backward(dst, src);# gives
-                                     * the solution #dst# of the linear system
-                                     * #(*this)dst=src#.
-                                     *
-                                     * Note that #src# and #(*this)# (i.e. the
-                                     * matrix itself) is changed in
-                                     * the process of the #householder(src)# function!!
-                                     */
-    void householder (dVector& src);
-
-                                    /**
-                                     * Least - Squares - Approximation by QR-factorization.
-                                     *
-                                     * Note that #src# and #(*this)# (i.e. the
-                                     * matrix itself) is changed in
-                                     * the process of this function!!
-                                     */
-    double least_squares (dVector& dst, dVector& src);
-
-                                    /**
-                                     *  A(i,i)+=B(i,1-n). Addition of complete
-                                     *  rows of B to diagonal-elements of this ; <p>
-                                     *  ( i = 1 ... m )
-                                     */
-    void add_diag (const double s, const dFMatrix& B);
-
-                                    /**
-                                     *  A(i,i)+=s  i=1-m.
-                                     * Add constant to diagonal elements of this
-                                     */
-    void diagadd (const double& src);
-
-                                    /**
-                                     *  w+=part(A)*v. Conditional partial
-                                     *  Matrix-vector-multiplication <p>
-                                     *  (used elements of v determined by x)
-                                     */
-    void gsmult (dVector& w, const dVector& v, const iVector& x) const;
-
-
-                                    /**
-                                     * Output of the matrix in user-defined format.
-                                     */
-    void print (ostream& s, int width=5, int precision=2) const;
-
-                                    /**
-                                     * Print the matrix in the usual format,
-                                     * i.e. as a matrix and not as a list of
-                                     * nonzero elements. For better
-                                     * readability, zero elements
-                                     * are displayed as empty space.
-                                     *
-                                     * Each entry is printed in scientific
-                                     * format, with one pre-comma digit and
-                                     * the number of digits given by
-                                     * #precision# after the comma, with one
-                                     * space following.
-                                     * The precision defaults to four, which
-                                     * suffices for most cases. The precision
-                                     * and output format are {\it not}
-                                     * properly reset to the old values
-                                     * when the function exits.
-                                     *
-                                     * You should be aware that this function
-                                     * may produce {\bf large} amounts of
-                                     * output if applied to a large matrix!
-                                     * Be careful with it.
-                                     */
-    void print_formatted (ostream &out,
-                         const unsigned int presicion=3) const;
-                                    //@}
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcInvalidIndex,
-                   int, int,
-                   << "The given index " << arg1
-                   << " should be less than " << arg2 << ".");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcDimensionMismatch,
-                   int, int,
-                   << "The two dimensions " << arg1 << " and " << arg2
-                   << " do not match here.");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNotQuadratic);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNotRegular);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInternalError);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException3 (ExcInvalidDestination,
-                   int, int, int,
-                   << "Target region not in matrix: size in this direction="
-                   << arg1 << ", size of new matrix=" << arg2
-                   << ", offset=" << arg3);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException1 (ExcNotImplemented,
-                   int,
-                   << "This function is not implemented for the given"
-                   << " matrix dimension " << arg1);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcIO);
-};
-
-
-
-
-
-/*-------------------------Inline functions -------------------------------*/
-
-
-inline
-double & dFMatrix::el (const unsigned int i, const unsigned int j)  {
-  return val[i*dim_range+j];
-};
-
-
-
-inline
-double dFMatrix::el (const unsigned int i, const unsigned int j) const {
-  return val[i*dim_range+j];
-};
-
-
-
-inline
-unsigned int dFMatrix::m() const {
-  return dim_image;
-};
-
-
-
-inline
-unsigned int dFMatrix::n() const {
-  return dim_range;
-};
-
-
-
-inline
-double dFMatrix::el (const unsigned int i) const {
-  return val[i];
-};
-
-
-
-inline
-double dFMatrix::operator() (const unsigned int i, const unsigned int j) const {  
-  Assert (i<dim_image, ExcInvalidIndex (i, dim_image));
-  Assert (j<dim_range, ExcInvalidIndex (i, dim_range));
-  return el(i,j);
-};
-
-
-
-inline
-double & dFMatrix::operator() (const unsigned int i, const unsigned int j) {
-  Assert (i<dim_image, ExcInvalidIndex (i, dim_image));
-  Assert (j<dim_range, ExcInvalidIndex (j, dim_range));
-  return el(i,j);
-}
-
-
-
-
-/*----------------------------   dfmatrix.h     ---------------------------*/
-/* end of #ifndef __dfmatrix_H */
-#endif
-/*----------------------------   dfmatrix.h     ---------------------------*/
diff --git a/deal.II/lac/include/lac/dsmatrix.h b/deal.II/lac/include/lac/dsmatrix.h
deleted file mode 100644 (file)
index 30484ce..0000000
+++ /dev/null
@@ -1,860 +0,0 @@
-/*----------------------------   dsmatrix.h     ---------------------------*/
-/*      $Id$                 */
-#ifndef __dsmatrix_H
-#define __dsmatrix_H
-/*----------------------------   dsmatrix.h     ---------------------------*/
-
-
-// This file is part of the DEAL Library
-// DEAL is Copyright(1995) by
-// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
-// Revised by Wolfgang Bangerth
-
-
-#include <base/exceptions.h>
-
-
-//forward declarations
-class dVector;
-class iVector;
-class ostream;
-
-
-
-/*
-CLASS
-   dSMatrixStruct
-
-   @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth
-   */
-class dSMatrixStruct
-{
-  private:
-                                    /**
-                                     * Copy constructor, made private in order to
-                                     * prevent copying such an object which does
-                                     * not make much sense because you can use
-                                     * a structure like this for more than one
-                                     * matrix.
-                                     *
-                                     * Because it is not needed, this function
-                                     * is not implemented.
-                                     */
-    dSMatrixStruct (const dSMatrixStruct &);
-    
-  public:
-                                    /**
-                                     * Initialize the matrix empty, i.e. with
-                                     * no memory allocated. This is useful if
-                                     * you want such objects as member
-                                     * variables in other classes. You can make
-                                     * the structure usable by calling the
-                                     * #reinit# function.
-                                     */
-    dSMatrixStruct ();
-    
-                                    /**
-                                     * Initialize a rectangular matrix with
-                                     * #m# rows and #n# columns,
-                                     * with at most #max_per_row#
-                                     * nonzero entries per row.
-                                     */
-    dSMatrixStruct (const unsigned int m,
-                   const unsigned int n,
-                   const unsigned int max_per_row);
-    
-                                    /**
-                                     * Initialize a square matrix of dimension
-                                     * #n# with at most #max_per_row#
-                                     * nonzero entries per row.
-                                     */
-    dSMatrixStruct (const unsigned int n,
-                   const unsigned int max_per_row);
-
-                                    /**
-                                     * Destructor.
-                                     */
-    ~dSMatrixStruct ();
-    
-                                    /**
-                                     * Reallocate memory and set up data
-                                     * structures for a new matrix with
-                                     * #m# rows and #n# columns,
-                                     * with at most #max_per_row#
-                                     * nonzero entries per row.
-                                     *
-                                     * If #m*n==0# all memory is freed,
-                                     * resulting in a total reinitialization
-                                     * of the object. If it is nonzero, new
-                                     * memory is only allocated if the new
-                                     * size extends the old one. This is done
-                                     * to save time and to avoid fragmentation
-                                     * of the heap.
-                                     */
-    void reinit (const unsigned int m,
-                const unsigned int n,
-                const unsigned int max_per_row);
-
-                                    /**
-                                     * This function compresses the sparsity
-                                     * structure that this object represents.
-                                     * It does so by eliminating unused
-                                     * entries and sorting the remaining
-                                     * ones to allow faster access by usage
-                                     * of binary search algorithms. A special
-                                     * sorting scheme is used for the diagonal
-                                     * entry of square matrices, which is
-                                     * always the first entry of each row.
-                                     *
-                                     * #dSMatrix# objects require the
-                                     * #dSMatrixStruct# objects they are
-                                     * initialized with to be compressed, to
-                                     * reduce memory requirements.
-                                     */
-    void compress ();
-
-                                    /**
-                                     * Return whether the object is empty. It
-                                     * is empty if no memory is allocated,
-                                     * which is the same as that both
-                                     * dimensions are zero.
-                                     */
-    bool empty () const;
-    
-
-                                    /**
-                                     * Return the index of the matrix
-                                     * element with row number #i# and
-                                     * column number #j#. If the matrix
-                                     * element is not a nonzero one,
-                                     * return -1.
-                                     *
-                                     * This function is usually called
-                                     * by the #operator()# of the
-                                     * #dSMatrix#. It shall only be
-                                     * called for compressed sparsity
-                                     * patterns, since in this case
-                                     * searching whether the entry
-                                     * exists can be done quite fast
-                                     * with a binary sort algorithm
-                                     * because the column numbers are
-                                     * sorted.
-                                     */
-    int operator() (const unsigned int i, const unsigned int j) const;
-
-                                    /**
-                                     * Add a nonzero entry to the matrix.
-                                     * This function may only be called
-                                     * for non-compressed sparsity patterns.
-                                     *
-                                     * If the entry already exists, nothing
-                                     * bad happens.
-                                     */
-    void add (const unsigned int i, const unsigned int j);
-    
-                                    /**
-                                     * This matrix adds a whole connectivity
-                                     * list to the sparsity structure
-                                     * respresented by this object. It assumes
-                                     * the #rowcols# array to be a list of
-                                     * indices which are all linked together,
-                                     * i.e. all entries
-                                     * #(rowcols[i], rowcols[j])# for all
-                                     * #i,j=0...n# for this sparsity pattern
-                                     * are created. #n# is assumed to be the
-                                     * number of elements pointed to by
-                                     * #rowcols#.
-                                     */
-    void add_matrix (const unsigned int n, const int* rowcols);
-
-                                    //////////
-    void add_matrix (const unsigned int m, const unsigned int n,
-                    const int* rows, const int* cols);
-                                    //////////
-    void add_matrix (const iVector& rowcols);
-                                    //////////
-    void add_matrix (const iVector& rows, const iVector& cols);
-
-                                    /**
-                                     * Print the sparsity of the matrix
-                                     * in a format that #gnuplot# understands
-                                     * and which can be used to plot the
-                                     * sparsity pattern in a graphical
-                                     * way. The format consists of pairs
-                                     * #i j# of nonzero elements, each
-                                     * representing one entry of this
-                                     * matrix, one per line of the output
-                                     * file. Indices are counted from
-                                     * zero on, as usual. Since sparsity
-                                     * patterns are printed in the same
-                                     * way as matrices are displayed, we
-                                     * print the negative of the column
-                                     * index, which means that the
-                                     * #(0,0)# element is in the top left
-                                     * rather than in the bottom left
-                                     * corner.
-                                     *
-                                     * Print the sparsity pattern in
-                                     * gnuplot by setting the data style
-                                     * to dots or points and use the
-                                     * #plot# command.
-                                     */
-    void print_gnuplot (ostream &out) const;
-
-                                    /**
-                                     * Return number of rows of this
-                                     * matrix, which equals the dimension
-                                     * of the image space.
-                                     */
-    unsigned int n_rows () const;
-
-                                    /**
-                                     * Return number of columns of this
-                                     * matrix, which equals the dimension
-                                     * of the range space.
-                                     */
-    unsigned int n_cols () const;
-
-                                    /**
-                                     * Compute the bandwidth of the matrix
-                                     * represented by this structure. The
-                                     * bandwidth is the maximum of
-                                     * $|i-j|$ for which the index pair
-                                     * $(i,j)$ represents a nonzero entry
-                                     * of the matrix.
-                                     */
-    unsigned int bandwidth () 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.
-                                     *
-                                     * This function may only be called if the
-                                     * matrix struct is compressed. It does not
-                                     * make too much sense otherwise anyway.
-                                     */
-    unsigned int n_nonzero_elements () const;
-
-                                    /**
-                                     * Return whether the structure is
-                                     * compressed or not.
-                                     */
-    bool is_compressed () const;
-    
-                                    /**
-                                     * This is kind of an expert mode: get
-                                     * access to the rowstart array, but
-                                     * readonly.
-                                     *
-                                     * Though the return value is declared
-                                     * #const#, you should be aware that it
-                                     * may change if you call any nonconstant
-                                     * function of objects which operate on
-                                     * it.
-                                     *
-                                     * 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!
-                                     */
-    const unsigned int * get_rowstart_indices () const;
-
-                                    /**
-                                     * This is kind of an expert mode: get
-                                     * access to the colnums array, but
-                                     * readonly.
-                                     *
-                                     * Though the return value is declared
-                                     * #const#, you shoudl be aware that it
-                                     * may change if you call any nonconstant
-                                     * function of objects which operate on
-                                     * it.
-                                     *
-                                     * 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!
-                                     */
-    const int * get_column_numbers () const;
-    
-    
-                                    /**
-                                     * Exception
-                                     */
-    DeclException1 (ExcInvalidNumber,
-                   int,
-                   << "The provided number is invalid here: " << arg1);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcInvalidIndex,
-                   int, int,
-                   << "The given index " << arg1
-                   << " should be less than " << arg2 << ".");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcNotEnoughSpace,
-                   int, int,
-                   << "Upon entering a new entry to row " << arg1
-                   << ": there was no free entry any more. " << endl
-                   << "(Maximum number of entries for this row: "
-                   << arg2 << "; maybe the matrix is already compressed?)");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNotCompressed);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcMatrixIsCompressed);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcEmptyObject);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcInternalError);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcIO);
-
-  private:
-    unsigned int max_dim;
-    unsigned int rows, cols;
-    unsigned int vec_len, max_vec_len;
-    unsigned int max_row_len;
-    unsigned int* rowstart;
-    int* colnums;
-
-                                    /**
-                                     * Store whether the #compress# function
-                                     * was called for this object.
-                                     */
-    bool compressed;
-    
-    friend class dSMatrix;
-};
-
-
-
-
-/*
-CLASS
-   dSMatrix
-
-   @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998
-   */
-class dSMatrix
-{
-  public:
-
-                                    /**
-                                     * 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
-                                     * #reinit(dSMatrixStruct)#.
-                                     */
-    dSMatrix ();
-    
-                                    /**
-                                     * Constructor. Takes the given matrix
-                                     * sparisty structure to represent the
-                                     * sparsity pattern of this matrix. You
-                                     * can change the sparsity pattern later
-                                     * on by calling the #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 #reinit# is not called with a
-                                     * new sparsity structure.
-                                     */
-    dSMatrix (const dSMatrixStruct &sparsity);
-    
-                                    /**
-                                     * Destructor. Free all memory, but do not
-                                     * release the memory of the sparsity
-                                     * structure.
-                                     */
-    virtual ~dSMatrix ();
-    
-
-                                    /**
-                                     * Reinitialize the object but keep to
-                                     * the sparsity pattern previously used.
-                                     * This may be necessary if you #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.
-                                     */
-    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 #reinit# is not called with a
-                                     * new sparsity structure.
-                                     */
-    virtual void reinit (const dSMatrixStruct &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 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;
-    
-                                    /**
-                                     * Set the element #(i,j)# to #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 double value);
-    
-                                    /**
-                                     * Add #value# to the element #(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 double value);
-
-                                    /**
-                                     * 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
-                                     * #operator =#, since this may lead
-                                     * to unwanted usage, e.g. in copy
-                                     * arguments to functions, which should
-                                     * really be arguments by reference.
-                                     *
-                                     * The function returns a reference to
-                                     * #this#.
-                                     */
-    dSMatrix & copy_from (const dSMatrix &);
-
-                                    /**
-                                     * Add #matrix# scaled by #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.
-                                     */
-    void add_scaled (const double factor, const dSMatrix &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 wanted
-                                     * element does not exist in the matrix.
-                                     */
-    double operator () (const unsigned int i, const unsigned int j) const;
-
-                                    /**
-                                     * Return the main diagonal element in
-                                     * the #i#th row. This function throws an
-                                     * error if the matrix is not square.
-                                     *
-                                     * This function is considerably faster
-                                     * than the #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.
-                                     */
-    double diag_element (const unsigned int i) const;
-
-                                    /**
-                                     * This is kind of an expert mode: get
-                                     * access to the #i#th element of this
-                                     * matrix. The elements are stored in
-                                     * a consecutive way, refer to the
-                                     * #dSMatrixStruct# 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!
-                                     */
-    double global_entry (const unsigned int i) const;
-
-                                    /**
-                                     * Same as above, but with write access.
-                                     * You certainly know what you do?
-                                     */
-    double & global_entry (const unsigned int i);
-
-                                    /**
-                                     * Matrix-vector multiplication: let
-                                     * #dst = M*src# with #M# being this matrix.
-                                     */
-    void vmult (dVector& dst, const dVector& src) const;
-    
-                                    /**
-                                     * Matrix-vector multiplication: let
-                                     * #dst = M^T*src# with #M# being this
-                                     * matrix. This function does the same as
-                                     * #vmult# but takes the transposed matrix.
-                                     */
-    void Tvmult (dVector& dst, const dVector& src) const;
-  
-
-                                    /**
-                                     * Return 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.
-                                     *
-                                     * Note the order in which the matrix
-                                     * appears. For non-symmetric matrices
-                                     * there is a difference whether the
-                                     * matrix operates on the first
-                                     * or on the second operand of the
-                                     * scalar product.
-                                     *
-                                     * Obviously, the matrix needs to be square
-                                     * for this operation.
-                                     */
-    double matrix_norm (const dVector &v) const;
-    
-                                    //
-    double residual (dVector& dst, const dVector& x,
-                    const dVector& b) const;
-                                    //
-    void precondition_Jacobi (dVector& dst, const dVector& src,
-                             const double om = 1.) const;
-                                    //
-    void precondition_SSOR (dVector& dst, const dVector& src,
-                           const double om = 1.) const;
-                                    //
-    void precondition_SOR (dVector& dst, const dVector& src,
-                          const double om = 1.) const;
-                                    //
-    void SSOR (dVector& dst, const double om = 1.) const;
-                                    //
-    void SOR (dVector& dst, const double om = 1.) const;
-                                    //
-    void precondition (dVector& dst, const dVector& src) const;
-
-                                    /**
-                                     * Return a (constant) reference to the
-                                     * underlying sparsity pattern of this
-                                     * matrix.
-                                     *
-                                     * Though the return value is declared
-                                     * #const#, you shoudl be aware that it
-                                     * may change if you call any nonconstant
-                                     * function of objects which operate on
-                                     * it.
-                                     */
-    const dSMatrixStruct & get_sparsity_pattern () const;
-
-                                    /**
-                                     * Print the matrix to the given stream,
-                                     * using the format
-                                     * #(line,col) value#, i.e. one
-                                     * nonzero entry of the matrix per line.
-                                     */
-    void print (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 explicitely
-                                     * set to zero are displayed as such.
-                                     *
-                                     * Each entry is printed in scientific
-                                     * format, with one pre-comma digit and
-                                     * the number of digits given by
-                                     * #precision# after the comma, with one
-                                     * space following.
-                                     * The precision defaults to four, which
-                                     * suffices for most cases. The precision
-                                     * and output format are {\it not}
-                                     * properly reset to the old values
-                                     * when the function exits.
-                                     *
-                                     * You should be aware that this function
-                                     * may produce {\bf large} amounts of
-                                     * output if applied to a large matrix!
-                                     * Be careful with it.
-                                     */
-    void print_formatted (ostream &out,
-                         const unsigned int presicion=3) const;
-    
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcNotCompressed);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcMatrixNotInitialized);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcDimensionsDontMatch,
-                   int, int,
-                   << "The dimensions " << arg1 << " and " << arg2
-                   << " do not match properly.");
-                                    /**
-                                     * 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 (ExcIO);
-    
-  private:
-    const dSMatrixStruct * cols;
-    double* val;
-    unsigned int max_len;
-};
-
-
-
-
-
-/*---------------------- Inline functions -----------------------------------*/
-
-inline
-unsigned int dSMatrixStruct::n_rows () const {
-  return rows;
-};
-
-
-
-inline
-unsigned int dSMatrixStruct::n_cols () const {
-  return cols;
-};
-
-
-
-inline
-bool dSMatrixStruct::is_compressed () const {
-  return compressed;
-};
-
-
-
-inline
-const unsigned int * dSMatrixStruct::get_rowstart_indices () const {
-  return rowstart;
-};
-
-
-
-inline
-const int * dSMatrixStruct::get_column_numbers () const {
-  return colnums;
-};
-
-
-
-inline
-unsigned int dSMatrix::m () const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  return cols->rows;
-};
-
-
-
-inline
-unsigned int dSMatrix::n () const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  return cols->cols;
-};
-
-
-
-inline
-void dSMatrix::set (const unsigned int i, const unsigned int j,
-                   const double value) {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert ((cols->operator()(i,j) != -1) || (value == 0.),
-         ExcInvalidIndex(i,j));
-
-  const int index = cols->operator()(i,j);
-
-  if (index >= 0) val[index] = value;
-};
-
-
-
-inline
-void dSMatrix::add (const unsigned int i, const unsigned int j,
-                   const double value) {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert ((cols->operator()(i,j) != -1) || (value == 0.),
-         ExcInvalidIndex(i,j));
-
-  const int index = cols->operator()(i,j);
-  
-  if (index >= 0) val[index] += value;
-};
-
-
-
-
-
-inline
-double dSMatrix::operator () (const unsigned int i, const unsigned int j) const {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (cols->operator()(i,j) != -1,
-         ExcInvalidIndex(i,j));
-  return val[cols->operator()(i,j)];
-};
-
-
-
-inline
-double dSMatrix::diag_element (const unsigned int i) const {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
-  Assert (i<max_len, ExcInvalidIndex1(i));
-  
-                                  // Use that the first element in each
-                                  // row of a square matrix is the main
-                                  // diagonal
-  return val[cols->rowstart[i]];
-};
-
-
-
-inline
-double dSMatrix::global_entry (const unsigned int j) const {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  return val[j];
-};
-
-
-
-inline
-double & dSMatrix::global_entry (const unsigned int j) {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  return val[j];
-};
-
-
-
-/*----------------------------   dsmatrix.h     ---------------------------*/
-/* end of #ifndef __dsmatrix_H */
-#endif
-/*----------------------------   dsmatrix.h     ---------------------------*/
-
-
diff --git a/deal.II/lac/include/lac/dvector.h b/deal.II/lac/include/lac/dvector.h
deleted file mode 100644 (file)
index 39f61ea..0000000
+++ /dev/null
@@ -1,483 +0,0 @@
-/*----------------------------   dvector.h     ---------------------------*/
-/*      $Id$                 */
-#ifndef __dvector_H
-#define __dvector_H
-/*----------------------------   dvector.h     ---------------------------*/
-
-// This file is part of the DEAL Library
-// DEAL is Copyright(1995) by
-// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
-// Revised by Wolfgang Bangerth
-
-#include <base/exceptions.h>
-#include <cstdio>
-
-
-
-/**
- *  Double precision Vector.
- *  Memory for Components is supplied explicitly <p>
- *  ( ! Amount of memory needs not to comply with actual dimension due to reinitializations ! ) <p>
- *  - all necessary methods for Vectors are supplied <p>
- *  - operators available are `=` , `*` and `( )` <p>
- *  CONVENTIONS for used `equations` : <p>
- *  - THIS vector is always named `U` <p>
- *  - vectors are always uppercase , scalars are lowercase
- *
- * @author Roland Becker, Guido Kanschat, Franz-Theo Suttmeier, revised and extended by Wolfgang Bangerth, documented by Klaus Mampel and Wolfgang Bangerth
- */
-class dVector {
-  friend class dFMatrix;
-
-  protected:
-
-                                    /**
-                                     * Dimension. Actual number of components
-                                     * contained in the vector.
-                                     * Get this number by calling #size()#.
-                                     */
-    unsigned int dim;
-
-                                    /**
-                                     * Amount of memory actually reserved for
-                                     * this vector. This number may be greater
-                                     * than #dim# if a #reinit# was called with
-                                     * less memory requirements than the vector
-                                     * needed last time. At present #reinit#
-                                     * does not free memory when the number of
-                                     * needed elements is reduced.
-                                     */
-    unsigned int maxdim;
-
-                                    /**
-                                     * Pointer to the array of components.
-                                     */
-    double *val;
-
-  public:
-
-                                    /**
-                                     * Declare iterator types just like those
-                                     * for the C++ standard library:
-                                     *
-                                     * Data type stored by this container.
-                                     */
-    typedef double value_type;
-
-                                    /**
-                                     * Declare standard types used in all
-                                     * containers.
-                                     */
-    typedef value_type* pointer;
-    typedef const value_type* const_pointer;
-    typedef value_type* iterator;
-    typedef const value_type* const_iterator;
-    typedef value_type& reference;
-    typedef const value_type& const_reference;
-    typedef size_t size_type;
-    
-
-                                    /**
-                                     * @name 1: Basic Object-handling 
-                                     */
-                                    //@{
-                                    /**
-                                     *  Dummy-Constructor. Dimension=0
-                                     */
-    dVector ();
-    
-                                    /**
-                                     *   Copy-Constructor. Dimension set to that of V , <p>
-                                     *                     all components are copied from V
-                                     */
-    dVector (const dVector& V);
-    
-                                    /**
-                                     * Constructor. Set dimension to #n# and
-                                     * initialize all elements with zero.
-                                     */
-    dVector (const unsigned int n);
-    
-                                    /**
-                                     *         Destructor. Clears memory
-                                     */
-    ~dVector ();
-
-                                    /**
-                                     * Set all entries to zero. Equivalent to
-                                     * #v = 0#, but more obvious and faster.
-                                     * Note that this function does not change
-                                     * the size of the vector, unlike the
-                                     * STL's #vector<>::clear# function.
-                                     */
-    void clear ();
-    
-                                    /**
-                                     *  U(0-N) = s       . Fill all components
-                                     */
-    dVector& operator= (const double s);
-    
-                                    /**
-                                     *  U = V            . Copy all components
-                                     */
-    dVector& operator= (const dVector& V);
-    
-                                    /**
-                                     *  U = U * V        . Scalar Produkt
-                                     */
-    double operator* (const dVector& V) const;
-
-                                    /**
-                                     * Return square of the l2-norm.
-                                     */
-    double norm_sqr () const;
-
-                                    /**
-                                     * Return the mean value of the elements of
-                                     * this vector.
-                                     */
-    double mean_value () const;
-
-                                    /**
-                                     * Return the l1-norm of the vector, i.e.
-                                     * the sum of the absolute values.
-                                     */
-    double l1_norm () const;
-
-                                    /**
-                                     * Return the l2-norm of the vector, i.e.
-                                     * the square root of the sum of the
-                                     * squares of the elements.
-                                     */
-    double l2_norm () const;
-
-                                    /**
-                                     * Return the maximum absolute value of the
-                                     * elements of this vector.
-                                     */
-    double linfty_norm () const;
-    
-    
-                                    /**
-                                     * Change the dimension of the vector to
-                                     * #N#. The reserved memory for this vector
-                                     * remains unchanged if possible, to make
-                                     * things faster, but this may waste some
-                                     * memory, so take this in the back of your
-                                     * head.
-                                     * However, if #N==0# all memory is freed,
-                                     * i.e. if you want to resize the vector
-                                     * and release the memory not needed, you
-                                     * have to first call #reinit(0)# and then
-                                     * #reinit(N)#. This cited behaviour is
-                                     * analogous to that of the STL containers.
-                                     *
-                                     * On #fast==false#, the vector is filled by
-                                     * zeros.
-                                     */ 
-    void reinit (const unsigned int N, const bool fast=false);
-    
-                                    /**
-                                     * Change the dimension to that of the
-                                     * vector #V#. The same applies as for
-                                     * the other #reinit# function.
-                                     *
-                                     * The elements of #V# are not copied, i.e.
-                                     * this function is the same as calling
-                                     * #reinit (V.size(), fast)#.
-                                     */
-    void reinit (const dVector& V, const bool fast=false);
-    
-                                    /**
-                                     * Return dimension of the vector. This
-                                     * function was formerly called #n()#, but
-                                     * was renamed to get the #dVector# class
-                                     * closer to the C++ standard library's
-                                     * #vector# container.
-                                     */
-    unsigned int size () const;
-
-                                    /**
-                                     * Return whether the vector contains only
-                                     * elements with value zero. This function
-                                     * is mainly for internal consistency
-                                     * check and should seldomly be used when
-                                     * not in debug mode since it uses quite
-                                     * some time.
-                                     */
-    bool all_zero () const;
-    
-                                    /**
-                                     * Make the #dVector# class a bit like the
-                                     * #vector<># class of the C++ standard
-                                     * library by returning iterators to
-                                     * the start and end of the elements of this
-                                     * vector.
-                                     */
-    iterator begin ();
-
-                                    /**
-                                     * Return constant iterator to the start of
-                                     * the vectors.
-                                     */
-    const_iterator begin () const;
-
-                                    /**
-                                     * Return an iterator pointing to the
-                                     * element past the end of the array.
-                                     */
-    iterator end ();
-
-                                    /**
-                                     * Return a constant iterator pointing to
-                                     * the element past the end of the array.
-                                     */
-    const_iterator end () const;  
-                                    //@}
-    
-    
-                                    /**
-                                     * @name 2: Data-Access
-                                     */
-                                    //@{
-                                    /**
-                                     *  Access Components. returns U(i) , 
-                                     *             INLINE
-                                     */
-    double operator() (const unsigned int i) const;
-    
-                                    /**
-                                     *  Access Components. returns U(i) , 
-                                     *             INLINE
-                                     */
-    double& operator() (const unsigned int i);
-                                    //@}
-    
-    
-                                    /**
-                                     * @name 3: Modification of vectors
-                                     */
-                                    //@{
-                                    /**
-                                     * Fast equivalent to #U.add(1, V)#.
-                                     */
-    dVector & operator += (const dVector &V);
-
-                                    /**
-                                     * Fast equivalent to #U.add(-1, V)#.
-                                     */
-    dVector & operator -= (const dVector &V);
-
-                                    /**
-                                     * U(0-DIM)+=s.
-                                     * Addition of #s# to all components. Note
-                                     * that #s# is a scalar and not a vector.
-                                     */
-    void add (const double s);
-    
-                                    /**
-                                     * U+=V.
-                                     * Simple vector addition, equal to the
-                                     * #operator +=#.
-                                     */
-    void add (const dVector& V);
-    
-                                    /**
-                                     * U+=a*V.
-                                     * Simple addition of a scaled vector.
-                                     */
-    void add (const double a, const dVector& V);
-    
-                                    /**
-                                     * U+=a*V+b*W.
-                                     * Multiple addition of scaled vectors.
-                                     */
-    void add (const double a, const dVector& V,
-             const double b, const dVector& W);
-    
-                                    /**
-                                     * U=s*U+V.
-                                     * Scaling and simple vector addition.
-                                     */
-    void sadd (const double s, const dVector& V);
-    
-                                    /**
-                                     * U=s*U+a*V.
-                                     * Scaling and simple addition.
-                                     */
-    void sadd (const double s, const double a, const dVector& V);
-    
-                                    /**
-                                     * U=s*U+a*V+b*W.
-                                     * Scaling and multiple addition.
-                                     */
-    void sadd (const double s, const double a,
-              const dVector& V, const double b, const dVector& W);
-    
-                                    /**
-                                     * U=s*U+a*V+b*W+c*X.
-                                     * Scaling and multiple addition.
-                                     */
-    void sadd (const double s, const double a,
-              const dVector& V, const double b, const dVector& W, 
-              const double c, const dVector& X);
-    
-                                    /**
-                                     * Scale each element of the vector by the
-                                     * given factor. This function was
-                                     * previously called #equ(double)#, which
-                                     * in my eyes is an extremely unintuitive
-                                     * naming and was thus replaced.
-                                     */
-    void scale (const double factor);
-    
-                                    /**
-                                     *  U=a*V. Replacing
-                                     */
-    void equ (const double a, const dVector& V);
-    
-                                    /**
-                                     * U=a*V+b*W.
-                                     * Replacing by sum.
-                                     */
-    void equ (const double a, const dVector& V,
-             const double b, const dVector& W);
-
-                                    /**
-                                     * Compute the elementwise ratio of the
-                                     * two given vectors, that is let
-                                     * #this[i] = a[i]/b[i]#. This is useful
-                                     * for example if you want to compute
-                                     * the cellwise ratio of true to estimated
-                                     * error.
-                                     *
-                                     * This vector is appropriately scaled to
-                                     * hold the result.
-                                     *
-                                     * If any of the #b[i]# is zero, the result
-                                     * is undefined. No attempt is made to
-                                     * catch such situations.
-                                     */
-    void ratio (const dVector &a, const dVector &b);
-                                    //@}
-    
-    
-                                    /**
-                                     * @name 5: Mixed stuff
-                                     */
-                                    //@{
-                                    /**
-                                     *  Output of vector in user-defined format.
-                                     */
-    void print (FILE* fp, const char* format = 0) const;
-    
-                                    /**
-                                     *  Output of vector in user-defined format.
-                                     */
-    void print (const char* format = 0) const;
-
-                                    /**
-                                     * Print to given stream, one element per line.
-                                     */
-    void print (ostream &) const;
-                                    //@}
-
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcDimensionsDontMatch,
-                   int, int,
-                   << "The dimensions " << arg1 << " and " << arg2
-                   << " do not match here.");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException2 (ExcInvalidIndex,
-                   int, int,
-                   << "The given index " << arg1
-                   << " should be less than " << arg2 << ".");
-                                    /**
-                                     * Exception
-                                     */
-    DeclException1 (ExcInvalidNumber,
-                   int,
-                   << "The provided number is invalid here: " << arg1);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcOutOfMemory);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcEmptyVector);
-                                    /**
-                                     * Exception
-                                     */
-    DeclException0 (ExcIO);
-};
-
-
-
-
-
-
-/*----------------------- Inline functions ----------------------------------*/
-
-
-inline unsigned int dVector::size () const
-{
-  return dim;
-}
-
-
-
-inline
-dVector::iterator dVector::begin () {
-  return &val[0];
-};
-
-
-
-inline
-dVector::const_iterator dVector::begin () const {
-  return &val[0];
-};
-
-
-
-inline
-dVector::iterator dVector::end () {
-  return &val[dim];
-};
-
-
-
-inline
-dVector::const_iterator dVector::end () const {
-  return &val[dim];
-};
-
-
-
-inline double dVector::operator() (const unsigned int i) const
-{
-  Assert (i<dim, ExcInvalidIndex(i,dim));
-  return val[i];
-}
-
-
-
-inline double& dVector::operator() (const unsigned int i)
-{
-  Assert (i<dim, ExcInvalidIndex(i,dim));
-  return val[i];
-}
-
-
-
-
-
-/*----------------------------   dvector.h     ---------------------------*/
-/* end of #ifndef __dvector_H */
-#endif
-/*----------------------------   dvector.h     ---------------------------*/
diff --git a/deal.II/lac/source/dblocksmatrix.cc b/deal.II/lac/source/dblocksmatrix.cc
deleted file mode 100644 (file)
index 4e11933..0000000
+++ /dev/null
@@ -1,175 +0,0 @@
-/*----------------------------   dblocksmatrix.cc     ---------------------------*/
-/*      $Id$                 */
-/*----------------------------   dblocksmatrix.cc     ---------------------------*/
-
-#include <lac/dblocksmatrix.h>
-#include <lac/dvector.h>
-
-
-dBlockSMatrix::dBlockSMatrix ():
-               blocksize(0) {};
-
-dBlockSMatrix::~dBlockSMatrix ()
-{
-  if (inverse.size()!=0)
-    inverse.erase(inverse.begin(), inverse.end());
-}
-
-
-void dBlockSMatrix::reinit ()
-{
-  if (inverse.size()!=0)
-    inverse.erase(inverse.begin(), inverse.end());
-  blocksize=0;
-  dSMatrix::reinit ();
-}
-
-
-void dBlockSMatrix::reinit (const dSMatrixStruct &sparsity)
-{
-  if (inverse.size()!=0)
-    inverse.erase(inverse.begin(), inverse.end());
-  blocksize=0;
-  dSMatrix::reinit (sparsity);
-}
-
-
-void dBlockSMatrix::clear ()
-{
-  dSMatrix::clear();
-  if (inverse.size()!=0)
-    inverse.erase(inverse.begin(), inverse.end());
-  blocksize=0;
-}
-
-
-void dBlockSMatrix::set_block_size(unsigned int bsize) {
-  blocksize=bsize;
-}
-
-
-
-unsigned int dBlockSMatrix::block_size() const {
-  return blocksize;
-}
-
-
-
-void dBlockSMatrix::precondition_BlockSOR (dVector &dst, const dVector &src) const
-{
-  Assert (m() == n(), ExcMatrixNotSquare());
-  Assert (blocksize!=0, ExcBlockSizeNotSet());
-  Assert (m()%blocksize==0, ExcWrongBlockSize(blocksize, m()));
-  unsigned int n_cells=m()/blocksize;
-  Assert (inverse.size()==0 || inverse.size()==n_cells,
-         ExcWrongInverses(inverse.size(), n_cells));
-
-  const dSMatrixStruct &spars=get_sparsity_pattern();
-  const unsigned int *rowstart = spars.get_rowstart_indices();
-  const int *columns = spars.get_column_numbers();
-
-  dVector b_cell(blocksize), x_cell(blocksize);
-
-                                      // cell_row, cell_column are the
-                                      // numbering of the blocks (cells).
-                                      // row_cell, column_cell are the local
-                                      // numbering of the unknowns in the
-                                      // blocks.
-                                      // row, column are the global numbering
-                                      // of the unkowns.
-  unsigned int row, column, row_cell, begin_diag_block=0;
-  double b_cell_row;
-
-  if (inverse.size()==0)
-    {
-      dFMatrix M_cell(blocksize);
-      for (unsigned int cell=0; cell<n_cells; ++cell)
-       {
-         for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
-           {
-             b_cell_row=src(row);
-             for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
-               if ((column=static_cast<unsigned int>(columns[j]))
-                   < begin_diag_block)
-                   b_cell_row -= global_entry(j) * dst(column);
-             b_cell(row_cell)=b_cell_row;
-             for (unsigned int column_cell=0, column=cell*blocksize;
-                  column_cell<blocksize; ++column_cell, ++column)
-                 M_cell(row_cell,column_cell)=(*this)(row,column);
-           }
-         M_cell.householder(b_cell);
-         M_cell.backward(x_cell,b_cell);
-                                          // distribute x_cell to dst
-         for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
-           dst(row)=x_cell(row_cell);
-         
-         begin_diag_block+=blocksize;
-       }      
-    }
-  else
-    for (unsigned int cell=0; cell<n_cells; ++cell)
-      {
-       for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
-         {
-           b_cell_row=src(row);
-           for (unsigned int j=rowstart[row]; j<rowstart[row+1]; ++j)
-             if ((column=static_cast<unsigned int>(columns[j])) < begin_diag_block)
-               {
-                 b_cell_row -= global_entry(j) * dst(column);
-               }
-           b_cell(row_cell)=b_cell_row;
-         }
-       inverse[cell].vmult(x_cell, b_cell);
-                                        // distribute x_cell to dst
-       for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
-         dst(row)=x_cell(row_cell);
-       
-       begin_diag_block+=blocksize;
-      }
-}
-
-
-void dBlockSMatrix::invert_diagblocks()
-{
-  Assert (m() == n(), ExcMatrixNotSquare());
-  Assert (inverse.size()==0, ExcInverseMatricesAlreadyExist());
-
-  Assert (blocksize!=0, ExcBlockSizeNotSet());
-  Assert (m()%blocksize==0, ExcWrongBlockSize(blocksize, m()));
-  unsigned int n_cells=m()/blocksize;
-
-  inverse.insert(inverse.begin(), n_cells, dFMatrix(blocksize));
-  
-                                  // cell_row, cell_column are the
-                                  // numbering of the blocks (cells).
-                                  // row_cell, column_cell are the local
-                                  // numbering of the unknowns in the
-                                  // blocks.
-                                  // row, column are the global numbering
-                                  // of the unkowns.
-  dFMatrix M_cell(blocksize);
-
-  for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
-    {
-      for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell, ++row)
-       for (unsigned int column_cell=0, column=cell*blocksize;
-            column_cell<blocksize; ++column_cell, ++column)
-           M_cell(row_cell,column_cell)=(*this)(row,column);
-
-                                      // perhaps #dFMatrix::invert# should
-                                      // be change such that it calls
-                                      // #gauss_jordan()# automatically
-                                      // if blocksize > 4
-      if (blocksize<=4)
-       inverse[cell].invert(M_cell);
-      else
-       {
-         M_cell.gauss_jordan();
-         inverse[cell]=M_cell;
-       }
-    }
-}
-
-
-
-/*----------------------------   dblocksmatrix.cc     ---------------------------*/
diff --git a/deal.II/lac/source/dfmatrix.cc b/deal.II/lac/source/dfmatrix.cc
deleted file mode 100644 (file)
index 040141a..0000000
+++ /dev/null
@@ -1,1333 +0,0 @@
-// $Id$
-
-#include <lac/dvector.h>
-#include <lac/ivector.h>
-#include <lac/dfmatrix.h>
-
-#include <cmath>
-#include <cstdlib>
-#include <cstdio>
-#include <iomanip>
-
-
-
-dFMatrix::dFMatrix (const unsigned int n) {
-  init (n,n);
-};
-
-
-
-dFMatrix::dFMatrix (const unsigned int m, const unsigned int n) {
-  init (m,n);
-};
-
-
-
-dFMatrix::dFMatrix (const dFMatrix &m) 
-{
-  init (m.dim_image, m.dim_range);
-  double       *       p = &val[0];
-  const double *      vp = &m.val[0];
-  const double * const e = &val[dim_image*dim_range];
-
-  while (p!=e)
-    *p++ = *vp++;
-};
-
-
-
-void dFMatrix::init (const unsigned int mm, const unsigned int nn)
-{
-  val_size = nn*mm;
-  val = new double[val_size];
-  dim_range = nn;
-  dim_image = mm;
-  clear ();
-};
-
-
-
-dFMatrix::~dFMatrix () {
-  delete[] val;
-};
-
-
-
-bool dFMatrix::all_zero () const {
-  const double *p = &val[0],
-              *e = &val[n()*m()];
-  while (p!=e)
-    if (*p++ != 0.0)
-      return false;
-
-  return true;
-};
-
-
-
-void dFMatrix::reinit (const unsigned int mm, const unsigned int nn)
-{
-  if (val_size<nn*mm)
-    {
-      delete[] val;
-      init(mm, nn);
-    }
-  else
-    {
-      dim_range = nn;
-      dim_image = mm;
-//      memset(val, 0, sizeof(double)*nn*mm);
-      clear ();
-    }
-}
-
-
-
-void dFMatrix::reinit (const unsigned int n) {
-  reinit (n, n);
-};
-
-
-
-void dFMatrix::reinit (const dFMatrix &B) {
-  reinit (B.m(), B.n());
-};
-
-
-
-void dFMatrix::vmult (dVector& dst, const dVector& src,
-                     const bool adding) const
-{
-  Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
-  Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
-
-  double s;
-  if ((n()==3) && (m()==3))
-  {
-    double s0,s1,s2;
-    s   = src(0);
-    s0  = s*val[0]; s1  = s*val[3]; s2  = s*val[6]; 
-    s   = src(1);
-    s0 += s*val[1]; s1 += s*val[4]; s2 += s*val[7];
-    s   = src(2);
-    s0 += s*val[2]; s1 += s*val[5]; s2 += s*val[8];
-
-    if (!adding)
-    {
-      dst(0) = s0;
-      dst(1) = s1;
-      dst(2) = s2;
-    }
-    else
-    {
-      dst(0) += s0;
-      dst(1) += s1;
-      dst(2) += s2;
-    }
-  }
-  else if ((n()==4) && (m()==4))
-  {
-    double s0,s1,s2,s3;
-    s = src(0);
-    s0  = s*val[0]; s1  = s*val[4]; s2  = s*val[8];  s3  = s*val[12];
-    s = src(1);
-    s0 += s*val[1]; s1 += s*val[5]; s2 += s*val[9];  s3 += s*val[13];
-    s = src(2);
-    s0 += s*val[2]; s1 += s*val[6]; s2 += s*val[10]; s3 += s*val[14];
-    s = src(3);
-    s0 += s*val[3]; s1 += s*val[7]; s2 += s*val[11]; s3 += s*val[15];
-    
-    if (!adding)
-    {
-      dst(0) = s0;
-      dst(1) = s1;
-      dst(2) = s2;
-      dst(3) = s3;
-    }
-    else
-    {
-      dst(0) += s0;
-      dst(1) += s1;
-      dst(2) += s2;
-      dst(3) += s3;
-    }
-  }
-  else if ((n()==8) && (m()==8))
-  {
-    double s0,s1,s2,s3,s4,s5,s6,s7;
-    s = src(0);
-    s0 = s*val[0]; s1 = s*val[8]; s2 = s*val[16]; s3 = s*val[24];
-    s4 = s*val[32]; s5 = s*val[40]; s6 = s*val[48]; s7 = s*val[56];
-    s = src(1);
-    s0 += s*val[1]; s1 += s*val[9]; s2 += s*val[17]; s3 += s*val[25];
-    s4 += s*val[33]; s5 += s*val[41]; s6 += s*val[49]; s7 += s*val[57];
-    s = src(2);
-    s0 += s*val[2]; s1 += s*val[10]; s2 += s*val[18]; s3 += s*val[26];
-    s4 += s*val[34]; s5 += s*val[42]; s6 += s*val[50]; s7 += s*val[58];
-    s = src(3);
-    s0 += s*val[3]; s1 += s*val[11]; s2 += s*val[19]; s3 += s*val[27];
-    s4 += s*val[35]; s5 += s*val[43]; s6 += s*val[51]; s7 += s*val[59];
-    s = src(4);
-    s0 += s*val[4]; s1 += s*val[12]; s2 += s*val[20]; s3 += s*val[28];
-    s4 += s*val[36]; s5 += s*val[44]; s6 += s*val[52]; s7 += s*val[60];
-    s = src(5);
-    s0 += s*val[5]; s1 += s*val[13]; s2 += s*val[21]; s3 += s*val[29];
-    s4 += s*val[37]; s5 += s*val[45]; s6 += s*val[53]; s7 += s*val[61];
-    s = src(6);
-    s0 += s*val[6]; s1 += s*val[14]; s2 += s*val[22]; s3 += s*val[30];
-    s4 += s*val[38]; s5 += s*val[46]; s6 += s*val[54]; s7 += s*val[62];
-    s = src(7);
-    s0 += s*val[7]; s1 += s*val[15]; s2 += s*val[23]; s3 += s*val[31];
-    s4 += s*val[39]; s5 += s*val[47]; s6 += s*val[55]; s7 += s*val[63];
-    
-    if (!adding)
-    {
-      dst(0) = s0;
-      dst(1) = s1;
-      dst(2) = s2;
-      dst(3) = s3;
-      dst(4) = s4;
-      dst(5) = s5;
-      dst(6) = s6;
-      dst(7) = s7;
-    }
-    else
-    {
-      dst(0) += s0;
-      dst(1) += s1;
-      dst(2) += s2;
-      dst(3) += s3;
-      dst(4) += s4;
-      dst(5) += s5;
-      dst(6) += s6;
-      dst(7) += s7;
-    }
-  }
-  else
-  {    
-    double* e = val;
-    const unsigned int size_m = m(),
-                      size_n = n();
-    for (unsigned int i=0; i<size_m; ++i)
-    {
-      s = 0.;
-      for (unsigned int j=0; j<size_n; ++j)
-       s += src(j) * *(e++);
-      if (!adding) dst(i) = s;
-      else dst(i) += s;
-    }
-  }
-}
-
-void dFMatrix::gsmult (dVector& dst, const dVector& src, const iVector& gl) const
-{
-  Assert(n() == m(), ExcNotQuadratic());
-  Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
-  Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
-  Assert(gl.n() == n(), ExcDimensionMismatch(gl.n(), n()));
-
-  double s;
-  if ((n()==3) && (m()==3))
-  {
-    double s0=0.,s1=0.,s2=0.;
-    s   = src(0);
-    if(gl(1)<gl(0)) s1  = s*val[3]; if(gl(2)<gl(0))  s2  = s*val[6]; 
-    s   = src(1);
-    if(gl(0)<gl(1)) s0 += s*val[1]; if(gl(2)<gl(1))  s2 += s*val[7];
-    s   = src(2);
-    if(gl(0)<gl(2)) s0 += s*val[2]; if(gl(1)<gl(2))  s1 += s*val[5]; 
-
-    dst(0) += s0;
-    dst(1) += s1;
-    dst(2) += s2;
-  }
-  else if ((n()==4) && (m()==4))
-  {
-    double s0=0.,s1=0.,s2=0.,s3=0.;
-    s = src(0);
-    if(gl(1)<gl(0)) s1 = s*val[4];  if(gl(2)<gl(0)) s2 = s*val[8]; if(gl(3)<gl(0)) s3 = s*val[12];
-    s = src(1);
-    if(gl(0)<gl(1)) s0 += s*val[1]; if(gl(2)<gl(1)) s2 += s*val[9]; if(gl(3)<gl(1)) s3 += s*val[13];
-    s = src(2);
-    if(gl(0)<gl(2)) s0 += s*val[2]; if(gl(1)<gl(2)) s1 += s*val[6]; if(gl(3)<gl(2)) s3 += s*val[14];
-    s = src(3);
-    if(gl(0)<gl(3)) s0 += s*val[3]; if(gl(1)<gl(3)) s1 += s*val[7]; if(gl(2)<gl(3)) s2 += s*val[11];
-
-    dst(0) += s0;
-    dst(1) += s1;
-    dst(2) += s2;
-    dst(3) += s3;
-  }
-  else if ((n()==8) && (m()==8))
-  {
-    double s0=0.,s1=0.,s2=0.,s3=0.,s4=0.,s5=0.,s6=0.,s7=0.;
-    s = src(0);
-    if(gl(1)<gl(0)) s1 = s*val[8]; 
-    if(gl(2)<gl(0)) s2 = s*val[16]; 
-    if(gl(3)<gl(0)) s3 = s*val[24];
-    if(gl(4)<gl(0)) s4 = s*val[32]; 
-    if(gl(5)<gl(0)) s5 = s*val[40]; 
-    if(gl(6)<gl(0)) s6 = s*val[48]; 
-    if(gl(7)<gl(0)) s7 = s*val[56];
-    s = src(1);
-    if(gl(0)<gl(1)) s0 += s*val[1]; 
-    if(gl(2)<gl(1)) s2 += s*val[17]; 
-    if(gl(3)<gl(1)) s3 += s*val[25];
-    if(gl(4)<gl(1)) s4 += s*val[33]; 
-    if(gl(5)<gl(1)) s5 += s*val[41]; 
-    if(gl(6)<gl(1)) s6 += s*val[49]; 
-    if(gl(7)<gl(1)) s7 += s*val[57];
-    s = src(2);
-    if(gl(0)<gl(2)) s0 += s*val[2]; 
-    if(gl(1)<gl(2)) s1 += s*val[10]; 
-    if(gl(3)<gl(2)) s3 += s*val[26];
-    if(gl(4)<gl(2)) s4 += s*val[34]; 
-    if(gl(5)<gl(2)) s5 += s*val[42]; 
-    if(gl(6)<gl(2)) s6 += s*val[50]; 
-    if(gl(7)<gl(2)) s7 += s*val[58];
-    s = src(3);
-    if(gl(0)<gl(3)) s0 += s*val[3]; 
-    if(gl(1)<gl(3)) s1 += s*val[11]; 
-    if(gl(2)<gl(3)) s2 += s*val[19]; 
-    if(gl(4)<gl(3)) s4 += s*val[35]; 
-    if(gl(5)<gl(3)) s5 += s*val[43]; 
-    if(gl(6)<gl(3)) s6 += s*val[51]; 
-    if(gl(7)<gl(3)) s7 += s*val[59];
-    s = src(4);
-    if(gl(0)<gl(4)) s0 += s*val[4]; 
-    if(gl(1)<gl(4)) s1 += s*val[12]; 
-    if(gl(2)<gl(4)) s2 += s*val[20]; 
-    if(gl(3)<gl(4)) s3 += s*val[28];
-    if(gl(5)<gl(4)) s5 += s*val[44]; 
-    if(gl(6)<gl(4)) s6 += s*val[52]; 
-    if(gl(7)<gl(4)) s7 += s*val[60];
-    s = src(5);
-    if(gl(0)<gl(5)) s0 += s*val[5]; 
-    if(gl(1)<gl(5)) s1 += s*val[13]; 
-    if(gl(2)<gl(5)) s2 += s*val[21]; 
-    if(gl(3)<gl(5)) s3 += s*val[29];
-    if(gl(4)<gl(5)) s4 += s*val[37]; 
-    if(gl(6)<gl(5)) s6 += s*val[53]; 
-    if(gl(7)<gl(5)) s7 += s*val[61];
-    s = src(6);
-    if(gl(0)<gl(6)) s0 += s*val[6]; 
-    if(gl(1)<gl(6)) s1 += s*val[14]; 
-    if(gl(2)<gl(6)) s2 += s*val[22]; 
-    if(gl(3)<gl(6)) s3 += s*val[30];
-    if(gl(4)<gl(6)) s4 += s*val[38]; 
-    if(gl(5)<gl(6)) s5 += s*val[46]; 
-    if(gl(7)<gl(6)) s7 += s*val[62];
-    s = src(7);
-    if(gl(0)<gl(7)) s0 += s*val[7]; 
-    if(gl(1)<gl(7)) s1 += s*val[15]; 
-    if(gl(2)<gl(7)) s2 += s*val[23]; 
-    if(gl(3)<gl(7)) s3 += s*val[31];
-    if(gl(4)<gl(7)) s4 += s*val[39]; 
-    if(gl(5)<gl(7)) s5 += s*val[47]; 
-    if(gl(6)<gl(7)) s6 += s*val[55]; 
-    
-    dst(0) += s0;
-    dst(1) += s1;
-    dst(2) += s2;
-    dst(3) += s3;
-    dst(4) += s4;
-    dst(5) += s5;
-    dst(6) += s6;
-    dst(7) += s7;
-  }
-  else
-  {    
-    double* e = val;
-    const unsigned int size_m = m(),
-                      size_n = n();
-    for (unsigned int i=0; i<size_m; ++i)
-    {
-      s = 0.;
-      for (unsigned int j=0; j<size_n; ++j)
-       if(gl(i)<gl(j)) s += src(j) * *(e++);
-      dst(i) += s;
-    }
-  }
-}
-
-void dFMatrix::Tvmult (dVector& dst, const dVector& src, const bool adding) const
-{
-  Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
-  Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m()));
-
-  unsigned int i,j;
-  double s;
-  const unsigned int size_m = m(),
-                    size_n = n();
-  for (i=0; i<size_m; ++i)
-  {
-    s = 0.;
-    for (j=0; j<size_n; ++j)
-      s += src(j) * el(j,i);
-    if(!adding) dst(i) = s;
-    else dst(i) += s;
-  }
-}
-
-double dFMatrix::residual (dVector& dst, const dVector& src,
-                          const dVector& right) const
-{
-  Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
-  Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
-  Assert(right.size() == m(), ExcDimensionMismatch(right.size(), m()));
-
-  unsigned int i,j;
-  double s, res = 0.;
-  const unsigned int size_m = m(),
-                    size_n = n();
-  for (i=0; i<size_n; ++i)
-    {
-      s = right(i);
-      for (j=0; j<size_m; ++j)
-       s -= src(j) * el(i,j);
-      dst(i) = s;
-      res += s*s;
-    }
-  return sqrt(res);
-}
-
-void dFMatrix::forward (dVector& dst, const dVector& src) const
-{
-  Assert(n() == m(), ExcNotQuadratic());
-  Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
-  Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
-
-  unsigned int i,j;
-  unsigned int nu = (m()<n() ? m() : n());
-  double s;
-  for (i=0; i<nu; ++i)
-    {
-      s = src(i);
-      for (j=0; j<i; ++j) s -= dst(j) * el(i,j);
-      dst(i) = s/el(i,i);
-    }
-}
-
-void dFMatrix::backward (dVector& dst, const dVector& src) const
-{
-  Assert(n() == m(), ExcNotQuadratic());
-  Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
-  Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
-
-  unsigned int j;
-  unsigned int nu = (m()<n() ? m() : n());
-  double s;
-  for (int i=nu-1; i>=0; --i)
-    {
-      s = src(i);
-      for (j=i+1; j<nu; ++j) s -= dst(j) * el(i,j);
-      dst(i) = s/el(i,i);
-    }
-}
-
-dFMatrix& dFMatrix::operator = (const dFMatrix& m) {
-  reinit(m);
-
-  double *             p = &val[0];
-  const double *      vp = &m.val[0];
-  const double * const e = &val[dim_image*dim_range];
-
-  while (p!=e)
-    *p++ = *vp++;
-
-  return *this;
-}
-
-void dFMatrix::fill (const dFMatrix& src,
-                    const unsigned int i, const unsigned int j)
-{
-  Assert (n() >= src.n() + j, ExcInvalidDestination(n(), src.n(), j));
-  Assert (m() >= src.m() + i, ExcInvalidDestination(m(), src.m(), i));
-
-  for (unsigned int ii=0; ii<src.m() ; ++ii)
-    for (unsigned int jj=0; jj<src.n() ; ++jj)
-      el(ii+i,jj+j) = src.el(ii,jj);
-}
-
-void dFMatrix::add_row (const unsigned int i,
-                       const double s, const unsigned int j)
-{
-  for (unsigned int k=0; k<m(); ++k)
-    el(i,k) += s*el(j,k);
-}
-
-void dFMatrix::add_row (const unsigned int i, const double s,
-                       const unsigned int j, const double t,
-                       const unsigned int k)
-{
-  const unsigned int size_m = m();
-  for (unsigned l=0; l<size_m; ++l)
-    el(i,l) += s*el(j,l) + t*el(k,l);
-}
-
-void dFMatrix::add_col (const unsigned int i, const double s,
-                       const unsigned int j)
-{
-  for (unsigned int k=0; k<n(); ++k)
-    el(k,i) += s*el(k,j);
-}
-
-void dFMatrix::add_col (const unsigned int i, const double s,
-                       const unsigned int j, const double t,
-                       const unsigned int k)
-{
-  for (unsigned int l=0; l<n(); ++l)
-    el(l,i) += s*el(l,j) + t*el(l,k);
-}
-
-void dFMatrix::swap_row (const unsigned int i, const unsigned int j)
-{
-  double s;
-  for (unsigned int k=0; k<m(); ++k)
-  {
-    s = el(i,k); el(i,k) = el(j,k); el(j,k) = s;
-  }
-}
-
-void dFMatrix::swap_col (const unsigned int i, const unsigned int j)
-{
-  double s;
-  for (unsigned int k=0; k<n(); ++k)
-  {
-    s = el(k,i); el(k,i) = el(k,j); el(k,j) = s;
-  }
-}
-
-void dFMatrix::diagadd (const double& src)
-{
-  Assert (m() == n(), ExcDimensionMismatch(m(),n()));
-  for (unsigned int i=0; i<n(); ++i)
-    el(i,i) += src;
-}
-
-void dFMatrix::mmult (dFMatrix& dst, const dFMatrix& src) const
-{
-  Assert (n() == src.m(), ExcDimensionMismatch(n(), src.m()));
-  unsigned int i,j,k;
-  double s = 1.;
-  dst.reinit(m(), src.n());
-
-  for (i=0;i<m();i++)
-    for (j=0; j<src.n(); ++j)
-      {
-       s = 0.;
-       for (k=0;k<n();k++) s+= el(i,k) * src.el(k,j);
-       dst.el(i,j) = s;
-      }
-}
-
-/*void dFMatrix::mmult (dFMatrix& dst, const dFMatrix& src) const
-{
-  Assert (m() == src.n(), ExcDimensionMismatch(m(), src.n()));
-
-  unsigned int i,j,k;
-  double s = 1.;
-
-  dst.reinit(n(), src.m());
-
-  for (i=0;i<n();i++)
-    for (j=0;j<src.m();j++)
-      {
-       s = 0.;
-       for (k=0;k<m();k++) s+= el(i,k) * src.el(k,j);
-       dst.el(i,j) = s;
-      }
-}*/
-
-void dFMatrix::Tmmult (dFMatrix& dst, const dFMatrix& src) const
-{
-  Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n()));
-
-  unsigned int i,j,k;
-  double s = 1.;
-  dst.reinit(m(), src.m());
-
-  for (i=0;i<m();i++)
-    for (j=0;j<src.m();j++)
-      {
-       s = 0.;
-       for (k=0;k<n();k++) s+= el(k,i) * src.el(k,j);
-       dst.el(i,j) = s;
-      }
-}
-
-/*void dFMatrix::Tmmult(dFMatrix& dst, const dFMatrix& src) const
-{
-  Assert (m() == src.n(), ExcDimensionMismatch(m(), src.n()));
-
-  unsigned int i,j,k;
-  double s = 1.;
-  
-  dst.reinit(n(), src.m());
-
-  for (i=0;i<n();i++)
-    for (j=0;j<src.m();j++)
-      {
-       s = 0.;
-       for (k=0;k<m();k++) s+= el(k,i) * src.el(k,j);
-       dst.el(i,j) = s;
-      }
-}*/
-
-
-
-double dFMatrix::matrix_norm (const dVector &v) const {
-  Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size()));
-  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
-
-  double sum = 0.;
-  const unsigned int n_rows = m();
-  const double *val_ptr = &val[0];
-  const double *v_ptr;
-  
-  for (unsigned int row=0; row<n_rows; ++row)
-    {
-      double s = 0.;
-      const double * const val_end_of_row = val_ptr+n_rows;
-      v_ptr = v.begin();
-      while (val_ptr != val_end_of_row)
-       s += *val_ptr++ * *v_ptr++;
-
-      sum += s* v(row);
-    };
-
-  return sum;
-};
-
-
-
-double dFMatrix::matrix_scalar_product (const dVector &u, const dVector &v) const {
-  Assert(m() == u.size(), ExcDimensionMismatch(m(),v.size()));
-  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
-
-  double sum = 0.;
-  const unsigned int n_rows = m();
-  const unsigned int n_cols = n();
-  const double *val_ptr = &val[0];
-  const double *v_ptr;
-  
-  for (unsigned int row=0; row<n_rows; ++row)
-    {
-      double s = 0.;
-      const double * const val_end_of_row = val_ptr+n_cols;
-      v_ptr = v.begin();
-      while (val_ptr != val_end_of_row)
-       s += *val_ptr++ * *v_ptr++;
-
-      sum += s* u(row);
-    };
-
-  return sum;
-};
-
-
-
-void dFMatrix::print (ostream& s, int w, int p) const
-{
-  unsigned int i,j;
-  for (i=0;i<m();i++)
-    {
-      for (j=0;j<n();j++) s << setw(w) << setprecision(p) << el(i,j);
-      s << endl;
-    }
-}
-
-void dFMatrix::add (const double s,const dFMatrix& src)
-{
-  Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m()));
-  Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n()));
-  if ((n()==3) && (m()==3))
-  {
-    val[0] += s * src.el(0);
-    val[1] += s * src.el(1);
-    val[2] += s * src.el(2);
-    val[3] += s * src.el(3);
-    val[4] += s * src.el(4);
-    val[5] += s * src.el(5);
-    val[6] += s * src.el(6);
-    val[7] += s * src.el(7);
-    val[8] += s * src.el(8);
-  }
-  else if ((n()==4) && (m()==4))
-  {
-    val[0] += s * src.el(0);
-    val[1] += s * src.el(1);
-    val[2] += s * src.el(2);
-    val[3] += s * src.el(3);
-    val[4] += s * src.el(4);
-    val[5] += s * src.el(5);
-    val[6] += s * src.el(6);
-    val[7] += s * src.el(7);
-    val[8] += s * src.el(8);
-    val[9] += s * src.el(9);
-    val[10] += s * src.el(10);
-    val[11] += s * src.el(11);
-    val[12] += s * src.el(12);
-    val[13] += s * src.el(13);
-    val[14] += s * src.el(14);
-    val[15] += s * src.el(15);
-  }
-  else if ((n()==8) && (m()==8))
-  {
-    val[0] += s * src.el(0);
-    val[1] += s * src.el(1);
-    val[2] += s * src.el(2);
-    val[3] += s * src.el(3);
-    val[4] += s * src.el(4);
-    val[5] += s * src.el(5);
-    val[6] += s * src.el(6);
-    val[7] += s * src.el(7);
-    val[8] += s * src.el(8);
-    val[9] += s * src.el(9);
-    val[10] += s * src.el(10);
-    val[11] += s * src.el(11);
-    val[12] += s * src.el(12);
-    val[13] += s * src.el(13);
-    val[14] += s * src.el(14);
-    val[15] += s * src.el(15);
-    val[16] += s * src.el(16);
-    val[17] += s * src.el(17);
-    val[18] += s * src.el(18);
-    val[19] += s * src.el(19);
-
-    val[20] += s * src.el(20);
-    val[21] += s * src.el(21);
-    val[22] += s * src.el(22);
-    val[23] += s * src.el(23);
-    val[24] += s * src.el(24);
-    val[25] += s * src.el(25);
-    val[26] += s * src.el(26);
-    val[27] += s * src.el(27);
-    val[28] += s * src.el(28);
-    val[29] += s * src.el(29);
-
-    val[30] += s * src.el(30);
-    val[31] += s * src.el(31);
-    val[32] += s * src.el(32);
-    val[33] += s * src.el(33);
-    val[34] += s * src.el(34);
-    val[35] += s * src.el(35);
-    val[36] += s * src.el(36);
-    val[37] += s * src.el(37);
-    val[38] += s * src.el(38);
-    val[39] += s * src.el(39);
-
-    val[40] += s * src.el(40);
-    val[41] += s * src.el(41);
-    val[42] += s * src.el(42);
-    val[43] += s * src.el(43);
-    val[44] += s * src.el(44);
-    val[45] += s * src.el(45);
-    val[46] += s * src.el(46);
-    val[47] += s * src.el(47);
-    val[48] += s * src.el(48);
-    val[49] += s * src.el(49);
-
-    val[50] += s * src.el(50);
-    val[51] += s * src.el(51);
-    val[52] += s * src.el(52);
-    val[53] += s * src.el(53);
-    val[54] += s * src.el(54);
-    val[55] += s * src.el(55);
-    val[56] += s * src.el(56);
-    val[57] += s * src.el(57);
-    val[58] += s * src.el(58);
-    val[59] += s * src.el(59);
-
-    val[60] += s * src.el(60);
-    val[61] += s * src.el(61);
-    val[62] += s * src.el(62);
-    val[63] += s * src.el(63);
-  }
-  else
-  {
-    const unsigned int size = n()*m();
-    for (unsigned int i=0; i<size; i++)
-      val[i] += s * src.el(i);
-  }
-}
-
-
-
-void dFMatrix::add_diag (const double s, const dFMatrix& src)
-{
-  Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m()));
-  Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n()));
-
-  if ((n()==3) && (m()==3))
-  {
-    val[0] += s * src.el(0);
-    val[0] += s * src.el(1);
-    val[0] += s * src.el(2);
-    val[3] += s * src.el(3);
-    val[3] += s * src.el(4);
-    val[3] += s * src.el(5);
-    val[6] += s * src.el(6);
-    val[6] += s * src.el(7);
-    val[6] += s * src.el(8);
-  }
-  else if ((n()==4) && (m()==4))
-  {
-    val[0] += s * src.el(0);
-    val[0] += s * src.el(1);
-    val[0] += s * src.el(2);
-    val[0] += s * src.el(3);
-    val[4] += s * src.el(4);
-    val[4] += s * src.el(5);
-    val[4] += s * src.el(6);
-    val[4] += s * src.el(7);
-    val[8] += s * src.el(8);
-    val[8] += s * src.el(9);
-    val[8] += s * src.el(10);
-    val[8] += s * src.el(11);
-    val[12] += s * src.el(12);
-    val[12] += s * src.el(13);
-    val[12] += s * src.el(14);
-    val[12] += s * src.el(15);
-  }
-  else if ((n()==8) && (m()==8))
-  {
-    val[0] += s * src.el(0);
-    val[0] += s * src.el(1);
-    val[0] += s * src.el(2);
-    val[0] += s * src.el(3);
-    val[0] += s * src.el(4);
-    val[0] += s * src.el(5);
-    val[0] += s * src.el(6);
-    val[0] += s * src.el(7);
-    val[8] += s * src.el(8);
-    val[8] += s * src.el(9);
-    val[8] += s * src.el(10);
-    val[8] += s * src.el(11);
-    val[8] += s * src.el(12);
-    val[8] += s * src.el(13);
-    val[8] += s * src.el(14);
-    val[8] += s * src.el(15);
-    val[16] += s * src.el(16);
-    val[16] += s * src.el(17);
-    val[16] += s * src.el(18);
-    val[16] += s * src.el(19);
-
-    val[16] += s * src.el(20);
-    val[16] += s * src.el(21);
-    val[16] += s * src.el(22);
-    val[16] += s * src.el(23);
-    val[24] += s * src.el(24);
-    val[24] += s * src.el(25);
-    val[24] += s * src.el(26);
-    val[24] += s * src.el(27);
-    val[24] += s * src.el(28);
-    val[24] += s * src.el(29);
-
-    val[24] += s * src.el(30);
-    val[24] += s * src.el(31);
-    val[32] += s * src.el(32);
-    val[32] += s * src.el(33);
-    val[32] += s * src.el(34);
-    val[32] += s * src.el(35);
-    val[32] += s * src.el(36);
-    val[32] += s * src.el(37);
-    val[32] += s * src.el(38);
-    val[32] += s * src.el(39);
-
-    val[40] += s * src.el(40);
-    val[40] += s * src.el(41);
-    val[40] += s * src.el(42);
-    val[40] += s * src.el(43);
-    val[40] += s * src.el(44);
-    val[40] += s * src.el(45);
-    val[40] += s * src.el(46);
-    val[40] += s * src.el(47);
-    val[48] += s * src.el(48);
-    val[48] += s * src.el(49);
-
-    val[48] += s * src.el(50);
-    val[48] += s * src.el(51);
-    val[48] += s * src.el(52);
-    val[48] += s * src.el(53);
-    val[48] += s * src.el(54);
-    val[48] += s * src.el(55);
-    val[56] += s * src.el(56);
-    val[56] += s * src.el(57);
-    val[56] += s * src.el(58);
-    val[56] += s * src.el(59);
-
-    val[56] += s * src.el(60);
-    val[56] += s * src.el(61);
-    val[56] += s * src.el(62);
-    val[56] += s * src.el(63);
-  }
-  else
-  {
-    const unsigned int size = n()*m();
-    for (unsigned int i=0; i<size; i++)
-      val[i] += s * src.el(i);
-  }
-}
-
-void dFMatrix::Tadd (const double s, const dFMatrix& src)
-{
-  Assert (m() == n(),     ExcNotQuadratic());
-  Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m()));
-  Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n()));
-
-  if ((n()==3) && (m()==3))
-  {
-    val[0] += s * src.el(0);
-    val[1] += s * src.el(3);
-    val[2] += s * src.el(6);
-
-    val[3] += s * src.el(1);
-    val[4] += s * src.el(4);
-    val[5] += s * src.el(7);
-
-    val[6] += s * src.el(2);
-    val[7] += s * src.el(5);
-    val[8] += s * src.el(8);
-  }
-  else if ((n()==4) && (m()==4))
-  {
-    val[0] += s * src.el(0);
-    val[1] += s * src.el(4);
-    val[2] += s * src.el(8);
-    val[3] += s * src.el(12);
-
-    val[4] += s * src.el(1);
-    val[5] += s * src.el(5);
-    val[6] += s * src.el(9);
-    val[7] += s * src.el(13);
-
-    val[8] += s * src.el(2);
-    val[9] += s * src.el(6);
-    val[10] += s * src.el(10);
-    val[11] += s * src.el(14);
-
-    val[12] += s * src.el(3);
-    val[13] += s * src.el(7);
-    val[14] += s * src.el(11);
-    val[15] += s * src.el(15);
-  }
-  else if ((n()==8) && (m()==8))
-  {
-    val[0] += s * src.el(0);
-    val[1] += s * src.el(8);
-    val[2] += s * src.el(16);
-    val[3] += s * src.el(24);
-    val[4] += s * src.el(32);
-    val[5] += s * src.el(40);
-    val[6] += s * src.el(48);
-    val[7] += s * src.el(56);
-
-    val[8] += s * src.el(1);
-    val[9] += s * src.el(9);
-    val[10] += s * src.el(17);
-    val[11] += s * src.el(25);
-    val[12] += s * src.el(33);
-    val[13] += s * src.el(41);
-    val[14] += s * src.el(49);
-    val[15] += s * src.el(57);
-
-    val[16] += s * src.el(2);
-    val[17] += s * src.el(10);
-    val[18] += s * src.el(18);
-    val[19] += s * src.el(26);
-    val[20] += s * src.el(34);
-    val[21] += s * src.el(42);
-    val[22] += s * src.el(50);
-    val[23] += s * src.el(58);
-
-    val[24] += s * src.el(3);
-    val[25] += s * src.el(11);
-    val[26] += s * src.el(19);
-    val[27] += s * src.el(27);
-    val[28] += s * src.el(35);
-    val[29] += s * src.el(43);
-    val[30] += s * src.el(51);
-    val[31] += s * src.el(59);
-
-    val[32] += s * src.el(4);
-    val[33] += s * src.el(12);
-    val[34] += s * src.el(20);
-    val[35] += s * src.el(28);
-    val[36] += s * src.el(36);
-    val[37] += s * src.el(44);
-    val[38] += s * src.el(52);
-    val[39] += s * src.el(60);
-
-    val[40] += s * src.el(5);
-    val[41] += s * src.el(13);
-    val[42] += s * src.el(21);
-    val[43] += s * src.el(29);
-    val[44] += s * src.el(37);
-    val[45] += s * src.el(45);
-    val[46] += s * src.el(53);
-    val[47] += s * src.el(61);
-
-    val[48] += s * src.el(6);
-    val[49] += s * src.el(14);
-    val[50] += s * src.el(22);
-    val[51] += s * src.el(30);
-    val[52] += s * src.el(38);
-    val[53] += s * src.el(46);
-    val[54] += s * src.el(54);
-    val[55] += s * src.el(62);
-
-    val[56] += s * src.el(7);
-    val[57] += s * src.el(15);
-    val[58] += s * src.el(23);
-    val[59] += s * src.el(31);
-    val[60] += s * src.el(39);
-    val[61] += s * src.el(47);
-    val[62] += s * src.el(55);
-    val[63] += s * src.el(63);
-  }
-  else
-    Assert (false, ExcInternalError());
-}
-
-
-bool
-dFMatrix::operator == (const dFMatrix &m) const
-{
-  bool q = (dim_range==m.dim_range) && (dim_image==m.dim_image);
-  if (!q) return false;
-
-  for (unsigned int i=0; i<dim_image; ++i)
-    for (unsigned int j=0; j<dim_range; ++j)
-      if (el(i,j) != m.el(i,j)) return false;
-  return true;
-};
-
-
-double dFMatrix::determinant () const {
-  Assert (dim_range == dim_image,
-         ExcDimensionMismatch(dim_range, dim_image));
-  Assert ((dim_range>=1) && (dim_range<=3), ExcNotImplemented(dim_range));
-  
-  switch (dim_range) 
-    {
-      case 1:
-            return el(0,0);
-      case 2:
-            return el(0,0)*el(1,1) - el(1,0)*el(0,1);
-      case 3:
-            return  (el(0,0)*el(1,1)*el(2,2)
-                     -el(0,0)*el(1,2)*el(2,1)
-                     -el(1,0)*el(0,1)*el(2,2)
-                     +el(1,0)*el(0,2)*el(2,1)
-                     +el(2,0)*el(0,1)*el(1,2)
-                     -el(2,0)*el(0,2)*el(1,1));
-      default:
-            return 0;
-    };
-};
-
-double dFMatrix::norm2 () const
-{
-  double s = 0.;
-  for (unsigned int i=0;i<dim_image*dim_range;++i)
-    s += val[i]*val[i];
-  return s;
-}
-
-
-void dFMatrix::clear () {
-  double       *val_ptr = &val[0];
-  const double *end_ptr = &val[n()*m()];
-  while (val_ptr != end_ptr)
-    *val_ptr++ = 0.;
-};
-
-
-
-void dFMatrix::invert (const dFMatrix &M) {
-  Assert (dim_range == dim_image, ExcNotQuadratic());
-  Assert ((dim_range>=1) && (dim_range<=4), ExcNotImplemented(dim_range));
-  Assert (dim_range == M.dim_range,
-          ExcDimensionMismatch(dim_range,M.dim_range));
-  Assert (dim_image == M.dim_image,
-         ExcDimensionMismatch(dim_image,M.dim_image));
-  
-  switch (dim_range) 
-    {
-      case 1:
-           val[0] = 1.0/M.val[0];
-           return;
-      case 2:
-                                            // this is Maple output,
-                                            // thus a bit unstructured
-      {
-           const double t4 = 1.0/(M.el(0,0)*M.el(1,1)-M.el(0,1)*M.el(1,0));
-           el(0,0) = M.el(1,1)*t4;
-           el(0,1) = -M.el(0,1)*t4;
-           el(1,0) = -M.el(1,0)*t4;
-           el(1,1) = M.el(0,0)*t4;
-           return;
-      };
-      
-      case 3:
-      {
-           const double t4 = M.el(0,0)*M.el(1,1),
-                        t6 = M.el(0,0)*M.el(1,2),
-                        t8 = M.el(0,1)*M.el(1,0),
-                       t00 = M.el(0,2)*M.el(1,0),
-                       t01 = M.el(0,1)*M.el(2,0),
-                       t04 = M.el(0,2)*M.el(2,0),
-                       t07 = 1.0/(t4*M.el(2,2)-t6*M.el(2,1)-t8*M.el(2,2)+
-                                  t00*M.el(2,1)+t01*M.el(1,2)-t04*M.el(1,1));
-           el(0,0) = (M.el(1,1)*M.el(2,2)-M.el(1,2)*M.el(2,1))*t07;
-           el(0,1) = -(M.el(0,1)*M.el(2,2)-M.el(0,2)*M.el(2,1))*t07;
-           el(0,2) = -(-M.el(0,1)*M.el(1,2)+M.el(0,2)*M.el(1,1))*t07;
-           el(1,0) = -(M.el(1,0)*M.el(2,2)-M.el(1,2)*M.el(2,0))*t07;
-           el(1,1) = (M.el(0,0)*M.el(2,2)-t04)*t07;
-           el(1,2) = -(t6-t00)*t07;
-           el(2,0) = -(-M.el(1,0)*M.el(2,1)+M.el(1,1)*M.el(2,0))*t07;
-           el(2,1) = -(M.el(0,0)*M.el(2,1)-t01)*t07;
-           el(2,2) = (t4-t8)*t07;
-           return;
-      };
-
-      case 4:
-      {
-                                        // with (linalg);
-                                        // a:=matrix(4,4);
-                                        // evalm(a);
-                                        // ai:=inverse(a);
-                                        // readlib(C);
-                                        // C(ai,optimized,filename=x4);
-
-       const double t14 = M.el(0,0)*M.el(1,1);
-       const double t15 = M.el(2,2)*M.el(3,3);
-       const double t17 = M.el(2,3)*M.el(3,2);
-       const double t19 = M.el(0,0)*M.el(2,1);
-       const double t20 = M.el(1,2)*M.el(3,3);
-       const double t22 = M.el(1,3)*M.el(3,2);
-       const double t24 = M.el(0,0)*M.el(3,1);
-       const double t25 = M.el(1,2)*M.el(2,3);
-       const double t27 = M.el(1,3)*M.el(2,2);
-       const double t29 = M.el(1,0)*M.el(0,1);
-       const double t32 = M.el(1,0)*M.el(2,1);
-       const double t33 = M.el(0,2)*M.el(3,3);
-       const double t35 = M.el(0,3)*M.el(3,2);
-       const double t37 = M.el(1,0)*M.el(3,1);
-       const double t38 = M.el(0,2)*M.el(2,3);
-       const double t40 = M.el(0,3)*M.el(2,2);
-       const double t42 = t14*t15-t14*t17-t19*t20+t19*t22+
-                          t24*t25-t24*t27-t29*t15+t29*t17+
-                          t32*t33-t32*t35-t37*t38+t37*t40;
-       const double t43 = M.el(2,0)*M.el(0,1);
-       const double t46 = M.el(2,0)*M.el(1,1);
-       const double t49 = M.el(2,0)*M.el(3,1);
-       const double t50 = M.el(0,2)*M.el(1,3);
-       const double t52 = M.el(0,3)*M.el(1,2);
-       const double t54 = M.el(3,0)*M.el(0,1);
-       const double t57 = M.el(3,0)*M.el(1,1);
-       const double t60 = M.el(3,0)*M.el(2,1);
-       const double t63 = t43*t20-t43*t22-t46*t33+t46*t35+
-                          t49*t50-t49*t52-t54*t25+t54*t27+
-                          t57*t38-t57*t40-t60*t50+t60*t52;
-       const double t65 = 1/(t42+t63);
-       const double t71 = M.el(0,2)*M.el(2,1);
-       const double t73 = M.el(0,3)*M.el(2,1);
-       const double t75 = M.el(0,2)*M.el(3,1);
-       const double t77 = M.el(0,3)*M.el(3,1);
-       const double t81 = M.el(0,1)*M.el(1,2);
-       const double t83 = M.el(0,1)*M.el(1,3);
-       const double t85 = M.el(0,2)*M.el(1,1);
-       const double t87 = M.el(0,3)*M.el(1,1);
-       const double t101 = M.el(1,0)*M.el(2,2);
-       const double t103 = M.el(1,0)*M.el(2,3);
-       const double t105 = M.el(2,0)*M.el(1,2);
-       const double t107 = M.el(2,0)*M.el(1,3);
-       const double t109 = M.el(3,0)*M.el(1,2);
-       const double t111 = M.el(3,0)*M.el(1,3);
-       const double t115 = M.el(0,0)*M.el(2,2);
-       const double t117 = M.el(0,0)*M.el(2,3);
-       const double t119 = M.el(2,0)*M.el(0,2);
-       const double t121 = M.el(2,0)*M.el(0,3);
-       const double t123 = M.el(3,0)*M.el(0,2);
-       const double t125 = M.el(3,0)*M.el(0,3);
-       const double t129 = M.el(0,0)*M.el(1,2);
-       const double t131 = M.el(0,0)*M.el(1,3);
-       const double t133 = M.el(1,0)*M.el(0,2);
-       const double t135 = M.el(1,0)*M.el(0,3);
-       el(0,0) = (M.el(1,1)*M.el(2,2)*M.el(3,3)-M.el(1,1)*M.el(2,3)*M.el(3,2)-
-                  M.el(2,1)*M.el(1,2)*M.el(3,3)+M.el(2,1)*M.el(1,3)*M.el(3,2)+
-                  M.el(3,1)*M.el(1,2)*M.el(2,3)-M.el(3,1)*M.el(1,3)*M.el(2,2))*t65;
-       el(0,1) = -(M.el(0,1)*M.el(2,2)*M.el(3,3)-M.el(0,1)*M.el(2,3)*M.el(3,2)-
-                   t71*M.el(3,3)+t73*M.el(3,2)+t75*M.el(2,3)-t77*M.el(2,2))*t65;
-       el(0,2) = (t81*M.el(3,3)-t83*M.el(3,2)-t85*M.el(3,3)+t87*M.el(3,2)+
-                  t75*M.el(1,3)-t77*M.el(1,2))*t65;
-       el(0,3) = -(t81*M.el(2,3)-t83*M.el(2,2)-t85*M.el(2,3)+t87*M.el(2,2)+
-                   t71*M.el(1,3)-t73*M.el(1,2))*t65;
-       el(1,0) = -(t101*M.el(3,3)-t103*M.el(3,2)-t105*M.el(3,3)+t107*M.el(3,2)+
-                   t109*M.el(2,3)-t111*M.el(2,2))*t65;
-       el(1,1) = (t115*M.el(3,3)-t117*M.el(3,2)-t119*M.el(3,3)+t121*M.el(3,2)+
-                  t123*M.el(2,3)-t125*M.el(2,2))*t65;
-       el(1,2) = -(t129*M.el(3,3)-t131*M.el(3,2)-t133*M.el(3,3)+t135*M.el(3,2)+
-                   t123*M.el(1,3)-t125*M.el(1,2))*t65;
-       el(1,3) = (t129*M.el(2,3)-t131*M.el(2,2)-t133*M.el(2,3)+t135*M.el(2,2)+
-                  t119*M.el(1,3)-t121*M.el(1,2))*t65;
-       el(2,0) = (t32*M.el(3,3)-t103*M.el(3,1)-t46*M.el(3,3)+t107*M.el(3,1)+
-                  t57*M.el(2,3)-t111*M.el(2,1))*t65;
-       el(2,1) = -(t19*M.el(3,3)-t117*M.el(3,1)-t43*M.el(3,3)+t121*M.el(3,1)+
-                   t54*M.el(2,3)-t125*M.el(2,1))*t65;
-       el(2,2) = (t14*M.el(3,3)-t131*M.el(3,1)-t29*M.el(3,3)+t135*M.el(3,1)+
-                  t54*M.el(1,3)-t125*M.el(1,1))*t65;
-       el(2,3) = -(t14*M.el(2,3)-t131*M.el(2,1)-t29*M.el(2,3)+t135*M.el(2,1)+
-                   t43*M.el(1,3)-t121*M.el(1,1))*t65;
-       el(3,0) = -(t32*M.el(3,2)-t101*M.el(3,1)-t46*M.el(3,2)+t105*M.el(3,1)+
-                   t57*M.el(2,2)-t109*M.el(2,1))*t65;
-       el(3,1) = (t19*M.el(3,2)-t115*M.el(3,1)-t43*M.el(3,2)+t119*M.el(3,1)+
-                  t54*M.el(2,2)-t123*M.el(2,1))*t65;
-       el(3,2) = -(t14*M.el(3,2)-t129*M.el(3,1)-t29*M.el(3,2)+t133*M.el(3,1)+
-                   t54*M.el(1,2)-t123*M.el(1,1))*t65;
-       el(3,3) = (t14*M.el(2,2)-t129*M.el(2,1)-t29*M.el(2,2)+t133*M.el(2,1)+
-                  t43*M.el(1,2)-t119*M.el(1,1))*t65;
-      }
-    };    
-};
-  
-
-
-void dFMatrix::print_formatted (ostream &out, const unsigned int precision) const {
-  out.precision (precision);
-  out.setf (ios::scientific, ios::floatfield);   // set output format
-  
-  for (unsigned int i=0; i<m(); ++i) 
-    {
-      for (unsigned int j=0; j<n(); ++j)
-       if (el(i,j) != 0)
-         out << setw(precision+7)
-             << el(i,j) << ' ';
-       else
-         out << setw(precision+8) << " ";
-      out << endl;
-    };
-
-  AssertThrow (out, ExcIO());
-  
-  out.setf (0, ios::floatfield);                 // reset output format
-};
-
-
-// Gauss-Jordan-Algorithmus
-// cf. Stoer I (4th Edition) p. 153
-
-void dFMatrix::gauss_jordan()
-{
-  Assert (dim_range == dim_image, ExcNotQuadratic());
-  iVector p(n());
-
-  unsigned int i,j,k,r;
-  double max, hr;
-
-  for (i=0; i<n(); ++i) p(i) = i;
-
-  for (j=0; j<n(); ++j)
-    {
-                                      // pivotsearch
-      max = fabs(el(j,j));
-      r = j;
-      for (i=j+1; i<n(); ++i)
-       {
-         if (fabs(el(i,j)) > max)
-           {
-             max = fabs(el(i,j));
-             r = i;
-           }
-       }
-      Assert(max>1.e-16, ExcNotRegular());
-                                      // rowinterchange
-      if (r>j)
-       {
-         for (k=0; k<n(); ++k)
-           {
-             hr = el(j,k) ; el(j,k) = el(r,k) ; el(r,k) = hr;
-           }
-         i = p(j) ; p(j) = p(r) ; p(r) = i;
-       }
-
-                                      // transformation
-      hr = 1./el(j,j);
-      el(j,j) = hr;
-      for (k=0; k<n(); ++k)
-       {
-         if (k==j) continue;
-         for (i=0; i<n(); ++i)
-           {
-             if (i==j) continue;
-             el(i,k) -= el(i,j)*el(j,k)*hr;
-           }
-       }
-      for (i=0; i<n(); ++i)
-       {
-         el(i,j) *= hr;
-         el(j,i) *= -hr;
-       }
-      el(j,j) = hr;
-    }
-                                  // columninterchange
-  dVector hv(n());
-  for (i=0; i<n(); ++i)
-    {
-      for (k=0; k<n(); ++k) hv(p(k)) = el(i,k);
-      for (k=0; k<n(); ++k) el(i,k) = hv(k);
-    }
-}
-
-// QR-transformation cf. Stoer 1 4.8.2 (p. 191)
-
-void dFMatrix::householder(dVector& src)
-{
-  // m > n, src.n() = m
-  Assert (dim_range <= dim_image, ExcDimensionMismatch(dim_range, dim_image));
-  Assert (src.size() == dim_range, ExcDimensionMismatch(src.size(), dim_range));
-
-  for (unsigned int j=0 ; j<n() ; ++j)
-  {
-    double sigma = 0;
-    unsigned int i;
-    for (i=j ; i<m() ; ++i) sigma += el(i,j)*el(i,j);
-    if (fabs(sigma) < 1.e-15) return;
-    double s = el(j,j);
-    s = (s<0) ? sqrt(sigma) : -sqrt(sigma);
-    double dj = s;
-
-    double beta = 1./(s*el(j,j)-sigma);
-    el(j,j) -= s;
-
-    for (unsigned int k=j+1 ; k<n() ; ++k)
-    {
-      double sum = 0.;
-      for (i=j ; i<m() ; ++i) sum += el(i,j)*el(i,k);
-      sum *= beta;
-
-      for (i=j ; i<m() ; ++i) el(i,k) += sum*el(i,j);
-    }
-
-    double sum = 0.;
-    for (i=j ; i<m() ; ++i) sum += el(i,j)*src(i);
-    sum *= beta;
-
-    for (i=j ; i<m() ; ++i) src(i) += sum*el(i,j);
-    el(j,j) = dj;
-  }
-}
-
-double dFMatrix::least_squares(dVector& dst, dVector& src)
-{
-  // m > n, m = src.n, n = dst.n
-
-  householder(src);
-  backward(dst, src);
-
-  double sum = 0.;
-  for (unsigned int i=n() ; i<m() ; ++i) sum += src(i) * src(i);
-  return sqrt(sum);
-}
diff --git a/deal.II/lac/source/dsmatrix.cc b/deal.II/lac/source/dsmatrix.cc
deleted file mode 100644 (file)
index 16ec70d..0000000
+++ /dev/null
@@ -1,783 +0,0 @@
-// $Id$
-
-// This file is part of the DEAL Library
-// DEAL is Copyright(1995) by
-// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
-
-#include <lac/dsmatrix.h>
-#include <lac/dvector.h>
-#include <lac/ivector.h>
-
-#include <iostream>
-#include <iomanip>
-#include <algorithm>
-
-
-
-dSMatrixStruct::dSMatrixStruct () :
-               max_dim(0),
-               max_vec_len(0),
-               rowstart(0),
-               colnums(0)
-{
-  reinit (0,0,0);
-};
-
-
-
-dSMatrixStruct::dSMatrixStruct (const unsigned int m, const unsigned int n,
-                               const unsigned int max_per_row) 
-               : max_dim(0),
-                 max_vec_len(0),
-                 rowstart(0),
-                 colnums(0)
-{
-  reinit (m,n,max_per_row);
-};
-
-
-
-dSMatrixStruct::dSMatrixStruct (const unsigned int n,
-                               const unsigned int max_per_row)
-               : max_dim(0),
-                 max_vec_len(0),
-                 rowstart(0),
-                 colnums(0)
-{
-  reinit (n,n,max_per_row);
-};
-
-
-
-dSMatrixStruct::~dSMatrixStruct ()
-{
-  if (rowstart != 0)  delete[] rowstart;
-  if (colnums != 0)   delete[] colnums;
-}
-
-
-
-
-void
-dSMatrixStruct::reinit (const unsigned int m, const unsigned int n,
-                       const unsigned int max_per_row)
-{
-  Assert ((max_per_row>0) || ((m==0) && (n==0)), ExcInvalidNumber(max_per_row));
-  rows = m;
-  cols = n;
-  vec_len = m * max_per_row;
-  max_row_len = max_per_row;
-
-                                  // delete empty matrices
-  if ((m==0) || (n==0))
-    {
-      if (rowstart)  delete[] rowstart;
-      if (colnums)   delete[] colnums;
-      rowstart = 0;
-      colnums = 0;
-      max_vec_len = vec_len = max_dim = rows = cols = 0;
-      compressed = false;
-      return;
-    };
-  
-  if (rows > max_dim)
-    {
-      if (rowstart) delete[] rowstart;
-      max_dim = rows;
-      rowstart = new unsigned int[max_dim+1];
-    };
-  
-  if (vec_len > max_vec_len)
-    {
-      if (colnums) delete[] colnums;
-      max_vec_len = vec_len;
-      colnums = new int[max_vec_len];
-    };
-  
-  for (unsigned int i=0; i<=rows; i++)
-    rowstart[i] = i * max_per_row;
-  fill_n (&colnums[0], vec_len, -1);
-
-  if (rows == cols)
-    for (unsigned int i=0;i<rows;i++)
-      colnums[rowstart[i]] = i;
-
-  compressed = false;
-}
-
-
-void
-dSMatrixStruct::compress ()
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
-  
-  if (compressed) return;
-  unsigned int next_free_entry = 0,
-               next_row_start = 0,
-                   row_length = 0;
-
-                                  // reserve temporary storage to
-                                  // store the entries of one wor
-  int *tmp_entries = new int[max_row_len];
-
-                                  // Traverse all rows
-  for (unsigned int line=0; line<rows; ++line)
-    {
-                                      // copy used entries, break if
-                                      // first unused entry is reached
-      row_length = 0;
-      for (unsigned int j=rowstart[line]; j<rowstart[line+1]; ++j,++row_length)
-       if (colnums[j] != -1)
-         tmp_entries[row_length] = colnums[j];
-       else
-         break;
-                                      // now #rowstart# is
-                                      // the number of entries in
-                                      // this line
-
-                                      // for square matrices, the
-                                      // first entry in each row
-                                      // is the diagonal one. In
-                                      // this case only sort the
-                                      // remaining entries, otherwise
-                                      // sort all
-      sort ((rows==cols) ? &tmp_entries[1] : &tmp_entries[0],
-           &tmp_entries[row_length]);
-
-                                      // Re-insert column numbers
-                                      // into the field
-      for (unsigned int j=0; j<row_length; ++j)
-       colnums[next_free_entry++] = tmp_entries[j];
-
-                                      // note new start of this and
-                                      // the next row
-      rowstart[line] = next_row_start;
-      next_row_start = next_free_entry;
-
-                                      // some internal checks
-      Assert ((rows!=cols) ||
-             (colnums[rowstart[line]] == static_cast<signed int>(line)),
-             ExcInternalError());
-                                      // assert that the first entry
-                                      // does not show up in
-                                      // the remaining ones and that
-                                      // the remaining ones are unique
-                                      // among themselves (this handles
-                                      // both cases, quadratic and
-                                      // rectangular matrices)
-      Assert (find (&colnums[rowstart[line]+1],
-                   &colnums[next_row_start],
-                   colnums[rowstart[line]]) ==
-             &colnums[next_row_start],
-             ExcInternalError());
-      Assert (adjacent_find(&colnums[rowstart[line]+1],
-                           &colnums[next_row_start]) ==
-             &colnums[next_row_start],
-             ExcInternalError());
-    };
-  
-  vec_len = rowstart[rows] = next_row_start;
-  compressed = true;
-
-  delete[] tmp_entries;
-};
-
-
-
-bool
-dSMatrixStruct::empty () const {
-                                  // let's try to be on the safe side of
-                                  // life by using multiple possibilities in
-                                  // the check for emptiness... (sorry for
-                                  // this kludge -- emptying matrices and
-                                  // freeing memory was not present in the
-                                  // original implementation and I don't
-                                  // know at how many places I missed
-                                  // something in adding it, so I try to
-                                  // be cautious. wb)
-  if ((rowstart==0) || (rows==0) || (cols==0))
-    {
-      Assert (rowstart==0, ExcInternalError());
-      Assert (rows==0, ExcInternalError());
-      Assert (cols==0, ExcInternalError());
-      Assert (colnums==0, ExcInternalError());
-      Assert (vec_len==0, ExcInternalError());
-      Assert (max_vec_len==0, ExcInternalError());
-      Assert (vec_len==0, ExcInternalError());
-
-      return true;
-    };
-  return false;
-};
-
-
-
-int
-dSMatrixStruct::operator () (const unsigned int i, const unsigned int j) const
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  Assert (i<rows, ExcInvalidIndex(i,rows));
-  Assert (j<cols, ExcInvalidIndex(j,cols));
-  Assert (compressed, ExcNotCompressed());
-
-                                  // check first entry separately, since
-                                  // for square matrices this is
-                                  // the diagonal entry (check only
-                                  // if a first entry exists)
-  if (rowstart[i] != rowstart[i+1]) 
-    {
-      if (static_cast<signed int>(j) == colnums[rowstart[i]])
-       return rowstart[i];
-    }
-  else
-                                    // no first entry exists for this
-                                    // line
-    return -1;
-
-                                  // all other entries are sorted, so
-                                  // we can use a binary seach algorithm
-  const int* p = lower_bound (&colnums[rowstart[i]+1],
-                             &colnums[rowstart[i+1]],
-                             static_cast<signed int>(j));
-  if ((*p == static_cast<signed int>(j)) &&
-      (p != &colnums[rowstart[i+1]]))
-    return (p - &colnums[0]);
-  else
-    return -1;
-}
-
-
-void
-dSMatrixStruct::add (const unsigned int i, const unsigned int j)
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  Assert (i<rows, ExcInvalidIndex(i,rows));
-  Assert (j<cols, ExcInvalidIndex(j,cols));
-  Assert (compressed==false, ExcMatrixIsCompressed());
-
-  for (unsigned int k=rowstart[i]; k<rowstart[i+1]; k++)
-    {
-                                      // entry already exists
-      if (colnums[k] == (signed int)j) return;
-                                      // empty entry found, put new
-                                      // entry here
-      if (colnums[k] == -1)
-       {
-         colnums[k] = j;
-         return;
-       };
-    };
-
-                                  // if we came thus far, something went
-                                  // wrong: there was not enough space
-                                  // in this line
-  Assert (false, ExcNotEnoughSpace(i, rowstart[i+1]-rowstart[i]));
-}
-
-
-
-void
-dSMatrixStruct::add_matrix (const unsigned int n, const int* rowcols)
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  for (unsigned int i=0; i<n; ++i)
-    for (unsigned int j=0; j<n; ++j)
-      add(rowcols[i], rowcols[j]);
-}
-
-
-
-void
-dSMatrixStruct::add_matrix (const unsigned int m, const unsigned int n,
-                           const int* rows, const int* cols)
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  for (unsigned i=0; i<m; ++i)
-    for (unsigned j=0; j<n; ++j)
-      add(rows[i], cols[j]);
-}
-
-
-
-void
-dSMatrixStruct::add_matrix (const iVector& rowcols)
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  unsigned int i,j;
-  for (i=0;i<rowcols.n();i++)
-    for (j=0;j<rowcols.n();j++)
-      add(rowcols(i), rowcols(j));
-}
-
-
-
-void
-dSMatrixStruct::add_matrix (const iVector& rows, const iVector& cols)
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  unsigned int i,j;
-  for (i=0;i<rows.n();i++)
-    for (j=0;j<cols.n();j++)
-      add(rows(i), cols(j));
-}
-
-
-
-void
-dSMatrixStruct::print_gnuplot (ostream &out) const
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  for (unsigned int i=0; i<rows; ++i)
-    for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
-      if (colnums[j]>=0)
-       out << i << " " << -colnums[j] << endl;
-
-  AssertThrow (out, ExcIO());
-}
-
-
-
-unsigned int
-dSMatrixStruct::bandwidth () const
-{
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  unsigned int b=0;
-  for (unsigned int i=0; i<rows; ++i)
-    for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
-      if (colnums[j]>=0) 
-       {
-         if (static_cast<unsigned int>(abs(static_cast<int>(i-colnums[j]))) > b)
-           b = abs(static_cast<int>(i-colnums[j]));
-       }
-      else
-                                        // leave if at the end of
-                                        // the entries of this line
-       break;
-  return b;
-};
-
-
-
-unsigned int
-dSMatrixStruct::n_nonzero_elements () const {
-  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
-  Assert (compressed, ExcNotCompressed());
-  return colnums[rows]-colnums[0];
-};
-
-
-
-
-
-/*-------------------------------------------------------------------------*/
-
-
-dSMatrix::dSMatrix () :
-               cols(0),
-               val(0),
-               max_len(0) {};
-
-
-
-dSMatrix::dSMatrix (const dSMatrixStruct &c)
-               : cols(&c), val(0), max_len(0)
-{
-  reinit();
-};
-
-
-
-dSMatrix::~dSMatrix ()
-{
-  delete[] val;
-};
-
-
-
-void
-dSMatrix::reinit ()
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (cols->compressed || cols->empty(), ExcNotCompressed());
-
-  if (cols->empty()) 
-    {
-      if (val) delete[] val;
-      val = 0;
-      max_len = 0;
-      return;
-    };
-        
-  if (max_len<cols->vec_len)
-    {
-      if (val) delete[] val;
-      val = new double[cols->vec_len];
-      max_len = cols->vec_len;
-    };
-
-  if (val)
-    fill_n (&val[0], cols->vec_len, 0);
-}
-
-
-
-void
-dSMatrix::reinit (const dSMatrixStruct &sparsity) {
-  cols = &sparsity;
-  reinit ();
-};
-
-
-
-void
-dSMatrix::clear () {
-  cols = 0;
-  if (val) delete[] val;
-  val = 0;
-  max_len = 0;
-};
-
-
-
-unsigned int
-dSMatrix::n_nonzero_elements () const {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  return cols->n_nonzero_elements ();
-};
-
-
-
-dSMatrix &
-dSMatrix::copy_from (const dSMatrix &matrix) {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
-
-  double             *val_ptr = &val[0];
-  const double    *matrix_ptr = &matrix.val[0];
-  const double *const end_ptr = &val[cols->vec_len];
-
-  while (val_ptr != end_ptr)
-    *val_ptr++ = *matrix_ptr++;
-  
-  return *this;
-};
-
-
-
-void
-dSMatrix::add_scaled (const double factor, const dSMatrix &matrix) {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
-
-  double             *val_ptr    = &val[0];
-  const double       *matrix_ptr = &matrix.val[0];
-  const double *const end_ptr    = &val[cols->vec_len];
-
-  while (val_ptr != end_ptr)
-    *val_ptr++ += factor * *matrix_ptr++;
-};
-
-
-
-void
-dSMatrix::vmult (dVector& dst, const dVector& src) const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
-  Assert(n() == src.size(), ExcDimensionsDontMatch(n(),src.size()));
-
-  const unsigned int n_rows = m();
-  const double *val_ptr = &val[cols->rowstart[0]];
-  const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
-  double       *dst_ptr = &dst(0);
-  for (unsigned int row=0; row<n_rows; ++row)
-    {
-      double s = 0.;
-      const double *const val_end_of_row = &val[cols->rowstart[row+1]];
-      while (val_ptr != val_end_of_row)
-       s += *val_ptr++ * src(*colnum_ptr++);
-      *dst_ptr++ = s;
-    };
-};
-
-
-void
-dSMatrix::Tvmult (dVector& dst, const dVector& src) const
-{
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert(n() == dst.size(), ExcDimensionsDontMatch(n(),dst.size()));
-  Assert(m() == src.size(), ExcDimensionsDontMatch(m(),src.size()));
-
-  dst.clear ();
-
-  for (unsigned int i=0;i<m();i++)
-    {
-      for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
-       {
-         int p = cols->colnums[j];
-         dst(p) += val[j] * src(i);
-       }
-    }
-}
-
-
-
-double
-dSMatrix::matrix_norm (const dVector& v) const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert(m() == v.size(), ExcDimensionsDontMatch(m(),v.size()));
-  Assert(n() == v.size(), ExcDimensionsDontMatch(n(),v.size()));
-
-  double sum = 0.;
-  const unsigned int n_rows = m();
-  const double *val_ptr = &val[cols->rowstart[0]];
-  const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
-  for (unsigned int row=0; row<n_rows; ++row)
-    {
-      double s = 0.;
-      const double *val_end_of_row = &val[cols->rowstart[row+1]];
-      while (val_ptr != val_end_of_row)
-       s += *val_ptr++ * v(*colnum_ptr++);
-
-      sum += s* v(row);
-    };
-
-  return sum;
-};
-
-
-
-double
-dSMatrix::residual (dVector& dst, const dVector& u, const dVector& b) const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
-  Assert(m() == b.size(), ExcDimensionsDontMatch(m(),b.size()));
-  Assert(n() == u.size(), ExcDimensionsDontMatch(n(),u.size()));
-
-  double s,norm=0.;   
-  
-  for (unsigned int i=0;i<m();i++)
-    {
-      s = b(i);
-      for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
-       {
-         int p = cols->colnums[j];
-         s -= val[j] * u(p);
-       }
-      dst(i) = s;
-      norm += dst(i)*dst(i);
-    }
-  return sqrt(norm);
-}
-
-void
-dSMatrix::precondition_Jacobi (dVector& dst, const dVector& src,
-                              const double om) const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
-
-  const unsigned int n = src.size();
-  double             *dst_ptr = dst.begin();
-  const double       *src_ptr = src.begin();
-  const unsigned int *rowstart_ptr = &cols->rowstart[0];
-  
-  for (unsigned int i=0; i<n; ++i, ++dst_ptr, ++src_ptr, ++rowstart_ptr)
-                                    // note that for square matrices,
-                                    // the diagonal entry is the first
-                                    // in each row, i.e. at index
-                                    // rowstart[i]
-    *dst_ptr = om * *src_ptr / val[*rowstart_ptr];
-};
-
-
-void
-dSMatrix::precondition_SSOR (dVector& dst, const dVector& src,
-                            const double om) const
-{
-                                  // to understand how this function works
-                                  // you may want to take a look at the CVS
-                                  // archives to see the original version
-                                  // which is much clearer...
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
-
-  const unsigned int  n            = src.size();
-  const unsigned int *rowstart_ptr = &cols->rowstart[0];
-  double             *dst_ptr      = &dst(0);
-  
-  for (unsigned int row=0; row<n; ++row, ++dst_ptr, ++rowstart_ptr)
-    {
-      *dst_ptr = src(row);
-                                      // find the first element in this line
-                                      // which is on the right of the diagonal.
-                                      // we need to precondition with the
-                                      // elements on the left only.
-                                      // note: the first entry in each
-                                      // line denotes the diagonal element,
-                                      // which we need not check.
-      const unsigned int first_right_of_diagonal_index
-       = (lower_bound (&cols->colnums[*rowstart_ptr+1],
-                       &cols->colnums[*(rowstart_ptr+1)],
-                       static_cast<signed int>(row)) -
-          &cols->colnums[0]);
-                                      
-      for (unsigned int j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
-       *dst_ptr -= om* val[j] * dst(cols->colnums[j]);
-      *dst_ptr /= val[*rowstart_ptr];
-    };
-  
-  rowstart_ptr = &cols->rowstart[0];
-  dst_ptr      = &dst(0);
-  for (unsigned int row=0; row<n; ++row, ++rowstart_ptr, ++dst_ptr)
-    *dst_ptr *= (2.-om)*val[*rowstart_ptr];
-
-  rowstart_ptr = &cols->rowstart[n-1];
-  dst_ptr      = &dst(n-1);
-  for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
-    {
-      const unsigned int first_right_of_diagonal_index
-       = (lower_bound (&cols->colnums[*rowstart_ptr+1],
-                       &cols->colnums[*(rowstart_ptr+1)],
-                       static_cast<signed int>(row)) -
-          &cols->colnums[0]);
-      for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j)
-       if (cols->colnums[j] > row)
-         *dst_ptr -= om* val[j] * dst(cols->colnums[j]);
-      
-      *dst_ptr /= val[*rowstart_ptr];
-    };
-}
-
-void
-dSMatrix::precondition_SOR (dVector& dst, const dVector& src,
-                           const double om) const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
-
-  dst = src;
-  SOR(dst,om);
-};
-
-
-void dSMatrix::precondition (dVector &dst, const dVector &src) const {
-  Assert (m() == n(), ExcMatrixNotSquare());
-  dst=src;
-};
-
-
-void
-dSMatrix::SOR (dVector& dst, const double om) const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
-  Assert (m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
-
-  for (unsigned int row=0; row<m(); ++row)
-    {
-      double s = dst(row);
-      for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
-       if ((unsigned int)cols->colnums[j] < row)
-         s -= val[j] * dst(cols->colnums[j]);
-
-      dst(row) = s * om / val[cols->rowstart[row]];
-    }
-}
-
-void
-dSMatrix::SSOR (dVector& dst, const double om) const
-{
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-
-  int p;
-  const unsigned int  n = dst.size();
-  unsigned int  j;
-  double s;
-  
-  for (unsigned int i=0; i<n; i++)
-    {
-      s = 0.;
-      for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
-       {
-         p = cols->colnums[j];
-         if (p>=0)
-           {
-             if (i>j) s += val[j] * dst(p);
-           }
-       }
-      dst(i) -= s * om;
-      dst(i) /= val[cols->rowstart[i]];
-    }
-
-  for (int i=n-1; i>=0; i--)  // this time, i is signed, but alsways positive!
-    {
-      s = 0.;
-      for (j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
-       {
-         p = cols->colnums[j];
-         if (p>=0)
-           {
-             if ((unsigned int)i<j) s += val[j] * dst(p);
-           }
-       }
-      dst(i) -= s * om / val[cols->rowstart[i]];
-    }
-}
-
-
-
-const dSMatrixStruct & dSMatrix::get_sparsity_pattern () const {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  return *cols;
-};
-
-
-
-void dSMatrix::print (ostream &out) const {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-
-  for (unsigned int i=0; i<cols->rows; ++i)
-    for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
-      out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << endl;
-
-  AssertThrow (out, ExcIO());
-};
-
-
-
-void dSMatrix::print_formatted (ostream &out, const unsigned int precision) const {
-  Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (val != 0, ExcMatrixNotInitialized());
-  out.precision (precision);
-  out.setf (ios::scientific, ios::floatfield);   // set output format
-  
-  for (unsigned int i=0; i<m(); ++i) 
-    {
-      for (unsigned int j=0; j<n(); ++j)
-       if ((*cols)(i,j) != -1)
-         out << setw(precision+7)
-             << val[cols->operator()(i,j)] << ' ';
-       else
-         out << setw(precision+8) << " ";
-      out << endl;
-    };
-  AssertThrow (out, ExcIO());
-
-  out.setf (0, ios::floatfield);                 // reset output format
-};
-
diff --git a/deal.II/lac/source/dvector.cc b/deal.II/lac/source/dvector.cc
deleted file mode 100644 (file)
index 481bf5f..0000000
+++ /dev/null
@@ -1,520 +0,0 @@
-// $Id$
-
-// This file is part of the DEAL Library
-// DEAL is Copyright(1995) by
-// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
-
-#include <lac/dvector.h>
-#include <cmath>
-#include <algorithm>
-
-
-static inline double sqr (const double x) {
-  return x*x;
-};
-
-
-
-dVector::dVector () :
-               dim(0),
-               maxdim(0),
-               val(0)
-{}
-
-
-dVector::dVector (const unsigned int n) :
-               dim(0),
-               maxdim(0),
-               val(0)
-{
-  reinit (n, false);
-}
-
-
-dVector::dVector (const dVector& v) :
-               dim(v.size()),
-               maxdim(v.size()),
-               val(0)
-{
-  if (dim)
-    {
-      val = new double[maxdim];
-      Assert (val != 0, ExcOutOfMemory());
-      copy (v.begin(), v.end(), begin());
-    }
-}
-
-
-
-void dVector::reinit (const unsigned int n, const bool fast) {
-  if (n==0) 
-    {
-      if (val) delete[] val;
-      val = 0;
-      maxdim = dim = 0;
-      return;
-    };
-  
-  if (n>maxdim)
-    {
-      if (val) delete[] val;
-      val = new double[n];
-      Assert (val != 0, ExcOutOfMemory());
-      maxdim = n;
-    };
-  dim = n;
-  if (fast == false)
-    clear ();
-}
-
-
-
-void dVector::reinit (const dVector& v, const bool fast) {
-  reinit (v.size(), fast);
-};
-
-
-
-
-dVector::~dVector ()
-{
-  if (val) delete[] val;
-}
-
-
-
-void dVector::clear () {
-  if (dim>0)
-    fill (begin(), end(), 0.);
-}
-
-
-
-bool dVector::all_zero () const {
-  Assert (dim!=0, ExcEmptyVector());
-  
-  const_iterator p = begin(),
-                e = end();
-  while (p!=e)
-    if (*p++ != 0.0)
-      return false;
-  return true;
-};
-
-
-
-double dVector::operator * (const dVector& v) const
-{
-  Assert (dim!=0, ExcEmptyVector());
-
-  if (&v == this)
-    return norm_sqr();
-  
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-  
-  double sum0 = 0,
-        sum1 = 0,
-        sum2 = 0,
-        sum3 = 0;
-
-                                  // use modern processors better by
-                                  // allowing pipelined commands to be
-                                  // executed in parallel
-  const_iterator ptr  = begin(),
-                vptr = v.begin(),
-                eptr = ptr + (dim/4)*4;
-  while (ptr!=eptr)
-    {
-      sum0 += (*ptr++ * *vptr++);
-      sum1 += (*ptr++ * *vptr++);
-      sum2 += (*ptr++ * *vptr++);
-      sum3 += (*ptr++ * *vptr++);
-    };
-                                  // add up remaining elements
-  while (ptr != end())
-    sum0 += *ptr++ * *vptr++;
-    
-  return sum0+sum1+sum2+sum3;
-}
-
-
-
-double dVector::norm_sqr () const
-{
-  Assert (dim!=0, ExcEmptyVector());
-
-  double sum0 = 0,
-        sum1 = 0,
-        sum2 = 0,
-        sum3 = 0;
-
-                                  // use modern processors better by
-                                  // allowing pipelined commands to be
-                                  // executed in parallel
-  const_iterator ptr  = begin(),
-                eptr = ptr + (dim/4)*4;
-  while (ptr!=eptr)
-    {
-      sum0 += sqr(*ptr++);
-      sum1 += sqr(*ptr++);
-      sum2 += sqr(*ptr++);
-      sum3 += sqr(*ptr++);
-    };
-                                  // add up remaining elements
-  while (ptr != end())
-    sum0 += sqr(*ptr++);
-  
-  return sum0+sum1+sum2+sum3;
-};
-
-
-
-double dVector::mean_value () const
-{
-  Assert (dim!=0, ExcEmptyVector());
-
-  double sum0 = 0,
-        sum1 = 0,
-        sum2 = 0,
-        sum3 = 0;
-
-                                  // use modern processors better by
-                                  // allowing pipelined commands to be
-                                  // executed in parallel
-  const_iterator ptr  = begin(),
-                eptr = ptr + (dim/4)*4;
-  while (ptr!=eptr)
-    {
-      sum0 += *ptr++;
-      sum1 += *ptr++;
-      sum2 += *ptr++;
-      sum3 += *ptr++;
-    };
-                                  // add up remaining elements
-  while (ptr != end())
-    sum0 += *ptr++;
-  
-  return (sum0+sum1+sum2+sum3)/size();
-};
-
-
-
-double dVector::l1_norm () const
-{
-  Assert (dim!=0, ExcEmptyVector());
-
-  double sum0 = 0,
-        sum1 = 0,
-        sum2 = 0,
-        sum3 = 0;
-
-                                  // use modern processors better by
-                                  // allowing pipelined commands to be
-                                  // executed in parallel
-  const_iterator ptr  = begin(),
-                eptr = ptr + (dim/4)*4;
-  while (ptr!=eptr)
-    {
-      sum0 += fabs(*ptr++);
-      sum1 += fabs(*ptr++);
-      sum2 += fabs(*ptr++);
-      sum3 += fabs(*ptr++);
-    };
-                                  // add up remaining elements
-  while (ptr != end())
-    sum0 += fabs(*ptr++);
-  
-  return sum0+sum1+sum2+sum3;
-};
-
-
-
-double dVector::l2_norm () const
-{
-  return sqrt(norm_sqr());
-};
-
-
-
-double dVector::linfty_norm () const {
-  Assert (dim!=0, ExcEmptyVector());
-
-  double max0=0.,
-        max1=0.,
-        max2=0.,
-        max3=0.;
-  for (unsigned int i=0; i<(dim/4); ++i) 
-    {
-      if (max0<fabs(val[4*i]))   max0=fabs(val[4*i]);
-      if (max1<fabs(val[4*i+1])) max1=fabs(val[4*i+1]);
-      if (max2<fabs(val[4*i+2])) max2=fabs(val[4*i+2]);
-      if (max3<fabs(val[4*i+3])) max3=fabs(val[4*i+3]);
-    };
-                                  // add up remaining elements
-  for (unsigned int i=(dim/4)*4; i<dim; ++i)
-    if (max0<fabs(val[i]))
-      max0 = fabs(val[i]);
-
-  return max (max(max0, max1),
-             max(max2, max3));
-};
-  
-
-
-
-
-dVector& dVector::operator += (const dVector& v)
-{
-  Assert (dim!=0, ExcEmptyVector());
-
-  add (v);
-  return *this;
-}
-
-
-
-dVector& dVector::operator -= (const dVector& v)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin();
-  while (i_ptr!=i_end)
-    *i_ptr++ -= *v_ptr++;
-
-  return *this;
-}
-
-
-
-void dVector::add (const double v)
-{
-  Assert (dim!=0, ExcEmptyVector());
-
-  iterator i_ptr = begin(),
-          i_end = end();
-  while (i_ptr!=i_end)
-    *i_ptr++ += v;
-}
-
-
-
-void dVector::add (const dVector& v)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin();
-  while (i_ptr!=i_end)
-    *i_ptr++ += *v_ptr++;
-}
-
-
-
-void dVector::add (const double a, const dVector& v)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin();
-  while (i_ptr!=i_end)
-    *i_ptr++ += a * *v_ptr++;
-}
-
-
-
-void dVector::add (const double a, const dVector& v,
-                  const double b, const dVector& w)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-  Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin(),
-                w_ptr = w.begin();
-  while (i_ptr!=i_end)
-    *i_ptr++ += a * *v_ptr++ + b * *w_ptr++;
-}
-
-
-
-void dVector::sadd (const double x, const dVector& v)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin();
-  for (; i_ptr!=i_end; ++i_ptr)
-    *i_ptr = x * *i_ptr  + *v_ptr++;
-}
-
-
-
-void dVector::sadd (const double x, const double a, const dVector& v)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin();
-  for (; i_ptr!=i_end; ++i_ptr)
-    *i_ptr = x * *i_ptr  +  a * *v_ptr++;
-}
-
-
-
-void dVector::sadd (const double x, const double a,
-                   const dVector& v, const double b, const dVector& w)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-  Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin(),
-                w_ptr = w.begin();
-  for (; i_ptr!=i_end; ++i_ptr)
-    *i_ptr = x * *i_ptr  +  a * *v_ptr++  + b * *w_ptr++;
-}
-
-
-
-void dVector::sadd (const double x, const double a,
-                   const dVector& v, const double b,
-                   const dVector& w, const double c, const dVector& y)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-  Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
-  Assert (dim == y.dim, ExcDimensionsDontMatch(dim, y.dim));
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator v_ptr = v.begin(),
-                w_ptr = w.begin(),
-                y_ptr = y.begin();
-  
-  for (; i_ptr!=i_end; ++i_ptr)
-    *i_ptr = (x * *i_ptr)  +  (a * *v_ptr++)  +  (b * *w_ptr++)  + (c * *y_ptr++);
-}
-
-
-
-void dVector::scale (const double factor)
-{
-  Assert (dim!=0, ExcEmptyVector());
-
-  iterator ptr=begin(), eptr=end();
-  while (ptr!=eptr)
-    *ptr++ *= factor;
-}
-
-
-
-void dVector::equ (const double a, const dVector& u,
-                  const double b, const dVector& v)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == u.dim, ExcDimensionsDontMatch(dim, u.dim));
-  Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator u_ptr = u.begin(),
-                v_ptr = v.begin();
-  while (i_ptr!=i_end)
-    *i_ptr++ = a * *u_ptr++  + b * *v_ptr++;
-}
-
-
-
-void dVector::equ (const double a, const dVector& u)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (dim == u.dim, ExcDimensionsDontMatch(dim, u.dim));
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator u_ptr = u.begin();
-  while (i_ptr!=i_end)
-    *i_ptr++ = a * *u_ptr++;
-}
-
-
-
-void dVector::ratio (const dVector &a, const dVector &b) {
-  Assert (dim!=0, ExcEmptyVector());
-  Assert (a.dim == b.dim, ExcDimensionsDontMatch (a.dim, b.dim));
-
-                                  // no need to reinit with zeros, since
-                                  // we overwrite them anyway
-  reinit (a.size(), true);
-  iterator i_ptr = begin(),
-          i_end = end();
-  const_iterator a_ptr = a.begin(),
-                b_ptr = b.begin();
-  while (i_ptr!=i_end)
-    *i_ptr++ = *a_ptr++ / *b_ptr++;
-};
-
-
-
-dVector& dVector::operator = (const double s)
-{
-  Assert (dim!=0, ExcEmptyVector());
-  fill (begin(), end(), s);
-  return *this;
-}
-
-
-
-dVector& dVector::operator = (const dVector& v)
-{
-  if (v.dim != dim)
-    reinit (v.dim, true);
-
-  if (dim!=0)
-    copy (v.begin(), v.end(), begin());
-  return *this;
-}
-
-
-
-void dVector::print (FILE* f, const char* format) const
-{
-  Assert (dim!=0, ExcEmptyVector());
-  if (!format) format = " %5.2f";
-  for (unsigned int j=0;j<size();j++)
-    fprintf(f, format, val[j]);
-  fputc('\n',f);
-}
-
-
-
-void dVector::print (const char* format) const
-{
-  Assert (dim!=0, ExcEmptyVector());
-  if (!format) format = " %5.2f";
-  for (unsigned int j=0;j<size();j++)
-    printf (format, val[j]);
-  printf ("\n");
-}
-
-
-
-void dVector::print (ostream &out) const {
-  Assert (dim!=0, ExcEmptyVector());
-  for (unsigned int i=0; i<size(); ++i)
-    out << val[i] << endl;
-
-  AssertThrow (out, ExcIO());
-};
-
-
-
-

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.