SparseMatrix_Eigen.hpp

/*
 * SparseMatrix_Eigen.hpp
 *
 *  Created on: Nov 27, 2015
 *      Author: i-bird
 */
 
#ifndef OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_EIGEN_HPP_
#define OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_EIGEN_HPP_
 
#include "Vector/map_vector.hpp"
#include <boost/mpl/int.hpp>
#include "VCluster/VCluster.hpp"
 
#define EIGEN_TRIPLET 1
 
template<typename T>
class triplet<T,EIGEN_TRIPLET>
{
    Eigen::Triplet<T,long int> triplet_ei;
 
public:
 
    long int & row()
    {
        size_t ptr = (size_t)&triplet_ei.row();
 
        return (long int &)*(long int *)ptr;
    }
 
    long int & col()
    {
        size_t ptr = (size_t)&triplet_ei.col();
 
        return (long int &)*(long int *)ptr;
    }
 
    T & value()
    {
        size_t ptr = (size_t)&triplet_ei.value();
 
        return (T &)*(T *)ptr;
    }
 
    triplet(long int i, long int j, T val)
    {
        row() = i;
        col() = j;
        value() = val;
    }
 
    // Default constructor
    triplet()   {};
};
 
/* ! \brief Sparse Matrix implementation, that map over Eigen
 *
 * \tparam T Type of the sparse Matrix store on each row,colums
 * \tparam id_t type of id
 * \tparam impl implementation
 *
 */
template<typename T, typename id_t>
class SparseMatrix<T,id_t,EIGEN_BASE>
{
public:
 
    typedef boost::mpl::int_<EIGEN_BASE> triplet_impl;
 
    typedef triplet<T,EIGEN_BASE> triplet_type;
 
private:
 
    Eigen::SparseMatrix<T,0,id_t> mat;
    openfpm::vector<triplet_type> trpl;
    openfpm::vector<triplet_type> trpl_recv;
 
    void assemble()
    {
        Vcluster<> & vcl = create_vcluster();
 
        // we assemble the Matrix from the collected data
        if (vcl.getProcessingUnits() != 1)
        {
            collect();
            // only master assemble the Matrix
            if (vcl.getProcessUnitID() == 0)
                mat.setFromTriplets(trpl_recv.begin(),trpl_recv.end());
        }
        else
            mat.setFromTriplets(trpl.begin(),trpl.end());
    }
 
    void collect()
    {
        Vcluster<> & vcl = create_vcluster();
 
        trpl_recv.clear();
 
        // here we collect all the triplet in one array on the root node
        vcl.SGather(trpl,trpl_recv,0);
    }
 
public:
 
 
 
    SparseMatrix(size_t N1, size_t N2)
    :mat(N1,N2)
    {
    }
 
 
    SparseMatrix()  {}
 
    openfpm::vector<triplet_type> & getMatrixTriplets()
    {
        return this->trpl;
    }
 
    const Eigen::SparseMatrix<T,0,id_t> & getMat() const
    {
        // Here we collect the information on master
        assemble();
 
        return mat;
    }
 
    Eigen::SparseMatrix<T,0,id_t> & getMat()
    {
        assemble();
 
        return mat;
    }
 
    void resize(size_t row, size_t col,size_t l_row, size_t l_col)
    {
        mat.resize(row,col);
    }
 
    T operator()(id_t i, id_t j)
    {
        return mat.coeffRef(i,j);
    }
 
    bool save(const std::string & file) const
    {
        size_t pap_prp;
 
        Packer<decltype(trpl),HeapMemory>::packRequest(trpl,pap_prp);
 
        // allocate the memory
        HeapMemory pmem;
        pmem.allocate(pap_prp);
        ExtPreAlloc<HeapMemory> mem(pap_prp,pmem);
 
        //Packing
 
        Pack_stat sts;
        Packer<decltype(trpl),HeapMemory>::pack(mem,trpl,sts);
 
        // Save into a binary file
        std::ofstream dump (file, std::ios::out | std::ios::binary);
        if (dump.is_open() == false)
            return false;
        dump.write ((const char *)pmem.getPointer(), pmem.size());
 
        return true;
    }
 
    bool load(const std::string & file)
    {
        std::ifstream fs (file, std::ios::in | std::ios::binary | std::ios::ate );
        if (fs.is_open() == false)
            return false;
 
        // take the size of the file
        size_t sz = fs.tellg();
 
        fs.close();
 
        // reopen the file without ios::ate to read
        std::ifstream input (file, std::ios::in | std::ios::binary );
        if (input.is_open() == false)
            return false;
 
        // Create the HeapMemory and the ExtPreAlloc memory
        HeapMemory pmem;
        ExtPreAlloc<HeapMemory> mem(sz,pmem);
 
        // read
        input.read((char *)pmem.getPointer(), sz);
 
        //close the file
        input.close();
 
        //Unpacking
        Unpack_stat ps;
 
        Unpacker<decltype(trpl),HeapMemory>::unpack(mem,trpl,ps);
 
        return true;
    }
 
    T getValue(size_t r, size_t c)
    {
        for (size_t i = 0 ; i < trpl.size() ; i++)
        {
            if (r == (size_t)trpl.get(i).row() && c == (size_t)trpl.get(i).col())
                return trpl.get(i).value();
        }
 
        return 0;
    }
};
 
 
 
#endif /* OPENFPM_NUMERICS_SRC_MATRIX_SPARSEMATRIX_EIGEN_HPP_ */