Qrd.h

Go to the documentation of this file.

// Honeycomb, Copyright (C) 2015 NewGamePlus Inc.  Distributed under the Boost Software License v1.0.
#pragma once

#include "Honey/Math/NumAnalysis/BackSub.h"

namespace honey
{


template<class Real>
class Qrd
{
    typedef typename Numeral<Real>::Real_   Real_;
    typedef typename Real_::DoubleType      Double_;
    typedef typename Double_::Real          Double;
    typedef Alge_<Real>                     Alge;
    typedef Alge_<Double>                   Alge_d;
    typedef BackSub<Real>                   BackSub;

public:
    typedef Matrix<matrix::dynamic, matrix::dynamic, Real> Matrix;

    enum class Mode
    {
        full,       
        reduced     
    };

    Qrd() {}
    template<class T>
    Qrd(const MatrixBase<T>& a, Mode mode = Mode::reduced)      { calc(a, mode); }

    template<class T>
    Qrd& calc(const MatrixBase<T>& a, Mode mode = Mode::reduced)
    {
        sdt m = a.rows(), n = a.cols();
        assert(m >= n);
        _h = a;
        _q.resize(m, mode == Mode::full ? m : n);
        _r.resize(n, n);
        householder(_h, _q, _r);
        return *this;
    }

    bool isFullRank() const                                     { return BackSub::isFullRank(_r); }

    template<class B, class X>
    void solve(const MatrixBase<B>& b, MatrixBase<X>& x)
    {
        sdt m = _q.rows(), n = _r.rows();
        bool full = m == _q.cols();
        if (full)   _q.block(0,0,m,n).transposeMul(b, _y);
        else        _q.transposeMul(b, _y);
        BackSub::solve(_r, _y, x);
    }

    const Matrix& q() const                                     { return _q; }
    const Matrix& r() const                                     { return _r; }
    const Matrix& h() const                                     { return _h; }

private:
    static void householder(Matrix& a, Matrix& q, Matrix& r);

    Matrix          _q;
    Matrix          _r;
    Matrix          _h;
    Matrix          _y;
};

extern template class Qrd<Float>;
extern template class Qrd<Double>;
extern template class Qrd<Quad>;

}