From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Sat, 12 Mar 2005 03:26:12 +0000 (+0000)
Subject: beginning TridiagonalMatrix
X-Git-Tag: v8.0.0~14424
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0dde13a5950412e4bee5e48e7af1ad530115b691;p=dealii.git

beginning TridiagonalMatrix


git-svn-id: https://svn.dealii.org/trunk@10109 0785d39b-7218-0410-832d-ea1e28bc413d
---

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 <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
index 0000000000..3c33373af0
--- /dev/null
+++ b/deal.II/lac/source/tridiagonal_matrix.cc
@@ -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>;
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 <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>();
+}