Svd.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"
#include "Honey/Math/Random/Chacha.h"

namespace honey
{


template<class Real>
class Svd
{
    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;
public:
    typedef Matrix<matrix::dynamic, matrix::dynamic, Real> Matrix;
    typedef Vec<matrix::dynamic, Double>    Vec_d;
    typedef Vec<matrix::dynamic, Real>      Vec;

    enum class Mode
    {
        full,       
        reduced     
    };

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

    template<class T>
    Svd& calc(const MatrixBase<T>& a, Mode mode = Mode::reduced)
    {
        sdt m = a.rows(), n = a.cols();
        _w.resize(Alge::min(m, n));
        _wd.resize(_w.size());

        if (m >= n)
        {
            _vt.resize(n, n);
            if (mode == Mode::full) a.transpose(_ut.resize(m,m).block(0,0,n,m));
            else                    a.transpose(_ut);
            jacobi(_ut, _w, _wd, _vt, _rand);
            _ut.transpose(_u);
        }
        else
        {
            //m < n
            _u.resize(m, m);
            if (mode == Mode::full) _vt.resize(n,n).block(0,0,m,n) = a;
            else                    _vt = a;
            jacobi(_vt, _w, _wd, _u, _rand);
            _u.transposeInPlace();
        }

        return *this;
    }

    template<class T>
    Svd& calcValues(const MatrixBase<T>& a)
    {
        sdt m = a.rows(), n = a.cols();
        _w.resize(Alge::min(m, n));
        _wd.resize(_w.size());

        if (m >= n)
        {
            a.transpose(_ut);
            jacobi(_ut, _w, _wd, optnull, _rand);
        }
        else
        {
            _vt = a;
            jacobi(_vt, _w, _wd, optnull, _rand);
        }

        return *this;
    }

    template<class B, class X>
    void solve(const MatrixBase<B>& b, MatrixBase<X>& x)        { _backSub.solve(_w, _u, _vt, b, x); }

    template<class T>
    void inverse(MatrixBase<T>& res)                            { _backSub.solve(_w, _u, _vt, res); }

    const Vec& w() const                                        { return _w; }
    const Matrix& u() const                                     { return _u; }
    const Matrix& vt() const                                    { return _vt; }

private:
    static void jacobi(Matrix& At, Vec& w, Vec_d& wd, optional<Matrix&> Vt, RandomGen& rand);

    Vec             _w;
    Vec_d           _wd;
    Matrix          _u;
    Matrix          _ut;
    Matrix          _vt;
    Chacha          _rand;
    BackSub<Real>   _backSub;
};

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

}