MetisDistribution.hpp

/*
 * MetisDistribution.hpp
 *
 *  Created on: Nov 19, 2015
 *      Author: Antonio Leo
 */
 
#ifndef SRC_DECOMPOSITION_METISDISTRIBUTION_HPP_
#define SRC_DECOMPOSITION_METISDISTRIBUTION_HPP_
 
#include "SubdomainGraphNodes.hpp"
#include "metis_util.hpp"
 
#define METIS_DISTRIBUTION_ERROR_OBJECT std::runtime_error("Metis runtime error");
 
struct met_sub_w
{
    size_t id;
 
    size_t w;
 
    static bool noPointers() {return true;}
};
 
template<unsigned int dim, typename T>
class MetisDistribution
{
    Vcluster<> & v_cl;
 
    grid_sm<dim, void> gr;
 
    Box<dim, T> domain;
 
    Graph_CSR<nm_v<dim>, nm_e> gp;
 
    bool testing = false;
 
    Metis<Graph_CSR<nm_v<dim>, nm_e>> metis_graph;
 
    std::unordered_map<size_t,size_t> owner_scs;
 
    openfpm::vector<met_sub_w> owner_cost_sub;
 
    openfpm::vector<met_sub_w> recv_ass;
 
    inline void check_overflow(size_t id)
    {
#ifdef SE_CLASS1
        if (id >= gp.getNVertex())
        {
            std::cerr << "Error " << __FILE__ ":" << __LINE__ << " such sub-sub-domain doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
            ACTION_ON_ERROR(METIS_DISTRIBUTION_ERROR_OBJECT)
        }
#endif
    }
 
    inline void check_overflowe(size_t id, size_t e)
    {
#ifdef SE_CLASS1
        if (e >= gp.getNChilds(id))
        {
            std::cerr << "Error " << __FILE__ ":" << __LINE__ << " for the sub-sub-domain " << id << " such neighborhood doesn't exist (e = " << e << ", " << "total size = " << gp.getNChilds(id) << ")\n";
            ACTION_ON_ERROR(METIS_DISTRIBUTION_ERROR_OBJECT)
        }
#endif
    }
 
public:
 
    static constexpr unsigned int computation = nm_v_computation;
 
    MetisDistribution(Vcluster<> & v_cl)
    :v_cl(v_cl),metis_graph(gp)
    {
#ifdef SE_CLASS2
            check_new(this,8,VECTOR_EVENT,1);
#endif
    }
 
    MetisDistribution(const MetisDistribution & mt)
    :v_cl(mt.v_cl)
    {
#ifdef SE_CLASS2
            check_valid(mt);
            check_new(this,8,VECTOR_EVENT,1);
#endif
        this->operator=(mt);
    }
 
    MetisDistribution(MetisDistribution && mt)
    {
#ifdef SE_CLASS2
            check_valid(mt);
            check_new(this,8,VECTOR_EVENT,1);
#endif
        this->operator=(mt);
    }
 
    ~MetisDistribution()
    {
#ifdef SE_CLASS2
        check_delete(this);
#endif
    }
 
 
    void createCartGraph(grid_sm<dim, void> & grid, Box<dim, T> dom)
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        // NON periodic boundary conditions
        size_t bc[dim];
 
        for (size_t i = 0 ; i < dim ; i++)
            bc[i] = NON_PERIODIC;
 
        // Set grid and domain
        gr = grid;
        domain = dom;
 
        // Create a cartesian grid graph
        CartesianGraphFactory<dim, Graph_CSR<nm_v<dim>, nm_e>> g_factory_part;
        gp = g_factory_part.template construct<NO_EDGE, nm_v_id, T, dim - 1, 0>(gr.getSize(), domain, bc);
 
        // Init to 0.0 axis z (to fix in graphFactory)
        if (dim < 3)
        {
            for (size_t i = 0; i < gp.getNVertex(); i++)
            {
                gp.vertex(i).template get<nm_v_x>()[2] = 0.0;
            }
        }
 
        for (size_t i = 0; i < gp.getNVertex(); i++)
            gp.vertex(i).template get<nm_v_global_id>() = i;
    }
 
    Graph_CSR<nm_v<dim>, nm_e> & getGraph()
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        return gp;
    }
 
    void decompose()
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
 
        // Gather the sub-domain weight in one processor
        recv_ass.clear();
        v_cl.SGather(owner_cost_sub,recv_ass,0);
 
        if (v_cl.getProcessUnitID() == 0)
        {
            if (recv_ass.size() != 0)
            {
                // we fill the assignment
                for (size_t i = 0 ; i < recv_ass.size() ; i++)
                    gp.template vertex_p<nm_v_computation>(recv_ass.get(i).id) = recv_ass.get(i).w;
 
                metis_graph.initMetisGraph(v_cl.getProcessingUnits(),true);
            }
            else
                metis_graph.initMetisGraph(v_cl.getProcessingUnits(),false);
            metis_graph.onTest(testing);
 
            // decompose
            metis_graph.template decompose<nm_v_proc_id>();
 
            if (recv_ass.size() != 0)
            {
                // we fill the assignment
                for (size_t i = 0 ; i < recv_ass.size() ; i++)
                    recv_ass.get(i).w = gp.template vertex_p<nm_v_proc_id>(recv_ass.get(i).id);
            }
            else
            {
                recv_ass.resize(gp.getNVertex());
 
                // we fill the assignment
                for (size_t i = 0 ; i < gp.getNVertex() ; i++)
                {
                    recv_ass.get(i).id = i;
                    recv_ass.get(i).w = gp.template vertex_p<nm_v_proc_id>(i);
                }
            }
        }
        else
        {
            metis_graph.inc_dec();
        }
 
        recv_ass.resize(gp.getNVertex());
 
        // broad cast the result
        v_cl.Bcast(recv_ass,0);
        v_cl.execute();
        owner_scs.clear();
        owner_cost_sub.clear();
 
        size_t j = 0;
 
        // Fill the metis graph
        for (size_t i = 0 ; i < recv_ass.size() ; i++)
        {
            gp.template vertex_p<nm_v_proc_id>(recv_ass.get(i).id) = recv_ass.get(i).w;
 
            if (recv_ass.get(i).w == v_cl.getProcessUnitID())
            {
                owner_scs[recv_ass.get(i).id] = j;
                j++;
                owner_cost_sub.add();
                owner_cost_sub.last().id = recv_ass.get(i).id;
                owner_cost_sub.last().w = 1;
            }
        }
    }
 
    void refine()
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
 
        decompose();
    }
 
    void redecompose()
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        decompose();
    }
 
 
    void getSSDomainPos(size_t id, T (&pos)[dim])
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        check_overflow(id);
 
        // Copy the geometrical informations inside the pos vector
        pos[0] = gp.vertex(id).template get<nm_v_x>()[0];
        pos[1] = gp.vertex(id).template get<nm_v_x>()[1];
        if (dim == 3)
        {pos[2] = gp.vertex(id).template get<nm_v_x>()[2];}
    }
 
    size_t getComputationalCost(size_t id)
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        check_overflow(id);
        return gp.vertex(id).template get<nm_v_computation>();
    }
 
 
    void setComputationCost(size_t id, size_t cost)
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
#ifdef SE_CLASS1
        check_overflow(id);
#endif
 
        auto fnd = owner_scs.find(id);
        if (fnd == owner_scs.end())
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Error you are setting a sub-sub-domain the processor does not own" << std::endl;
        }
        else
        {
            size_t id = fnd->second;
            owner_cost_sub.get(id).w = cost;
        }
    }
 
    void setMigrationCost(size_t id, size_t cost)
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
#ifdef SE_CLASS1
        check_overflow(id);
#endif
 
        gp.vertex(id).template get<nm_v_migration>() = cost;
    }
 
    void setCommunicationCost(size_t id, size_t e, size_t cost)
    {
#ifdef SE_CLASS2
        check_valid(this,8);
#endif
#ifdef SE_CLASS1
        check_overflow(id);
        check_overflowe(id,e);
#endif
 
        gp.getChildEdge(id, e).template get<nm_e::communication>() = cost;
    }
 
    size_t getNSubSubDomains()
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        return gp.getNVertex();
    }
 
    size_t getNSubSubDomainNeighbors(size_t id)
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
#ifdef SE_CLASS1
        check_overflow(id);
#endif
 
        return gp.getNChilds(id);
    }
 
    float getUnbalance()
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        size_t load_p = getProcessorLoad();
 
        float load_avg = load_p;
        v_cl.sum(load_avg);
        v_cl.execute();
 
        if (load_avg == 0)
        {
            // count the number if sub-sub-domain assigned
            load_avg = owner_cost_sub.size();
 
            v_cl.sum(load_avg);
            v_cl.execute();
        }
 
        load_avg /= v_cl.getProcessingUnits();
 
        return ((float)load_p - load_avg) / load_avg;
    }
 
    size_t getNOwnerSubSubDomains() const
    {
        return owner_cost_sub.size();
    }
 
    size_t getOwnerSubSubDomain(size_t id) const
    {
        return owner_cost_sub.get(id).id;
    }
 
    void onTest()
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        testing = true;
    }
 
    void write(std::string out)
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
 
        VTKWriter<Graph_CSR<nm_v<dim>, nm_e>, VTK_GRAPH> gv2(gp);
        gv2.write(std::to_string(v_cl.getProcessUnitID()) + "_" + out + ".vtk");
 
    }
 
    size_t getProcessorLoad()
    {
#ifdef SE_CLASS2
            check_valid(this,8);
#endif
        openfpm::vector<size_t> loads(v_cl.getProcessingUnits());
 
        size_t load = 0;
 
        if (v_cl.getProcessUnitID() == 0)
        {
            for (size_t i = 0; i < gp.getNVertex(); i++)
            {loads.get(gp.template vertex_p<nm_v_proc_id>(i)) += gp.template vertex_p<nm_v_computation>(i);}
 
            for (size_t i = 0 ; i < v_cl.getProcessingUnits() ; i++)
            {
                v_cl.send(i,1234,&loads.get(i),sizeof(size_t));
            }
        }
        v_cl.recv(0,1234,&load,sizeof(size_t));
        v_cl.execute();
 
        return load;
    }
 
    MetisDistribution & operator=(const MetisDistribution & mt)
    {
#ifdef SE_CLASS2
            check_valid(&mt,8);
            check_valid(this,8);
#endif
        this->gr = mt.gr;
        this->domain = mt.domain;
        this->gp = mt.gp;
        this->owner_cost_sub = mt.owner_cost_sub;
        this->owner_scs = mt.owner_scs;
        return *this;
    }
 
    MetisDistribution & operator=(MetisDistribution && mt)
    {
#ifdef SE_CLASS2
            check_valid(mt);
            check_valid(this,8);
#endif
        this->gr = mt.gr;
        this->domain = mt.domain;
        this->gp.swap(mt.gp);
        this->owner_cost_sub.swap(mt.owner_cost_sub);
        this->owner_scs.swap(mt.owner_scs);
        return *this;
    }
 
    inline bool operator==(const MetisDistribution & mt)
    {
#ifdef SE_CLASS2
            check_valid(&mt,8);
            check_valid(this,8);
#endif
        bool ret = true;
 
        ret &= (this->gr == mt.gr);
        ret &= (this->domain == mt.domain);
        ret &= (this->gp == mt.gp);
 
        return ret;
    }
 
    void setDistTol(double tol)
    {
        metis_graph.setDistTol(tol);
    }
 
 
 
    size_t getSubSubDomainComputationCost(size_t id)
    {
#ifdef SE_CLASS1
        if (id >= gp.getNVertex())
            std::cerr << __FILE__ << ":" << __LINE__ << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
#endif
 
        auto fnd = owner_scs.find(id);
        if (fnd == owner_scs.end())
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Error you are setting a sub-sub-domain that the processor does not own" << std::endl;
            return 0;
        }
 
        size_t ids = fnd->second;
        return owner_cost_sub.get(ids).w;
    }
 
    size_t get_ndec()
    {
        return metis_graph.get_ndec();
    }
};
 
#endif /* SRC_DECOMPOSITION_METISDISTRIBUTION_HPP_ */