--- /dev/null
+/*---------------------------- 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 ---------------------------*/
+
+
--- /dev/null
+// $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
+};
+
--- /dev/null
+// $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;