metis_util.hpp

/*
 * metis_util.hpp
 *
 *  Created on: Nov 21, 2014
 *      Author: Pietro Incardona
 */
 
#ifndef METIS_UTIL_HPP
#define METIS_UTIL_HPP
 
#include <iostream>
#include "metis.h"
#include "SubdomainGraphNodes.hpp"
#include "VTKWriter/VTKWriter.hpp"
 
struct Metis_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;
};
 
#define BALANCE_CC 1
#define BALANCE_CCM 2
#define BALANCE_CC_O(c) c+1
 
template<typename Graph>
class Metis
{
    Metis_graph Mg;
 
    Graph & g;
 
    size_t n_dec;
 
    real_t dist_tol = 1.05;
 
    void constructAdjList(Graph & g)
    {
        // create xadj and adjlist
 
        // create xadj, adjlist, vwgt, adjwgt and vsize
        Mg.xadj = new idx_t[g.getNVertex() + 1];
        Mg.adjncy = new idx_t[g.getNEdge()];
 
        size_t prev = 0;
 
        // actual position
        size_t id = 0;
 
        // for each vertex calculate the position of the starting point in the adjacency list
        for (size_t i = 0; i < g.getNVertex(); i++)
        {
            // Calculate the starting point in the adjacency list
            Mg.xadj[id] = prev;
 
            // Create the adjacency list
            for (size_t s = 0; s < g.getNChilds(i); s++)
            {
                Mg.adjncy[prev + s] = g.getChild(i, s);
            }
 
            // update the position for the next vertex
            prev += g.getNChilds(i);
 
            id++;
        }
 
        // Fill the last point
        Mg.xadj[id] = prev;
    }
 
    void constructAdjListWithWeights(Graph & g)
    {
        // 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[g.getNVertex() + 1];
        Mg.adjncy = new idx_t[g.getNEdge()];
        Mg.vwgt = new idx_t[g.getNVertex()];
        Mg.adjwgt = new idx_t[g.getNEdge()];
        Mg.vsize = new idx_t[g.getNVertex()];
 
        size_t prev = 0;
 
        // actual position
        size_t id = 0;
 
        // for each vertex calculate the position of the starting point in the adjacency list
        for (size_t i = 0; i < g.getNVertex(); i++)
        {
            // Add weight to vertex and migration cost
            Mg.vwgt[i] = g.vertex(i).template get<nm_v_computation>();
            Mg.vwgt[i] = (Mg.vwgt[i] == 0)?1:Mg.vwgt[i];
            Mg.vsize[i] = g.vertex(i).template get<nm_v_migration>();
            Mg.vsize[i] = (Mg.vsize[i] == 0)?1:Mg.vsize[i];
 
            // Calculate the starting point in the adjacency list
            Mg.xadj[id] = prev;
 
            // Create the adjacency list
            for (size_t s = 0; s < g.getNChilds(i); s++)
            {
                Mg.adjncy[prev + s] = g.getChild(i, s);
 
                // zero values on Metis are dangerous
                Mg.adjwgt[prev + s] = g.getChildEdge(i, s).template get<nm_e::communication>();
                Mg.adjwgt[prev + s] = (Mg.adjwgt[prev + s] == 0)?1:Mg.adjwgt[prev + s];
            }
 
            // update the position for the next vertex
            prev += g.getNChilds(i);
 
            id++;
        }
 
        // Fill the last point
        Mg.xadj[id] = prev;
    }
 
 
    void reset_pointer()
    {
        Mg.nvtxs = NULL;
        Mg.ncon = NULL;
        Mg.xadj = NULL;
        Mg.adjncy = NULL;
        Mg.vwgt = NULL;
        Mg.adjwgt = NULL;
        Mg.nparts = NULL;
        Mg.tpwgts = NULL;
        Mg.ubvec = NULL;
        Mg.options = NULL;
        Mg.objval = NULL;
        Mg.part = NULL;
        Mg.vsize = NULL;
    }
 
public:
 
    Metis(Graph & g, size_t nc, bool useWeights)
    :g(g),n_dec(0)
    {
        reset_pointer();
        initMetisGraph(nc,useWeights);
    }
 
    Metis(Graph & g, size_t nc)
    :g(g),n_dec(0)
    {
        reset_pointer();
        initMetisGraph(nc,false);
    }
 
    Metis(Graph & g)
    :g(g),n_dec(0)
    {
        reset_pointer();
    }
 
 
    void initMetisGraph(int nc, bool useWeights)
    {
 
        // Get the number of vertex
 
        if (Mg.nvtxs != NULL)
        {delete[] Mg.nvtxs;}
        Mg.nvtxs = new idx_t[1];
        Mg.nvtxs[0] = g.getNVertex();
 
        // 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
 
        if (Mg.vwgt != NULL)
        {delete[] Mg.vwgt;}
        Mg.vwgt = NULL;
 
        // Put the total communication size to NULL
 
        if (Mg.vsize != NULL)
        {delete[] Mg.vsize;}
        Mg.vsize = NULL;
 
        // Set to null the weight of the edge
 
        if (Mg.adjwgt != NULL)
        {delete[] Mg.adjwgt;}
        Mg.adjwgt = NULL;
 
        // construct the adjacency list
        if (useWeights)
            constructAdjListWithWeights(g);
        else
            constructAdjList(g);
 
        // Set the total number of partitions
 
        if (Mg.nparts != NULL)
        {delete[] Mg.nparts;}
        Mg.nparts = new idx_t[1];
        Mg.nparts[0] = nc;
 
        // Set to null the desired weight for each partition (one for each constrain)
 
        if (Mg.tpwgts != NULL)
        {delete[] Mg.tpwgts;}
        Mg.tpwgts = NULL;
 
        if (Mg.ubvec != NULL)
        {delete[] Mg.ubvec;}
        Mg.ubvec = NULL;
 
 
        if (Mg.options != NULL)
        {delete[] Mg.options;}
        Mg.options = NULL;
 
        if (Mg.objval != NULL)
        {delete[] Mg.objval;}
        Mg.objval = new idx_t[1];
 
        if (Mg.part != NULL)
        {delete[] Mg.part;}
        Mg.part = new idx_t[g.getNVertex()];
 
        for (size_t i = 0; i < g.getNVertex(); i++)
            Mg.part[i] = 0;
    }
 
    ~Metis()
    {
        // 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.objval != NULL)
        {
            delete[] Mg.objval;
        }
        if (Mg.part != NULL)
        {
            delete[] Mg.part;
        }
        if (Mg.vsize != NULL)
        {
            delete[] Mg.vsize;
        }
    }
 
    template<unsigned int i>
    void decompose()
    {
        if (Mg.nparts[0] != 1)
        {
            // Decompose
            METIS_PartGraphRecursive(Mg.nvtxs, Mg.ncon, Mg.xadj, Mg.adjncy, Mg.vwgt, Mg.vsize, Mg.adjwgt, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.options, Mg.objval, Mg.part);
 
            // vertex id
 
            size_t id = 0;
 
            // For each vertex store the processor that contain the data
 
            auto it = g.getVertexIterator();
 
            while (it.isNext())
            {
                g.vertex(it).template get<i>() = Mg.part[id];
 
                ++id;
                ++it;
            }
        }
        else
        {
            // Trivially assign all the domains to the processor 0
 
            auto it = g.getVertexIterator();
 
            while (it.isNext())
            {
                g.vertex(it).template get<i>() = 0;
 
                ++it;
            }
        }
 
        n_dec++;
    }
 
    template<unsigned int i, typename Graph_part>
    void decompose(Graph_part & gp)
    {
        // Decompose
        METIS_PartGraphRecursive(Mg.nvtxs, Mg.ncon, Mg.xadj, Mg.adjncy, Mg.vwgt, Mg.vsize, Mg.adjwgt, Mg.nparts, Mg.tpwgts, Mg.ubvec, Mg.options, Mg.objval, Mg.part);
 
        // vertex id
 
        size_t id = 0;
 
        // For each vertex store the processor that contain the data
 
        auto it = gp.getVertexIterator();
 
        while (it.isNext())
        {
            gp.vertex(it).template get<i>() = Mg.part[id];
 
            ++id;
            ++it;
        }
 
        n_dec++;
    }
 
    void onTest(bool testing)
    {
        if (testing == false)
            return;
 
        if (Mg.options == NULL)
        {
            // allocate
            Mg.options = new idx_t[METIS_NOPTIONS];
 
            // set default options
            METIS_SetDefaultOptions(Mg.options);
        }
 
        Mg.options[METIS_OPTION_SEED] = 0;
    }
 
    void setDistTol(real_t tol)
    {
        dist_tol = tol;
    }
 
    size_t get_ndec()
    {
        return n_dec;
    }
 
    void inc_dec()
    {
        n_dec++;
    }
};
 
#endif