]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add templates sparse matrices; they will eventually replace the old dSMatrix.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Feb 1999 14:36:30 +0000 (14:36 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 23 Feb 1999 14:36:30 +0000 (14:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@879 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_matrix.h [new file with mode: 0644]
deal.II/lac/include/lac/sparse_matrix.templates.h [new file with mode: 0644]
deal.II/lac/source/sparse_matrix.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h
new file mode 100644 (file)
index 0000000..10a72bc
--- /dev/null
@@ -0,0 +1,896 @@
+/*----------------------------   sparsematrix.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __sparsematrix_H
+#define __sparsematrix_H
+/*----------------------------   sparsematrix.h     ---------------------------*/
+
+
+// This file was once part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+// Revised, modified and extended by Wolfgang Bangerth
+
+
+#include <base/exceptions.h>
+
+
+//forward declarations
+template <typename number> class Vector;
+template <typename number> class SparseMatrix;
+
+class iVector;
+class ostream;
+
+
+
+/*
+
+   @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth
+   */
+class SparseMatrixStruct
+{
+  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.
+                                     */
+    SparseMatrixStruct (const SparseMatrixStruct &);
+    
+  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.
+                                     */
+    SparseMatrixStruct ();
+    
+                                    /**
+                                     * Initialize a rectangular matrix with
+                                     * #m# rows and #n# columns,
+                                     * with at most #max_per_row#
+                                     * nonzero entries per row.
+                                     */
+    SparseMatrixStruct (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.
+                                     */
+    SparseMatrixStruct (const unsigned int n,
+                       const unsigned int max_per_row);
+
+                                    /**
+                                     * Destructor.
+                                     */
+    ~SparseMatrixStruct ();
+    
+                                    /**
+                                     * 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.
+                                     *
+                                     * #SparseMatrix# objects require the
+                                     * #SparseMatrixStruct# 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
+                                     * #SparseMatrix#. 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;
+    
+    template <typename number> friend class SparseMatrix;
+};
+
+
+
+
+/*
+CLASS
+   SparseMatrix
+
+   @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998
+   */
+template <typename number>
+class SparseMatrix
+{
+  public:
+
+                                    /**
+                                     * 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(SparseMatrixStruct)#.
+                                     */
+    SparseMatrix ();
+    
+                                    /**
+                                     * 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.
+                                     */
+    SparseMatrix (const SparseMatrixStruct &sparsity);
+    
+                                    /**
+                                     * Destructor. Free all memory, but do not
+                                     * release the memory of the sparsity
+                                     * structure.
+                                     */
+    virtual ~SparseMatrix ();
+    
+
+                                    /**
+                                     * 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 SparseMatrixStruct &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 number 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 number 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 source matrix may be a matrix
+                                     * of arbitrary type, as long as its
+                                     * data type is convertible to the
+                                     * data type of this matrix.
+                                     *
+                                     * The function returns a reference to
+                                     * #this#.
+                                     */
+    template <typename somenumber>
+    SparseMatrix<number> & copy_from (const SparseMatrix<somenumber> &);
+
+                                    /**
+                                     * 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.
+                                     *
+                                     * The source matrix may be a matrix
+                                     * of arbitrary type, as long as its
+                                     * data type is convertible to the
+                                     * data type of this matrix.
+                                     */
+    template <typename somenumber>
+    void add_scaled (const number factor, const SparseMatrix<somenumber> &matrix);
+    
+                                    /**
+                                     * Return the value of the entry (i,j).
+                                     * This may be an expensive operation
+                                     * and you should always take care
+                                     * where to call this function.
+                                     * In order to avoid abuse, this function
+                                     * throws an exception if the wanted
+                                     * element does not exist in the matrix.
+                                     */
+    number 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.
+                                     */
+    number 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
+                                     * #SparseMatrixStruct# class for more details.
+                                     *
+                                     * You should use this interface very
+                                     * carefully and only if you are absolutely
+                                     * sure to know what you do. You should
+                                     * also note that the structure of these
+                                     * arrays may change over time.
+                                     * If you change the layout yourself, you
+                                     * should also rename this function to
+                                     * avoid programs relying on outdated
+                                     * information!
+                                     */
+    number global_entry (const unsigned int i) const;
+
+                                    /**
+                                     * Same as above, but with write access.
+                                     * You certainly know what you do?
+                                     */
+    number & global_entry (const unsigned int i);
+
+                                    /**
+                                     * Matrix-vector multiplication: let
+                                     * #dst = M*src# with #M# being this matrix.
+                                     */
+    template <typename somenumber>
+    void vmult (Vector<somenumber>& dst, const Vector<somenumber>& src) const;
+    
+                                    /**
+                                     * Matrix-vector multiplication: let
+                                     * #dst = M^T*src# with #M# being this
+                                     * matrix. This function does the same as
+                                     * #vmult# but takes the transposed matrix.
+                                     */
+    template <typename somenumber>
+    void Tvmult (Vector<somenumber>& dst, const Vector<somenumber>& 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.
+                                     */
+    template <typename somenumber>
+    double matrix_norm (const Vector<somenumber> &v) const;
+    
+                                    //
+    template <typename somenumber>
+    double residual (Vector<somenumber>& dst, const Vector<somenumber>& x,
+                    const Vector<somenumber>& b) const;
+                                    //
+    template <typename somenumber>
+    void precondition_Jacobi (Vector<somenumber>& dst, const Vector<somenumber>& src,
+                             const number om = 1.) const;
+                                    //
+    template <typename somenumber>
+    void precondition_SSOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
+                           const number om = 1.) const;
+                                    //
+    template <typename somenumber>
+    void precondition_SOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
+                          const number om = 1.) const;
+                                    //
+    template <typename somenumber>
+    void SSOR (Vector<somenumber>& dst, const number om = 1.) const;
+                                    //
+    template <typename somenumber>
+    void SOR (Vector<somenumber>& dst, const number om = 1.) const;
+                                    //
+    template <typename somenumber>
+    void precondition (Vector<somenumber>& dst, const Vector<somenumber>& 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 SparseMatrixStruct & 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 SparseMatrixStruct * cols;
+    number* val;
+    unsigned int max_len;
+
+                                    // make all other sparse matrices
+                                    // friends
+    template <typename somenumber> friend class SparseMatrix<somenumber>;
+};
+
+
+
+
+
+/*---------------------- Inline functions -----------------------------------*/
+
+
+inline
+unsigned int SparseMatrixStruct::n_rows () const {
+  return rows;
+};
+
+
+
+inline
+unsigned int SparseMatrixStruct::n_cols () const {
+  return cols;
+};
+
+
+
+inline
+bool SparseMatrixStruct::is_compressed () const {
+  return compressed;
+};
+
+
+
+inline
+const unsigned int * SparseMatrixStruct::get_rowstart_indices () const {
+  return rowstart;
+};
+
+
+
+inline
+const int * SparseMatrixStruct::get_column_numbers () const {
+  return colnums;
+};
+
+
+
+template <typename number>
+inline
+unsigned int SparseMatrix<number>::m () const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  return cols->rows;
+};
+
+
+
+template <typename number>
+inline
+unsigned int SparseMatrix<number>::n () const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  return cols->cols;
+};
+
+
+
+template <typename number>
+inline
+void SparseMatrix<number>::set (const unsigned int i, const unsigned int j,
+                               const number value) {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert ((cols->operator()(i,j) != -1) || (value == 0.),
+         ExcInvalidIndex(i,j));
+
+  const int index = cols->operator()(i,j);
+
+  if (index >= 0) val[index] = value;
+};
+
+
+
+template <typename number>
+inline
+void SparseMatrix<number>::add (const unsigned int i, const unsigned int j,
+                               const number value) {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert ((cols->operator()(i,j) != -1) || (value == 0.),
+         ExcInvalidIndex(i,j));
+
+  const int index = cols->operator()(i,j);
+  
+  if (index >= 0) val[index] += value;
+};
+
+
+
+
+
+template <typename number>
+inline
+number SparseMatrix<number>::operator () (const unsigned int i, const unsigned int j) const {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (cols->operator()(i,j) != -1,
+         ExcInvalidIndex(i,j));
+  return val[cols->operator()(i,j)];
+};
+
+
+
+template <typename number>
+inline
+number SparseMatrix<number>::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]];
+};
+
+
+
+template <typename number>
+inline
+number SparseMatrix<number>::global_entry (const unsigned int j) const {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  return val[j];
+};
+
+
+
+template <typename number>
+inline
+number & SparseMatrix<number>::global_entry (const unsigned int j) {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  return val[j];
+};
+
+
+
+/*----------------------------   sparsematrix.h     ---------------------------*/
+/* end of #ifndef __sparsematrix_H */
+#endif
+/*----------------------------   sparsematrix.h     ---------------------------*/
+
+
diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h
new file mode 100644 (file)
index 0000000..0e741ca
--- /dev/null
@@ -0,0 +1,472 @@
+// $Id$
+
+// This file was once part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+// Revised, modified and extended by Wolfgang Bangerth, 1998, 1999
+
+
+#include <lac/sparsematrix.h>
+#include <lac/vector.h>
+
+
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+
+
+
+
+template <typename number>
+SparseMatrix<number>::SparseMatrix () :
+               cols(0),
+               val(0),
+               max_len(0) {};
+
+
+
+template <typename number>
+SparseMatrix<number>::SparseMatrix (const SparseMatrixStruct &c)
+               : cols(&c), val(0), max_len(0)
+{
+  reinit();
+};
+
+
+
+template <typename number>
+SparseMatrix<number>::~SparseMatrix ()
+{
+  delete[] val;
+};
+
+
+
+template <typename number>
+void
+SparseMatrix<number>::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 number[cols->vec_len];
+      max_len = cols->vec_len;
+    };
+
+  if (val)
+    fill_n (&val[0], cols->vec_len, 0);
+}
+
+
+
+template <typename number>
+void
+SparseMatrix<number>::reinit (const SparseMatrixStruct &sparsity) {
+  cols = &sparsity;
+  reinit ();
+};
+
+
+
+template <typename number>
+void
+SparseMatrix<number>::clear () {
+  cols = 0;
+  if (val) delete[] val;
+  val = 0;
+  max_len = 0;
+};
+
+
+
+template <typename number>
+unsigned int
+SparseMatrix<number>::n_nonzero_elements () const {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  return cols->n_nonzero_elements ();
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+SparseMatrix<number> &
+SparseMatrix<number>::copy_from (const SparseMatrix<somenumber> &matrix) {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
+  Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
+
+  number                 *val_ptr = &val[0];
+  const somenumber    *matrix_ptr = &matrix.val[0];
+  const number     *const end_ptr = &val[cols->vec_len];
+
+  while (val_ptr != end_ptr)
+    *val_ptr++ = *matrix_ptr++;
+  
+  return *this;
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::add_scaled (const number factor,
+                                 const SparseMatrix<somenumber> &matrix) {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
+  Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
+
+  number             *val_ptr    = &val[0];
+  const somenumber   *matrix_ptr = &matrix.val[0];
+  const number *const end_ptr    = &val[cols->vec_len];
+
+  while (val_ptr != end_ptr)
+    *val_ptr++ += factor * *matrix_ptr++;
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::vmult (Vector<somenumber>& dst, const Vector<somenumber>& 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 number *val_ptr = &val[cols->rowstart[0]];
+  const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
+  somenumber   *dst_ptr = &dst(0);
+  for (unsigned int row=0; row<n_rows; ++row)
+    {
+      double s = 0.;
+      const number *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;
+    };
+};
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::Tvmult (Vector<somenumber>& dst, const Vector<somenumber>& 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);
+       }
+    }
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+double
+SparseMatrix<number>::matrix_norm (const Vector<somenumber>& 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 number *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 number *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;
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+double
+SparseMatrix<number>::residual (Vector<somenumber>& dst, const Vector<somenumber>& u, const Vector<somenumber>& 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);
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>& dst, const Vector<somenumber>& src,
+                                          const number om) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
+  Assert (m() == n(), ExcMatrixNotSquare());
+
+  const unsigned int n = src.size();
+  somenumber              *dst_ptr = dst.begin();
+  const somenumber        *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];
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
+                                        const number 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];
+  somenumber         *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];
+    };
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::precondition_SOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
+                           const number om) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
+  Assert (m() == n(), ExcMatrixNotSquare());
+
+  dst = src;
+  SOR(dst,om);
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void SparseMatrix<number>::precondition (Vector<somenumber> &dst, const Vector<somenumber> &src) const {
+  Assert (m() == n(), ExcMatrixNotSquare());
+  dst=src;
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::SOR (Vector<somenumber>& dst, const number 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)
+    {
+      somenumber 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]];
+    }
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::SSOR (Vector<somenumber>& dst, const number 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]];
+    }
+}
+
+
+
+template <typename number>
+const SparseMatrixStruct & SparseMatrix<number>::get_sparsity_pattern () const {
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  return *cols;
+};
+
+
+
+template <typename number>
+void SparseMatrix<number>::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());
+};
+
+
+
+template <typename number>
+void SparseMatrix<number>::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/sparse_matrix.cc b/deal.II/lac/source/sparse_matrix.cc
new file mode 100644 (file)
index 0000000..70d892a
--- /dev/null
@@ -0,0 +1,499 @@
+// $Id$
+
+// This file was once part of the DEAL Library
+// DEAL is Copyright(1995) by
+// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier
+// Revised, modified and extended by Wolfgang Bangerth, 1998, 1999
+
+
+#include <lac/sparsematrix.h>
+#include <lac/sparsematrix.templates.h>
+#include <lac/ivector.h>
+
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+
+
+
+SparseMatrixStruct::SparseMatrixStruct () :
+               max_dim(0),
+               max_vec_len(0),
+               rowstart(0),
+               colnums(0)
+{
+  reinit (0,0,0);
+};
+
+
+
+SparseMatrixStruct::SparseMatrixStruct (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);
+};
+
+
+
+SparseMatrixStruct::SparseMatrixStruct (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);
+};
+
+
+
+SparseMatrixStruct::~SparseMatrixStruct ()
+{
+  if (rowstart != 0)  delete[] rowstart;
+  if (colnums != 0)   delete[] colnums;
+}
+
+
+
+
+void
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::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
+SparseMatrixStruct::n_nonzero_elements () const {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
+  Assert (compressed, ExcNotCompressed());
+  return colnums[rows]-colnums[0];
+};
+
+
+
+
+
+
+
+
+// explicit instantiations for "float"
+template class SparseMatrix<float>;
+
+template SparseMatrix<float> & SparseMatrix<float>::copy_from (const SparseMatrix<float> &);
+template SparseMatrix<float> & SparseMatrix<float>::copy_from (const SparseMatrix<double> &);
+
+template void SparseMatrix<float>::add_scaled (const float, const SparseMatrix<float> &);
+template void SparseMatrix<float>::add_scaled (const float, const SparseMatrix<double> &);
+
+template void SparseMatrix<float>::vmult (Vector<float> &, const Vector<float> &) const;
+template void SparseMatrix<float>::vmult (Vector<double> &, const Vector<double> &) const;
+
+template void SparseMatrix<float>::Tvmult (Vector<float> &, const Vector<float> &) const;
+template void SparseMatrix<float>::Tvmult (Vector<double> &, const Vector<double> &) const;
+
+template double SparseMatrix<float>::matrix_norm (const Vector<float> &) const;
+template double SparseMatrix<float>::matrix_norm (const Vector<double> &) const;
+
+template double SparseMatrix<float>::residual (Vector<float> &,
+                                              const Vector<float> &,
+                                              const Vector<float> &) const;
+template double SparseMatrix<float>::residual (Vector<double> &,
+                                              const Vector<double> &,
+                                              const Vector<double> &) const;
+
+template void SparseMatrix<float>::precondition (Vector<float> &,
+                                                const Vector<float> &) const;
+template void SparseMatrix<float>::precondition (Vector<double> &,
+                                                const Vector<double> &) const;
+
+template void SparseMatrix<float>::precondition_SSOR (Vector<float> &,
+                                                     const Vector<float> &,
+                                                     float) const;
+template void SparseMatrix<float>::precondition_SSOR (Vector<double> &,
+                                                     const Vector<double> &,
+                                                     float) const;
+
+template void SparseMatrix<float>::precondition_SOR (Vector<float> &,
+                                                    const Vector<float> &,
+                                                    float) const;
+template void SparseMatrix<float>::precondition_SOR (Vector<double> &,
+                                                    const Vector<double> &,
+                                                    float) const;
+
+template void SparseMatrix<float>::precondition_Jacobi (Vector<float> &,
+                                                       const Vector<float> &,
+                                                       float) const;
+template void SparseMatrix<float>::precondition_Jacobi (Vector<double> &,
+                                                       const Vector<double> &,
+                                                       float) const;
+
+template void SparseMatrix<float>::SOR (Vector<float> &,
+                                       float) const;
+template void SparseMatrix<float>::SOR (Vector<double> &,
+                                       float) const;
+
+template void SparseMatrix<float>::SSOR (Vector<float> &,
+                                        float) const;
+template void SparseMatrix<float>::SSOR (Vector<double> &,
+                                        float) const;
+
+
+
+// explicit instantiations for "double"
+template class SparseMatrix<double>;
+
+template SparseMatrix<double> & SparseMatrix<double>::copy_from (const SparseMatrix<float> &);
+template SparseMatrix<double> & SparseMatrix<double>::copy_from (const SparseMatrix<double> &);
+
+template void SparseMatrix<double>::add_scaled (const double, const SparseMatrix<float> &);
+template void SparseMatrix<double>::add_scaled (const double, const SparseMatrix<double> &);
+
+template void SparseMatrix<double>::vmult (Vector<float> &, const Vector<float> &) const;
+template void SparseMatrix<double>::vmult (Vector<double> &, const Vector<double> &) const;
+
+template void SparseMatrix<double>::Tvmult (Vector<float> &, const Vector<float> &) const;
+template void SparseMatrix<double>::Tvmult (Vector<double> &, const Vector<double> &) const;
+
+template double SparseMatrix<double>::matrix_norm (const Vector<float> &) const;
+template double SparseMatrix<double>::matrix_norm (const Vector<double> &) const;
+
+template double SparseMatrix<double>::residual (Vector<float> &,
+                                               const Vector<float> &,
+                                               const Vector<float> &) const;
+template double SparseMatrix<double>::residual (Vector<double> &,
+                                               const Vector<double> &,
+                                               const Vector<double> &) const;
+
+template void SparseMatrix<double>::precondition (Vector<float> &,
+                                                 const Vector<float> &) const;
+template void SparseMatrix<double>::precondition (Vector<double> &,
+                                                 const Vector<double> &) const;
+
+template void SparseMatrix<double>::precondition_SSOR (Vector<float> &,
+                                                      const Vector<float> &,
+                                                      double) const;
+template void SparseMatrix<double>::precondition_SSOR (Vector<double> &,
+                                                      const Vector<double> &,
+                                                      double) const;
+
+template void SparseMatrix<double>::precondition_SOR (Vector<float> &,
+                                                     const Vector<float> &,
+                                                     double) const;
+template void SparseMatrix<double>::precondition_SOR (Vector<double> &,
+                                                     const Vector<double> &,
+                                                     double) const;
+
+template void SparseMatrix<double>::precondition_Jacobi (Vector<float> &,
+                                                        const Vector<float> &,
+                                                        double) const;
+template void SparseMatrix<double>::precondition_Jacobi (Vector<double> &,
+                                                        const Vector<double> &,
+                                                        double) const;
+
+template void SparseMatrix<double>::SOR (Vector<float> &,
+                                        double) const;
+template void SparseMatrix<double>::SOR (Vector<double> &,
+                                        double) const;
+
+template void SparseMatrix<double>::SSOR (Vector<float> &,
+                                         double) const;
+template void SparseMatrix<double>::SSOR (Vector<double> &,
+                                         double) const;

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.