From: bangerth Date: Wed, 6 Aug 2008 05:00:56 +0000 (+0000) Subject: Add the Chunk version of the SparseMatrix class. Pretty much all preconditioners... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0859bd7455c7d14dd24e8d36df44c20c20cb1035;p=dealii-svn.git Add the Chunk version of the SparseMatrix class. Pretty much all preconditioners aren't currently implemented, but will hopefully soon. git-svn-id: https://svn.dealii.org/trunk@16496 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/chunk_sparse_matrix.h b/deal.II/lac/include/lac/chunk_sparse_matrix.h new file mode 100644 index 0000000000..986c53b0ae --- /dev/null +++ b/deal.II/lac/include/lac/chunk_sparse_matrix.h @@ -0,0 +1,1467 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name: $ +// +// Copyright (C) 2008 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. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__chunk_sparse_matrix_h +#define __deal2__chunk_sparse_matrix_h + + +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template class Vector; +template class FullMatrix; + +/*! @addtogroup Matrix1 + *@{ + */ + + +/** + * Sparse matrix. This class implements the function to store values + * in the locations of a sparse matrix denoted by a + * SparsityPattern. The separation of sparsity pattern and values is + * done since one can store data elements of different type in these + * locations without the SparsityPattern having to know this, and more + * importantly one can associate more than one matrix with the same + * sparsity pattern. + * + * @note Instantiations for this template are provided for @ and + * @; others can be generated in application programs (see the + * section on @ref Instantiations in the manual). + * + * @author Wolfgang Bangerth, 2008 + */ +template +class ChunkSparseMatrix : public virtual Subscriptor +{ + public: + /** + * Type of matrix entries. In analogy to + * the STL container classes. + */ + typedef number value_type; + + /** + * Declare a type that has holds + * real-valued numbers with the + * same precision as the template + * argument to this class. If the + * template argument of this + * class is a real data type, + * then real_type equals the + * template argument. If the + * template argument is a + * std::complex type then + * real_type equals the type + * underlying the complex + * numbers. + * + * This typedef is used to + * represent the return type of + * norms. + */ + typedef typename numbers::NumberTraits::real_type real_type; + + /** + * A structure that describes some of the + * traits of this class in terms of its + * run-time behavior. Some other classes + * (such as the block matrix classes) + * that take one or other of the matrix + * classes as its template parameters can + * tune their behavior based on the + * variables in this class. + */ + struct Traits + { + /** + * It is safe to elide additions of + * zeros to individual elements of + * this matrix. + */ + static const bool zero_addition_can_be_elided = true; + }; + +/** + * @name Constructors and initalization. + */ +//@{ + /** + * 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(const ChunkSparsityPattern&). + */ + ChunkSparseMatrix (); + + /** + * 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 + * ChunkSparsityPattern, see there for the + * details. + * + * If you really want to copy a whole + * matrix, you can do so by using the + * copy_from() function. + */ + ChunkSparseMatrix (const ChunkSparseMatrix &); + + /** + * 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 reinit(const + * ChunkSparsityPattern&) 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(const + * ChunkSparsityPattern&) is not + * called with a new sparsity + * pattern. + * + * 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 ChunkSparseMatrix (const ChunkSparsityPattern &sparsity); + + /** + * Copy constructor: initialize + * the matrix with the identity + * matrix. This constructor will + * throw an exception if the + * sizes of the sparsity pattern + * and the identity matrix do not + * coincide, or if the sparsity + * pattern does not provide for + * nonzero entries on the entire + * diagonal. + */ + ChunkSparseMatrix (const ChunkSparsityPattern &sparsity, + const IdentityMatrix &id); + + /** + * Destructor. Free all memory, but do not + * release the memory of the sparsity + * structure. + */ + virtual ~ChunkSparseMatrix (); + + /** + * Copy operator. Since copying + * entire sparse matrices is a + * very expensive operation, we + * disallow doing so except for + * the special case of empty + * matrices of size zero. This + * doesn't seem particularly + * useful, but is exactly what + * one needs if one wanted to + * have a + * std::vector@ + * @>: in that case, one + * can create a vector (which + * needs the ability to copy + * objects) of empty matrices + * that are then later filled + * with something useful. + */ + ChunkSparseMatrix& operator = (const ChunkSparseMatrix &); + + /** + * Copy operator: initialize + * the matrix with the identity + * matrix. This operator will + * throw an exception if the + * sizes of the sparsity pattern + * and the identity matrix do not + * coincide, or if the sparsity + * pattern does not provide for + * nonzero entries on the entire + * diagonal. + */ + ChunkSparseMatrix & + operator= (const IdentityMatrix &id); + + /** + * This operator assigns a scalar to + * a matrix. Since this does usually + * not make much sense (should we set + * all matrix entries to this value? + * Only the nonzero entries of the + * sparsity pattern?), this operation + * is only allowed if the actual + * value to be assigned is zero. This + * operator only exists to allow for + * the obvious notation + * matrix=0, which sets all + * elements of the matrix to zero, + * but keep the sparsity pattern + * previously used. + */ + ChunkSparseMatrix & operator = (const double d); + + /** + * 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(const + * ChunkSparsityPattern &) is not + * called with a new sparsity + * structure. + * + * The elements of the matrix are + * set to zero by this function. + */ + virtual void reinit (const ChunkSparsityPattern &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 (); +//@} +/** + * @name Information on the matrix + */ +//@{ + /** + * Return whether the object is + * empty. It is empty if either + * both dimensions are zero or no + * ChunkSparsityPattern 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 + * 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; + + /** + * Return a (constant) reference + * to the underlying sparsity + * pattern of this matrix. + * + * 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. + */ + const ChunkSparsityPattern & get_sparsity_pattern () const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. See + * MemoryConsumption. + */ + unsigned int memory_consumption () const; + +//@} +/** + * @name Modifying entries + */ +//@{ + /** + * Set the element (i,j) + * to value. Throws an + * error if the entry does not + * exist or if value is + * not a finite number. 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 or if value is + * not a finite number. 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); + + /** + * Multiply the entire matrix by a + * fixed factor. + */ + ChunkSparseMatrix & operator *= (const number factor); + + /** + * Divide the entire matrix by a + * fixed factor. + */ + ChunkSparseMatrix & operator /= (const number factor); + + /** + * 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 + * ChunkSparsityPattern::symmetrize(). + */ + 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 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 + ChunkSparseMatrix & + copy_from (const ChunkSparseMatrix &source); + + /** + * This function is complete + * analogous to the + * ChunkSparsityPattern::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 std::pair, where + * value needs to be + * convertible to the element + * type of this class, as + * specified by the + * number template + * argument. + * + * Previous content of the matrix + * is overwritten. Note that the + * entries specified by the input + * parameters need not + * necessarily cover all elements + * of the matrix. Elements not + * covered remain untouched. + */ + template + 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 matrix scaled by + * factor to this matrix, + * i.e. the matrix factor*matrix + * is added to this. This + * 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 sparse + * matrix over an arbitrary underlying + * scalar type, as long as its data type + * is convertible to the data type of + * this matrix. + */ + template + void add (const number factor, + const ChunkSparseMatrix &matrix); + +//@} +/** + * @name Entry Access + */ +//@{ + + /** + * 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 el() + * function. + * + * If you are looping over all elements, + * consider using one of the iterator + * classes instead, since they are + * tailored better to a sparse matrix + * structure. + */ + number operator () (const unsigned int i, + const unsigned int j) const; + + /** + * This function is mostly like + * operator()() in that it + * returns the value of the + * matrix entry (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. + * + * If you are looping over all elements, + * consider using one of the iterator + * classes instead, since they are + * tailored better to a sparse matrix + * structure. + */ + number el (const unsigned int i, + const unsigned int j) const; + + /** + * Return the main diagonal + * element in the ith + * row. This function throws an + * error if the matrix is not + * quadratic (see + * ChunkSparsityPattern::optimize_diagonal()). + * + * This function is considerably + * faster than the operator()(), + * since for quadratic matrices, the + * diagonal entry may be 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); + +//@} +/** + * @name Matrix vector multiplications + */ +//@{ + /** + * Matrix-vector multiplication: + * let dst = M*src with + * M being this matrix. + * + * Note that while this function can + * operate on all vectors that offer + * iterator classes, it is only really + * effective for objects of type @ref + * Vector. For all classes for which + * iterating over elements, or random + * member access is expensive, this + * function is not efficient. In + * particular, if you want to multiply + * with BlockVector objects, you should + * consider using a BlockChunkSparseMatrix as + * well. + * + * Source and destination must + * not be the same vector. + */ + template + void vmult (OutVector& dst, + const InVector& src) const; + + /** + * Matrix-vector multiplication: + * let dst = MT*src with + * M being this + * matrix. This function does the + * same as vmult() but takes + * the transposed matrix. + * + * Note that while this function can + * operate on all vectors that offer + * iterator classes, it is only really + * effective for objects of type @ref + * Vector. For all classes for which + * iterating over elements, or random + * member access is expensive, this + * function is not efficient. In + * particular, if you want to multiply + * with BlockVector objects, you should + * consider using a BlockChunkSparseMatrix as + * well. + * + * Source and destination must + * not be the same vector. + */ + template + void Tvmult (OutVector& dst, + const InVector& src) const; + + /** + * Adding Matrix-vector + * multiplication. Add + * M*src on dst + * with M being this + * matrix. + * + * Note that while this function can + * operate on all vectors that offer + * iterator classes, it is only really + * effective for objects of type @ref + * Vector. For all classes for which + * iterating over elements, or random + * member access is expensive, this + * function is not efficient. In + * particular, if you want to multiply + * with BlockVector objects, you should + * consider using a BlockChunkSparseMatrix as + * well. + * + * Source and destination must + * not be the same vector. + */ + template + void vmult_add (OutVector& dst, + const InVector& src) const; + + /** + * Adding Matrix-vector + * multiplication. Add + * MT*src to + * dst with M being + * this matrix. This function + * does the same as vmult_add() + * but takes the transposed + * matrix. + * + * Note that while this function can + * operate on all vectors that offer + * iterator classes, it is only really + * effective for objects of type @ref + * Vector. For all classes for which + * iterating over elements, or random + * member access is expensive, this + * function is not efficient. In + * particular, if you want to multiply + * with BlockVector objects, you should + * consider using a BlockChunkSparseMatrix as + * well. + * + * Source and destination must + * not be the same vector. + */ + template + void Tvmult_add (OutVector& dst, + const InVector& 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 + * quadratic for this operation, and for + * the result to actually be a norm it + * also needs to be either real symmetric + * or complex hermitian. + * + * The underlying template types of both + * this matrix and the given vector + * should either both be real or + * complex-valued, but not mixed, for + * this function to make sense. + */ + 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; + /** + * Compute the residual of an + * equation Mx=b, where + * the residual is defined to be + * r=b-Mx. Write the + * residual into + * dst. The + * l2 norm of + * the residual vector is + * returned. + * + * Source x and destination + * dst must not be the same + * vector. + */ + template + somenumber residual (Vector &dst, + const Vector &x, + const Vector &b) const; + +//@} +/** + * @name Matrix norms + */ +//@{ + + /** + * 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) + */ + real_type 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) + */ + real_type linfty_norm () const; + + /** + * Return the frobenius norm of the + * matrix, i.e. the square root of the + * sum of squares of all entries in the + * matrix. + */ + real_type frobenius_norm () const; +//@} +/** + * @name Preconditioning methods + */ +//@{ + + /** + * Apply the Jacobi + * preconditioner, which + * multiplies every element of + * the src vector by the + * inverse of the respective + * diagonal element and + * multiplies the result with the + * relaxation factor omega. + */ + template + void precondition_Jacobi (Vector &dst, + const Vector &src, + const number omega = 1.) const; + + /** + * Apply SSOR preconditioning to + * src. + */ + template + void precondition_SSOR (Vector &dst, + const Vector &src, + const number om = 1.) const; + + /** + * Apply SOR preconditioning + * matrix to src. + */ + template + void precondition_SOR (Vector &dst, + const Vector &src, + const number om = 1.) const; + + /** + * Apply transpose SOR + * preconditioning matrix to + * 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. + * omega is the relaxation + * parameter. + */ + template + void SSOR (Vector &v, + const number omega = 1.) const; + + /** + * Perform an SOR preconditioning + * in-place. omega is + * the relaxation parameter. + */ + template + void SOR (Vector &v, + const number om = 1.) const; + + /** + * Perform a transpose SOR + * preconditioning in-place. + * omega is the + * relaxation parameter. + */ + template + void TSOR (Vector &v, + const number om = 1.) const; + + /** + * Perform a permuted SOR + * preconditioning in-place. + * + * The standard SOR method is + * applied in the order + * prescribed by permutation, + * that is, first the row + * permutation[0], then + * permutation[1] and so + * on. For efficiency reasons, + * the permutation as well as its + * inverse are required. + * + * omega is the + * relaxation parameter. + */ + template + void PSOR (Vector &v, + const std::vector& permutation, + const std::vector& inverse_permutation, + const number om = 1.) const; + + /** + * Perform a transposed permuted SOR + * preconditioning in-place. + * + * The transposed SOR method is + * applied in the order + * prescribed by + * permutation, that is, + * first the row + * permutation[m()-1], + * then + * permutation[m()-2] + * and so on. For efficiency + * reasons, the permutation as + * well as its inverse are + * required. + * + * omega is the + * relaxation parameter. + */ + template + void TPSOR (Vector &v, + const std::vector& permutation, + const std::vector& inverse_permutation, + const number om = 1.) const; + + /** + * Do one SOR step on v. + * Performs a direct SOR step + * with right hand side + * b. + */ + template + void SOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; + + /** + * Do one adjoint SOR step on + * v. Performs a direct + * TSOR step with right hand side + * b. + */ + template + void TSOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; + + /** + * Do one SSOR step on + * v. Performs a direct + * SSOR step with right hand side + * b by performing TSOR + * after SOR. + */ + template + void SSOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; +//@} +/** + * @name Input/Output + */ +//@{ + + /** + * 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 (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: precision and + * scientific are used + * to determine the number + * format, where scientific = + * false means fixed point + * notation. A zero entry for + * 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. + * + * @attention This function may + * produce 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; + + /** + * Print the actual pattern of + * the matrix. For each entry + * with an absolute value larger + * than threshold, a '*' is + * printed, a ':' for every value + * smaller and a '.' for every + * entry not allocated. + */ + void print_pattern(std::ostream& out, + const double threshold = 0.) 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 block_write() + * 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 ChunkSparsityPattern + * 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 matrix stored + * bitwise to a file that wasn't + * actually created that way, but + * not more. + */ + void block_read (std::istream &in); +//@} + /** @addtogroup Exceptions + * @{ */ + + /** + * 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 (ExcDifferentChunkSparsityPatterns); + /** + * Exception + */ + DeclException2 (ExcIteratorRange, + int, int, + << "The iterators denote a range of " << arg1 + << " elements, but the given number of rows was " << arg2); + /** + * Exception + */ + DeclException0 (ExcSourceEqualsDestination); + //@} + 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 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 #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 reinit() + * function. + */ + unsigned int max_len; + + /** + * Return the location of entry + * $(i,j)$ within the val array. + */ + unsigned int compute_location (const unsigned int i, + const unsigned int j) const; + + // make all other sparse matrices + // friends + template friend class ChunkSparseMatrix; +}; + +/*@}*/ + +#ifndef DOXYGEN +/*---------------------- Inline functions -----------------------------------*/ + + + +template +inline +unsigned int ChunkSparseMatrix::m () const +{ + Assert (cols != 0, ExcNotInitialized()); + return cols->rows; +} + + +template +inline +unsigned int ChunkSparseMatrix::n () const +{ + Assert (cols != 0, ExcNotInitialized()); + return cols->cols; +} + + + +template +inline +unsigned int +ChunkSparseMatrix::compute_location (const unsigned int i, + const unsigned int j) const +{ + const unsigned int chunk_size = cols->get_chunk_size(); + const unsigned int chunk_index + = cols->sparsity_pattern(i/chunk_size, j/chunk_size); + + if (chunk_index == ChunkSparsityPattern::invalid_entry) + return ChunkSparsityPattern::invalid_entry; + else + { + return (chunk_index * chunk_size * chunk_size + + + (i % chunk_size) * chunk_size + + + (j % chunk_size)); + } +} + + +template +inline +void ChunkSparseMatrix::set (const unsigned int i, + const unsigned int j, + const number value) +{ + + Assert (numbers::is_finite(value), + ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); + + Assert (cols != 0, ExcNotInitialized()); + // 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 = compute_location(i,j); + Assert ((index != SparsityPattern::invalid_entry) || + (value == 0.), + ExcInvalidIndex(i,j)); + + if (index != SparsityPattern::invalid_entry) + val[index] = value; +} + + + +template +inline +void ChunkSparseMatrix::add (const unsigned int i, + const unsigned int j, + const number value) +{ + + Assert (numbers::is_finite(value), + ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); + + Assert (cols != 0, ExcNotInitialized()); + + const unsigned int index = compute_location(i,j); + Assert ((index != ChunkSparsityPattern::invalid_entry) || + (value == 0.), + ExcInvalidIndex(i,j)); + + if (value != 0.) + val[index] += value; +} + + + +template +inline +ChunkSparseMatrix & +ChunkSparseMatrix::operator *= (const number factor) +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + + const unsigned int chunk_size = cols->get_chunk_size(); + + // multiply all elements of the matrix with + // the given factor. this includes the + // padding elements in chunks that overlap + // the boundaries of the actual matrix -- + // but since multiplication with a number + // does not violate the invariant of + // keeping these elements at zero nothing + // can happen + number *val_ptr = val; + const number *const end_ptr = val + + cols->sparsity_pattern.n_nonzero_elements() + * + chunk_size * chunk_size; + while (val_ptr != end_ptr) + *val_ptr++ *= factor; + + return *this; +} + + + +template +inline +ChunkSparseMatrix & +ChunkSparseMatrix::operator /= (const number factor) +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (factor !=0, ExcDivideByZero()); + + const number factor_inv = 1. / factor; + + const unsigned int chunk_size = cols->get_chunk_size(); + + // multiply all elements of the matrix with + // the given factor. this includes the + // padding elements in chunks that overlap + // the boundaries of the actual matrix -- + // but since multiplication with a number + // does not violate the invariant of + // keeping these elements at zero nothing + // can happen + number *val_ptr = val; + const number *const end_ptr = val + + cols->sparsity_pattern.n_nonzero_elements() + * + chunk_size * chunk_size; + + while (val_ptr != end_ptr) + *val_ptr++ *= factor_inv; + + return *this; +} + + + +template +inline +number ChunkSparseMatrix::operator () (const unsigned int i, + const unsigned int j) const +{ + Assert (cols != 0, ExcNotInitialized()); + AssertThrow (compute_location(i,j) != SparsityPattern::invalid_entry, + ExcInvalidIndex(i,j)); + return val[compute_location(i,j)]; +} + + + +template +inline +number ChunkSparseMatrix::el (const unsigned int i, + const unsigned int j) const +{ + Assert (cols != 0, ExcNotInitialized()); + const unsigned int index = compute_location(i,j); + + if (index != ChunkSparsityPattern::invalid_entry) + return val[index]; + else + return 0; +} + + + +template +inline +number ChunkSparseMatrix::diag_element (const unsigned int i) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), ExcNotQuadratic()); + Assert (iget_chunk_size(); + return val[cols->sparsity_pattern.rowstart[i/chunk_size] + * + chunk_size * chunk_size + + + (i % chunk_size) * chunk_size + + + (i % chunk_size)]; +} + + + +template +inline +number & ChunkSparseMatrix::diag_element (const unsigned int i) +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), ExcNotQuadratic()); + Assert (iget_chunk_size(); + return val[cols->sparsity_pattern.rowstart[i/chunk_size] + * + chunk_size * chunk_size + + + (i % chunk_size) * chunk_size + + + (i % chunk_size)]; +} + + + +template +template +void +ChunkSparseMatrix::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); + } +} + + +#endif // DOXYGEN + + +/*---------------------------- chunk_sparse_matrix.h ---------------------------*/ + +DEAL_II_NAMESPACE_CLOSE + +#endif +/*---------------------------- chunk_sparse_matrix.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/chunk_sparse_matrix.templates.h b/deal.II/lac/include/lac/chunk_sparse_matrix.templates.h new file mode 100644 index 0000000000..b0d63b6275 --- /dev/null +++ b/deal.II/lac/include/lac/chunk_sparse_matrix.templates.h @@ -0,0 +1,1591 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2008 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. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__chunk_sparse_matrix_templates_h +#define __deal2__chunk_sparse_matrix_templates_h + + +#include +#include +#include +#include + + +// we only need output streams, but older compilers did not provide +// them in a separate include file +#ifdef HAVE_STD_OSTREAM_HEADER +# include +#else +# include +#endif + +#include +#include +#include +#include + +#include +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace internal +{ +//TODO: the goal of the ChunkSparseMatrix class is to stream data and use +// the vectorization features of modern processors. to make this happen, +// we will have to vectorize the functions in the following namespace, either +// by hand or by using, for example, optimized BLAS versions for them. + namespace ChunkSparseMatrix + { + /** + * Add the result of multiplying a chunk + * of size chunk_size times chunk_size by + * a source vector fragment of size + * chunk_size to the destination vector + * fragment. + */ + template + inline + void + chunk_vmult_add (const unsigned int chunk_size, + const MatrixIterator matrix, + const SrcIterator src, + DstIterator dst) + { + MatrixIterator matrix_row = matrix; + + for (unsigned int i=0; i::value_type + sum = 0; + + for (unsigned int j=0; j + inline + void + chunk_vmult_subtract (const unsigned int chunk_size, + const MatrixIterator matrix, + const SrcIterator src, + DstIterator dst) + { + MatrixIterator matrix_row = matrix; + + for (unsigned int i=0; i::value_type + sum = 0; + + for (unsigned int j=0; j + inline + void + chunk_Tvmult_add (const unsigned int chunk_size, + const MatrixIterator matrix, + const SrcIterator src, + DstIterator dst) + { + for (unsigned int i=0; i::value_type + sum = 0; + + for (unsigned int j=0; j + inline + result_type + chunk_matrix_scalar_product (const unsigned int chunk_size, + const MatrixIterator matrix, + const SrcIterator1 u, + const SrcIterator2 v) + { + result_type result = 0; + + MatrixIterator matrix_row = matrix; + + for (unsigned int i=0; i::value_type + sum = 0; + + for (unsigned int j=0; j +ChunkSparseMatrix::ChunkSparseMatrix () + : + cols(0, "ChunkSparseMatrix"), + val(0), + max_len(0) +{} + + + +template +ChunkSparseMatrix::ChunkSparseMatrix (const ChunkSparseMatrix &m) + : + Subscriptor (m), + cols(0, "ChunkSparseMatrix"), + val(0), + max_len(0) +{ + Assert (m.cols==0, ExcInvalidConstructorCall()); + Assert (m.val==0, ExcInvalidConstructorCall()); + Assert (m.max_len==0, ExcInvalidConstructorCall()); +} + + + +template +ChunkSparseMatrix& +ChunkSparseMatrix::operator = (const ChunkSparseMatrix &m) +{ + Assert (m.cols==0, ExcInvalidConstructorCall()); + Assert (m.val==0, ExcInvalidConstructorCall()); + Assert (m.max_len==0, ExcInvalidConstructorCall()); + + return *this; +} + + + +template +ChunkSparseMatrix::ChunkSparseMatrix (const ChunkSparsityPattern &c) + : + cols(0, "ChunkSparseMatrix"), + val(0), + max_len(0) +{ + reinit (c); +} + + + +template +ChunkSparseMatrix::ChunkSparseMatrix (const ChunkSparsityPattern &c, + const IdentityMatrix &id) + : + cols(0, "ChunkSparseMatrix"), + val(0), + max_len(0) +{ + Assert (c.n_rows() == id.m(), ExcDimensionMismatch (c.n_rows(), id.m())); + Assert (c.n_cols() == id.n(), ExcDimensionMismatch (c.n_cols(), id.n())); + + reinit (c); + for (unsigned int i=0; iset(i,i,1.); +} + + + +template +ChunkSparseMatrix::~ChunkSparseMatrix () +{ + cols = 0; + + if (val != 0) + delete[] val; +} + + + +template +ChunkSparseMatrix & +ChunkSparseMatrix::operator = (const double d) +{ + Assert (d==0, ExcScalarAssignmentOnlyForZeroValue()); + + Assert (cols != 0, ExcNotInitialized()); + Assert (cols->sparsity_pattern.compressed || cols->empty(), + ChunkSparsityPattern::ExcNotCompressed()); + + if (val) + { + const unsigned int chunk_size = cols->get_chunk_size(); + std::fill_n (val, + val + + cols->sparsity_pattern.n_nonzero_elements() * + chunk_size * chunk_size, + 0.); + } + + return *this; +} + + + +template +ChunkSparseMatrix & +ChunkSparseMatrix::operator= (const IdentityMatrix &id) +{ + Assert (cols->n_rows() == id.m(), + ExcDimensionMismatch (cols->n_rows(), id.m())); + Assert (cols->n_cols() == id.n(), + ExcDimensionMismatch (cols->n_cols(), id.n())); + + *this = 0; + for (unsigned int i=0; iset(i,i,1.); + + return *this; +} + + + +template +void +ChunkSparseMatrix::reinit (const ChunkSparsityPattern &sparsity) +{ + cols = &sparsity; + + if (cols->empty()) + { + if (val != 0) + delete[] val; + val = 0; + max_len = 0; + return; + } + + // allocate not just m() * n() elements but + // enough so that we can store full + // chunks. this entails some padding + // elements + const unsigned int chunk_size = cols->get_chunk_size(); + const unsigned int N = cols->sparsity_pattern.n_nonzero_elements() * + chunk_size * chunk_size; + if (N > max_len || max_len == 0) + { + if (val != 0) + delete[] val; + val = new number[N]; + max_len = N; + } + + // fill with zeros. do not just fill N + // elements but all that we allocated to + // ensure that also the padding elements + // are zero and not left at previous values + if (val != 0) + std::fill_n (&val[0], max_len, 0); +} + + + +template +void +ChunkSparseMatrix::clear () +{ + cols = 0; + if (val) delete[] val; + val = 0; + max_len = 0; +} + + + +template +bool +ChunkSparseMatrix::empty () const +{ + if (cols == 0) + return true; + else + return cols->empty(); +} + + + +template +unsigned int +ChunkSparseMatrix::n_nonzero_elements () const +{ + Assert (cols != 0, ExcNotInitialized()); + return cols->n_nonzero_elements (); +} + + + +template +unsigned int +ChunkSparseMatrix::n_actually_nonzero_elements () const +{ + Assert (cols != 0, ExcNotInitialized()); + + // count those elements that are nonzero, + // even if they lie in the padding around + // the matrix. since we have the invariant + // that padding elements are zero, nothing + // bad can happen here + const unsigned int chunk_size = cols->get_chunk_size(); + return std::count_if(&val[0], + &val[cols->sparsity_pattern.n_nonzero_elements () * + chunk_size * chunk_size], + std::bind2nd(std::not_equal_to(), 0)); +} + + + +template +void +ChunkSparseMatrix::symmetrize () +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (cols->rows == cols->cols, ExcNotQuadratic()); + + Assert (false, ExcNotImplemented()); +} + + + +template +template +ChunkSparseMatrix & +ChunkSparseMatrix::copy_from (const ChunkSparseMatrix &matrix) +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns()); + + // copy everything, including padding + // elements + const unsigned int chunk_size = cols->get_chunk_size(); + std::copy (&matrix.val[0], + &matrix.val[cols->sparsity_pattern.n_nonzero_elements() + * chunk_size * chunk_size], + &val[0]); + + return *this; +} + + + +template +template +void +ChunkSparseMatrix::copy_from (const FullMatrix &matrix) +{ + // first delete previous content + *this = 0; + + // then copy old matrix + for (unsigned int row=0; row +template +void +ChunkSparseMatrix::add (const number factor, + const ChunkSparseMatrix &matrix) +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns()); + + // add everything, including padding + // elements + const unsigned int chunk_size = cols->get_chunk_size(); + number *val_ptr = &val[0]; + const somenumber *matrix_ptr = &matrix.val[0]; + const number *const end_ptr = &val[cols->sparsity_pattern.n_nonzero_elements() + * chunk_size * chunk_size]; + + while (val_ptr != end_ptr) + *val_ptr++ += factor * *matrix_ptr++; +} + + +template +template +void +ChunkSparseMatrix::vmult (OutVector& dst, + const InVector& src) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); + + Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination()); + + // set the output vector to zero and then + // add to it the contributions of vmults + // from individual chunks. this is what + // vmult_add does + dst = 0; + vmult_add (dst, src); +} + + + +template +template +void +ChunkSparseMatrix::Tvmult (OutVector& dst, + const InVector& src) const +{ + Assert (val != 0, ExcNotInitialized()); + Assert (cols != 0, ExcNotInitialized()); + Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size())); + Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size())); + + Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination()); + + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); + + Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination()); + + // set the output vector to zero and then + // add to it the contributions of vmults + // from individual chunks. this is what + // vmult_add does + dst = 0; + Tvmult_add (dst, src); +} + + + +template +template +void +ChunkSparseMatrix::vmult_add (OutVector& dst, + const InVector& src) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); + + Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination()); + + const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); + + // loop over all chunks. note that we need + // to treat the last chunk row and column + // differently if they have padding + // elements + const bool rows_have_padding = (m() % cols->chunk_size != 0), + cols_have_padding = (n() % cols->chunk_size != 0); + + const unsigned int n_regular_chunk_rows + = (rows_have_padding ? + n_chunk_rows-1 : + n_chunk_rows); + + const number *val_ptr = val; + const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums; + typename OutVector::iterator dst_ptr = dst.begin(); + + for (unsigned int chunk_row=0; chunk_rowsparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + internal::ChunkSparseMatrix::chunk_vmult_add + (cols->chunk_size, + val_ptr, + src.begin() + *colnum_ptr * cols->chunk_size, + dst_ptr); + else + // we're at a chunk column that + // has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(chunk_row * cols->chunk_size + r) + += (val_ptr[r*cols->chunk_size + c] * + src(*colnum_ptr * cols->chunk_size + c)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + + + dst_ptr += cols->chunk_size; + } + + // now deal with last chunk row if + // necessary + if (rows_have_padding) + { + const unsigned int chunk_row = n_chunk_rows - 1; + + const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + { + // we're at a chunk row but not + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(chunk_row * cols->chunk_size + r) + += (val_ptr[r*cols->chunk_size + c] * + src(*colnum_ptr * cols->chunk_size + c)); + } + else + // we're at a chunk row and + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(chunk_row * cols->chunk_size + r) + += (val_ptr[r*cols->chunk_size + c] * + src(*colnum_ptr * cols->chunk_size + c)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + } +} + + +template +template +void +ChunkSparseMatrix::Tvmult_add (OutVector& dst, + const InVector& src) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); + + Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination()); + + const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); + + // loop over all chunks. note that we need + // to treat the last chunk row and column + // differently if they have padding + // elements + const bool rows_have_padding = (m() % cols->chunk_size != 0), + cols_have_padding = (n() % cols->chunk_size != 0); + + const unsigned int n_regular_chunk_rows + = (rows_have_padding ? + n_chunk_rows-1 : + n_chunk_rows); + + // like in vmult_add, but don't keep an + // iterator into dst around since we're not + // traversing it sequentially this time + const number *val_ptr = val; + const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums; + + for (unsigned int chunk_row=0; chunk_rowsparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + internal::ChunkSparseMatrix::chunk_Tvmult_add + (cols->chunk_size, + val_ptr, + src.begin() + chunk_row * cols->chunk_size, + dst.begin() + *colnum_ptr * cols->chunk_size); + else + // we're at a chunk column that + // has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(*colnum_ptr * cols->chunk_size + c) + += (val_ptr[r*cols->chunk_size + c] * + src(chunk_row * cols->chunk_size + r)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + } + + // now deal with last chunk row if + // necessary + if (rows_have_padding) + { + const unsigned int chunk_row = n_chunk_rows - 1; + + const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + { + // we're at a chunk row but not + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(*colnum_ptr * cols->chunk_size + c) + += (val_ptr[r*cols->chunk_size + c] * + src(chunk_row * cols->chunk_size + r)); + } + else + // we're at a chunk row and + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(*colnum_ptr * cols->chunk_size + c) + += (val_ptr[r*cols->chunk_size + c] * + src(chunk_row * cols->chunk_size + r)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + } +} + + +template +template +somenumber +ChunkSparseMatrix::matrix_norm_square (const Vector& v) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size())); + Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size())); + + somenumber result = 0; + + //////////////// + // like matrix_scalar_product, except that + // the two vectors are now the same + + const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); + + // loop over all chunks. note that we need + // to treat the last chunk row and column + // differently if they have padding + // elements + const bool rows_have_padding = (m() % cols->chunk_size != 0), + cols_have_padding = (n() % cols->chunk_size != 0); + + const unsigned int n_regular_chunk_rows + = (rows_have_padding ? + n_chunk_rows-1 : + n_chunk_rows); + + const number *val_ptr = val; + const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums; + typename Vector::const_iterator v_ptr = v.begin(); + + for (unsigned int chunk_row=0; chunk_rowsparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + result += + internal::ChunkSparseMatrix:: + chunk_matrix_scalar_product + (cols->chunk_size, + val_ptr, + v_ptr, + v.begin() + *colnum_ptr * cols->chunk_size); + else + // we're at a chunk column that + // has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + result + += + v(chunk_row * cols->chunk_size + r) + * (val_ptr[r*cols->chunk_size + c] * + v(*colnum_ptr * cols->chunk_size + c)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + + + v_ptr += cols->chunk_size; + } + + // now deal with last chunk row if + // necessary + if (rows_have_padding) + { + const unsigned int chunk_row = n_chunk_rows - 1; + + const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + { + // we're at a chunk row but not + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + result + += + v(chunk_row * cols->chunk_size + r) + * (val_ptr[r*cols->chunk_size + c] * + v(*colnum_ptr * cols->chunk_size + c)); + } + else + // we're at a chunk row and + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + result + += + v(chunk_row * cols->chunk_size + r) + * (val_ptr[r*cols->chunk_size + c] * + v(*colnum_ptr * cols->chunk_size + c)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + } + + return result; +} + + + +template +template +somenumber +ChunkSparseMatrix::matrix_scalar_product (const Vector& u, + const Vector& v) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert(m() == u.size(), ExcDimensionMismatch(m(),u.size())); + Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size())); + + // the following works like the vmult_add + // function + somenumber result = 0; + + const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); + + // loop over all chunks. note that we need + // to treat the last chunk row and column + // differently if they have padding + // elements + const bool rows_have_padding = (m() % cols->chunk_size != 0), + cols_have_padding = (n() % cols->chunk_size != 0); + + const unsigned int n_regular_chunk_rows + = (rows_have_padding ? + n_chunk_rows-1 : + n_chunk_rows); + + const number *val_ptr = val; + const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums; + typename Vector::const_iterator u_ptr = u.begin(); + + for (unsigned int chunk_row=0; chunk_rowsparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + result += + internal::ChunkSparseMatrix:: + chunk_matrix_scalar_product + (cols->chunk_size, + val_ptr, + u_ptr, + v.begin() + *colnum_ptr * cols->chunk_size); + else + // we're at a chunk column that + // has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + result + += + u(chunk_row * cols->chunk_size + r) + * (val_ptr[r*cols->chunk_size + c] * + v(*colnum_ptr * cols->chunk_size + c)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + + + u_ptr += cols->chunk_size; + } + + // now deal with last chunk row if + // necessary + if (rows_have_padding) + { + const unsigned int chunk_row = n_chunk_rows - 1; + + const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + { + // we're at a chunk row but not + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + result + += + u(chunk_row * cols->chunk_size + r) + * (val_ptr[r*cols->chunk_size + c] * + v(*colnum_ptr * cols->chunk_size + c)); + } + else + // we're at a chunk row and + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + result + += + u(chunk_row * cols->chunk_size + r) + * (val_ptr[r*cols->chunk_size + c] * + v(*colnum_ptr * cols->chunk_size + c)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + } + + return result; +} + + + +template +typename ChunkSparseMatrix::real_type +ChunkSparseMatrix::l1_norm () const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + + const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); + + // loop over all rows and columns; it is + // safe to also loop over the padding + // elements (they are zero) if we make sure + // that the vector into which we sum column + // sums is large enough + Vector column_sums(cols->sparsity_pattern.n_cols() * + cols->chunk_size); + + for (unsigned int chunk_row=0; chunk_rowsparsity_pattern.rowstart[chunk_row]; + jsparsity_pattern.rowstart[chunk_row+1] ; ++j) + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int s=0; schunk_size; ++s) + column_sums(cols->sparsity_pattern.colnums[j] * + cols->chunk_size + s) += + numbers::NumberTraits::abs(val[j * cols->chunk_size * + cols->chunk_size + + r * cols->chunk_size + + s]); + + return column_sums.linfty_norm(); +} + + + +template +typename ChunkSparseMatrix::real_type +ChunkSparseMatrix::linfty_norm () const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + + // this function works like l1_norm(). it + // can be made more efficient (without + // allocating a temporary vector) as is + // done in the SparseMatrix class but since + // it is rarely called in time critical + // places it is probably not worth it + const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); + + // loop over all rows and columns; it is + // safe to also loop over the padding + // elements (they are zero) if we make sure + // that the vector into which we sum column + // sums is large enough + Vector row_sums(cols->sparsity_pattern.n_rows() * + cols->chunk_size); + + for (unsigned int chunk_row=0; chunk_rowsparsity_pattern.rowstart[chunk_row]; + jsparsity_pattern.rowstart[chunk_row+1] ; ++j) + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int s=0; schunk_size; ++s) + row_sums(chunk_row * cols->chunk_size + r) += + numbers::NumberTraits::abs(val[j * cols->chunk_size * + cols->chunk_size + + r * cols->chunk_size + + s]); + + return row_sums.linfty_norm(); +} + + + +template +typename ChunkSparseMatrix::real_type +ChunkSparseMatrix::frobenius_norm () const +{ + // simply add up all entries in the + // sparsity pattern, without taking any + // reference to rows or columns + // + // padding elements are zero, so we can add + // them up as well + real_type norm_sqr = 0; + const unsigned int n_rows = m(); + for (const number *ptr = &val[0]; ptr != &val[max_len]; ++ptr) + norm_sqr += numbers::NumberTraits::abs_square(*ptr); + + return std::sqrt (norm_sqr); +} + + + +template +template +somenumber +ChunkSparseMatrix::residual (Vector &dst, + const Vector &u, + const Vector &b) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert(m() == b.size(), ExcDimensionMismatch(m(),b.size())); + Assert(n() == u.size(), ExcDimensionMismatch(n(),u.size())); + + Assert (&u != &dst, ExcSourceEqualsDestination()); + + // set dst=b, then subtract the result of + // A*u from it. since the purpose of the + // current class is to promote streaming of + // data rather than more random access + // patterns, breaking things up into two + // loops may be reasonable + dst = b; + + ///////// + // the rest of this function is like + // vmult_add, except that we subtract + // rather than add A*u + ///////// + const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); + + // loop over all chunks. note that we need + // to treat the last chunk row and column + // differently if they have padding + // elements + const bool rows_have_padding = (m() % cols->chunk_size != 0), + cols_have_padding = (n() % cols->chunk_size != 0); + + const unsigned int n_regular_chunk_rows + = (rows_have_padding ? + n_chunk_rows-1 : + n_chunk_rows); + + const number *val_ptr = val; + const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums; + typename Vector::iterator dst_ptr = dst.begin(); + + for (unsigned int chunk_row=0; chunk_rowsparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + internal::ChunkSparseMatrix::chunk_vmult_subtract + (cols->chunk_size, + val_ptr, + u.begin() + *colnum_ptr * cols->chunk_size, + dst_ptr); + else + // we're at a chunk column that + // has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(chunk_row * cols->chunk_size + r) + -= (val_ptr[r*cols->chunk_size + c] * + u(*colnum_ptr * cols->chunk_size + c)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + + + dst_ptr += cols->chunk_size; + } + + // now deal with last chunk row if + // necessary + if (rows_have_padding) + { + const unsigned int chunk_row = n_chunk_rows - 1; + + const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1] + * cols->chunk_size + * cols->chunk_size]; + while (val_ptr != val_end_of_row) + { + if ((cols_have_padding == false) + || + (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) + { + // we're at a chunk row but not + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(chunk_row * cols->chunk_size + r) + -= (val_ptr[r*cols->chunk_size + c] * + u(*colnum_ptr * cols->chunk_size + c)); + } + else + // we're at a chunk row and + // column that has padding + for (unsigned int r=0; rchunk_size; ++r) + for (unsigned int c=0; cchunk_size; ++c) + dst(chunk_row * cols->chunk_size + r) + -= (val_ptr[r*cols->chunk_size + c] * + u(*colnum_ptr * cols->chunk_size + c)); + + ++colnum_ptr; + val_ptr += cols->chunk_size * cols->chunk_size; + } + + + dst_ptr += cols->chunk_size; + } + + // finally compute the norm + return dst.l2_norm(); +} + + + +template +template +void +ChunkSparseMatrix::precondition_Jacobi (Vector &dst, + const Vector &src, + const number /*om*/) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n())); + Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n())); + + Assert (false, ExcNotImplemented()); +} + + + +template +template +void +ChunkSparseMatrix::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, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n())); + Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n())); + + Assert (false, ExcNotImplemented()); +} + + +template +template +void +ChunkSparseMatrix::precondition_SOR (Vector& dst, + const Vector& src, + const number om) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + + dst = src; + SOR(dst,om); +} + + +template +template +void +ChunkSparseMatrix::precondition_TSOR (Vector& dst, + const Vector& src, + const number om) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + + dst = src; + TSOR(dst,om); +} + + +template +template +void +ChunkSparseMatrix::SOR (Vector& dst, + const number /*om*/) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + + Assert (false, ExcNotImplemented()); +} + + +template +template +void +ChunkSparseMatrix::TSOR (Vector& dst, + const number /*om*/) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + + Assert (false, ExcNotImplemented()); +} + + +template +template +void +ChunkSparseMatrix::PSOR (Vector& dst, + const std::vector& permutation, + const std::vector& inverse_permutation, + const number /*om*/) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert (m() == permutation.size(), + ExcDimensionMismatch(m(), permutation.size())); + Assert (m() == inverse_permutation.size(), + ExcDimensionMismatch(m(), inverse_permutation.size())); + + Assert (false, ExcNotImplemented()); +} + + +template +template +void +ChunkSparseMatrix::TPSOR (Vector& dst, + const std::vector& permutation, + const std::vector& inverse_permutation, + const number /*om*/) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + Assert (m() == permutation.size(), + ExcDimensionMismatch(m(), permutation.size())); + Assert (m() == inverse_permutation.size(), + ExcDimensionMismatch(m(), inverse_permutation.size())); + + Assert (false, ExcNotImplemented()); +} + + + +template +template +void +ChunkSparseMatrix::SOR_step (Vector &v, + const Vector &b, + const number /*om*/) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size())); + Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size())); + + Assert (false, ExcNotImplemented()); +} + + + +template +template +void +ChunkSparseMatrix::TSOR_step (Vector &v, + const Vector &b, + const number /*om*/) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size())); + Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size())); + + Assert (false, ExcNotImplemented()); +} + + + +template +template +void +ChunkSparseMatrix::SSOR_step (Vector &v, + const Vector &b, + const number om) const +{ + SOR_step(v,b,om); + TSOR_step(v,b,om); +} + + + +template +template +void +ChunkSparseMatrix::SSOR (Vector& dst, + const number /*om*/) const +{ + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + Assert (cols->optimize_diagonal(), + typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + + Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + + Assert (false, ExcNotImplemented()); +} + + + +template +const ChunkSparsityPattern & +ChunkSparseMatrix::get_sparsity_pattern () const +{ + Assert (cols != 0, ExcNotInitialized()); + return *cols; +} + + + +template +void ChunkSparseMatrix::print (std::ostream &out) const +{ + AssertThrow (out, ExcIO()); + + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + + Assert (false, ExcNotImplemented()); + + AssertThrow (out, ExcIO()); +} + + +template +void ChunkSparseMatrix::print_formatted (std::ostream &out, + const unsigned int precision, + const bool scientific, + const unsigned int width_, + const char* zero_string, + const double denominator) const +{ + AssertThrow (out, ExcIO()); + + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + + unsigned int width = width_; + + Assert (false, ExcNotImplemented()); + + std::ios::fmtflags old_flags = out.flags(); + unsigned int old_precision = out.precision (precision); + + if (scientific) + { + out.setf (std::ios::scientific, std::ios::floatfield); + if (!width) + width = precision+7; + } else { + out.setf (std::ios::fixed, std::ios::floatfield); + if (!width) + width = precision+2; + } + + for (unsigned int i=0; isparsity_pattern(i,j) != SparsityPattern::invalid_entry) + out << std::setw(width) + << val[cols->sparsity_pattern(i,j)] * denominator << ' '; + else + out << std::setw(width) << zero_string << ' '; + out << std::endl; + }; + AssertThrow (out, ExcIO()); + + // reset output format + out.precision(old_precision); + out.flags (old_flags); +} + + + +template +void ChunkSparseMatrix::print_pattern (std::ostream &out, + const double threshold) const +{ + AssertThrow (out, ExcIO()); + + Assert (cols != 0, ExcNotInitialized()); + Assert (val != 0, ExcNotInitialized()); + + const unsigned int chunk_size = cols->get_chunk_size(); + + // loop over all chunk rows and columns, + // and each time we find something repeat + // it chunk_size times in both directions + for (unsigned int i=0; isparsity_pattern.n_rows(); ++i) + { + for (unsigned int d=0; dsparsity_pattern.n_cols(); ++j) + if (cols->sparsity_pattern(i,j) == SparsityPattern::invalid_entry) + { + for (unsigned int e=0; esparsity_pattern(i,j)]) > threshold) + { + for (unsigned int e=0; e +void +ChunkSparseMatrix::block_write (std::ostream &out) const +{ + AssertThrow (out, ExcIO()); + + // first the simple objects, + // bracketed in [...] + out << '[' << max_len << "]["; + // then write out real data + out.write (reinterpret_cast(&val[0]), + reinterpret_cast(&val[max_len]) + - reinterpret_cast(&val[0])); + out << ']'; + + AssertThrow (out, ExcIO()); +} + + + +template +void +ChunkSparseMatrix::block_read (std::istream &in) +{ + AssertThrow (in, ExcIO()); + + char c; + + // first read in simple data + in >> c; + AssertThrow (c == '[', ExcIO()); + in >> max_len; + + in >> c; + AssertThrow (c == ']', ExcIO()); + in >> c; + AssertThrow (c == '[', ExcIO()); + + // reallocate space + delete[] val; + val = new number[max_len]; + + // then read data + in.read (reinterpret_cast(&val[0]), + reinterpret_cast(&val[max_len]) + - reinterpret_cast(&val[0])); + in >> c; + AssertThrow (c == ']', ExcIO()); +} + + + +template +unsigned int +ChunkSparseMatrix::memory_consumption () const +{ + return sizeof(*this) + max_len*sizeof(number); +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/lac/source/chunk_sparse_matrix.cc b/deal.II/lac/source/chunk_sparse_matrix.cc new file mode 100644 index 0000000000..cc2825941d --- /dev/null +++ b/deal.II/lac/source/chunk_sparse_matrix.cc @@ -0,0 +1,20 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2008 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. +// +//--------------------------------------------------------------------------- + + +#include +#include + +DEAL_II_NAMESPACE_OPEN +#include "chunk_sparse_matrix.inst" +DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/lac/source/chunk_sparse_matrix.inst.in b/deal.II/lac/source/chunk_sparse_matrix.inst.in new file mode 100644 index 0000000000..f592b880bb --- /dev/null +++ b/deal.II/lac/source/chunk_sparse_matrix.inst.in @@ -0,0 +1,225 @@ +//---------------------------- sparse_matrix.inst.in --------------------------- +// $Id: sparse_matrix_matrix.in.h 15011 2007-08-22 16:59:41Z kanschat $ +// Version: $Name$ +// +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2008 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.inst.in --------------------------- + + +// real instantiations + +for (S : REAL_SCALARS) + { + template class ChunkSparseMatrix; + } + + + +for (S1, S2 : REAL_SCALARS) + { + template ChunkSparseMatrix & + ChunkSparseMatrix::copy_from (const ChunkSparseMatrix &); + + template + void ChunkSparseMatrix::copy_from (const FullMatrix &); + + template void ChunkSparseMatrix::add (const S1, + const ChunkSparseMatrix &); + } + + +for (S1, S2 : REAL_SCALARS) + { + template S2 + ChunkSparseMatrix:: + matrix_norm_square (const Vector &) const; + + template S2 + ChunkSparseMatrix:: + matrix_scalar_product (const Vector &, + const Vector &) const; + + template S2 ChunkSparseMatrix:: + residual (Vector &, + const Vector &, + const Vector &) const; + + template void ChunkSparseMatrix:: + precondition_SSOR (Vector &, + const Vector &, + const S1) const; + + template void ChunkSparseMatrix:: + precondition_SOR (Vector &, + const Vector &, + const S1) const; + + template void ChunkSparseMatrix:: + precondition_TSOR (Vector &, + const Vector &, + const S1) const; + + template void ChunkSparseMatrix:: + precondition_Jacobi (Vector &, + const Vector &, + const S1) const; + + template void ChunkSparseMatrix:: + SOR (Vector &, + const S1) const; + template void ChunkSparseMatrix:: + TSOR (Vector &, + const S1) const; + template void ChunkSparseMatrix:: + SSOR (Vector &, + const S1) const; + template void ChunkSparseMatrix:: + PSOR (Vector &, + const std::vector&, + const std::vector&, + const S1) const; + template void ChunkSparseMatrix:: + TPSOR (Vector &, + const std::vector&, + const std::vector&, + const S1) const; + template void ChunkSparseMatrix:: + SOR_step (Vector &, + const Vector &, + const S1) const; + template void ChunkSparseMatrix:: + TSOR_step (Vector &, + const Vector &, + const S1) const; + template void ChunkSparseMatrix:: + SSOR_step (Vector &, + const Vector &, + const S1) const; + } + + +for (S1, S2, S3 : REAL_SCALARS; + V1, V2 : DEAL_II_VEC_TEMPLATES) + { + template void ChunkSparseMatrix:: + vmult (V1 &, const V2 &) const; + template void ChunkSparseMatrix:: + Tvmult (V1 &, const V2 &) const; + template void ChunkSparseMatrix:: + vmult_add (V1 &, const V2 &) const; + template void ChunkSparseMatrix:: + Tvmult_add (V1 &, const V2 &) const; + } + + + +// complex instantiations + +// for (S : COMPLEX_SCALARS) +// { +// template class ChunkSparseMatrix; +// } + + + +// for (S1, S2 : COMPLEX_SCALARS) +// { +// template ChunkSparseMatrix & +// ChunkSparseMatrix::copy_from (const ChunkSparseMatrix &); + +// template +// void ChunkSparseMatrix::copy_from (const FullMatrix &); + +// template void ChunkSparseMatrix::add (const S1, +// const ChunkSparseMatrix &); +// } + + +// for (S1, S2 : COMPLEX_SCALARS) +// { +// template S2 +// ChunkSparseMatrix:: +// matrix_norm_square (const Vector &) const; + +// template S2 +// ChunkSparseMatrix:: +// matrix_scalar_product (const Vector &, +// const Vector &) const; + +// template S2 ChunkSparseMatrix:: +// residual (Vector &, +// const Vector &, +// const Vector &) const; + +// template void ChunkSparseMatrix:: +// precondition_SSOR (Vector &, +// const Vector &, +// const S1) const; + +// template void ChunkSparseMatrix:: +// precondition_SOR (Vector &, +// const Vector &, +// const S1) const; + +// template void ChunkSparseMatrix:: +// precondition_TSOR (Vector &, +// const Vector &, +// const S1) const; + +// template void ChunkSparseMatrix:: +// precondition_Jacobi (Vector &, +// const Vector &, +// const S1) const; + +// template void ChunkSparseMatrix:: +// SOR (Vector &, +// const S1) const; +// template void ChunkSparseMatrix:: +// TSOR (Vector &, +// const S1) const; +// template void ChunkSparseMatrix:: +// SSOR (Vector &, +// const S1) const; +// template void ChunkSparseMatrix:: +// PSOR (Vector &, +// const std::vector&, +// const std::vector&, +// const S1) const; +// template void ChunkSparseMatrix:: +// TPSOR (Vector &, +// const std::vector&, +// const std::vector&, +// const S1) const; +// template void ChunkSparseMatrix:: +// SOR_step (Vector &, +// const Vector &, +// const S1) const; +// template void ChunkSparseMatrix:: +// TSOR_step (Vector &, +// const Vector &, +// const S1) const; +// template void ChunkSparseMatrix:: +// SSOR_step (Vector &, +// const Vector &, +// const S1) const; +// } + + +// for (S1, S2, S3 : COMPLEX_SCALARS; +// V1, V2 : DEAL_II_VEC_TEMPLATES) +// { +// template void ChunkSparseMatrix:: +// vmult (V1 &, const V2 &) const; +// template void ChunkSparseMatrix:: +// Tvmult (V1 &, const V2 &) const; +// template void ChunkSparseMatrix:: +// vmult_add (V1 &, const V2 &) const; +// template void ChunkSparseMatrix:: +// Tvmult_add (V1 &, const V2 &) const; +// } diff --git a/tests/bits/Makefile b/tests/bits/Makefile index 77f9e512ac..d35b167863 100644 --- a/tests/bits/Makefile +++ b/tests/bits/Makefile @@ -1,6 +1,6 @@ ############################################################ # Makefile,v 1.15 2002/06/13 12:51:13 hartmann Exp -# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors ############################################################ ############################################################ @@ -63,6 +63,7 @@ tests_x = grid_generator_?? \ find_cell_* \ sparsity_pattern_* \ sparse_matrix_* \ + chunk_sparse_matrix_* \ full_matrix_vector_* \ solver_* \ deal_solver_* \ diff --git a/tests/bits/chunk_sparse_matrix_00.cc b/tests/bits/chunk_sparse_matrix_00.cc new file mode 100644 index 0000000000..ade67fa117 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_00.cc @@ -0,0 +1,95 @@ +//---------------------------- chunk_sparse_matrix_00.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2008 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. +// +//---------------------------- chunk_sparse_matrix_00.cc --------------------------- + + +// set a few elements in a chunk sparse matrix and output them again. should +// yield the same result for all chunk sizes, of course + +#include "../tests.h" +#include +#include +#include + + +void test (const unsigned int chunk_size) +{ + deallog << "Chunk size = " << chunk_size << std::endl; + + ChunkSparsityPattern sp (5,5,3,chunk_size,false); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i +#include +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i +#include +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size,false); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i1, then this isn't + // necessarily true + bool exc_thrown = false; + double d; + try + { + d = m(i,j); + } + catch (const std::exception &) + { + exc_thrown = true; + } + Assert ((exc_thrown == true) || (chunk_size > 1), + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_01a/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 }; + for (unsigned int i=0; + i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + Assert (m.m() == 5, ExcInternalError()); + Assert (m.n() == 5, ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_04/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 }; + for (unsigned int i=0; + i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries. count how many + // entries we have. note that for square + // matrices we also always store the + // diagonal element, so add one per row, + // but don't count it when traversing the + // row + unsigned int counter = 0; + for (unsigned int i=0; i= counter, + ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_05/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 }; + for (unsigned int i=0; + i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size,false); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries. count how many + // entries we have. note that for square + // matrices we also always store the + // diagonal element, except when as above + // we set the special flag for the matrix + // sparsity pattern + unsigned int counter = 0; + for (unsigned int i=0; i= counter, + ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_05a/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 }; + for (unsigned int i=0; + i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries. count how many + // entries we have + for (unsigned int i=0; i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries. count how many + // entries we have + for (unsigned int i=0; i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries. count how many + // entries we have + double norm = 0; + for (unsigned int i=0; i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i +#include + + +void test (const unsigned int chunk_size) +{ + ChunkSparsityPattern sp (5,5,3,chunk_size); + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if ((i+2*j+1) % 3 == 0) + sp.add (i,j); + sp.compress (); + + ChunkSparseMatrix m(sp); + + // first set a few entries + for (unsigned int i=0; i +#include +#include +#include +#include + + +void test (const unsigned int chunk_size, + Vector &v, + Vector &w) +{ + // set some entries in the + // matrix. actually, set them all + ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size); + for (unsigned int i=0; i m(sp); + for (unsigned int i=0; i v (100); + Vector w (100); + test (chunk_sizes[i], v,w); + } + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_vector_01/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_01/cmp/generic new file mode 100644 index 0000000000..c9bdf33b16 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_01/cmp/generic @@ -0,0 +1,6 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_vector_02.cc b/tests/bits/chunk_sparse_matrix_vector_02.cc new file mode 100644 index 0000000000..318b33d289 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_02.cc @@ -0,0 +1,110 @@ +//---------------------------- chunk_sparse_matrix_vector_02.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2008 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. +// +//---------------------------- chunk_sparse_matrix_vector_02.cc --------------------------- + + +// check ChunkSparseMatrix::Tvmult + +#include "../tests.h" +#include +#include +#include +#include +#include + + +void test (const unsigned int chunk_size, + Vector &v, + Vector &w) +{ + // set some entries in the + // matrix. actually, set them all + ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size); + for (unsigned int i=0; i m(sp); + for (unsigned int i=0; i v (100); + Vector w (100); + test (chunk_sizes[i], v,w); + } + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_vector_02/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_02/cmp/generic new file mode 100644 index 0000000000..c9bdf33b16 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_02/cmp/generic @@ -0,0 +1,6 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_vector_03.cc b/tests/bits/chunk_sparse_matrix_vector_03.cc new file mode 100644 index 0000000000..59e19be090 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_03.cc @@ -0,0 +1,113 @@ +//---------------------------- chunk_sparse_matrix_vector_03.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2008 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. +// +//---------------------------- chunk_sparse_matrix_vector_03.cc --------------------------- + + +// check ChunkSparseMatrix::vmult_add + +#include "../tests.h" +#include +#include +#include +#include +#include + + +void test (const unsigned int chunk_size, + Vector &v, + Vector &w) +{ + // set some entries in the + // matrix. actually, set them all + ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size); + for (unsigned int i=0; i m(sp); + for (unsigned int i=0; i v (100); + Vector w (100); + test (chunk_sizes[i], v,w); + } + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_vector_03/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_03/cmp/generic new file mode 100644 index 0000000000..c9bdf33b16 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_03/cmp/generic @@ -0,0 +1,6 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_vector_04.cc b/tests/bits/chunk_sparse_matrix_vector_04.cc new file mode 100644 index 0000000000..9f85cb0977 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_04.cc @@ -0,0 +1,113 @@ +//---------------------------- chunk_sparse_matrix_vector_04.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2008 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. +// +//---------------------------- chunk_sparse_matrix_vector_04.cc --------------------------- + + +// check ChunkSparseMatrix::Tvmult_add + +#include "../tests.h" +#include +#include +#include +#include +#include + + +void test (const unsigned int chunk_size, + Vector &v, + Vector &w) +{ + // set some entries in the + // matrix. actually, set them all + ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size); + for (unsigned int i=0; i m(sp); + for (unsigned int i=0; i v (100); + Vector w (100); + test (chunk_sizes[i], v,w); + } + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_vector_04/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_04/cmp/generic new file mode 100644 index 0000000000..c9bdf33b16 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_04/cmp/generic @@ -0,0 +1,6 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_vector_05.cc b/tests/bits/chunk_sparse_matrix_vector_05.cc new file mode 100644 index 0000000000..4f2efc4205 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_05.cc @@ -0,0 +1,116 @@ +//---------------------------- chunk_sparse_matrix_vector_05.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2008 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. +// +//---------------------------- chunk_sparse_matrix_vector_05.cc --------------------------- + + +// check ChunkSparseMatrix::matrix_scalar_product + +#include "../tests.h" +#include +#include +#include +#include +#include + + +void test (const unsigned int chunk_size, + Vector &v, + Vector &w) +{ + // set some entries in the + // matrix. actually, set them all + ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size); + for (unsigned int i=0; i m(sp); + for (unsigned int i=0; i + const double s = m.matrix_scalar_product (w,v); + + // make sure we get the expected result + for (unsigned int i=0; i v (100); + Vector w (100); + test (chunk_sizes[i],v,w); + } + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_vector_05/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_05/cmp/generic new file mode 100644 index 0000000000..c9bdf33b16 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_05/cmp/generic @@ -0,0 +1,6 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_vector_06.cc b/tests/bits/chunk_sparse_matrix_vector_06.cc new file mode 100644 index 0000000000..493dd4537c --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_06.cc @@ -0,0 +1,107 @@ +//---------------------------- chunk_sparse_matrix_vector_06.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2008 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. +// +//---------------------------- chunk_sparse_matrix_vector_06.cc --------------------------- + + +// check ChunkSparseMatrix::matrix_norm_square + +#include "../tests.h" +#include +#include +#include +#include +#include + + +void test (const unsigned int chunk_size, + Vector &v) +{ + // set some entries in the + // matrix. actually, set them all + ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size); + for (unsigned int i=0; i m(sp); + for (unsigned int i=0; i + const double s = m.matrix_norm_square (v); + + // make sure we get the expected result + for (unsigned int i=0; i v (100); + test (chunk_sizes[i], v); + } + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_vector_06/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_06/cmp/generic new file mode 100644 index 0000000000..c9bdf33b16 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_06/cmp/generic @@ -0,0 +1,6 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_vector_07.cc b/tests/bits/chunk_sparse_matrix_vector_07.cc new file mode 100644 index 0000000000..a6f536d21d --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_07.cc @@ -0,0 +1,119 @@ +//---------------------------- chunk_sparse_matrix_vector_07.cc --------------------------- +// chunk_sparse_matrix_vector_07.cc,v 1.5 2004/02/26 17:25:45 wolf Exp +// Version: +// +// Copyright (C) 2004, 2005, 2008 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. +// +//---------------------------- chunk_sparse_matrix_vector_07.cc --------------------------- + + +// check ChunkSparseMatrix::residual + +#include "../tests.h" +#include +#include +#include +#include +#include + + +void test (const unsigned int chunk_size, + Vector &v, + Vector &w, + Vector &x) +{ + // set some entries in the + // matrix. actually, set them all + ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size); + for (unsigned int i=0; i m(sp); + for (unsigned int i=0; i v (100); + Vector w (100); + Vector x (100); + test (chunk_sizes[i],v,w,x); + } + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_vector_07/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_07/cmp/generic new file mode 100644 index 0000000000..c9bdf33b16 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_vector_07/cmp/generic @@ -0,0 +1,6 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK