parmetis_util.hpp

/*
 * parmetis_util.hpp
 *
 *  Created on: Oct 07, 2015
 *      Author: Antonio Leo
 */
 
#ifndef PARMETIS_UTIL_HPP
#define PARMETIS_UTIL_HPP
 
#include <iostream>
#include "parmetis.h"
#include "VTKWriter/VTKWriter.hpp"
#include "VCluster/VCluster.hpp"
#include "Graph/ids.hpp"
 
#define PARMETIS_ERROR_OBJECT std::runtime_error("Runtime Parmetis error");
 
 
struct Parmetis_graph
{
    idx_t * nvtxs;
 
    idx_t * ncon;
 
    idx_t * xadj;
 
    idx_t * adjncy;
 
    idx_t * vwgt;
 
    idx_t * vsize;
 
    idx_t * adjwgt;
 
    idx_t * nparts;
 
    real_t * tpwgts;
 
    real_t * ubvec;
 
    idx_t * options;
 
    idx_t * objval;
 
    idx_t * part;
 
    idx_t * edgecut;
 
    real_t * itr;
 
    idx_t * numflag;
 
    // 0 No weights (vwgt and adjwgt are both NULL).
    // 1 Weights on the edges only (vwgt is NULL).
    // 2 Weights on the vertices only (adjwgt is NULL).
    // 3 Weights on both the vertices and edges.
    idx_t * wgtflag;
};
 
#define BALANCE_CC 1
#define BALANCE_CCM 2
#define BALANCE_CC_O(c) c+1
 
template<typename Graph>
class Parmetis
{
    Parmetis_graph Mg;
 
    // Original graph
    //  Graph & g;
 
    MPI_Comm comm = (MPI_Comm)NULL;
 
    Vcluster<> & v_cl;
 
    int p_id = 0;
 
    size_t nc = 0;
 
    rid first;
 
    rid last;
 
    size_t nvertex;
 
    size_t n_dec;
 
    real_t dist_tol = 1.05;
 
    void constructAdjList(Graph &g, const std::unordered_map<rid, gid> & m2g)
    {
        // init basic graph informations and part vector
        // Put the total communication size to NULL
 
        Mg.nvtxs[0] = nvertex;
        if (Mg.part != NULL)
        {delete[] Mg.part;}
        Mg.part = new idx_t[nvertex];
 
        size_t nedge = 0;
        size_t i = 0;
        for (rid j = first; i < nvertex; i++, ++j)
        {
            Mg.part[i] = p_id;
            nedge += g.getNChilds(m2g.find(j)->second.id);
        }
 
        // create xadj, adjlist, vwgt, adjwgt and vsize
        if (Mg.xadj != NULL)
        {delete[] Mg.xadj;}
        if (Mg.adjncy != NULL)
        {delete[] Mg.adjncy;}
        if (Mg.vwgt != NULL)
        {delete[] Mg.vwgt;}
        if (Mg.adjwgt != NULL)
        {delete[] Mg.adjwgt;}
        if (Mg.vsize != NULL)
        {delete[] Mg.vsize;}
        Mg.xadj = new idx_t[nvertex + 1];
        Mg.adjncy = new idx_t[nedge];
        Mg.vwgt = new idx_t[nvertex];
        Mg.adjwgt = new idx_t[nedge];
        Mg.vsize = new idx_t[nvertex];
 
        size_t prev = 0;
 
        // actual position
        size_t id = 0;
 
        size_t j = 0;
 
        // for each vertex calculate the position of the starting point in the adjacency list
        for (rid i = first; i <= last; ++i, j++)
        {
            gid idx = m2g.find(i)->second;
 
            // Add weight to vertex and migration cost
            Mg.vwgt[j] = g.vertex(idx.id).template get<nm_v_computation>();
            Mg.vsize[j] = g.vertex(idx.id).template get<nm_v_migration>();
 
            // Calculate the starting point in the adjacency list
            Mg.xadj[id] = prev;
 
            // Create the adjacency list and the weights for edges
            for (size_t s = 0; s < g.getNChilds(idx.id); s++)
            {
 
                size_t child = g.getChild(idx.id, s);
 
                Mg.adjncy[prev + s] = g.vertex(child).template get<nm_v_id>();
                Mg.adjwgt[prev + s] = g.getChildEdge(idx.id, s).template get<nm_e::communication>();
            }
 
            // update the position for the next vertex
            prev += g.getNChilds(idx.id);
 
            id++;
        }
 
        // Fill the last point
        Mg.xadj[id] = prev;
    }
 
public:
 
    Parmetis(Vcluster<> & v_cl, size_t nc)
    :v_cl(v_cl), nc(nc),n_dec(0)
    {
#ifdef SE_CLASS1
 
        if (sizeof(idx_t) != 8)
        {
            std::cerr << __FILE__ << ":" << __LINE__ << " Error detected invalid installation of Parmetis. OpenFPM support Parmetis/Metis version with 64 bit idx_t" << std::endl;
            ACTION_ON_ERROR(PARMETIS_ERROR_OBJECT);
        }
 
#endif
 
        // TODO Move into VCluster
        MPI_Comm_dup(MPI_COMM_WORLD, &comm);
 
        // Nullify Mg
        Mg.nvtxs = NULL;
        Mg.ncon = NULL;
        Mg.xadj = NULL;
        Mg.adjncy = NULL;
        Mg.vwgt = NULL;
        Mg.vsize = NULL;
        Mg.adjwgt = NULL;
        Mg.nparts = NULL;
        Mg.tpwgts = NULL;
        Mg.ubvec = NULL;
        Mg.options = NULL;
        Mg.objval = NULL;
        Mg.part = NULL;
        Mg.edgecut = NULL;
        Mg.itr = NULL;
        Mg.numflag = NULL;
        Mg.wgtflag = NULL;
        first.id = 0;
        last.id = 0;
        nvertex = 0;
    }
 
    ~Parmetis()
    {
        // Deallocate the Mg structure
        if (Mg.nvtxs != NULL)
        {
            delete[] Mg.nvtxs;
        }
 
        if (Mg.ncon != NULL)
        {
            delete[] Mg.ncon;
        }
 
        if (Mg.xadj != NULL)
        {
            delete[] Mg.xadj;
        }
 
        if (Mg.adjncy != NULL)
        {
            delete[] Mg.adjncy;
        }
 
        if (Mg.vwgt != NULL)
        {
            delete[] Mg.vwgt;
        }
 
        if (Mg.adjwgt != NULL)
        {
            delete[] Mg.adjwgt;
        }
 
        if (Mg.nparts != NULL)
        {
            delete[] Mg.nparts;
        }
 
        if (Mg.tpwgts != NULL)
        {
            delete[] Mg.tpwgts;
        }
 
        if (Mg.ubvec != NULL)
        {
            delete[] Mg.ubvec;
        }
 
        if (Mg.options != NULL)
        {
            delete[] Mg.options;
        }
 
        if (Mg.part != NULL)
        {
            delete[] Mg.part;
        }
 
        if (Mg.edgecut != NULL)
        {
            delete[] Mg.edgecut;
        }
 
        if (Mg.numflag != NULL)
        {
            delete[] Mg.numflag;
        }
 
        if (Mg.wgtflag != NULL)
        {
            delete[] Mg.wgtflag;
        }
 
        if (Mg.itr != NULL)
        {
            delete[] Mg.itr;
        }
 
        if (Mg.vsize != NULL)
        {
            delete[] Mg.vsize;
        }
 
        if (Mg.objval != NULL)
        {
            delete[] Mg.objval;
        }
 
        if (is_openfpm_init() == true)
        {MPI_Comm_free(&comm);}
    }
 
    void initSubGraph(Graph & g,
                      const openfpm::vector<rid> & vtxdist,
                      const std::unordered_map<rid, gid> & m2g,
                      bool w)
    {
        p_id = v_cl.getProcessUnitID();
 
        first = vtxdist.get(p_id);
        last = vtxdist.get(p_id + 1) - 1;
        nvertex = last.id - first.id + 1;
 
        setDefaultParameters(w);
 
        // construct the adjacency list
        constructAdjList(g, m2g);
    }
 
    void decompose(const openfpm::vector<rid> & vtxdist)
    {
        // Decompose
 
        ParMETIS_V3_PartKway((idx_t *) vtxdist.getPointer(), Mg.xadj, Mg.adjncy, Mg.vwgt, Mg.adjwgt, Mg.wgtflag, Mg.numflag, Mg.ncon, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.options, Mg.edgecut, Mg.part, &comm);
 
        n_dec++;
    }
 
    void refine(openfpm::vector<rid> & vtxdist)
    {
        // Refine
 
        ParMETIS_V3_RefineKway((idx_t *) vtxdist.getPointer(), Mg.xadj, Mg.adjncy, Mg.vwgt, Mg.adjwgt, Mg.wgtflag, Mg.numflag, Mg.ncon, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.options, Mg.edgecut, Mg.part, &comm);
 
        n_dec++;
    }
 
    void redecompose(openfpm::vector<rid> & vtxdist)
    {
        ParMETIS_V3_AdaptiveRepart((idx_t *)vtxdist.getPointer(), Mg.xadj, Mg.adjncy, Mg.vwgt, Mg.vsize, Mg.adjwgt, Mg.wgtflag, Mg.numflag, Mg.ncon, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.itr, Mg.options, Mg.edgecut, Mg.part, &comm);
 
        n_dec++;
    }
 
    idx_t* getPartition()
    {
        return Mg.part;
    }
 
    void reset(Graph & g, const openfpm::vector<rid> & vtxdist, const std::unordered_map<rid, gid> & m2g, bool vgw)
    {
        first = vtxdist.get(p_id);
        last = vtxdist.get(p_id + 1) - 1;
        nvertex = last.id - first.id + 1;
 
        // Deallocate the graph structures
 
        setDefaultParameters(vgw);
 
        // construct the adjacency list
        constructAdjList(g, m2g);
    }
 
    void setDefaultParameters(bool w)
    {
        if (Mg.nvtxs != NULL)
        {delete[] Mg.nvtxs;}
        Mg.nvtxs = new idx_t[1];
 
        // Set the number of constrains
        if (Mg.ncon != NULL)
        {delete[] Mg.ncon;}
        Mg.ncon = new idx_t[1];
        Mg.ncon[0] = 1;
 
        // Set to null the weight of the vertex (init after in constructAdjList) (can be removed)
        if (Mg.vwgt != NULL)
        {delete[] Mg.vwgt;}
        Mg.vwgt = NULL;
 
        // Set to null the weight of the edge (init after in constructAdjList) (can be removed)
        if (Mg.adjwgt != NULL)
        {delete[] Mg.adjwgt;}
        Mg.adjwgt = NULL;
 
        // Set the total number of partitions
        if (Mg.nparts != NULL)
        {delete[] Mg.nparts;}
        Mg.nparts = new idx_t[1];
        Mg.nparts[0] = nc;
 
 
        if (Mg.options != NULL)
        {delete[] Mg.options;}
        Mg.options = new idx_t[4];
        Mg.options[0] = 0;
        Mg.options[1] = 0;
        Mg.options[2] = 0;
        Mg.options[3] = 0;
 
 
        if(Mg.itr != NULL)
        {delete[] Mg.itr;}
        Mg.itr = new real_t[1];
        Mg.itr[0] = 1000.0;
 
        if (Mg.objval != NULL)
        {delete[] Mg.objval;}
        Mg.objval = new idx_t[1];
 
 
        if (Mg.tpwgts != NULL)
        {delete[] Mg.tpwgts;}
        if (Mg.ubvec != NULL)
        {delete[] Mg.ubvec;}
        Mg.tpwgts = new real_t[Mg.nparts[0]];
        Mg.ubvec = new real_t[Mg.nparts[0]];
 
        for (size_t s = 0; s < (size_t) Mg.nparts[0]; s++)
        {
            Mg.tpwgts[s] = 1.0 / Mg.nparts[0];
            Mg.ubvec[s] = dist_tol;
        }
 
        if (Mg.edgecut != NULL)
        {delete[] Mg.edgecut;}
        Mg.edgecut = new idx_t[1];
        Mg.edgecut[0] = 0;
 
        if (Mg.numflag != NULL)
        {delete[] Mg.numflag;}
        Mg.numflag = new idx_t[1];
        Mg.numflag[0] = 0;
 
        if (Mg.wgtflag != NULL)
        {delete[] Mg.wgtflag;}
        Mg.wgtflag = new idx_t[1];
 
        if (w)
            Mg.wgtflag[0] = 3;
        else
            Mg.wgtflag[0] = 0;
    }
 
    const Parmetis<Graph> & operator=(const Parmetis<Graph> & pm)
    {
        if (comm != MPI_COMM_NULL){MPI_Comm_free(&comm);}
        MPI_Comm_dup(pm.comm, &comm);
        p_id = pm.p_id;
        nc = pm.nc;
        n_dec = pm.n_dec;
        dist_tol = pm.dist_tol;
 
        setDefaultParameters(pm.Mg.wgtflag[0] == 3);
 
        return *this;
    }
 
    const Parmetis<Graph> & operator=(Parmetis<Graph> && pm)
    {
        // TODO Move into VCluster
        if (comm != MPI_COMM_NULL){MPI_Comm_free(&comm);}
        MPI_Comm_dup(pm.comm, &comm);
        p_id = pm.p_id;
        nc = pm.nc;
        n_dec = pm.n_dec;
        dist_tol = pm.dist_tol;
 
        setDefaultParameters(pm.Mg.wgtflag[0] == 3);
 
        return *this;
    }
 
    size_t get_ndec()
    {
        return n_dec;
    }
 
    void setDistTol(real_t tol)
    {
        dist_tol = tol;
    }
};
 
#endif