]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new ProductSparseMatrix
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sat, 12 Mar 2005 02:51:48 +0000 (02:51 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sat, 12 Mar 2005 02:51:48 +0000 (02:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@10106 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.html
deal.II/lac/include/lac/matrix_lib.h
deal.II/lac/source/matrix_lib.cc
tests/lac/Makefile
tests/lac/matrix_lib.cc [new file with mode: 0644]

index b351e2cc7ef78771760c28c39513e9267f7a7aa8..b0607df23e4f4c4ee8f2d82008cca183896ce548 100644 (file)
@@ -161,6 +161,14 @@ inconvenience this causes.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>
+       New: The <code class="class">ProductSparseMatrix</code>
+       implements the product of two rectangular sparse matrices with
+       the same <code class="member">value_type</code>
+       <br> 
+       (GK, 2005/03/11)
+       </p>
+
   <li> <p>
        New: The <code class="class">PreconditionRichardson</code>
        implements a Richardson preconditioner.
index f073c1a7b4012aab78f4a172c9e5498192c402d5..e8bfd742453c20d7451b607aa45738573f6d4786 100644 (file)
@@ -9,6 +9,7 @@
 
 template<typename number> class Vector;
 template<typename number> class BlockVector;
+template<typename number> class SparseMatrix;
 
 /*! @addtogroup Matrix2
  *@{
@@ -96,6 +97,92 @@ class ProductMatrix : public PointerMatrixBase<VECTOR>
 };
 
 
+/**
+ * Poor man's matrix product of two sparse matrices. Stores two
+ * matrices #m1 and #m2 of arbitrary type SparseMatrix and implements
+ * matrix-vector multiplications for the product
+ * <i>M<sub>1</sub>M<sub>2</sub></i> by performing multiplication with
+ * both factors consecutively.
+ *
+ * The documentation of @ProductMatrix applies with exception that
+ * these matrices here may be rectangular.
+ *
+ * @author Guido Kanschat, 2000, 2001, 2002, 2005
+ */
+template<typename number, typename vector_number>
+class ProductSparseMatrix : public PointerMatrixBase<Vector<vector_number> >
+{
+  public:
+                                    /**
+                                     * Define the type of matrices used.
+                                     */
+    typedef SparseMatrix<number> MatrixType;
+    
+                                    /**
+                                     * Define the type of vectors we
+                                     * plly this matrix to.
+                                     */
+    typedef Vector<vector_number> VectorType;
+    
+                                    /**
+                                     * Constructor.  Additionally to
+                                     * the two constituting matrices, a
+                                     * memory pool for the auxiliary
+                                     * vector must be provided.
+                                     */
+    ProductSparseMatrix(const MatrixType& m1,
+                       const MatrixType& m2,
+                       VectorMemory<VectorType>& mem);
+
+                                    /**
+                                     * Matrix-vector product <i>w =
+                                     * m1 * m2 * v</i>.
+                                     */
+    virtual void vmult (VectorType&       w,
+                       const VectorType& v) const;
+    
+                                    /**
+                                     * Tranposed matrix-vector
+                                     * product <i>w = m2<sup>T</sup>
+                                     * * m1<sup>T</sup> * v</i>.
+                                     */
+    virtual void Tvmult (VectorType&       w,
+                        const VectorType& v) const;
+    
+                                    /**
+                                     * Adding matrix-vector product
+                                     * <i>w += m1 * m2 * v</i>
+                                     */
+    virtual void vmult_add (VectorType&       w,
+                           const VectorType& v) const;
+    
+                                    /**
+                                     * Adding, tranposed
+                                     * matrix-vector product <i>w +=
+                                     * m2<sup>T</sup> *
+                                     * m1<sup>T</sup> * v</i>.
+                                     */
+    virtual void Tvmult_add (VectorType&       w,
+                            const VectorType& v) const;
+    
+  private:
+                                    /**
+                                     * The left matrix of the product.
+                                     */
+    SmartPointer<const MatrixType> m1;
+    
+                                    /**
+                                     * The right matrix of the product.
+                                     */
+    SmartPointer<const MatrixType> m2;
+    
+                                    /**
+                                     * Memory for auxiliary vector.
+                                     */
+    SmartPointer<VectorMemory<VectorType> > mem;
+};
+
+
 /**
  * Mean value filter.  The @p vmult functions of this matrix filter
  * out mean values of the vector.  If the vector is of type
index 3cf1345ed2bea9ba4fe6fe83d82e71aa1c86707a..5e16eb1e50c16b92ad8684cc94b784ee79d45547 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
+//    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
 
 
 #include <lac/matrix_lib.templates.h>
+#include <lac/sparse_matrix.h>
 
 MeanValueFilter::MeanValueFilter(unsigned int component)
                :
                component(component)
 {}
 
+
+template<typename number, typename vnumber>
+ProductSparseMatrix<number, vnumber>::ProductSparseMatrix(
+  const MatrixType& mat1,
+  const MatrixType& mat2,
+  VectorMemory<VectorType>& mem)
+               :
+               m1(&mat1), m2(&mat2),
+               mem(&mem)
+{
+  Assert(mat1.n() == mat2.m(), ExcDimensionMismatch(mat1.n(),mat2.m()));
+}
+
+
+template<typename number, typename vnumber>
+void
+ProductSparseMatrix<number, vnumber>::vmult (VectorType& dst, const VectorType& src) const
+{
+  VectorType* v = mem->alloc();
+  v->reinit(m1->n());
+  m2->vmult (*v, src);
+  m1->vmult (dst, *v);
+  mem->free(v);
+}
+
+
+template<typename number, typename vnumber>
+void
+ProductSparseMatrix<number, vnumber>::vmult_add (VectorType& dst, const VectorType& src) const
+{
+  VectorType* v = mem->alloc();
+  v->reinit(m1->n());
+  m2->vmult (*v, src);
+  m1->vmult_add (dst, *v);
+  mem->free(v);
+}
+
+
+template<typename number, typename vnumber>
+void
+ProductSparseMatrix<number, vnumber>::Tvmult (VectorType& dst, const VectorType& src) const
+{
+  VectorType* v = mem->alloc();
+  v->reinit(m1->n());
+  m1->Tvmult (*v, src);
+  m2->Tvmult (dst, *v);
+  mem->free(v);
+}
+
+
+template<typename number, typename vnumber>
+void
+ProductSparseMatrix<number, vnumber>::Tvmult_add (VectorType& dst, const VectorType& src) const
+{
+  VectorType* v = mem->alloc();
+  v->reinit(m1->n());
+  m1->Tvmult (*v, src);
+  m2->Tvmult_add (dst, *v);
+  mem->free(v);
+}
+
+
+template class ProductSparseMatrix<double, double>;
+template class ProductSparseMatrix<double, float>;
+template class ProductSparseMatrix<float, double>;
+template class ProductSparseMatrix<float, float>;
+
 template void MeanValueFilter::filter(Vector<float>&) const;
 template void MeanValueFilter::filter(Vector<double>&) const;
 template void MeanValueFilter::filter(BlockVector<float>&) const;
index 4a885957a2a4b2a4cab7e9cbe0f247b68c1de4d5..533f16558192932ba131a52df01d690dd5a1da97 100644 (file)
@@ -1,6 +1,6 @@
 ############################################################
 # $Id$
-# Copyright (C) 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
 ############################################################
 
 ############################################################
@@ -20,18 +20,13 @@ default: run-tests
 
 ############################################################
 
-tests_x = $(sort \
-            sparse_matrices \
-            block_matrices \
-            block_vector \
-            block_vector_iterator \
-            full_matrix \
-            matrix_out \
-            solver \
-            eigen \
-            sparsity_pattern \
-            sparse_ilu \
-            vector-vector)
+tests_x = vector-vector \
+       full_matrix sparsity_pattern sparse_matrices \
+       block_vector block_vector_iterator block_matrices \
+       matrix_lib matrix_out \
+       solver eigen \
+       sparse_ilu
+
 
 # from above list of regular expressions, generate the real set of
 # tests
diff --git a/tests/lac/matrix_lib.cc b/tests/lac/matrix_lib.cc
new file mode 100644 (file)
index 0000000..5f9e584
--- /dev/null
@@ -0,0 +1,89 @@
+//---------------------------------------------------------------------------
+//    $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/matrix_lib.h>
+#include <lac/sparse_matrix.h>
+#include <lac/vector.h>
+#include <lac/vector_memory.h>
+
+#include <fstream>
+#include <iostream>
+#include <cmath>
+
+template<typename number>
+void check_sparse_product(const SparseMatrix<number>& m1, SparseMatrix<number>& m2)
+{
+  Vector<double> v(m2.n());
+  Vector<double> w(m1.m());
+  deallog << "Sizes\t" << v.size() << '\t' << w.size() << std::endl;
+  
+  for (unsigned int i=0;i<v.size();++i)
+    v(i) = i+1.;
+  
+  for (unsigned int i=0;i<v.size();++i)
+    deallog << ' ' << v(i);
+  deallog << std::endl;
+  
+  GrowingVectorMemory<Vector<double> > mem;
+  
+  ProductSparseMatrix<number, number> product(m1, m2, mem);
+  product.vmult(w,v);
+
+  for (unsigned int i=0;i<w.size();++i)
+    deallog << ' ' << w(i);
+  deallog << std::endl;
+
+  product.Tvmult(v,w);
+  
+  for (unsigned int i=0;i<v.size();++i)
+    deallog << ' ' << v(i);
+  deallog << std::endl;
+}
+
+
+int main()
+{
+  std::ofstream logfile("matrix_lib.output");
+  logfile.setf(std::ios::fixed);
+  logfile.precision(0);
+  deallog.attach(logfile);
+  deallog.depth_console(10);
+  
+  SparsityPattern sparsity1(2,3,3);
+  SparsityPattern sparsity2(3,4,4);
+
+  for (unsigned int i=0;i<2;++i)
+    for (unsigned int j=0;j<3;++j)
+      {
+       sparsity1.add(i,j);
+       sparsity2.add(j,i);
+       sparsity2.add(j,2+i);
+      }
+  sparsity1.compress();
+  sparsity2.compress();
+
+  SparseMatrix<double> m1(sparsity1);
+  SparseMatrix<double> m2(sparsity2);
+
+  for (unsigned int i=0;i<2;++i)
+    for (unsigned int j=0;j<3;++j)
+      {
+       m1.set(i, j, 1.*i-j);
+       m2.set(j, i, 1.*i-j);
+       m2.set(j, 2+i, 1.*j-i);
+      }
+  check_sparse_product(m1, m2);
+}

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.