]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new TransposeMatrix
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 7 Jul 2006 09:39:32 +0000 (09:39 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 7 Jul 2006 09:39:32 +0000 (09:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@13349 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.html
deal.II/lac/include/lac/transpose_matrix.h [new file with mode: 0644]
tests/lac/Makefile
tests/lac/pointer_matrix.cc [new file with mode: 0644]
tests/lac/pointer_matrix/.cvsignore [new file with mode: 0644]
tests/lac/pointer_matrix/cmp/generic [new file with mode: 0644]

index 1e6420d8c014a898b32b8d9dc1d93a4ed16d327b..c25d9e5bc036d0bf2a731d80ea0887cac39d24d4 100644 (file)
@@ -529,7 +529,15 @@ inconvenience this causes.
 <h3>lac</h3>
 
 <ol>
-  
+
+  <li> <p>New: The class <code class="class">TransposeMatrix</code>
+  modeled after <code class="class">PointerMatrix</code> swaps the
+  <code class="member">vmult</code> functions such that its effect is
+  the transpose of the matrix it points to.
+  <br>
+  (GK 2006/07/07)
+  </p>  
+
   <li> <p>New: There is now a function <code
        class="member">FullMatrix::trace</code> that does what its name suggests
        it does.
diff --git a/deal.II/lac/include/lac/transpose_matrix.h b/deal.II/lac/include/lac/transpose_matrix.h
new file mode 100644 (file)
index 0000000..28fa8e0
--- /dev/null
@@ -0,0 +1,228 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2002, 2003, 2004, 2005, 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__transpose_matrix_h
+#define __deal2__transpose_matrix_h
+
+
+#include <base/subscriptor.h>
+#include <lac/pointer_matrix.h>
+#include <boost/shared_ptr.hpp>
+
+
+/**
+ * The transpose of a given matrix.  This auxiliary class swaps the
+ * effect ov vmult() and Tvmult() as well as vmult_add() and
+ * Tvmult_add().
+ *
+ * The implementation is analogous to the class PointerMatrix.
+ *
+ * @note The transposed matrix is never actually assembled. Instead,
+ * only the matrix vector multiplication is performed in a transposed
+ * way.
+ *
+ * @ingroup Matrix2
+ * @author Guido Kanschat, 2006
+ */
+
+template<class MATRIX, class VECTOR>
+class
+TransposeMatrix : public PointerMatrixBase<VECTOR>
+{
+  public:
+                                    /**
+                                     * Constructor.  The pointer in the
+                                     * argument is stored in this
+                                     * class. As usual, the lifetime of
+                                     * <tt>*M</tt> must be longer than the
+                                     * one of the PointerMatrix.
+                                     *
+                                     * If <tt>M</tt> is zero, no
+                                     * matrix is stored.
+                                     */
+    TransposeMatrix (const MATRIX* M=0);
+
+                                    /**
+                                     * Constructor. The name argument
+                                     * is used to identify the
+                                     * SmartPointer for this object.
+                                     */
+    TransposeMatrix(const char* name);
+    
+                                    /**
+                                     * Constructor. <tt>M</tt> points
+                                     * to a matrix which must live
+                                     * longer than the
+                                     * TransposeMatrix. The name
+                                     * argument is used to identify
+                                     * the SmartPointer for this
+                                     * object.
+                                     */
+    TransposeMatrix(const MATRIX* M,
+                 const char* name);
+
+                                    // Use doc from base class
+    virtual void clear();
+    
+                                    /**
+                                     * Return whether the object is
+                                     * empty. 
+                                     */
+    bool empty () const;
+    
+                                    /**
+                                     * Assign a new matrix
+                                     * pointer. Deletes the old pointer
+                                     * and releases its matrix.
+                                     * @see SmartPointer
+                                     */
+    const TransposeMatrix& operator= (const MATRIX* M);
+    
+                                    /**
+                                     * Matrix-vector product.
+                                     */
+    virtual void vmult (VECTOR& dst,
+                       const VECTOR& src) const;
+    
+                                    /**
+                                     * Tranposed matrix-vector product.
+                                     */
+    virtual void Tvmult (VECTOR& dst,
+                        const VECTOR& src) const;
+    
+                                    /**
+                                     * Matrix-vector product, adding to
+                                     * <tt>dst</tt>.
+                                     */
+    virtual void vmult_add (VECTOR& dst,
+                           const VECTOR& src) const;
+    
+                                    /**
+                                     * Tranposed matrix-vector product,
+                                     * adding to <tt>dst</tt>.
+                                     */
+    virtual void Tvmult_add (VECTOR& dst,
+                            const VECTOR& src) const;
+    
+  private:
+                                    /**
+                                     * Return the address of the
+                                     * matrix for comparison.
+                                     */
+    virtual const void* get() const;
+
+                                    /**
+                                     * The pointer to the actual matrix.
+                                     */
+    SmartPointer<const MATRIX> m;
+};
+
+
+//----------------------------------------------------------------------//
+
+
+template<class MATRIX, class VECTOR>
+TransposeMatrix<MATRIX, VECTOR>::TransposeMatrix (const MATRIX* M)
+  : m(M, typeid(*this).name())
+{}
+
+
+template<class MATRIX, class VECTOR>
+TransposeMatrix<MATRIX, VECTOR>::TransposeMatrix (const char* name)
+  : m(0, name)
+{}
+
+
+template<class MATRIX, class VECTOR>
+TransposeMatrix<MATRIX, VECTOR>::TransposeMatrix (
+  const MATRIX* M,
+  const char* name)
+  : m(M, name)
+{}
+
+
+template<class MATRIX, class VECTOR>
+inline void
+TransposeMatrix<MATRIX, VECTOR>::clear ()
+{
+  m = 0;
+}
+
+
+template<class MATRIX, class VECTOR>
+inline const TransposeMatrix<MATRIX, VECTOR>&
+TransposeMatrix<MATRIX, VECTOR>::operator= (const MATRIX* M)
+{
+  m = M;
+  return *this;
+}
+
+
+template<class MATRIX, class VECTOR>
+inline bool
+TransposeMatrix<MATRIX, VECTOR>::empty () const
+{
+  if (m == 0)
+    return true;
+  return m->empty();
+}
+
+template<class MATRIX, class VECTOR>
+inline void
+TransposeMatrix<MATRIX, VECTOR>::vmult (VECTOR& dst,
+                                     const VECTOR& src) const
+{
+  Assert (m != 0, ExcNotInitialized());
+  m->Tvmult (dst, src);
+}
+
+
+template<class MATRIX, class VECTOR>
+inline void
+TransposeMatrix<MATRIX, VECTOR>::Tvmult (VECTOR& dst,
+                                      const VECTOR& src) const
+{
+  Assert (m != 0, ExcNotInitialized());
+  m->vmult (dst, src);
+}
+
+
+template<class MATRIX, class VECTOR>
+inline void
+TransposeMatrix<MATRIX, VECTOR>::vmult_add (VECTOR& dst,
+                                         const VECTOR& src) const
+{
+  Assert (m != 0, ExcNotInitialized());
+  m->Tvmult_add (dst, src);
+}
+
+
+template<class MATRIX, class VECTOR>
+inline void
+TransposeMatrix<MATRIX, VECTOR>::Tvmult_add (VECTOR& dst,
+                                          const VECTOR& src) const
+{
+  Assert (m != 0, ExcNotInitialized());
+  m->vmult_add (dst, src);
+}
+
+
+template<class MATRIX, class VECTOR>
+inline const void*
+TransposeMatrix<MATRIX, VECTOR>::get () const
+{
+  return m;
+}
+
+
+
+#endif
index 8365ca85b167a225d811ac5cf8929d87e3524126..1dd33e50ae39887eaf48d15bfeb51422f299b04c 100644 (file)
@@ -25,7 +25,7 @@ tests_x = vector-vector \
        sparsity_pattern sparse_matrices matrices \
        block_sparsity_pattern_01 \
        block_vector block_vector_* block_matrices \
-       matrix_lib matrix_out \
+       matrix_lib matrix_out pointer_matrix \
        solver eigen \
        sparse_ilu sparse_ilu_t lapack lapack_fill block_minres \
        identity_matrix* trace
diff --git a/tests/lac/pointer_matrix.cc b/tests/lac/pointer_matrix.cc
new file mode 100644 (file)
index 0000000..4936e95
--- /dev/null
@@ -0,0 +1,70 @@
+//----------------------------------------------------------------------
+//    $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.
+//
+//----------------------------------------------------------------------
+
+// Test vmult and Tvmult of PointerMatrix and TransposeMatrix
+
+#include <base/logstream.h>
+#include <lac/pointer_matrix.h>
+#include <lac/transpose_matrix.h>
+#include <lac/full_matrix.h>
+#include <lac/vector.h>
+
+#include <fstream>
+
+int main()
+{
+  std::ofstream logfile("pointer_matrix/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  FullMatrix<double> A(5,5);
+  unsigned int k=0;
+  for (unsigned int i=0;i<A.m();++i)
+    for (unsigned int j=0;j<A.n();++j)
+      A(i,j) = ++k;
+
+  PointerMatrix<FullMatrix<double>, Vector<double> > P(&A, "P");
+  TransposeMatrix<FullMatrix<double>, Vector<double> > T(&A, "T");
+  
+  Vector<double> x(5);
+  Vector<double> y(5);
+  Vector<double> y2(5);
+  Vector<double> diff(5);
+  
+  for (unsigned int i=0;i<x.size();++i)
+    {
+      x = 0.;
+      x(i) = 1.;
+      A.vmult(y,x);
+      P.vmult(y2,x);
+      diff = y;
+      diff.add(-1., y2);
+      deallog << "P.vmult:  diff " << diff.l2_norm() << std::endl;
+      T.Tvmult(y2,x);
+      diff = y;
+      diff.add(-1., y2);
+      deallog << "T.Tvmult: diff " << diff.l2_norm() << std::endl;
+      
+      A.Tvmult(y,x);
+      P.Tvmult(y2,x);
+      diff = y;
+      diff.add(-1., y2);
+      deallog << "P.Tvmult: diff " << diff.l2_norm() << std::endl;
+      T.vmult(y2,x);
+      diff = y;
+      diff.add(-1., y2);
+      deallog << "T.vmult:  diff " << diff.l2_norm() << std::endl;
+      
+      
+    }
+}
diff --git a/tests/lac/pointer_matrix/.cvsignore b/tests/lac/pointer_matrix/.cvsignore
new file mode 100644 (file)
index 0000000..d196dfb
--- /dev/null
@@ -0,0 +1 @@
+obj.* exe OK output
diff --git a/tests/lac/pointer_matrix/cmp/generic b/tests/lac/pointer_matrix/cmp/generic
new file mode 100644 (file)
index 0000000..0e3ab3c
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::P.vmult:  diff 0
+DEAL::T.Tvmult: diff 0
+DEAL::P.Tvmult: diff 0
+DEAL::T.vmult:  diff 0
+DEAL::P.vmult:  diff 0
+DEAL::T.Tvmult: diff 0
+DEAL::P.Tvmult: diff 0
+DEAL::T.vmult:  diff 0
+DEAL::P.vmult:  diff 0
+DEAL::T.Tvmult: diff 0
+DEAL::P.Tvmult: diff 0
+DEAL::T.vmult:  diff 0
+DEAL::P.vmult:  diff 0
+DEAL::T.Tvmult: diff 0
+DEAL::P.Tvmult: diff 0
+DEAL::T.vmult:  diff 0
+DEAL::P.vmult:  diff 0
+DEAL::T.Tvmult: diff 0
+DEAL::P.Tvmult: diff 0
+DEAL::T.vmult:  diff 0

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.