algebra/qs_matrix.h

// -*- c++ -*-

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

#ifndef _WAVELETTL_QS_MATRIX_H
#define _WAVELETTL_QS_MATRIX_H

#include <cmath>
#include <map>
#include <algebra/matrix.h>
#include <algebra/sparse_matrix.h>

namespace MathTL
{
  template <class C>
  class QuasiStationaryMatrix
  {
  public:
    typedef typename Vector<C>::size_type size_type;

    QuasiStationaryMatrix();

    QuasiStationaryMatrix(const int j0,
                          const size_type mj0, const size_type nj0,
                          const Matrix<C>& ML, const Matrix<C>& MR,
                          const Vector<C>& bandL, const Vector<C>& bandR,
                          const size_type offsetL, const size_type offsetR,
                          const double factor = 1.0);

    QuasiStationaryMatrix<C>& operator = (const QuasiStationaryMatrix<C>& M);
    
    const size_type row_dimension() const;

    const size_type column_dimension() const;

    const size_type ML_column_dimension() const { return ML_.column_dimension(); }

    const size_type MR_column_dimension() const { return MR_.column_dimension(); }

    void set_level(const int j) const;

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

    void to_sparse(SparseMatrix<C>& S) const;

    template <class VECTOR>
    void apply(const VECTOR& x, VECTOR& Mx,
               const size_type x_offset = 0,
               const size_type Mx_offset = 0,
               const bool add_to = false) const;

    void apply(const std::map<size_type, C>& x, std::map<size_type, C>& Mx,
               const size_type x_offset = 0,
               const size_type Mx_offset = 0,
               const bool add_to = false) const;
    
    template <class VECTOR>
    void apply_transposed(const VECTOR& x, VECTOR& Mtx,
                          const size_type x_offset = 0,
                          const size_type Mtx_offset = 0,
                          const bool add_to = false) const;
    
    void apply_transposed(const std::map<size_type, C>& x, std::map<size_type, C>& Mtx,
                          const size_type x_offset = 0,
                          const size_type Mtx_offset = 0,
                          const bool add_to = false) const;
    
    void apply_central_block
    (const std::map<size_type, C>& x, std::map<size_type, C>& Mx,
     const size_type x_offset = 0,
     const size_type Mx_offset = 0,
     const bool add_to = false) const;

    void apply_central_block_transposed
    (const std::map<size_type, C>& x, std::map<size_type, C>& Mtx,
     const size_type x_offset = 0,
     const size_type Mtx_offset = 0,
     const bool add_to = false) const;

    void apply_central_columns
    (const std::map<size_type, C>& x, std::map<size_type, C>& Mx,
     const size_type x_offset = 0,
     const size_type Mx_offset = 0,
     const bool add_to = false) const;

    void apply_central_columns_transposed
    (const std::map<size_type, C>& x, std::map<size_type, C>& Mtx,
     const size_type x_offset = 0,
     const size_type Mtx_offset = 0,
     const bool add_to = false) const;
    
    void print(std::ostream& os,
               const unsigned int tabwidth = 8,
               const unsigned int precision = 3) const;

  protected:
    int j0_;
    mutable int j_;
    size_type mj0_, nj0_;
    Matrix<C> ML_, MR_;
    Vector<C> bandL_,bandR_;
    size_type offsetL_, offsetR_;
    double factor_;
  };

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


  template <class C>
  class PeriodicQuasiStationaryMatrix
  {
  public:
    typedef typename Vector<C>::size_type size_type;

    PeriodicQuasiStationaryMatrix(const int j0,
                                  const int offset,
                                  const Array1D<C>& band,
                                  const double factor = 1.0);

    PeriodicQuasiStationaryMatrix<C>& operator = (const PeriodicQuasiStationaryMatrix<C>& M);
    
    const size_type row_dimension() const { return 1<<(j_+1); }

    const size_type column_dimension() const { return 1<<j_; }

    void set_level(const int j) const;

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

    void to_sparse(SparseMatrix<C>& S) const;

    template <class VECTOR>
    void apply(const VECTOR& x, VECTOR& Mx,
               const size_type x_offset = 0,
               const size_type Mx_offset = 0,
               const bool add_to = false) const;

    template <class VECTOR>
    void apply_transposed(const VECTOR& x, VECTOR& Mtx,
                          const size_type x_offset = 0,
                          const size_type Mtx_offset = 0,
                          const bool add_to = false) const;
    
    void print(std::ostream& os,
               const unsigned int tabwidth = 8,
               const unsigned int precision = 3) const;

  protected:
    int j0_, offset_;
    mutable int j_;
    Array1D<C> band_;
    double factor_;
  };

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

}

#include <algebra/qs_matrix.cpp>

#endif