From 632ab0afbb334419ece6a26e40051d0db215ab5d Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 23 Feb 1999 14:36:30 +0000 Subject: [PATCH] Add templates sparse matrices; they will eventually replace the old dSMatrix. git-svn-id: https://svn.dealii.org/trunk@879 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_matrix.h | 896 ++++++++++++++++++ .../lac/include/lac/sparse_matrix.templates.h | 472 +++++++++ deal.II/lac/source/sparse_matrix.cc | 499 ++++++++++ 3 files changed, 1867 insertions(+) create mode 100644 deal.II/lac/include/lac/sparse_matrix.h create mode 100644 deal.II/lac/include/lac/sparse_matrix.templates.h create mode 100644 deal.II/lac/source/sparse_matrix.cc diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h new file mode 100644 index 0000000000..10a72bc5e1 --- /dev/null +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -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 + + +//forward declarations +template class Vector; +template 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 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 +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 + SparseMatrix & copy_from (const SparseMatrix &); + + /** + * 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 + void add_scaled (const number factor, const SparseMatrix &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 + void vmult (Vector& dst, const Vector& 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 + void Tvmult (Vector& dst, const Vector& src) const; + + + /** + * Return the norm of the vector #v# with + * respect to the norm induced by this + * matrix, i.e. $\left$. 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 + double matrix_norm (const Vector &v) const; + + // + template + double residual (Vector& dst, const Vector& x, + const Vector& b) const; + // + template + void precondition_Jacobi (Vector& dst, const Vector& src, + const number om = 1.) const; + // + template + void precondition_SSOR (Vector& dst, const Vector& src, + const number om = 1.) const; + // + template + void precondition_SOR (Vector& dst, const Vector& src, + const number om = 1.) const; + // + template + void SSOR (Vector& dst, const number om = 1.) const; + // + template + void SOR (Vector& dst, const number om = 1.) const; + // + template + void precondition (Vector& dst, const Vector& 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 friend class SparseMatrix; +}; + + + + + +/*---------------------- 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 +inline +unsigned int SparseMatrix::m () const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + return cols->rows; +}; + + + +template +inline +unsigned int SparseMatrix::n () const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + return cols->cols; +}; + + + +template +inline +void SparseMatrix::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 +inline +void SparseMatrix::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 +inline +number SparseMatrix::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 +inline +number SparseMatrix::diag_element (const unsigned int i) const { + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (irowstart[i]]; +}; + + + +template +inline +number SparseMatrix::global_entry (const unsigned int j) const { + Assert (cols != 0, ExcMatrixNotInitialized()); + return val[j]; +}; + + + +template +inline +number & SparseMatrix::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 index 0000000000..0e741caa99 --- /dev/null +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -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 +#include + + +#include +#include +#include + + + + +template +SparseMatrix::SparseMatrix () : + cols(0), + val(0), + max_len(0) {}; + + + +template +SparseMatrix::SparseMatrix (const SparseMatrixStruct &c) + : cols(&c), val(0), max_len(0) +{ + reinit(); +}; + + + +template +SparseMatrix::~SparseMatrix () +{ + delete[] val; +}; + + + +template +void +SparseMatrix::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_lenvec_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 +void +SparseMatrix::reinit (const SparseMatrixStruct &sparsity) { + cols = &sparsity; + reinit (); +}; + + + +template +void +SparseMatrix::clear () { + cols = 0; + if (val) delete[] val; + val = 0; + max_len = 0; +}; + + + +template +unsigned int +SparseMatrix::n_nonzero_elements () const { + Assert (cols != 0, ExcMatrixNotInitialized()); + return cols->n_nonzero_elements (); +}; + + + +template +template +SparseMatrix & +SparseMatrix::copy_from (const SparseMatrix &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 +template +void +SparseMatrix::add_scaled (const number factor, + const SparseMatrix &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 +template +void +SparseMatrix::vmult (Vector& dst, const Vector& 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; rowrowstart[row+1]]; + while (val_ptr != val_end_of_row) + s += *val_ptr++ * src(*colnum_ptr++); + *dst_ptr++ = s; + }; +}; + + +template +template +void +SparseMatrix::Tvmult (Vector& dst, const Vector& 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;irowstart[i]; jrowstart[i+1] ;j++) + { + int p = cols->colnums[j]; + dst(p) += val[j] * src(i); + } + } +} + + + +template +template +double +SparseMatrix::matrix_norm (const Vector& 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; rowrowstart[row+1]]; + while (val_ptr != val_end_of_row) + s += *val_ptr++ * v(*colnum_ptr++); + + sum += s* v(row); + }; + + return sum; +}; + + + +template +template +double +SparseMatrix::residual (Vector& dst, const Vector& u, const Vector& 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;irowstart[i]; jrowstart[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 +template +void +SparseMatrix::precondition_Jacobi (Vector& dst, const Vector& 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 +template +void +SparseMatrix::precondition_SSOR (Vector& dst, const Vector& 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; rowcolnums[*rowstart_ptr+1], + &cols->colnums[*(rowstart_ptr+1)], + static_cast(row)) - + &cols->colnums[0]); + + for (unsigned int j=(*rowstart_ptr)+1; jcolnums[j]); + *dst_ptr /= val[*rowstart_ptr]; + }; + + rowstart_ptr = &cols->rowstart[0]; + dst_ptr = &dst(0); + for (unsigned int row=0; rowrowstart[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(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 +template +void +SparseMatrix::precondition_SOR (Vector& dst, const Vector& src, + const number om) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (val != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + + dst = src; + SOR(dst,om); +}; + + + +template +template +void SparseMatrix::precondition (Vector &dst, const Vector &src) const { + Assert (m() == n(), ExcMatrixNotSquare()); + dst=src; +}; + + + +template +template +void +SparseMatrix::SOR (Vector& 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; rowrowstart[row]; jrowstart[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 +template +void +SparseMatrix::SSOR (Vector& 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; irowstart[i]; jrowstart[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]; jrowstart[i+1] ;j++) + { + p = cols->colnums[j]; + if (p>=0) + { + if ((unsigned int)irowstart[i]]; + } +} + + + +template +const SparseMatrixStruct & SparseMatrix::get_sparsity_pattern () const { + Assert (cols != 0, ExcMatrixNotInitialized()); + return *cols; +}; + + + +template +void SparseMatrix::print (ostream &out) const { + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (val != 0, ExcMatrixNotInitialized()); + + for (unsigned int i=0; irows; ++i) + for (unsigned int j=cols->rowstart[i]; jrowstart[i+1]; ++j) + out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << endl; + + AssertThrow (out, ExcIO()); +}; + + + +template +void SparseMatrix::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; ioperator()(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 index 0000000000..70d892a245 --- /dev/null +++ b/deal.II/lac/source/sparse_matrix.cc @@ -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 +#include +#include + +#include +#include +#include + + + +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(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(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(j)); + if ((*p == static_cast(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=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=0) + { + if (static_cast(abs(static_cast(i-colnums[j]))) > b) + b = abs(static_cast(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; + +template SparseMatrix & SparseMatrix::copy_from (const SparseMatrix &); +template SparseMatrix & SparseMatrix::copy_from (const SparseMatrix &); + +template void SparseMatrix::add_scaled (const float, const SparseMatrix &); +template void SparseMatrix::add_scaled (const float, const SparseMatrix &); + +template void SparseMatrix::vmult (Vector &, const Vector &) const; +template void SparseMatrix::vmult (Vector &, const Vector &) const; + +template void SparseMatrix::Tvmult (Vector &, const Vector &) const; +template void SparseMatrix::Tvmult (Vector &, const Vector &) const; + +template double SparseMatrix::matrix_norm (const Vector &) const; +template double SparseMatrix::matrix_norm (const Vector &) const; + +template double SparseMatrix::residual (Vector &, + const Vector &, + const Vector &) const; +template double SparseMatrix::residual (Vector &, + const Vector &, + const Vector &) const; + +template void SparseMatrix::precondition (Vector &, + const Vector &) const; +template void SparseMatrix::precondition (Vector &, + const Vector &) const; + +template void SparseMatrix::precondition_SSOR (Vector &, + const Vector &, + float) const; +template void SparseMatrix::precondition_SSOR (Vector &, + const Vector &, + float) const; + +template void SparseMatrix::precondition_SOR (Vector &, + const Vector &, + float) const; +template void SparseMatrix::precondition_SOR (Vector &, + const Vector &, + float) const; + +template void SparseMatrix::precondition_Jacobi (Vector &, + const Vector &, + float) const; +template void SparseMatrix::precondition_Jacobi (Vector &, + const Vector &, + float) const; + +template void SparseMatrix::SOR (Vector &, + float) const; +template void SparseMatrix::SOR (Vector &, + float) const; + +template void SparseMatrix::SSOR (Vector &, + float) const; +template void SparseMatrix::SSOR (Vector &, + float) const; + + + +// explicit instantiations for "double" +template class SparseMatrix; + +template SparseMatrix & SparseMatrix::copy_from (const SparseMatrix &); +template SparseMatrix & SparseMatrix::copy_from (const SparseMatrix &); + +template void SparseMatrix::add_scaled (const double, const SparseMatrix &); +template void SparseMatrix::add_scaled (const double, const SparseMatrix &); + +template void SparseMatrix::vmult (Vector &, const Vector &) const; +template void SparseMatrix::vmult (Vector &, const Vector &) const; + +template void SparseMatrix::Tvmult (Vector &, const Vector &) const; +template void SparseMatrix::Tvmult (Vector &, const Vector &) const; + +template double SparseMatrix::matrix_norm (const Vector &) const; +template double SparseMatrix::matrix_norm (const Vector &) const; + +template double SparseMatrix::residual (Vector &, + const Vector &, + const Vector &) const; +template double SparseMatrix::residual (Vector &, + const Vector &, + const Vector &) const; + +template void SparseMatrix::precondition (Vector &, + const Vector &) const; +template void SparseMatrix::precondition (Vector &, + const Vector &) const; + +template void SparseMatrix::precondition_SSOR (Vector &, + const Vector &, + double) const; +template void SparseMatrix::precondition_SSOR (Vector &, + const Vector &, + double) const; + +template void SparseMatrix::precondition_SOR (Vector &, + const Vector &, + double) const; +template void SparseMatrix::precondition_SOR (Vector &, + const Vector &, + double) const; + +template void SparseMatrix::precondition_Jacobi (Vector &, + const Vector &, + double) const; +template void SparseMatrix::precondition_Jacobi (Vector &, + const Vector &, + double) const; + +template void SparseMatrix::SOR (Vector &, + double) const; +template void SparseMatrix::SOR (Vector &, + double) const; + +template void SparseMatrix::SSOR (Vector &, + double) const; +template void SparseMatrix::SSOR (Vector &, + double) const; -- 2.39.5