From 0dde13a5950412e4bee5e48e7af1ad530115b691 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Sat, 12 Mar 2005 03:26:12 +0000 Subject: [PATCH] beginning TridiagonalMatrix git-svn-id: https://svn.dealii.org/trunk@10109 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/tridiagonal_matrix.h | 405 +++++++++++++++++++ deal.II/lac/source/tridiagonal_matrix.cc | 145 +++++++ tests/lac/Makefile | 3 +- tests/lac/tridiagonal_matrix.cc | 72 ++++ 4 files changed, 624 insertions(+), 1 deletion(-) create mode 100644 deal.II/lac/include/lac/tridiagonal_matrix.h create mode 100644 deal.II/lac/source/tridiagonal_matrix.cc create mode 100644 tests/lac/tridiagonal_matrix.cc diff --git a/deal.II/lac/include/lac/tridiagonal_matrix.h b/deal.II/lac/include/lac/tridiagonal_matrix.h new file mode 100644 index 0000000000..64a7fd9c75 --- /dev/null +++ b/deal.II/lac/include/lac/tridiagonal_matrix.h @@ -0,0 +1,405 @@ +//------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 1998 - 2005 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__tridiagonal_matrix_h +#define __deal2__tridiagonal_matrix_h + +#include +#include +#include + +#include + +// forward declarations +template class Vector; + + +/*! @addtogroup Matrix1 + *@{ + */ + + +/** + * A quadratic tridiagonal matrix. That is, a matrix where all entries + * are zero, except the diagonal and the entries left and right of it. + * + * @note Only data management and entry functions are implemented + * directly. All more complex functions require LAPACK support. + */ +template +class TridiagonalMatrix +{ + public: +/** + * @name Constructors and initalization. + */ + /** + * Constructor generating an + * empty matrix of dimension + * n. + */ + TridiagonalMatrix(unsigned int n = 0); + +//@} +///@name Non-modifying operators +//@{ + + /** + * Number of rows of this matrix. + * To remember: this matrix is an + * m x m-matrix. + */ + unsigned int m () const; + + /** + * Number of columns of this matrix. + * To remember: this matrix is an + * n x n-matrix. + */ + unsigned int n () const; + + /** + * Return whether the matrix + * contains only elements with + * value zero. This function is + * mainly for internal + * consistency checks and should + * seldomly be used when not in + * debug mode since it uses quite + * some time. + */ + bool all_zero () const; + + + +//@} +///@name Element access +//@{ + /** + * Read-only access to a + * value. This is restricted to + * the case where |i-j| <= + * 1. + */ + number operator()(unsigned int i, unsigned int j) const; + + /** + * Read-write access to a + * value. This is restricted to + * the case where |i-j| <= + * 1. + */ + number& operator()(unsigned int i, unsigned int j); + +//@} +///@name Multiplications with vectors +//@{ + + /** + * Matrix-vector-multiplication. Multiplies + * v from the right and + * stores the result in + * w. + * + * If the optional parameter + * adding is true, the + * result is added to w. + * + * Source and destination must + * not be the same vector. + */ + void vmult (Vector &w, + const Vector &v, + const bool adding=false) const; + + /** + * Adding + * Matrix-vector-multiplication. Same + * as vmult() with parameter + * adding=true, but + * widely used in + * deal.II classes. + * + * Source and destination must + * not be the same vector. + */ + void vmult_add (Vector &w, + const Vector &v) const; + + /** + * Transpose + * matrix-vector-multiplication. + * Multiplies + * vT from + * the left and stores the result + * in w. + * + * If the optional parameter + * adding is true, the + * result is added to w. + * + * Source and destination must + * not be the same vector. + */ + void Tvmult (Vector &w, + const Vector &v, + const bool adding=false) const; + + /** + * Adding transpose + * matrix-vector-multiplication. Same + * as Tvmult() with parameter + * adding=true, but + * widely used in + * deal.II classes. + * + * Source and destination must + * not be the same vector. + */ + void Tvmult_add (Vector &w, + const Vector &v) const; + + /** + * Build the matrix scalar product + * u^T M v. This function is mostly + * useful when building the cellwise + * scalar product of two functions in + * the finite element context. + */ + number matrix_scalar_product (const Vector &u, + const Vector &v) const; + + /** + * Return the square of the norm + * of the vector v with + * respect to the norm induced by + * this matrix, + * i.e. (v,Mv). This is + * useful, e.g. in the finite + * element context, where the + * L2 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. + */ + number matrix_norm_square (const Vector &v) const; + +//@} +///@name Matrixnorms +//@{ + + /** + * Return the $l_1$-norm of the matrix, i.e. + * $|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 $l_1$-norm for vectors, i.e. + * $|Mv|_1\leq |M|_1 |v|_1$. + * (cf. Rannacher Numerik0) + */ + number l1_norm () const; + + /** + * Return the $l_\infty$-norm of the + * matrix, i.e. + * $|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 $l_\infty$-norm of vectors, i.e. + * $|Mv|_\infty \leq |M|_\infty |v|_\infty$. + */ + number linfty_norm () const; + + /** + * The Frobenius norm of the matrix. + * Return value is the root of the square + * sum of all matrix entries. + */ + number frobenius_norm () const; + + /** + * Compute the relative norm of + * the skew-symmetric part. The + * return value is the Frobenius + * norm of the skew-symmetric + * part of the matrix divided by + * that of the matrix. + * + * Main purpose of this function + * is to check, if a matrix is + * symmetric within a certain + * accuracy, or not. + */ + number relative_symmetry_norm2 () const; +//@} +///@name Miscellanea +//@{ + /** + * Output of the matrix in + * user-defined format. + */ + void print (std::ostream &s, + const unsigned int width=5, + const unsigned int precision=2) 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. + */ + void print_formatted (std::ostream &out, + const unsigned int presicion=3, + const bool scientific = true, + const unsigned int width = 0, + const char *zero_string = " ", + const double denominator = 1.) const; + + /** + * Determine an estimate for the + * memory consumption (in bytes) + * of this object. + */ + unsigned int memory_consumption () const; +//@} + + private: + /** + * The diagonal entries. + */ + std::vector diagonal; + /** + * The entries left of the + * diagonal. The entry with index + * zero is always zero, since the + * first row has no entry left of + * the diagonal. Therefore, the + * length of this vector is the + * same as that of #diagonal. + */ + std::vector left; + /** + * The entries right of the + * diagonal. The last entry is + * always zero, since the last + * row has no entry right of the + * diagonal. Therefore, the + * length of this vector is the + * same as that of #diagonal. + */ + std::vector right; +}; + +//----------------------------------------------------------------------// +///@if NoDoc + +template +unsigned int +TridiagonalMatrix::m() const +{ + return diagonal.size(); +} + + + +template +unsigned int +TridiagonalMatrix::n() const +{ + return diagonal.size(); +} + + +template +inline +number +TridiagonalMatrix::operator()(unsigned int i, unsigned int j) const +{ + Assert(i +inline +number& +TridiagonalMatrix::operator()(unsigned int i, unsigned int j) +{ + Assert(i +#include + +template +TridiagonalMatrix::TridiagonalMatrix(unsigned int size) + : + diagonal(size, 0.), + left(size, 0.), + right(size, 0.) +{} + + +template +bool +TridiagonalMatrix::all_zero() const +{ + typename std::vector::const_iterator i; + typename std::vector::const_iterator e; + + e = diagonal.end(); + for (i=diagonal.begin() ; i != e ; ++i) + if (*i != 0.) return false; + + e = left.end(); + for (i=left.begin() ; i != e ; ++i) + if (*i != 0.) return false; + + e = right.end(); + for (i=right.begin() ; i != e ; ++i) + if (*i != 0.) return false; + return true; +} + + +template +void +TridiagonalMatrix::vmult ( + Vector &w, + const Vector &v, + const bool adding) const +{ + Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n())); + Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n())); + + if (n()==0) return; + + if (adding) + { + w(0) += diagonal[0]*v(0) + right[0]*v(1); + const unsigned int e=n()-1; + for (unsigned int i=1;i +void +TridiagonalMatrix::vmult_add ( + Vector &w, + const Vector &v) const +{ + vmult(w, v, true); +} + + +template +void +TridiagonalMatrix::Tvmult ( + Vector &w, + const Vector &v, + const bool adding) const +{ + Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n())); + Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n())); + + if (n()==0) return; +//TODO:[GK] Check this!!! + if (adding) + { + w(0) += diagonal[0]*v(0) + left[1]*v(1); + const unsigned int e=n()-1; + for (unsigned int i=1;i +void +TridiagonalMatrix::Tvmult_add ( + Vector &w, + const Vector &v) const +{ + Tvmult(w, v, true); +} + +/* +template +TridiagonalMatrix:: +{ +} + + +template +TridiagonalMatrix:: +{ +} + + +*/ + +template TridiagonalMatrix; +template TridiagonalMatrix; diff --git a/tests/lac/Makefile b/tests/lac/Makefile index 533f165581..6884ed1878 100644 --- a/tests/lac/Makefile +++ b/tests/lac/Makefile @@ -21,7 +21,8 @@ default: run-tests ############################################################ tests_x = vector-vector \ - full_matrix sparsity_pattern sparse_matrices \ + full_matrix tridiagonal_matrix \ + sparsity_pattern sparse_matrices \ block_vector block_vector_iterator block_matrices \ matrix_lib matrix_out \ solver eigen \ diff --git a/tests/lac/tridiagonal_matrix.cc b/tests/lac/tridiagonal_matrix.cc new file mode 100644 index 0000000000..4a28aac561 --- /dev/null +++ b/tests/lac/tridiagonal_matrix.cc @@ -0,0 +1,72 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005 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 "../tests.h" +#include +#include +#include + +#include +#include +#include + +template +void +check_vmult() +{ + TridiagonalMatrix M(4); + Vector u(4); + Vector v(4); + + for (unsigned int i=0;i<4;++i) + { + u(i) = i+1; + M(i,i) = i+1; + if (i>0) + M(i,i-1) = 0.-i; + if (i<3) + M(i,i+1) = 4.-i; + } + M.vmult(v,u); + for (unsigned int i=0;i(); +} -- 2.39.5