<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.
template<typename number> class Vector;
template<typename number> class BlockVector;
+template<typename number> class SparseMatrix;
/*! @addtogroup Matrix2
*@{
};
+/**
+ * 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
// $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;
############################################################
# $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
############################################################
############################################################
############################################################
-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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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);
+}
--- /dev/null
+
+DEAL::Sizes 4 2
+DEAL:: 1. 2. 3. 4.
+DEAL:: -14. -8.
+DEAL:: -86. -44. 86. 44.
+DEAL::GrowingVectorMemory:Overall allocated vectors: 2
+DEAL::GrowingVectorMemory:Maximum allocated vectors: 1