*/
-/**
- * Poor man's matrix product of two quadratic matrices. Stores two quadratic
- * matrices #m1 and #m2 of arbitrary types 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. Because the
- * types of the matrices are opaque (i.e., they can be arbitrary), you can
- * stack products of three or more matrices by making one of the two matrices
- * an object of the current type handles be a ProductMatrix itself.
- *
- * Here is an example multiplying two different FullMatrix objects:
- * @include product_matrix.cc
- *
- * @deprecated If deal.II was configured with C++11 support, use the
- * LinearOperator class instead, see the module on
- * @ref LAOperators "linear operators"
- * for further details.
- *
- * @author Guido Kanschat, 2000, 2001, 2002, 2005
- */
-template<typename VectorType>
-class ProductMatrix : public PointerMatrixBase<VectorType>
-{
-public:
- /**
- * Standard constructor. Matrices and the memory pool must be added later
- * using initialize().
- */
- ProductMatrix();
-
- /**
- * Constructor only assigning the memory pool. Matrices must be added by
- * reinit() later.
- */
- ProductMatrix(VectorMemory<VectorType> &mem);
-
- /**
- * Constructor. Additionally to the two constituting matrices, a memory
- * pool for the auxiliary vector must be provided.
- */
- template <typename MatrixType1, typename MatrixType2>
- ProductMatrix (const MatrixType1 &m1,
- const MatrixType2 &m2,
- VectorMemory<VectorType> &mem);
-
- /**
- * Destructor.
- */
- ~ProductMatrix();
-
- /**
- * Change the matrices.
- */
- template <typename MatrixType1, typename MatrixType2>
- void reinit (const MatrixType1 &m1, const MatrixType2 &m2);
-
- /**
- * Change the matrices and memory pool.
- */
- template <typename MatrixType1, typename MatrixType2>
- void initialize (const MatrixType1 &m1,
- const MatrixType2 &m2,
- VectorMemory<VectorType> &mem);
-
- // Doc in PointerMatrixBase
- void clear();
-
- /**
- * Matrix-vector product <i>w = m1 * m2 * v</i>.
- */
- virtual void vmult (VectorType &w,
- const VectorType &v) const;
-
- /**
- * Transposed 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, transposed 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.
- */
- PointerMatrixBase<VectorType> *m1;
-
- /**
- * The right matrix of the product.
- */
- PointerMatrixBase<VectorType> *m2;
-
- /**
- * Memory for auxiliary vector.
- */
- SmartPointer<VectorMemory<VectorType>,ProductMatrix<VectorType> > mem;
-};
-
-
/**
* A matrix that is the multiple of another matrix.
*
}
-//---------------------------------------------------------------------------
-
-template<typename VectorType>
-ProductMatrix<VectorType>::ProductMatrix ()
- : m1(0), m2(0), mem(0)
-{}
-
-
-template<typename VectorType>
-ProductMatrix<VectorType>::ProductMatrix (VectorMemory<VectorType> &m)
- : m1(0), m2(0), mem(&m)
-{}
-
-
-template<typename VectorType>
-template<typename MatrixType1, typename MatrixType2>
-ProductMatrix<VectorType>::ProductMatrix (const MatrixType1 &mat1,
- const MatrixType2 &mat2,
- VectorMemory<VectorType> &m)
- : mem(&m)
-{
- m1 = new PointerMatrix<MatrixType1, VectorType>(&mat1, typeid(*this).name());
- m2 = new PointerMatrix<MatrixType2, VectorType>(&mat2, typeid(*this).name());
-}
-
-
-template<typename VectorType>
-template<typename MatrixType1, typename MatrixType2>
-void
-ProductMatrix<VectorType>::reinit (const MatrixType1 &mat1, const MatrixType2 &mat2)
-{
- if (m1) delete m1;
- if (m2) delete m2;
- m1 = new PointerMatrix<MatrixType1, VectorType>(&mat1, typeid(*this).name());
- m2 = new PointerMatrix<MatrixType2, VectorType>(&mat2, typeid(*this).name());
-}
-
-
-template<typename VectorType>
-template<typename MatrixType1, typename MatrixType2>
-void
-ProductMatrix<VectorType>::initialize (const MatrixType1 &mat1,
- const MatrixType2 &mat2,
- VectorMemory<VectorType> &memory)
-{
- mem = &memory;
- if (m1) delete m1;
- if (m2) delete m2;
- m1 = new PointerMatrix<MatrixType1, VectorType>(&mat1, typeid(*this).name());
- m2 = new PointerMatrix<MatrixType2, VectorType>(&mat2, typeid(*this).name());
-}
-
-
-template<typename VectorType>
-ProductMatrix<VectorType>::~ProductMatrix ()
-{
- if (m1) delete m1;
- if (m2) delete m2;
-}
-
-
-template<typename VectorType>
-void
-ProductMatrix<VectorType>::clear ()
-{
- if (m1) delete m1;
- m1 = 0;
- if (m2) delete m2;
- m2 = 0;
-}
-
-
-template<typename VectorType>
-void
-ProductMatrix<VectorType>::vmult (VectorType &dst, const VectorType &src) const
-{
- Assert (mem != 0, ExcNotInitialized());
- Assert (m1 != 0, ExcNotInitialized());
- Assert (m2 != 0, ExcNotInitialized());
-
- VectorType *v = mem->alloc();
- v->reinit(dst);
- m2->vmult (*v, src);
- m1->vmult (dst, *v);
- mem->free(v);
-}
-
-
-template<typename VectorType>
-void
-ProductMatrix<VectorType>::vmult_add (VectorType &dst, const VectorType &src) const
-{
- Assert (mem != 0, ExcNotInitialized());
- Assert (m1 != 0, ExcNotInitialized());
- Assert (m2 != 0, ExcNotInitialized());
-
- VectorType *v = mem->alloc();
- v->reinit(dst);
- m2->vmult (*v, src);
- m1->vmult_add (dst, *v);
- mem->free(v);
-}
-
-
-template<typename VectorType>
-void
-ProductMatrix<VectorType>::Tvmult (VectorType &dst, const VectorType &src) const
-{
- Assert (mem != 0, ExcNotInitialized());
- Assert (m1 != 0, ExcNotInitialized());
- Assert (m2 != 0, ExcNotInitialized());
-
- VectorType *v = mem->alloc();
- v->reinit(dst);
- m1->Tvmult (*v, src);
- m2->Tvmult (dst, *v);
- mem->free(v);
-}
-
-
-template<typename VectorType>
-void
-ProductMatrix<VectorType>::Tvmult_add (VectorType &dst, const VectorType &src) const
-{
- Assert (mem != 0, ExcNotInitialized());
- Assert (m1 != 0, ExcNotInitialized());
- Assert (m2 != 0, ExcNotInitialized());
-
- VectorType *v = mem->alloc();
- v->reinit(dst);
- m1->Tvmult (*v, src);
- m2->Tvmult_add (dst, *v);
- mem->free(v);
-}
-
-
-
//---------------------------------------------------------------------------
template <typename VectorType>