]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
beginning TridiagonalMatrix
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 12 Mar 2005 03:26:12 +0000 (03:26 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 12 Mar 2005 03:26:12 +0000 (03:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@10109 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/tridiagonal_matrix.h [new file with mode: 0644]
deal.II/lac/source/tridiagonal_matrix.cc [new file with mode: 0644]
tests/lac/Makefile
tests/lac/tridiagonal_matrix.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/tridiagonal_matrix.h b/deal.II/lac/include/lac/tridiagonal_matrix.h
new file mode 100644 (file)
index 0000000..64a7fd9
--- /dev/null
@@ -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 <base/config.h>
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+
+#include <vector>
+
+// forward declarations
+template<typename number> 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<typename number>
+class TridiagonalMatrix
+{
+  public:
+/**
+ * @name Constructors and initalization.
+ */
+                                    /**
+                                     * Constructor generating an
+                                     * empty matrix of dimension
+                                     * <tt>n</tt>.
+                                     */
+    TridiagonalMatrix(unsigned int n = 0);
+
+//@}
+///@name Non-modifying operators
+//@{
+
+                                    /**
+                                     * Number of rows of this matrix.
+                                     * To remember: this matrix is an
+                                     * <i>m x m</i>-matrix.
+                                     */
+    unsigned int m () const;
+    
+                                    /**
+                                     * Number of columns of this matrix.
+                                     * To remember: this matrix is an
+                                     * <i>n x n</i>-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>|i-j| <=
+                                     * 1</i>.
+                                     */
+    number operator()(unsigned int i, unsigned int j) const;
+    
+                                    /**
+                                     * Read-write access to a
+                                     * value. This is restricted to
+                                     * the case where <i>|i-j| <=
+                                     * 1</i>.
+                                     */
+    number& operator()(unsigned int i, unsigned int j);
+    
+//@}
+///@name Multiplications with vectors
+//@{    
+
+                                    /**
+                                     * Matrix-vector-multiplication. Multiplies
+                                     * <tt>v</tt> from the right and
+                                     * stores the result in
+                                     * <tt>w</tt>.
+                                     *
+                                     * If the optional parameter
+                                     * <tt>adding</tt> is <tt>true</tt>, the
+                                     * result is added to <tt>w</tt>.
+                                      *
+                                      * Source and destination must
+                                      * not be the same vector.
+                                     */
+    void vmult (Vector<number>       &w,
+               const Vector<number> &v,
+               const bool            adding=false) const;
+    
+                                    /**
+                                     * Adding
+                                     * Matrix-vector-multiplication. Same
+                                     * as vmult() with parameter
+                                     * <tt>adding=true</tt>, but
+                                     * widely used in
+                                     * <tt>deal.II</tt> classes.
+                                      *
+                                      * Source and destination must
+                                      * not be the same vector.
+                                     */
+    void vmult_add (Vector<number>       &w,
+                   const Vector<number> &v) const;
+    
+                                    /**
+                                     * Transpose
+                                     * matrix-vector-multiplication.
+                                     * Multiplies
+                                     * <tt>v<sup>T</sup></tt> from
+                                     * the left and stores the result
+                                     * in <tt>w</tt>.
+                                      *
+                                     * If the optional parameter
+                                     * <tt>adding</tt> is <tt>true</tt>, the
+                                     * result is added to <tt>w</tt>.
+                                      *
+                                      * Source and destination must
+                                      * not be the same vector.
+                                     */
+    void Tvmult (Vector<number>       &w,
+                const Vector<number> &v,
+                const bool            adding=false) const;
+
+                                    /**
+                                     * Adding transpose
+                                     * matrix-vector-multiplication. Same
+                                     * as Tvmult() with parameter
+                                     * <tt>adding=true</tt>, but
+                                     * widely used in
+                                     * <tt>deal.II</tt> classes.
+                                      *
+                                      * Source and destination must
+                                      * not be the same vector.
+                                     */
+    void Tvmult_add (Vector<number>       &w,
+                    const Vector<number> &v) const;
+
+                                    /**
+                                     * Build the matrix scalar product
+                                     * <tt>u^T M v</tt>. 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<number> &u,
+                                 const Vector<number> &v) const;
+
+                                    /**
+                                     * Return the square of the norm
+                                     * of the vector <tt>v</tt> with
+                                     * respect to the norm induced by
+                                     * this matrix,
+                                     * i.e. <i>(v,Mv)</i>. This is
+                                     * useful, e.g. in the finite
+                                     * element context, where the
+                                     * <i>L<sup>2</sup></i> 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<number> &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: <tt>precision</tt> and
+                                     * <tt>scientific</tt> are used to
+                                     * determine the number format,
+                                     * where <tt>scientific</tt> = <tt>false</tt>
+                                     * means fixed point notation.  A
+                                     * zero entry for <tt>width</tt> 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<number> 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<number> 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<number> right;
+};
+
+//----------------------------------------------------------------------//
+///@if NoDoc
+
+template<typename number>
+unsigned int
+TridiagonalMatrix<number>::m() const
+{
+  return diagonal.size();
+}
+
+
+
+template<typename number>
+unsigned int
+TridiagonalMatrix<number>::n() const
+{
+  return diagonal.size();
+}
+
+
+template<typename number>
+inline
+number
+TridiagonalMatrix<number>::operator()(unsigned int i, unsigned int j) const
+{
+  Assert(i<n(), ExcIndexRange(i,0,n()));
+  Assert(j<n(), ExcIndexRange(j,0,n()));
+  Assert (i<=j+1, ExcIndexRange(i,j-1,j+2));
+  Assert (j<=i+1, ExcIndexRange(j,i-1,i+2));
+  
+  if (j==i)
+    return diagonal[i];
+  if (j==i-1)
+    return left[i];
+  if (j==i+1)
+    return right[i];
+  AssertThrow(false, ExcInternalError());
+  return 0;
+}
+
+
+template<typename number>
+inline
+number&
+TridiagonalMatrix<number>::operator()(unsigned int i, unsigned int j)
+{
+  Assert(i<n(), ExcIndexRange(i,0,n()));
+  Assert(j<n(), ExcIndexRange(j,0,n()));
+  Assert (i<=j+1, ExcIndexRange(i,j-1,j+2));
+  Assert (j<=i+1, ExcIndexRange(j,i-1,i+2));
+  
+  if (j==i)
+    return diagonal[i];
+  if (j==i-1)
+    return left[i];
+  if (j==i+1)
+    return right[i];
+  AssertThrow(false, ExcInternalError());
+  return diagonal[0];
+}
+
+
+///@endif
+
+#endif
+
diff --git a/deal.II/lac/source/tridiagonal_matrix.cc b/deal.II/lac/source/tridiagonal_matrix.cc
new file mode 100644 (file)
index 0000000..3c33373
--- /dev/null
@@ -0,0 +1,145 @@
+//-------------------------------------------------------------------
+//    $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.
+//
+//-------------------------------------------------------------------
+
+#include <lac/tridiagonal_matrix.h>
+#include <lac/vector.h>
+
+template<typename number>
+TridiagonalMatrix<number>::TridiagonalMatrix(unsigned int size)
+               :
+               diagonal(size, 0.),
+               left(size, 0.),
+               right(size, 0.)
+{}
+
+
+template<typename number>
+bool
+TridiagonalMatrix<number>::all_zero() const
+{
+  typename std::vector<number>::const_iterator i;
+  typename std::vector<number>::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<typename number>
+void
+TridiagonalMatrix<number>::vmult (
+  Vector<number>       &w,
+  const Vector<number> &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<e;++i)
+       w(i) += left[i]*v(i-1)+diagonal[i]*v(i)+right[i]*v(i+1);
+      w(e) += left[e]*v(e-1)+diagonal[e]*v(e);
+    }
+  else
+    {
+      w(0) = diagonal[0]*v(0) + right[0]*v(1);
+      const unsigned int e=n()-1;
+      for (unsigned int i=1;i<e;++i)
+       w(i) = left[i]*v(i-1)+diagonal[i]*v(i)+right[i]*v(i+1);
+      w(e) = left[e]*v(e-1)+diagonal[e]*v(e);
+    }
+}
+
+
+template<typename number>
+void
+TridiagonalMatrix<number>::vmult_add (
+  Vector<number>       &w,
+  const Vector<number> &v) const
+{
+  vmult(w, v, true);
+}
+
+
+template<typename number>
+void
+TridiagonalMatrix<number>::Tvmult (
+  Vector<number>       &w,
+  const Vector<number> &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<e;++i)
+       w(i) += left[i+1]*v(i+1)+diagonal[i]*v(i)+right[i-1]*v(i-1);
+      w(e) += right[e-1]*v(e-1)+diagonal[e]*v(e);
+    }
+  else
+    {
+      w(0) = diagonal[0]*v(0) + left[1]*v(1);
+      const unsigned int e=n()-1;
+      for (unsigned int i=1;i<e;++i)
+       w(i) = left[i+1]*v(i+1)+diagonal[i]*v(i)+right[i-1]*v(i-1);
+      w(e) = right[e-1]*v(e-1)+diagonal[e]*v(e);
+    }
+}
+
+
+template<typename number>
+void
+TridiagonalMatrix<number>::Tvmult_add (
+  Vector<number>       &w,
+  const Vector<number> &v) const
+{
+  Tvmult(w, v, true);
+}
+
+/*
+template<typename number>
+TridiagonalMatrix<number>::
+{
+}
+
+
+template<typename number>
+TridiagonalMatrix<number>::
+{
+}
+
+
+*/
+
+template TridiagonalMatrix<float>;
+template TridiagonalMatrix<double>;
index 533f16558192932ba131a52df01d690dd5a1da97..6884ed1878f30ea68362921d414abddd645a97f8 100644 (file)
@@ -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 (file)
index 0000000..4a28aac
--- /dev/null
@@ -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 <base/logstream.h>
+#include <lac/tridiagonal_matrix.h>
+#include <lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+#include <cmath>
+
+template<typename number>
+void
+check_vmult()
+{
+  TridiagonalMatrix<number> M(4);
+  Vector<number> u(4);
+  Vector<number> 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<v.size();++i)
+    deallog << ' ' << v(i);
+  deallog << std::endl;
+
+  M.vmult_add(v,u);
+  for (unsigned int i=0;i<v.size();++i)
+    deallog << ' ' << v(i);
+  deallog << std::endl;
+
+  M.Tvmult(v,u);
+  for (unsigned int i=0;i<v.size();++i)
+    deallog << ' ' << v(i);
+  deallog << std::endl;
+
+  M.Tvmult_add(v,u);
+  for (unsigned int i=0;i<v.size();++i)
+    deallog << ' ' << v(i);
+  deallog << std::endl;
+}
+
+
+int main()
+{
+  std::ofstream logfile("tridiagonal_matrix.output");
+  logfile.setf(std::ios::fixed);
+  logfile.precision(0);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  
+  check_vmult<double>();
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.