From: wolf Date: Mon, 2 Dec 2002 17:29:16 +0000 (+0000) Subject: These files should not have been picked up. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0d5ed92119500a10d7e85085e4cdd677fed61692;p=dealii-svn.git These files should not have been picked up. git-svn-id: https://svn.dealii.org/trunk@6793 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h.x b/deal.II/lac/include/lac/sparse_ilu.templates.h.x deleted file mode 100644 index 608d889a0c..0000000000 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h.x +++ /dev/null @@ -1,353 +0,0 @@ -//---------------------------- sparse_ilu.templates.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_ilu.templates.h --------------------------- -#ifndef __deal2__sparse_ilu_templates_h -#define __deal2__sparse_ilu_templates_h - - - -#include -#include -#include - -#include -#include - - -template -SparseILU::SparseILU () -{}; - - - -template -SparseILU::SparseILU (const SparsityPattern &sparsity) : - SparseMatrix (sparsity) -{}; - - - -template -void SparseILU::reinit () -{ - SparseMatrix::reinit (); -}; - - - -template -void SparseILU::reinit (const SparsityPattern &sparsity) -{ - SparseMatrix::reinit (sparsity); -}; - - - -template -template -void SparseILU::decompose (const SparseMatrix &matrix, - const double strengthen_diagonal) -{ - Assert (matrix.m()==matrix.n(), ExcMatrixNotSquare ()); - Assert (this->m()==this->n(), ExcMatrixNotSquare ()); - Assert (matrix.m()==this->m(), ExcSizeMismatch(matrix.m(), this->m())); - - Assert (strengthen_diagonal>=0, ExcInvalidStrengthening (strengthen_diagonal)); - - - // first thing: copy over all elements - // of @p{matrix} to the present object - // - // note that some elements in this - // matrix may not be in @p{matrix}, - // so we need to preset our matrix - // by zeroes. - if (true) - { - // preset the elements - std::fill_n (&this->global_entry(0), - this->n_nonzero_elements(), - 0); - - // note: pointers to the sparsity - // pattern of the old matrix! - const unsigned int * const rowstart_indices - = matrix.get_sparsity_pattern().get_rowstart_indices(); - const unsigned int * const column_numbers - = matrix.get_sparsity_pattern().get_column_numbers(); - - for (unsigned int row=0; rowm(); ++row) - for (const unsigned int * col = &column_numbers[rowstart_indices[row]]; - col != &column_numbers[rowstart_indices[row+1]]; ++col) - set (row, *col, matrix.global_entry(col-column_numbers)); - }; - - if (strengthen_diagonal > 0) - for (unsigned int row=0; rowm(); ++row) - { - // get the length of the row - // (without the diagonal element) - const unsigned int - rowlength = (this->get_sparsity_pattern().get_rowstart_indices()[row+1] - - - this->get_sparsity_pattern().get_rowstart_indices()[row] - - - 1); - - // get the global index of the first - // non-diagonal element in this row - const unsigned int rowstart - = this->get_sparsity_pattern().get_rowstart_indices()[row] + 1; - number * const diagonal_element = &this->global_entry(rowstart-1); - - number rowsum = 0; - for (unsigned int global_index=rowstart; - global_indexglobal_entry(global_index)); - - *diagonal_element += strengthen_diagonal * rowsum; - }; - - - // now work only on this - // matrix - const SparsityPattern &sparsity = this->get_sparsity_pattern(); - const unsigned int * const rowstart_indices = sparsity.get_rowstart_indices(); - const unsigned int * const column_numbers = sparsity.get_column_numbers(); - -/* - PSEUDO-ALGORITHM - (indices=0..N-1) - - for i=1..N-1 - a[i-1,i-1] = a[i-1,i-1]^{-1} - - for k=0..i-1 - a[i,k] = a[i,k] * a[k,k] - - for j=k+1..N-1 - if (a[i,j] exists & a[k,j] exists) - a[i,j] -= a[i,k] * a[k,j] -*/ - - - // i := row - for (unsigned int row=1; rowm(); ++row) - { - // invert diagonal element of the - // previous row. this is a hack, - // which is possible since this - // element is not needed any more - // in the process of decomposition - // and since it makes the backward - // step when applying the decomposition - // significantly faster - AssertThrow((this->global_entry(rowstart_indices[row-1]) !=0), - ExcDivideByZero()); - - this->global_entry (rowstart_indices[row-1]) - = 1./this->global_entry (rowstart_indices[row-1]); - - // let k run over all lower-left - // elements of row i; skip - // diagonal element at start - const unsigned int * first_of_row - = &column_numbers[rowstart_indices[row]+1]; - const unsigned int * first_after_diagonal - = std::lower_bound (&column_numbers[rowstart_indices[row]+1], - &column_numbers[rowstart_indices[row+1]], - row); - - // k := *col_ptr - for (const unsigned int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++col_ptr) - { - const unsigned int global_index_ik = col_ptr-column_numbers; - this->global_entry(global_index_ik) *= this->diag_element(*col_ptr); - - // now do the inner loop over - // j. note that we need to do - // it in the right order, i.e. - // taking into account that the - // columns are sorted within each - // row correctly, but excluding - // the main diagonal entry - const int global_index_ki = sparsity(*col_ptr,row); - - if (global_index_ki != -1) - this->diag_element(row) -= this->global_entry(global_index_ik) * - this->global_entry(global_index_ki); - - for (const unsigned int * j = col_ptr+1; - j<&column_numbers[rowstart_indices[row+1]]; - ++j) - { -//TODO:[WB] make code faster by using the following comment - // note: this inner loop could - // be made considerable faster - // if we consulted the row - // with number *col_ptr, - // instead of always asking - // sparsity(*col_ptr,*j), - // since we traverse this - // row linearly. I just didn't - // have the time to figure out - // the details. - const int global_index_ij = j - &column_numbers[0], - global_index_kj = sparsity(*col_ptr,*j); - if ((global_index_ij != -1) && - (global_index_kj != -1)) - this->global_entry(global_index_ij) -= this->global_entry(global_index_ik) * - this->global_entry(global_index_kj); - }; - }; - }; - - // Here the very last diagonal - // element still has to be inverted - // because the for-loop doesn't do - // it... - this->diag_element(this->m()-1) = 1./this->diag_element(this->m()-1); - -/* - OLD CODE, rather crude first implementation with an algorithm taken - from 'W. Hackbusch, G. Wittum: Incomplete Decompositions (ILU)- - Algorithms, Theory, and Applications', page 6. - - for (unsigned int k=0; k -template -void SparseILU::apply_decomposition (Vector &dst, - const Vector &src) const -{ - Assert (dst.size() == src.size(), ExcSizeMismatch(dst.size(), src.size())); - Assert (dst.size() == this->m(), ExcSizeMismatch(dst.size(), this->m())); - - const unsigned int N=dst.size(); - const unsigned int * const rowstart_indices - = this->get_sparsity_pattern().get_rowstart_indices(); - const unsigned int * const column_numbers - = this->get_sparsity_pattern().get_column_numbers(); - // solve LUx=b in two steps: - // first Ly = b, then - // Ux = y - // - // first a forward solve. since - // the diagonal values of L are - // one, there holds - // y_i = b_i - // - sum_{j=0}^{i-1} L_{ij}y_j - // we split the y_i = b_i off and - // perform it at the outset of the - // loop - dst = src; - for (unsigned int row=0; rowglobal_entry (col-column_numbers) * dst(*col); - }; - - // now the backward solve. same - // procedure, but we need not set - // dst before, since this is already - // done. - // - // note that we need to scale now, - // since the diagonal is not zero - // now - for (int row=N-1; row>=0; --row) - { - // get end of this row - const unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]]; - // find the position where the part - // right of the diagonal starts - const unsigned int * const first_after_diagonal - = std::lower_bound (&column_numbers[rowstart_indices[row]+1], - &column_numbers[rowstart_indices[row+1]], - static_cast(row)); - - for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col) - dst(row) -= this->global_entry (col-column_numbers) * dst(*col); - - // scale by the diagonal element. - // note that the diagonal element - // was stored inverted - dst(row) *= this->diag_element(row); - }; -}; - - - -template -unsigned int -SparseILU::memory_consumption () const -{ - return SparseMatrix::memory_consumption (); -}; - - - -/*---------------------------- sparse_ilu.templates.h ---------------------------*/ - -#endif -/*---------------------------- sparse_ilu.templates.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h deleted file mode 100644 index 29c71926bc..0000000000 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ /dev/null @@ -1,1491 +0,0 @@ -//---------------------------- 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 - - -// explicit instantiations -template class SparseILU; -template void SparseILU::decompose (const SparseMatrix &, - const double); -template void SparseILU::apply_decomposition (Vector &, - const Vector &) const; -template void SparseILU::decompose (const SparseMatrix &, - const double); -template void SparseILU::apply_decomposition (Vector &, - const Vector &) const; - - -template class SparseILU; -template void SparseILU::decompose (const SparseMatrix &, - const double); -template void SparseILU::apply_decomposition (Vector &, - const Vector &) const; -template void SparseILU::decompose (const SparseMatrix &, - const double); -template void SparseILU::apply_decomposition (Vector &, - const Vector &) const;