]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
matrix library started
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 13 Aug 2002 16:01:43 +0000 (16:01 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 13 Aug 2002 16:01:43 +0000 (16:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@6323 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/matrix_lib.h [new file with mode: 0644]
deal.II/lac/include/lac/matrix_lib.templates.h [new file with mode: 0644]
deal.II/lac/include/lac/pointer_matrix.h
deal.II/lac/source/matrix_lib.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/matrix_lib.h b/deal.II/lac/include/lac/matrix_lib.h
new file mode 100644 (file)
index 0000000..4aa5d48
--- /dev/null
@@ -0,0 +1,245 @@
+// $Id$
+
+#ifndef __deal2__matrix_lib_h
+#define __deal2__matrix_lib_h
+
+#include <base/subscriptor.h>
+#include <lac/vector_memory.h>
+#include <lac/pointer_matrix.h>
+
+template<typename number> class Vector;
+template<typename number> class BlockVector;
+
+/**
+ * Poor man's matrix product. Stores two matrices $m_1$ and $m_2$ and
+ * implements matrix-vector multiplications for the product $m_1 m_2$
+ * by performing multiplication with both factors consecutively.
+ *
+ * @author Guido Kanschat, 2000, 2001, 2002
+ */
+template<class VECTOR>
+class ProductMatrix : public PointerMatrixBase<VECTOR>
+{
+  public:
+                                    /**
+                                     * Constructor.  Additionally to
+                                     * the two constituting matrices, a
+                                     * memory pool for the auxiliary
+                                     * vector must be provided.
+                                     */
+    template <class MATRIX1, class MATRIX2>
+    ProductMatrix(const MATRIX1& m1,
+                 const MATRIX2& m2,
+                 VectorMemory<VECTOR>& mem);
+
+                                    /**
+                                     * Destructor.
+                                     */
+    ~ProductMatrix();
+    
+                                    /**
+                                     * 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
+                                     * @p{dst}.
+                                     */
+    virtual void vmult_add (VECTOR& dst,
+                           const VECTOR& src) const;
+    
+                                    /**
+                                     * Tranposed matrix-vector product,
+                                     * adding to @p{dst}.
+                                     */
+    virtual void Tvmult_add (VECTOR& dst,
+                            const VECTOR& src) const;
+    
+  private:
+                                    /**
+                                     * The left matrix of the product.
+                                     */
+    PointerMatrixBase<VECTOR>* m1;
+    
+                                    /**
+                                     * The right matrix of the product.
+                                     */
+    PointerMatrixBase<VECTOR>* m2;
+    
+                                    /**
+                                     * Memory for auxiliary vector.
+                                     */
+    SmartPointer<VectorMemory<VECTOR> > mem;
+};
+
+
+/**
+ * Mean value filter.  The @p{vmult} functions of this matrix filter
+ * out mean values of the vector.  If the vector is of type
+ * @p{BlockVector}, then an additional parameter selects a single
+ * component for this operation.
+ *
+ * @author Guido Kanschat, 2002
+ */
+class MeanValueFilter : public Subscriptor
+{
+  public:
+                                    /**
+                                     * Constructor, optionally
+                                     * selecting a component.
+                                     */
+    MeanValueFilter(unsigned int component = static_cast<unsigned int>(-1));
+
+                                    /**
+                                     * Return the source vector with
+                                     * subtracted mean value.
+                                     */
+    template <typename number>
+    void vmult (Vector<number>& dst,
+               const Vector<number>& src) const;
+    
+                                    /**
+                                     * Add source vector with
+                                     * subtracted mean value to dest.
+                                     */
+    template <typename number>
+    void vmult_add (Vector<number>& dst,
+                   const Vector<number>& src) const;
+    
+                                    /**
+                                     * Return the source vector with
+                                     * subtracted mean value in
+                                     * selected component.
+                                     */
+    template <typename number>
+    void vmult (BlockVector<number>& dst,
+               const BlockVector<number>& src) const;
+    
+                                    /**
+                                     * Add a soruce to dest, where
+                                     * the mean value in the selected
+                                     * component is subtracted.
+                                     */
+    template <typename number>
+    void vmult_add (BlockVector<number>& dst,
+                   const BlockVector<number>& src) const;
+
+
+                                    /**
+                                     * Not implemented.
+                                     */
+    template <typename VECTOR>
+    void Tvmult(VECTOR&, const VECTOR&) const;
+    
+                                    /**
+                                     * Not implemented.
+                                     */
+    template <typename VECTOR>
+    void Tvmult_add(VECTOR&, const VECTOR&) const;
+    
+  private:
+                                    /**
+                                     * Component for filtering block vectors.
+                                     */
+    unsigned int component;
+};
+
+
+//----------------------------------------------------------------------//
+
+template<class VECTOR>
+template<class MATRIX1, class MATRIX2>
+ProductMatrix<VECTOR>::ProductMatrix (
+  const MATRIX1& mat1,
+  const MATRIX2& mat2,
+  VectorMemory<VECTOR>& m)
+  : mem(&m)
+{
+  m1 = new PointerMatrix<MATRIX1, VECTOR>(&mat1);
+  m2 = new PointerMatrix<MATRIX2, VECTOR>(&mat2);
+}
+
+
+template<class VECTOR>
+ProductMatrix<VECTOR>::~ProductMatrix ()
+{
+  delete m1;
+  delete m2;
+}
+
+
+template<class VECTOR>
+void
+ProductMatrix<VECTOR>::vmult (VECTOR& dst, const VECTOR& src) const
+{
+  VECTOR* v = mem->alloc();
+  v->reinit(dst);
+  m2->vmult (*v, src);
+  m1->vmult (dst, *v);
+  mem->free(v);
+}
+
+
+template<class VECTOR>
+void
+ProductMatrix<VECTOR>::vmult_add (VECTOR& dst, const VECTOR& src) const
+{
+  VECTOR* v = mem->alloc();
+  v->reinit(dst);
+  m2->vmult (*v, src);
+  m1->vmult_add (dst, *v);
+  mem->free(v);
+}
+
+
+template<class VECTOR>
+void
+ProductMatrix<VECTOR>::Tvmult (VECTOR& dst, const VECTOR& src) const
+{
+  VECTOR* v = mem->alloc();
+  v->reinit(dst);
+  m1->Tvmult (*v, src);
+  m2->Tvmult (dst, *v);
+  mem->free(v);
+}
+
+
+template<class VECTOR>
+void
+ProductMatrix<VECTOR>::Tvmult_add (VECTOR& dst, const VECTOR& src) const
+{
+  VECTOR* v = mem->alloc();
+  v->reinit(dst);
+  m1->Tvmult (*v, src);
+  m2->Tvmult_add (dst, *v);
+  mem->free(v);
+}
+
+
+//----------------------------------------------------------------------//
+
+template <class VECTOR>
+inline void
+MeanValueFilter::Tvmult(VECTOR&, const VECTOR&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template <class VECTOR>
+inline void
+MeanValueFilter::Tvmult_add(VECTOR&, const VECTOR&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+#endif
diff --git a/deal.II/lac/include/lac/matrix_lib.templates.h b/deal.II/lac/include/lac/matrix_lib.templates.h
new file mode 100644 (file)
index 0000000..9dd11e2
--- /dev/null
@@ -0,0 +1,91 @@
+//----------------------------  vector.templates.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 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.
+//
+//----------------------------  vector.templates.h  ---------------------------
+#ifndef __deal2__matrix_lib_templates_h
+#define __deal2__matrix_lib_templates_h
+
+#include <lac/matrix_lib.h>
+#include <lac/vector.h>
+#include <lac/block_vector.h>
+
+
+template <typename number>
+void
+MeanValueFilter::vmult(Vector<number>& dst,
+                      const Vector<number>& src) const
+{
+  Assert (dst.size() == src.size(),
+         ExcDimensionMismatch(dst.size(), src.size()));
+  
+  number mean = src.mean_value();
+
+  for (unsigned int i=0;i<dst.size();++i)
+    dst(i) = src(i) - mean;
+}
+
+
+
+template <typename number>
+void
+MeanValueFilter::vmult_add(Vector<number>& dst,
+                          const Vector<number>& src) const
+{
+  Assert (dst.size() == src.size(),
+         ExcDimensionMismatch(dst.size(), src.size()));
+  
+  number mean = src.mean_value();
+  
+  for (unsigned int i=0;i<dst.size();++i)
+    dst(i) += src(i) - mean;
+}
+
+
+
+template <typename number>
+void
+MeanValueFilter::vmult(BlockVector<number>& dst,
+                          const BlockVector<number>& src) const
+{
+  Assert (component != static_cast<unsigned int>(-1),
+         ExcNotInitialized());
+  
+  Assert (dst.n_blocks() == src.n_blocks(),
+         ExcDimensionMismatch(dst.n_blocks(), src.n_blocks()));
+  
+  for (unsigned int i=0;i<dst.n_blocks();++i)
+    if (i == component)
+      vmult(dst.block(i), src.block(i));
+    else
+      dst.block(i) = src.block(i);
+}
+
+
+
+template <typename number>
+void
+MeanValueFilter::vmult_add(BlockVector<number>& dst,
+                          const BlockVector<number>& src) const
+{
+  Assert (component != static_cast<unsigned int>(-1),
+         ExcNotInitialized());
+  
+  Assert (dst.n_blocks() == src.n_blocks(),
+         ExcDimensionMismatch(dst.n_blocks(), src.n_blocks()));
+  
+  for (unsigned int i=0;i<dst.n_blocks();++i)
+    if (i == component)
+      vmult_add(dst.block(i), src.block(i));
+    else
+      dst.block(i).add(src.block(i));
+}
+
+#endif
index 5eddb1e271e6755e16065dd71abaca2bdb421131..bd89093c7570ecf604618b512966b64e2ef873c8 100644 (file)
@@ -4,6 +4,7 @@
 #define __deal2__pointer_matrix_h
 
 #include <base/subscriptor.h>
+#include <base/smartpointer.h>
 #include <lac/vector_memory.h>
 
 /**
diff --git a/deal.II/lac/source/matrix_lib.cc b/deal.II/lac/source/matrix_lib.cc
new file mode 100644 (file)
index 0000000..0422d0d
--- /dev/null
@@ -0,0 +1,39 @@
+//----------------------------  vector.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 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.
+//
+//----------------------------  vector.cc  ---------------------------
+
+
+#include <lac/matrix_lib.templates.h>
+
+MeanValueFilter::MeanValueFilter(unsigned int component)
+               :
+               component(component)
+{}
+
+template void MeanValueFilter::vmult(Vector<float>&,
+                                    const Vector<float>&) const;
+template void MeanValueFilter::vmult(Vector<double>&,
+                                    const Vector<double>&) const;
+template void MeanValueFilter::vmult(BlockVector<float>&,
+                                    const BlockVector<float>&) const;
+template void MeanValueFilter::vmult(BlockVector<double>&,
+                                    const BlockVector<double>&) const;
+
+template void MeanValueFilter::vmult_add(Vector<float>&,
+                                        const Vector<float>&) const;
+template void MeanValueFilter::vmult_add(Vector<double>&,
+                                        const Vector<double>&) const;
+template void MeanValueFilter::vmult_add(BlockVector<float>&,
+                                        const BlockVector<float>&) const;
+template void MeanValueFilter::vmult_add(BlockVector<double>&,
+                                        const BlockVector<double>&) const;
+

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.