From: wolf Date: Mon, 2 Dec 2002 17:32:38 +0000 (+0000) Subject: Should not have deleted this file. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ce46688feba5df24c69ab1dc93c0513f21cbbdf2;p=dealii-svn.git Should not have deleted this file. git-svn-id: https://svn.dealii.org/trunk@6794 0785d39b-7218-0410-832d-ea1e28bc413d --- 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..29c71926bc --- /dev/null +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -0,0 +1,1491 @@ +//---------------------------- sparse_matrix.h --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- sparse_matrix.h --------------------------- +#ifndef __deal2__sparse_matrix_h +#define __deal2__sparse_matrix_h + + +#include +#include +#include +#include +#include + +template class Vector; +template class FullMatrix; + +/** + * Sparse matrix. + * + * + * @sect2{On template instantiations} + * + * Member functions of this class are either implemented in this file + * or in a file of the same name with suffix ``.templates.h''. For the + * most common combinations of the template parameters, instantiations + * of this class are provided in a file with suffix ``.cc'' in the + * ``source'' directory. If you need an instantiation that is not + * listed there, you have to include this file along with the + * corresponding ``.templates.h'' file and instantiate the respective + * class yourself. + * + * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998 + */ +template +class SparseMatrix : public Subscriptor +{ + public: + /** + * Accessor class for iterators + */ + class Accessor + { + public: + /** + * Constructor. Since we use + * accessors only for read + * access, a const matrix + * pointer is sufficient. + */ + Accessor (const SparseMatrix*, + unsigned int row, + unsigned short index); + + /** + * Row number of the element + * represented by this + * object. + */ + unsigned int row() const; + + /** + * Index in row of the element + * represented by this + * object. + */ + unsigned short index() const; + + /** + * Column number of the + * element represented by + * this object. + */ + unsigned int column() const; + + /** + * Value of this matrix entry. + */ + number value() const; + + protected: + /** + * The matrix accessed. + */ + const SparseMatrix* matrix; + + /** + * Current row number. + */ + unsigned int a_row; + + /** + * Current index in row. + */ + unsigned short a_index; + }; + + /** + * STL conforming iterator. + */ + class const_iterator : private Accessor + { + public: + /** + * Constructor. + */ + const_iterator(const SparseMatrix*, + unsigned int row, + unsigned short index); + + /** + * Prefix increment. + */ + const_iterator& operator++ (); + + /** + * Postfix increment. + */ + const_iterator& operator++ (int); + + /** + * Dereferencing operator. + */ + const Accessor& operator* () const; + + /** + * Dereferencing operator. + */ + const Accessor* operator-> () const; + + /** + * Comparison. True, if + * both iterators point to + * the same matrix + * position. + */ + bool operator == (const const_iterator&) const; + /** + * Inverse of @p{==}. + */ + bool operator != (const const_iterator&) const; + + /** + * Comparison + * operator. Result is true + * if either the first row + * number is smaller or if + * the row numbers are + * equal and the first + * index is smaller. + */ + bool operator < (const const_iterator&) const; + }; + + /** + * Type of matrix entries. In analogy to + * the STL container classes. + */ + typedef number value_type; + + /** + * Constructor; initializes the matrix to + * be empty, without any structure, i.e. + * the matrix is not usable at all. This + * constructor is therefore only useful + * for matrices which are members of a + * class. All other matrices should be + * created at a point in the data flow + * where all necessary information is + * available. + * + * You have to initialize + * the matrix before usage with + * @p{reinit(SparsityPattern)}. + */ + SparseMatrix (); + + /** + * Copy constructor. This constructor is + * only allowed to be called if the matrix + * to be copied is empty. This is for the + * same reason as for the + * @p{SparsityPattern}, see there for the + * details. + * + * If you really want to copy a whole + * matrix, you can do so by using the + * @p{copy_from} function. + */ + SparseMatrix (const SparseMatrix &); + + /** + * Constructor. Takes the given + * matrix sparsity structure to + * represent the sparsity pattern + * of this matrix. You can change + * the sparsity pattern later on + * by calling the @p{reinit} + * function. + * + * You have to make sure that the + * lifetime of the sparsity + * structure is at least as long + * as that of this matrix or as + * long as @p{reinit} is not + * called with a new sparsity + * structure. + * + * The constructor is marked + * explicit so as to disallow + * that someone passes a sparsity + * pattern in place of a sparse + * matrix to some function, where + * an empty matrix would be + * generated then. + */ + explicit SparseMatrix (const SparsityPattern &sparsity); + + /** + * Destructor. Free all memory, but do not + * release the memory of the sparsity + * structure. + */ + virtual ~SparseMatrix (); + + /** + * Pseudo operator only copying + * empty objects. + */ + SparseMatrix& operator = (const SparseMatrix &); + + /** + * Reinitialize the object but + * keep to the sparsity pattern + * previously used. This may be + * necessary if you @p{reinit}'d + * the sparsity structure and + * want to update the size of the + * matrix. + * + * Note that memory is only + * reallocated if the new size + * exceeds the old size. If that + * is not the case, the allocated + * memory is not reduced. However, + * if the sparsity structure is + * empty (i.e. the dimensions are + * zero), then all memory is + * freed. + * + * If the sparsity pattern has + * not changed, then the effect + * of this function is simply to + * reset all matrix entries to + * zero. + */ + virtual void reinit (); + + /** + * Reinitialize the sparse matrix + * with the given sparsity + * pattern. The latter tells the + * matrix how many nonzero + * elements there need to be + * reserved. + * + * Regarding memory allocation, + * the same applies as said + * above. + * + * You have to make sure that the + * lifetime of the sparsity + * structure is at least as long + * as that of this matrix or as + * long as @p{reinit} is not called + * with a new sparsity structure. + * + * The elements of the matrix are + * set to zero by this function. + */ + virtual void reinit (const SparsityPattern &sparsity); + + /** + * Release all memory and return + * to a state just like after + * having called the default + * constructor. It also forgets + * the sparsity pattern it was + * previously tied to. + */ + virtual void clear (); + + /** + * Return whether the object is + * empty. It is empty if either + * both dimensions are zero or no + * @p{SparsityPattern} is + * associated. + */ + bool empty () const; + + /** + * Return the dimension of the + * image space. To remember: the + * matrix is of dimension + * $m \times n$. + */ + unsigned int m () const; + + /** + * Return the dimension of the + * range space. To remember: the + * matrix is of dimension + * $m \times n$. + */ + unsigned int n () const; + + /** + * Return the number of nonzero + * elements of this + * matrix. Actually, it returns + * the number of entries in the + * sparsity pattern; if any of + * the entries should happen to + * be zero, it is counted anyway. + */ + unsigned int n_nonzero_elements () const; + + /** + * Return the number of actually + * nonzero elements of this + * matrix. + * + * Note, that this function does + * (in contrary to the + * @p{n_nonzero_elements}) NOT + * count all entries of the + * sparsity pattern but only the + * ones that are nonzero. + */ + unsigned int n_actually_nonzero_elements () const; + + /** + * Set the element @p{(i,j)} to @p{value}. + * Throws an error if the entry does + * not exist. Still, it is allowed to store + * zero values in non-existent fields. + */ + void set (const unsigned int i, const unsigned int j, + const number value); + + /** + * Add @p{value} to the element + * @p{(i,j)}. Throws an error if + * the entry does not + * exist. Still, it is allowed to + * store zero values in + * non-existent fields. + */ + void add (const unsigned int i, const unsigned int j, + const number value); + + /** + * Symmetrize the matrix by + * forming the mean value between + * the existing matrix and its + * transpose, $A = \frac 12(A+A^T)$. + * + * This operation assumes that + * the underlying sparsity + * pattern represents a symmetric + * object. If this is not the + * case, then the result of this + * operation will not be a + * symmetric matrix, since it + * only explicitly symmetrizes + * by looping over the lower left + * triangular part for efficiency + * reasons; if there are entries + * in the upper right triangle, + * then these elements are missed + * in the + * symmetrization. Symmetrization + * of the sparsity pattern can be + * obtain by the + * @ref{SparsityPattern}@p{::symmetrize} + * function. + */ + void symmetrize (); + + /** + * Copy the given matrix to this + * one. The operation throws an + * error if the sparsity patterns + * of the two involved matrices + * do not point to the same + * object, since in this case the + * copy operation is + * cheaper. Since this operation + * is notheless not for free, we + * do not make it available + * through @p{operator =}, since + * this may lead to unwanted + * usage, e.g. in copy arguments + * to functions, which should + * really be arguments by + * reference. + * + * The source matrix may be a matrix + * of arbitrary type, as long as its + * data type is convertible to the + * data type of this matrix. + * + * The function returns a reference to + * @p{this}. + */ + template + SparseMatrix & + copy_from (const SparseMatrix &source); + + /** + * This function is complete + * analogous to the + * @ref{SparsityPattern}@p{::copy_from} + * function in that it allows to + * initialize a whole matrix in + * one step. See there for more + * information on argument types + * and their meaning. You can + * also find a small example on + * how to use this function + * there. + * + * The only difference to the + * cited function is that the + * objects which the inner + * iterator points to need to be + * of type @p{std::pair + void copy_from (const ForwardIterator begin, + const ForwardIterator end); + + /** + * Copy the nonzero entries of a + * full matrix into this + * object. Previous content is + * deleted. Note that the + * underlying sparsity pattern + * must be appropriate to hold + * the nonzero entries of the + * full matrix. + */ + template + void copy_from (const FullMatrix &matrix); + + /** + * Add @p{matrix} scaled by + * @p{factor} to this matrix. The + * function throws an error if + * the sparsity patterns of the + * two involved matrices do not + * point to the same object, + * since in this case the + * operation is cheaper. + * + * The source matrix may be a matrix + * of arbitrary type, as long as its + * data type is convertible to the + * data type of this matrix. + */ + template + 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 required element does + * not exist in the matrix. + * + * In case you want a function + * that returns zero instead (for + * entries that are not in the + * sparsity pattern of the + * matrix), use the @p{el} + * function. + */ + number operator () (const unsigned int i, + const unsigned int j) const; + + /** + * This function is mostly like + * @p{operator()} in that it + * returns the value of the + * matrix entry @p{(i,j)}. The only + * difference is that if this + * entry does not exist in the + * sparsity pattern, then instead + * of raising an exception, zero + * is returned. While this may be + * convenient in some cases, note + * that it is simple to write + * algorithms that are slow + * compared to an optimal + * solution, since the sparsity + * of the matrix is not used. + */ + number el (const unsigned int i, + const unsigned int j) const; + + /** + * Return the main diagonal element in + * the @p{i}th row. This function throws an + * error if the matrix is not square. + * + * This function is considerably + * faster than the @p{operator()}, + * since for square matrices, the + * diagonal entry is always the + * first to be stored in each row + * and access therefore does not + * involve searching for the + * right column number. + */ + number diag_element (const unsigned int i) const; + + /** + * Same as above, but return a + * writeable reference. You're + * sure you know what you do? + */ + number & diag_element (const unsigned int i); + + /** + * Access to values in internal + * mode. Returns the value of + * the @p{index}th entry in + * @p{row}. Here, @p{index} refers to + * the internal representation of + * the matrix, not the column. Be + * sure to understand what you are + * doing here. + */ + number raw_entry (const unsigned int row, + const unsigned int index) const; + + /** + * This is for hackers. Get + * access to the @p{i}th element of + * this matrix. The elements are + * stored in a consecutive way, + * refer to the @p{SparsityPattern} + * class for more details. + * + * You should use this interface + * very carefully and only if you + * are absolutely sure to know + * what you do. You should also + * note that the structure of + * these arrays may change over + * time. If you change the + * layout yourself, you should + * also rename this function to + * avoid programs relying on + * outdated information! + */ + number global_entry (const unsigned int i) const; + + /** + * Same as above, but with write + * access. You certainly know + * what you do? + */ + number & global_entry (const unsigned int i); + + /** + * Matrix-vector multiplication: + * let $dst = M*src$ with $M$ + * being this matrix. + */ + template + 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 + * @p{vmult} but takes the + * transposed matrix. + */ + template + void Tvmult (Vector &dst, + const Vector &src) const; + + /** + * Adding Matrix-vector + * multiplication. Add $M*src$ on + * $dst$ with $M$ being this + * matrix. + */ + template + void vmult_add (Vector &dst, + const Vector &src) const; + + /** + * Adding Matrix-vector + * multiplication. Add $M^T*src$ + * to $dst$ with $M$ being this + * matrix. This function does the + * same as @p{vmult_add} but takes + * the transposed matrix. + */ + template + void Tvmult_add (Vector &dst, + const Vector &src) const; + + /** + * Return the square of the norm + * of the vector $v$ with respect + * to the norm induced by this + * matrix, + * i.e. $\left(v,Mv\right)$. This + * is useful, e.g. in the finite + * element context, where the + * $L_2$ norm of a function + * equals the matrix norm with + * respect to the mass matrix of + * the vector representing the + * nodal values of the finite + * element function. + * + * Obviously, the matrix needs to + * be square for this operation. + */ + template + somenumber matrix_norm_square (const Vector &v) const; + + /** + * Compute the matrix scalar + * product $\left(u,Mv\right)$. + */ + template + somenumber matrix_scalar_product (const Vector &u, + const Vector &v) const; + + /** + * Return the l1-norm of the matrix, that is + * $|M|_1=max_{all columns j}\sum_{all + * rows i} |M_ij|$, + * (max. sum of columns). + * This is the + * natural matrix norm that is compatible + * to the l1-norm for vectors, i.e. + * $|Mv|_1\leq |M|_1 |v|_1$. + * (cf. Haemmerlin-Hoffmann : Numerische Mathematik) + */ + number l1_norm () const; + + /** + * Return the linfty-norm of the + * matrix, that is + * $|M|_infty=max_{all rows i}\sum_{all + * columns j} |M_ij|$, + * (max. sum of rows). + * This is the + * natural matrix norm that is compatible + * to the linfty-norm of vectors, i.e. + * $|Mv|_infty \leq |M|_infty |v|_infty$. + * (cf. Haemmerlin-Hoffmann : Numerische Mathematik) + */ + number linfty_norm () const; + + /** + * Compute the residual of an + * equation @p{Mx=b}, where the + * residual is defined to be + * @p{r=b-Mx} with @p{x} typically + * being an approximate of the + * true solution of the + * equation. Write the residual + * into @p{dst}. The l2 norm of + * the residual vector is + * returned. + */ + template + somenumber residual (Vector &dst, + const Vector &x, + const Vector &b) const; + + /** + * Apply the Jacobi + * preconditioner, which + * multiplies every element of + * the @p{src} vector by the + * inverse of the respective + * diagonal element and + * multiplies the result with the + * damping factor @p{omega}. + */ + template + void precondition_Jacobi (Vector &dst, + const Vector &src, + const number omega = 1.) const; + + /** + * Apply SSOR preconditioning to + * @p{src}. + */ + template + void precondition_SSOR (Vector &dst, + const Vector &src, + const number om = 1.) const; + + /** + * Apply SOR preconditioning matrix to @p{src}. + * The result of this method is + * $dst = (om D - L)^{-1} src$. + */ + template + void precondition_SOR (Vector &dst, + const Vector &src, + const number om = 1.) const; + + /** + * Apply transpose SOR preconditioning matrix to @p{src}. + * The result of this method is + * $dst = (om D - U)^{-1} src$. + */ + template + void precondition_TSOR (Vector &dst, + const Vector &src, + const number om = 1.) const; + + /** + * Perform SSOR preconditioning + * in-place. Apply the + * preconditioner matrix without + * copying to a second vector. + * @p{omega} is the relaxation + * parameter. + */ + template + void SSOR (Vector &v, + const number omega = 1.) const; + + /** + * Perform an SOR preconditioning in-place. + * The result is $v = (\omega D - L)^{-1} v$. + * @p{omega} is the damping parameter. + */ + template + void SOR (Vector &v, + const number om = 1.) const; + + /** + * Perform a transpose SOR preconditioning in-place. + * The result is $v = (\omega D - L)^{-1} v$. + * @p{omega} is the damping parameter. + */ + template + void TSOR (Vector &v, + const number om = 1.) const; + + /** + * Do one SOR step on @p{v}. + * Performs a direct SOR step + * with right hand side @p{b}. + */ + template + void SOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; + + /** + * Do one adjoint SOR step on + * @p{v}. Performs a direct TSOR + * step with right hand side @p{b}. + */ + template + void TSOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; + + /** + * Do one adjoint SSOR step on + * @p{v}. Performs a direct SSOR + * step with right hand side @p{b} + * by performing TSOR after SOR. + */ + template + void SSOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; + + /** + * Return a (constant) reference + * to the underlying sparsity + * pattern of this matrix. + * + * Though the return value is + * declared @p{const}, you should + * be aware that it may change if + * you call any nonconstant + * function of objects which + * operate on it. + */ + const SparsityPattern & get_sparsity_pattern () const; + + /** + * STL-like iterator with the + * first entry. + */ + const_iterator begin () const; + + /** + * Final iterator. + */ + const_iterator end () const; + + /** + * STL-like iterator with the + * first entry of row @p{r}. + */ + const_iterator begin (unsigned int r) const; + + /** + * Final iterator of row @p{r}. + */ + const_iterator end (unsigned int r) const; + + /** + * Print the matrix to the given + * stream, using the format + * @p{(line,col) value}, i.e. one + * nonzero entry of the matrix + * per line. + */ + void print (std::ostream &out) const; + + /** + * Print the matrix in the usual + * format, i.e. as a matrix and + * not as a list of nonzero + * elements. For better + * readability, elements not in + * the matrix are displayed as + * empty space, while matrix + * elements which are explicitly + * set to zero are displayed as + * such. + * + * The parameters allow for a + * flexible setting of the output + * format: @p{precision} and + * @p{scientific} are used to + * determine the number format, + * where @p{scientific} = @p{false} + * means fixed point notation. A + * zero entry for @p{width} makes + * the function compute a width, + * but it may be changed to a + * positive value, if output is + * crude. + * + * Additionally, a character for + * an empty value may be + * specified. + * + * Finally, the whole matrix can + * be multiplied with a common + * denominator to produce more + * readable output, even + * integers. + * + * This function + * may produce @em{large} amounts of + * output if applied to a large matrix! + */ + void print_formatted (std::ostream &out, + const unsigned int precision = 3, + const bool scientific = true, + const unsigned int width = 0, + const char *zero_string = " ", + const double denominator = 1.) const; + + /** + * Write the data of this object + * en bloc to a file. This is + * done in a binary mode, so the + * output is neither readable by + * humans nor (probably) by other + * computers using a different + * operating system of number + * format. + * + * The purpose of this function + * is that you can swap out + * matrices and sparsity pattern + * if you are short of memory, + * want to communicate between + * different programs, or allow + * objects to be persistent + * across different runs of the + * program. + */ + void block_write (std::ostream &out) const; + + /** + * Read data that has previously + * been written by + * @p{block_write} en block from + * a file. This is done using the + * inverse operations to the + * above function, so it is + * reasonably fast because the + * bitstream is not interpreted + * except for a few numbers up + * front. + * + * The object is resized on this + * operation, and all previous + * contents are lost. Note, + * however, that no checks are + * performed whether new data and + * the underlying + * @ref{SparsityPattern} object + * fit together. It is your + * responsibility to make sure + * that the sparsity pattern and + * the data to be read match. + * + * A primitive form of error + * checking is performed which + * will recognize the bluntest + * attempts to interpret some + * data as a vector stored + * bitwise to a file, but not + * more. + */ + void block_read (std::istream &in); + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + */ + unsigned int memory_consumption () const; + + /** + * Exception + */ + DeclException0 (ExcNotCompressed); + /** + * Exception + */ + DeclException0 (ExcMatrixNotInitialized); + /** + * Exception + */ + DeclException2 (ExcInvalidIndex, + int, int, + << "The entry with index <" << arg1 << ',' << arg2 + << "> does not exist."); + /** + * Exception + */ + DeclException1 (ExcInvalidIndex1, + int, + << "The index " << arg1 << " is not in the allowed range."); + /** + * Exception + */ + DeclException0 (ExcMatrixNotSquare); + /** + * Exception + */ + DeclException0 (ExcDifferentSparsityPatterns); + /** + * Exception + */ + DeclException0 (ExcInvalidConstructorCall); + /** + * Exception + */ + DeclException2 (ExcIteratorRange, + int, int, + << "The iterators denote a range of " << arg1 + << " elements, but the given number of rows was " << arg2); + + private: + /** + * Pointer to the sparsity + * pattern used for this + * matrix. In order to guarantee + * that it is not deleted while + * still in use, we subscribe to + * it using the @p{SmartPointer} + * class. + */ + SmartPointer cols; + + /** + * Array of values for all the + * nonzero entries. The position + * within the matrix, i.e. the + * row and column number for a + * given entry can only be + * deduced using the sparsity + * pattern. The same holds for + * the more common operation of + * finding an entry by its + * coordinates. + */ + number *val; + + /** + * Allocated size of + * @p{val}. This can be larger + * than the actually used part if + * the size of the matrix was + * reduced somewhen in the past + * by associating a sparsity + * pattern with a smaller size to + * this object, using the + * @p{reinit} function. + */ + unsigned int max_len; + + /** + * Version of @p{vmult} which only + * performs its actions on the + * region defined by + * @p{[begin_row,end_row)}. This + * function is called by @p{vmult} + * in the case of enabled + * multithreading. + */ + template + void threaded_vmult (Vector &dst, + const Vector &src, + const unsigned int begin_row, + const unsigned int end_row) const; + + /** + * Version of + * @p{matrix_norm_square} which + * only performs its actions on + * the region defined by + * @p{[begin_row,end_row)}. This + * function is called by + * @p{matrix_norm_square} in the + * case of enabled + * multithreading. + */ + template + void threaded_matrix_norm_square (const Vector &v, + const unsigned int begin_row, + const unsigned int end_row, + somenumber *partial_sum) const; + + /** + * Version of + * @p{matrix_scalar_product} which + * only performs its actions on + * the region defined by + * @p{[begin_row,end_row)}. This + * function is called by + * @p{matrix_scalar_product} in the + * case of enabled + * multithreading. + */ + template + void threaded_matrix_scalar_product (const Vector &u, + const Vector &v, + const unsigned int begin_row, + const unsigned int end_row, + somenumber *partial_sum) const; + + /** + * Version of @p{residual} which + * only performs its actions on + * the region defined by + * @p{[begin_row,end_row)} (these + * numbers are the components of + * @p{interval}). This function is + * called by @p{residual} in the + * case of enabled + * multithreading. + */ + template + void threaded_residual (Vector &dst, + const Vector &u, + const Vector &b, + const std::pair interval, + somenumber *partial_norm) const; + + // make all other sparse matrices + // friends + template friend class SparseMatrix; +}; + + +/*---------------------- Inline functions -----------------------------------*/ + + + +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()); + // it is allowed to set elements of + // the matrix that are not part of + // the sparsity pattern, if the + // value to which we set it is zero + const unsigned int index = cols->operator()(i,j); + Assert ((index != SparsityPattern::invalid_entry) || + (value == 0.), + ExcInvalidIndex(i,j)); + + if (index != SparsityPattern::invalid_entry) + val[index] = value; +}; + + + +template +inline +void SparseMatrix::add (const unsigned int i, + const unsigned int j, + const number value) +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + + const unsigned int index = cols->operator()(i,j); + Assert ((index != SparsityPattern::invalid_entry) || + (value == 0.), + ExcInvalidIndex(i,j)); + + if (value != 0.) + val[index] += value; +}; + + + +template +inline +number SparseMatrix::operator () (const unsigned int i, + const unsigned int j) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (cols->operator()(i,j) != SparsityPattern::invalid_entry, + ExcInvalidIndex(i,j)); + return val[cols->operator()(i,j)]; +}; + + + +template +inline +number SparseMatrix::el (const unsigned int i, + const unsigned int j) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + const unsigned int index = cols->operator()(i,j); + + if (index != SparsityPattern::invalid_entry) + return val[index]; + else + return 0; +}; + + + +template +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::diag_element (const unsigned int i) +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (irowstart[i]]; +}; + + + +template +inline +number +SparseMatrix::raw_entry (const unsigned int row, + const unsigned int index) const +{ + Assert(rowrows, ExcIndexRange(row,0,cols->rows)); + Assert(indexrow_length(row), + ExcIndexRange(index,0,cols->row_length(row))); + + return val[cols->rowstart[row]+index]; +}; + + + +template +inline +number SparseMatrix::global_entry (const unsigned int j) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (j < cols->n_nonzero_elements(), + ExcIndexRange (j, 0, cols->n_nonzero_elements())); + + return val[j]; +}; + + + +template +inline +number & SparseMatrix::global_entry (const unsigned int j) +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (j < cols->n_nonzero_elements(), + ExcIndexRange (j, 0, cols->n_nonzero_elements())); + + return val[j]; +}; + + + +template +template +void +SparseMatrix::copy_from (const ForwardIterator begin, + const ForwardIterator end) +{ + Assert (static_cast(std::distance (begin, end)) == m(), + ExcIteratorRange (std::distance (begin, end), m())); + + // for use in the inner loop, we + // define a typedef to the type of + // the inner iterators + typedef typename std::iterator_traits::value_type::const_iterator inner_iterator; + unsigned int row=0; + for (ForwardIterator i=begin; i!=end; ++i, ++row) + { + const inner_iterator end_of_row = i->end(); + for (inner_iterator j=i->begin(); j!=end_of_row; ++j) + // write entries + set (row, j->first, j->second); + }; +}; + + +//----------------------------------------------------------------------// + +template +inline +SparseMatrix::Accessor::Accessor ( + const SparseMatrix* matrix, + unsigned int r, + unsigned short i) + : + matrix(matrix), + a_row(r), + a_index(i) +{} + + +template +inline +unsigned int +SparseMatrix::Accessor::row() const +{ + return a_row; +} + + +template +inline +unsigned int +SparseMatrix::Accessor::column() const +{ + const SparsityPattern& pat = matrix->get_sparsity_pattern(); + return pat.get_column_numbers()[pat.get_rowstart_indices()[a_row]+a_index]; +} + + +template +inline +unsigned short +SparseMatrix::Accessor::index() const +{ + return a_index; +} + + + +template +inline +number +SparseMatrix::Accessor::value() const +{ + return matrix->raw_entry(a_row, a_index); +} + + +template +inline +SparseMatrix::const_iterator::const_iterator( + const SparseMatrix* matrix, + unsigned int r, + unsigned short i) + : + Accessor(matrix, r, i) +{} + + +template +inline +typename SparseMatrix::const_iterator& +SparseMatrix::const_iterator::operator++ () +{ + Assert (a_row < matrix->m(), ExcIteratorPastEnd()); + + ++a_index; + if (a_index >= matrix->get_sparsity_pattern().row_length(a_row)) + { + a_index = 0; + a_row++; + } + return *this; +} + + +template +inline +const typename SparseMatrix::Accessor& +SparseMatrix::const_iterator::operator* () const +{ + return *this; +} + + +template +inline +const typename SparseMatrix::Accessor* +SparseMatrix::const_iterator::operator-> () const +{ + return this; +} + + +template +inline +bool +SparseMatrix::const_iterator::operator == ( + const const_iterator& other) const +{ + return (row() == other->row() && index() == other->index()); +} + + +template +inline +bool +SparseMatrix::const_iterator::operator != ( + const const_iterator& other) const +{ + return ! (*this == other); +} + + +template +inline +bool +SparseMatrix::const_iterator::operator < ( + const const_iterator& other) const +{ + return (row() < other->row() || + (row() == other->row() && index() < other->index())); +} + + +template +inline +typename SparseMatrix::const_iterator +SparseMatrix::begin () const +{ + return const_iterator(this, 0, 0); +} + +template +inline +typename SparseMatrix::const_iterator +SparseMatrix::end () const +{ + return const_iterator(this, m(), 0); +} + +template +inline +typename SparseMatrix::const_iterator +SparseMatrix::begin (unsigned int r) const +{ + Assert (r +inline +typename SparseMatrix::const_iterator +SparseMatrix::end (unsigned int r) const +{ + Assert (r