numerics/decomposable_matrix.h

// -*- c++ -*-

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

#ifndef _MATHTL_DECOMPOSABLE_MATRIX_H
#define _MATHTL_DECOMPOSABLE_MATRIX_H

#include <iostream>
#include <algebra/vector.h>
#include <algebra/matrix.h>
#include <utils/array1d.h>

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

namespace MathTL
{
  template <class C>
  class DecomposableMatrix
    : protected Matrix<C>
  {
  public:
    typedef typename Vector<C>::size_type size_type;

    enum DecompositionType {
      none,
      LU,
      QU
    };

    explicit DecomposableMatrix(const size_type n = 0);

    DecomposableMatrix(const DecomposableMatrix<C>& M);

    DecomposableMatrix(const Matrix<C>& M);

    DecomposableMatrix<C>& operator = (const DecomposableMatrix<C>& M);

    const size_type row_dimension() const;

    const size_type column_dimension() const;

    const size_type size() const;

    void decompose(DecompositionType d = LU);

    void solve(const Vector<C>& b, Vector<C>& x,
               DecompositionType d = LU);

    void print(std::ostream& os,
               const unsigned int tabwidth = 10,
               const unsigned int precision = 3) const;

 protected:
    // decomposition type
    DecompositionType decomposition;

    // storage for row equilibration matrix entries (-> LU factorization)
    Vector<double> D;

    // storage for permuation matrix (-> LU factorization)
    Array1D<size_type> P;

    // LU decomposition (row equilibration, row pivoting, in place (up to D and P))
    void LU_decomposition();

    // QR decomposition (in place)
    void QU_decomposition();

    // revert decomposition
    void revert_decomposition();
  };

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

// include implementation of inline functions
#include <numerics/decomposable_matrix.cpp>

#endif