matrix.hpp

#ifndef MPQC_MATRIX_HPP
#define MPQC_MATRIX_HPP
 
#include "mpqc/range.hpp"
#include "math/scmat/matrix.h"
 
#define EIGEN_DONT_PARALLELIZE
 
#include <Eigen/Eigen>
// #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
// #include <Eigen/Sparse>
 
namespace mpqc {
 
    
 
    template<typename T, int Order = Eigen::ColMajor>
    struct matrix :
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Order>
    {
        typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Order> EigenType;
 
    matrix() : EigenType() {}
 
    matrix(size_t m, size_t n) : EigenType(m,n) {}
 
        template<class A>
    matrix(const Eigen::EigenBase<A> &a) : EigenType(a) {}
 
        matrix(sc::RefSCMatrix a) {
            this->resize(a.nrow(), a.ncol());
            apply(assign(), this->rows(), this->cols(), *this, a);
        }
 
        matrix(const sc::RefSymmSCMatrix &a) {
            this->resize(a.n(), a.n());
            apply(assign(), this->rows(), this->cols(), *this, a);
        }
 
        struct Assignable {
            virtual ~Assignable() {}
            virtual void assign_to(matrix &m) const = 0;
        };
 
        matrix(const Assignable &a) {
            a.assign_to(*this);
        }
 
#ifdef DOXYGEN
 
    Type operator()(i, j);
 
#else // DOXYGEN
 
        using EigenType::operator();
 
        template<typename T_, typename U_>
        struct is_index : boost::mpl::and_<
            boost::is_integral<T_>,
            boost::is_integral<U_> > {};
 
        template<class Ri, class Rj>
        typename boost::disable_if<
            is_index<Ri, Rj>,
            Eigen::Block<EigenType>
            >::type 
        operator()(const Ri &i, const Rj &j) {
            range ri = range_cast(i);
            range rj = range_cast(j);
            // printf("i=(%i,%i),j=(%i:%i)\n",
            //        *ri.begin(), *ri.end(),
            //        *rj.begin(), *rj.end());
        return this->block(*ri.begin(), *rj.begin(), ri.size(), rj.size());
        }
 
        template<class Ri, class Rj>
        typename boost::disable_if<
            is_index<Ri, Rj>,
            Eigen::Block<const EigenType>
            >::type 
        operator()(const Ri &i, const Rj &j) const {
            range ri = range_cast(i);
            range rj = range_cast(j);
            // printf("i=(%i,%i),j=(%i:%i)\n",
            //        *i.begin(), *i.end(),
            //        *j.begin(), *j.end());
        return this->block(*ri.begin(), *rj.begin(), ri.size(), rj.size());
        }
 
#endif // DOXYGEN
 
        void reshape(int m, int n);
 
    private:
        template<class F, class A, class B>
        static void apply(const F &f, size_t m, size_t n, A &a, B &b) {
        for (int j = 0; j < n; ++j) {
        for (int i = 0; i < m; ++i) {
            f(a(i, j), b(i, j));
        }
        }
        }
        struct assign {
            template<class A, class B>
            void operator()(A &a, const B &b) const {
        a = b;
            }
        };
    };
 
    template<typename T>
    struct vector: Eigen::Matrix<T, Eigen::Dynamic, 1> {
 
        typedef Eigen::Matrix<T, Eigen::Dynamic, 1> EigenType;
 
    explicit vector(size_t m = 0) : EigenType(m) {}
 
        template<class A>
    vector(const Eigen::EigenBase<A> &a) : EigenType(a) {}
 
        vector(const T *begin, const T *end) {
        EigenType::resize(end - begin, 1);
        std::copy(begin, end, EigenType::data());
        }
 
        using EigenType::operator();
 
        Eigen::Block<EigenType> operator()(range i) {
        return this->block(*i.begin(), 0, i.size(), 1);
        }
 
        Eigen::Block<const EigenType> operator()(range i) const {
        return this->block(*i.begin(), 0, i.size(), 1);
        }
 
    };
 
 
    typedef matrix<double> Matrix;
    typedef vector<double> Vector;
 
    // typedef Eigen::SparseMatrix<double> Sparse;
 
    template<class E>
    double absmax(const E &e) {
        return std::max(fabs(e.maxCoeff()), fabs(e.minCoeff()));
    }
 
    template<class T>
    T dot(const matrix<T> &a, const matrix<T> &b) {
    T q = 0;
    MPQC_ASSERT(a.cols() == b.cols());
    for (size_t j = 0; j < a.cols(); ++j) {
        q += a.col(j).dot(b.col(j));
    }
    return q;
    }
 
    template<class T>
    Eigen::SelfAdjointEigenSolver<Matrix::EigenType> symmetric(const matrix<T> &a) {
    Eigen::SelfAdjointEigenSolver<Matrix::EigenType> es(a);
    if (es.info() != Eigen::Success)
        throw std::runtime_error("Eigen solver failed");
    return es;
    }
 
    template<class T>
    T norm(const matrix<T> &a) {
    return a.norm();
    }
 
    template<class T>
    void normalize(matrix<T> &a) {
    a *= 1/a.norm();
    }
 
    template<class T>
    void orthonormalize(matrix<T> &d, const matrix<T> &b) {
    T db = dot(d, b);
    T bb = 1; //dot(b,b); // should be normalized already
    d = d - (db/bb)*b;
    normalize(d);
        //d /= sqrt(d.norm());
    }
 
 
} // namespace mpqc
 
namespace sc {
    using mpqc::Vector;
    using mpqc::Matrix;
    using mpqc::dot;
}
 
#endif /* MPQC_MATRIX_HPP */