From 1e0f838ea81d782dfa06abfb3ae015b525fe5fbe Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 24 Apr 2006 20:24:27 +0000 Subject: [PATCH] Add identity matrix class. git-svn-id: https://svn.dealii.org/trunk@12867 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/identity_matrix.h | 275 ++++++++++++++++++++++ 1 file changed, 275 insertions(+) create mode 100644 deal.II/lac/include/lac/identity_matrix.h diff --git a/deal.II/lac/include/lac/identity_matrix.h b/deal.II/lac/include/lac/identity_matrix.h new file mode 100644 index 0000000000..a244633cd4 --- /dev/null +++ b/deal.II/lac/include/lac/identity_matrix.h @@ -0,0 +1,275 @@ +//--------------------------------------------------------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2006 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__identity_matrix_h +#define __deal2__identity_matrix_h + + +#include +#include + + +/*! @addtogroup Matrix1 + *@{ + */ + + +/** + * Implementation of a simple class representing the identity matrix + * of a given size, i.e. a matrix with entries + * $A_{ij}=\delta_{ij}$. While it has the most important ingredients + * of a matrix, in particular that one can ask for its size and + * perform matrix-vector products with it, a matrix of this type is + * really only useful in two contexts: preconditioning and + * initializing other matrices. + * + + *

Initialization

+ * + * The main usefulness of this class lies in its ability to initialize + * other matrix, like this: + @verbatim + FullMatrix identity (IdentityMatrix(10)); + @endverbatim + * + * This creates a $10\times 10$ matrix with ones on the diagonal and + * zeros everywhere else. All matrix types have conversion + * constructors and can therefore be filled rather easily with + * identity matrices. + * + * + *

Preconditioning

+ * + * No preconditioning at all is equivalent to preconditioning with + * preconditioning with the identity matrix. deal.II has a specialized + * class for this purpose, PreconditionIdentity, than can be used in a + * context as shown in the documentation of that class. The present + * class can be used in much the same way, although without any + * additional benefit: + @verbatim + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + cg.solve (system_matrix, solution, system_rhs, + IdentityMatrix(solution.size())); + @endverbatim + * + * + * @author Wolfgang Bangerth, 2006 + */ +class IdentityMatrix +{ + public: + /** + * Default constructor. Creates a + * zero-sized matrix that should + * be resized later on using the + * reinit() function. + */ + IdentityMatrix (); + + /** + * Constructor. Creates a + * identity matrix of size #n. + */ + IdentityMatrix (const unsigned int n); + + /** + * Resize the matrix to be of + * size #n by #n. + */ + void reinit (const unsigned int n); + + /** + * Number of rows of this + * matrix. For the present + * matrix, the number of rows and + * columns are equal, of course. + */ + unsigned int m () const; + + /** + * Number of columns of this + * matrix. For the present + * matrix, the number of rows and + * columns are equal, of course. + */ + unsigned int n () const; + + /** + * Matrix-vector + * multiplication. For the + * present case, this of course + * amounts to simply copying the + * input vector to the output + * vector. + */ + template + void vmult (VECTOR1 &out, + const VECTOR2 &in) const; + + /** + * Matrix-vector multiplication + * with addition to the output + * vector. For the present case, + * this of course amounts to + * simply adding the input + * vector to the output vector. + */ + template + void vmult_add (VECTOR1 &out, + const VECTOR2 &in) const; + + /** + * Matrix-vector multiplication + * with the transpose matrix. For + * the present case, this of + * course amounts to simply + * copying the input vector to + * the output vector. + */ + template + void Tvmult (VECTOR1 &out, + const VECTOR2 &in) const; + + + /** + * Matrix-vector multiplication + * with the transpose matrix, + * with addition to the output + * vector. For the present case, + * this of course amounts to + * simply adding the input vector + * to the output vector. + */ + template + void Tvmult_add (VECTOR1 &out, + const VECTOR2 &in) const; + private: + + /** + * Number of rows and columns of + * this matrix. + */ + unsigned int size; +}; + + + + +// ------------------------- inline and template functions ------------- +#ifndef DOXYGEN + + +inline +IdentityMatrix::IdentityMatrix () + : + size (0) +{} + + + +inline +IdentityMatrix::IdentityMatrix (const unsigned int n) + : + size (n) +{} + + + +inline +void +IdentityMatrix::reinit (const unsigned int n) +{ + size = n; +} + + + +inline +unsigned int +IdentityMatrix::m () const +{ + return size; +} + + + +inline +unsigned int +IdentityMatrix::n () const +{ + return size; +} + + + +template +inline +void +IdentityMatrix::vmult (VECTOR1 &out, + const VECTOR2 &in) const +{ + Assert (out.size() == size, ExcDimensionMismatch (out.size(), size)); + Assert (in.size() == size, ExcDimensionMismatch (in.size(), size)); + + out = in; +} + + + +template +inline +void +IdentityMatrix::vmult_add (VECTOR1 &out, + const VECTOR2 &in) const +{ + Assert (out.size() == size, ExcDimensionMismatch (out.size(), size)); + Assert (in.size() == size, ExcDimensionMismatch (in.size(), size)); + + out += in; +} + + + +template +inline +void +IdentityMatrix::Tvmult (VECTOR1 &out, + const VECTOR2 &in) const +{ + Assert (out.size() == size, ExcDimensionMismatch (out.size(), size)); + Assert (in.size() == size, ExcDimensionMismatch (in.size(), size)); + + out = in; +} + + + +template +inline +void +IdentityMatrix::Tvmult_add (VECTOR1 &out, + const VECTOR2 &in) const +{ + Assert (out.size() == size, ExcDimensionMismatch (out.size(), size)); + Assert (in.size() == size, ExcDimensionMismatch (in.size(), size)); + + out += in; +} + + +#endif + +/**@}*/ + +#endif + -- 2.39.5