algebra/kronecker_matrix.h

// -*- c++ -*-

// +--------------------------------------------------------------------+
// | This file is part of MathTL - the Mathematical Template Library    |
// |                                                                    |
// | Copyright (c) 2002-2009                                            |
// | Thorsten Raasch, Manuel Werner                                     |
// +--------------------------------------------------------------------+

#ifndef _MATHTL_KRONECKER_MATRIX_H
#define _MATHTL_KRONECKER_MATRIX_H

#include <iostream>
#include <algebra/vector.h>
#include <algebra/matrix_block.h>

// matrix norms, for convenience
#include <algebra/matrix_norms.h>

namespace MathTL
{
  template <class C, class MATRIX1, class MATRIX2>
  class KroneckerMatrix
    : public MatrixBlock<C>
  {
  public:
    typedef typename Vector<C>::size_type size_type;
    
    explicit KroneckerMatrix(const MATRIX1& A, const MATRIX2& B,
                             const double factor = 1.0);
    
    KroneckerMatrix(const KroneckerMatrix<C,MATRIX1,MATRIX2>& M);

    MatrixBlock<C>* clone() const;

    MatrixBlock<C>* clone_transposed() const;

    const size_type row_dimension() const;

    const size_type column_dimension() const;

    bool empty() const;

    const C operator () (const size_type row, const size_type column) const;

    const C get_entry(const size_type row, const size_type column) const;

    template <class VECTOR>
    void apply(const VECTOR& x, VECTOR& Mx) const;

    void apply(const Vector<C>& x, Vector<C>& Mx) const;

    template <class VECTOR>
    void apply_transposed(const VECTOR& x, VECTOR& Mtx) const;

    void apply_transposed(const Vector<C>& x, Vector<C>& Mtx) const;
    
    void print(std::ostream& os,
               const unsigned int tabwidth = 10,
               const unsigned int precision = 3) const;

  protected:
    MATRIX1 A;
    MATRIX2 B;
    C factor_;
  };

  template <class C, class MATRIX1, class MATRIX2>
  std::ostream& operator << (std::ostream& os, const KroneckerMatrix<C,MATRIX1,MATRIX2>& M);

  
  template <class C, class MATRIX1, class MATRIX2>
  class KroneckerHelper
  {
  public:
    typedef typename Vector<C>::size_type size_type;
    
    explicit KroneckerHelper(const MATRIX1& A, const MATRIX2& B,
                             const double factor = 1.0);
    
    KroneckerHelper(const KroneckerHelper<C,MATRIX1,MATRIX2>& M);

    const size_type row_dimension() const;

    const size_type column_dimension() const;

    void apply(const Vector<C>& x, Vector<C>& Mx,
               const size_type x_offset,
               const size_type Mx_offset,
               const bool add_to) const;
    
    void apply_transposed(const Vector<C>& x, Vector<C>& Mtx,
                          const size_type x_offset,
                          const size_type Mtx_offset,
                          const bool add_to) const;

  protected:
    MATRIX1 A;
    MATRIX2 B;
    C factor_;
  };  
}

// include implementation of inline functions
#include <algebra/kronecker_matrix.cpp>

#endif