Eigen.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/Misc/Lazy.h"

namespace honey
{


template<class Real>
class Eigen
{
    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, Real>      Vec;

    Eigen()                                                     { init(); }
    template<class T>
    Eigen(const MatrixBase<T>& a)                               { init(); calc(a); }

    template<class T>
    Eigen& calc(const MatrixBase<T>& a)
    {
        assert(a.rows() == a.cols(), "Matrix must be square and symmetric");
        sdt n = a.rows();
        _w.resize(n); _indR.resize(n); _indC.resize(n);
        _v.resize(n, n);
        _a = a;
        jacobi(_a, _w, _v, _indR, _indC);
        _vt.setDirty(true);
        return *this;
    }

    template<class T>
    Eigen& calcValues(const MatrixBase<T>& a)
    {
        assert(a.rows() == a.cols(), "Matrix must be square and symmetric");
        sdt n = a.rows();
        _w.resize(n); _indR.resize(n); _indC.resize(n);
        _a = a;
        jacobi(_a, _w, optnull, _indR, _indC);
        return *this;
    }

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

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

    const Vec& w() const                                        { return _w; }
    const Matrix& v() const                                     { return _v; }
    const Matrix& vt() const                                    { return _vt; }

private:
    void init()                                                 { _vt.setEval([&](Matrix& val) { _v.transpose(val); }); }

    static void jacobi(Matrix& A, Vec& W, optional<Matrix&> V, vector<sdt>& indR, vector<sdt>& indC);

    Vec             _w;
    Matrix          _v;
    lazy<Matrix>    _vt;
    Matrix          _a;
    vector<sdt>    _indR;
    vector<sdt>    _indC;
    BackSub<Real>   _backSub;
};

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

}