#include <iostream>
#include <fstream>
#include <queue>
#include <string>
#include <cstring>
#include <omp.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
using namespace std;

#define T_num 2

struct AdjListNode
{
    int dest;
    AdjListNode* next;
};

struct AdjList
{    
   AdjListNode *head;
};

class Graph
{
    private:
        int V;
        AdjList* array;
    
    public:
    	Graph()
    	{
    	    //always call setV if a Graph object is created this way
    	}
        Graph(int V)
        {
            this->V = V;
            array = new AdjList [V];
            for (int i = 0; i < V; ++i)
               array[i].head = NULL;
        }
	void setV(int V)
	{
	    this->V = V;
            array = new AdjList [V];
            for (int i = 0; i < V; ++i)
               array[i].head = NULL;
	}
        AdjListNode* newAdjListNode(int dest)
        {
            AdjListNode* newNode = new AdjListNode;
            newNode->dest = dest;
            newNode->next = NULL;
            return newNode;
        }
        void addEdge(int src, int dest)
        {
            AdjListNode* newNode = newAdjListNode(dest);
            newNode->next = array[src].head;
            array[src].head = newNode;
        }
        int findDegree(int row)
        {
            int count = -1;
            AdjListNode* pCrawl = array[row].head;
            while (pCrawl)
            {
                pCrawl = pCrawl->next;
                count++;
            }
            return count;
        }
        queue<int> getNeighbours(int v,int m,int n)
        {
           //cout<<"entered getNeighbours\n";
           queue<int> neighbours;
           bool row = true;
           if(v>=m)
              row = false;
           if(row)
           {
              //cout<<"accesing array for: "<<v<<endl;
              AdjListNode* pCrawl = array[v].head;
              //cout<<"Values:";
              while (pCrawl)
              {
                neighbours.push(pCrawl->dest);
                //cout<<pCrawl->dest<<" ";
                pCrawl = pCrawl->next;
              }
              //cout<<endl;
           }
           else
           {
             //cout<<"for node: "<<v<<endl;
             for( int i=0;i<m;i++)
             {
                //cout<<"accesing array for: "<<i<<endl;
                AdjListNode* pCrawl = array[i].head;
              	//cout<<"Found:";
              	while (pCrawl)
                {
                   if(pCrawl->dest == v)
                   {
                      neighbours.push(i);
                      //cout<<"yes\n";
                      break;
                   }
                   pCrawl = pCrawl->next;
                }
             }
          }
          return neighbours;
      }	
      void printGraph()
      {
         int v;
         for (v = 0; v < V; ++v)
         {
             AdjListNode* pCrawl = array[v].head;
             while (pCrawl)
             {
                 //cout<<"-> "<<pCrawl->dest;
                 pCrawl = pCrawl->next;
             }
             //cout<<endl;
         }
      }
};

class graphProcessor
{
    int m,n;
    int **g;
    int *degree;
    int *matching;
    int *visited;
    int *label;

    void MatchAndUpdate(int v);
    void globalRelabeling();
    
public:
    Graph graph;
    graphProcessor(){}
    void doSetup();
    graphProcessor(string filename);
    void doInitialMatching();
    int* pushRelabel(int breakpoint,int RFVrow);
    void pushrelableSeq(int RFVrow);
    int getN();
    int getM();
    void setM(int M){m = M;}
    void setN(int N){n = N;}
    ~graphProcessor();
    
    int getSuitableVertex(int v);
};

inline int fetch_and_add(int &var, int value)
{
	asm volatile("lock; xaddl %%eax, %2;"
			:"=a" (value)             //Output
			:"a" (value), "m" (var)  //Input
			:"memory");
	return value;
}

inline bool compare_and_swap(bool &var, bool value)
{
	asm volatile("lock; xaddl %%eax, %2;"
			:"=a" (value)             //Output
			:"a" (value), "m" (var)  //Input
			:"memory");
	return value;
}

void initialize_in_parallel(int *arr, int init_value, int size)
{
        #pragma omp parallel for shared(arr)
        for (int i = 0; i < size; i++) 
	{
		arr[i] = init_value; 	
	}
}

void get_large_input(graphProcessor& g) 
{	
    std::ifstream fin("input.mtx");

    int M, N, L;

    while (fin.peek() == '%') fin.ignore(2048, '\n');

    fin >> M >> N >> L;
    g.setM(M);
    
    g.setN(N);
    
    g.graph.setV(M);
    
    for (int l = 0; l < L; l++)
    {
	int m, n;
	float blah;
	fin >> m >> n >>blah;
	g.graph.addEdge(m-1,M+n-1);
    }
    fin.close();
}

void initialize_in_parallel(bool *arr, bool init_value, int size)
{
        #pragma omp parallel for default(shared) 
	for (int i = 0; i < size; i++) 
	     arr[i] = init_value; 	
}

void graphProcessor::MatchAndUpdate(int v)
{
    if(fetch_and_add(visited[v],1) == 0)
    {
        queue<int> neighbours;
        neighbours = graph.getNeighbours(v,m,n);
        
        while(!neighbours.empty())
        {
            int u = neighbours.front();
            neighbours.pop();
            if(fetch_and_add(visited[u],1) == 0)
            {
            	#pragma omp critical
                matching[v] = u;
                queue<int>new_neighbours = graph.getNeighbours(u,m,n);
                while(!new_neighbours.empty())
                {
                    int w = new_neighbours.front();
                    new_neighbours.pop();
                    
                    if(fetch_and_add(degree[w],-1) == 2)
                    {
                        MatchAndUpdate(w);
                    }
                }
                break;
            }
       }
    }
}

void graphProcessor::globalRelabeling()
{

    clock_t t1,t2;
    t1=clock();

    queue<int> Q;
    int dist = 1;
    
    for(int i=m ; i<m+n; i++)
    {
        label[i] = m+n;
    }
    
    for(int i=0;i<m;i++)
    {
        if(matching[i]!=-1)
            label[i]=m+n;
        else
        {
            Q.push(i);
            label[i]=0;
        }
    }
    
   while(!Q.empty())
    {
        int q_size = Q.size();
        queue<int> Q_dash;
        
        #pragma omp parallel for shared(Q,Q_dash)
        for(int i = 0; i < q_size; i++)
        {
            int u;
    	    #pragma omp critical
            {
        	u = Q.front();
	       	Q.pop();
            }
            
            queue<int> neighbours = graph.getNeighbours(u,m,n);
            while(!neighbours.empty())
            {
            	int v;
    	   	#pragma omp critical
            	{
        		v = neighbours.front();
	       		neighbours.pop();
            	}
            	
                if(label[v]==m+n)
                {
                    #pragma omp critical
                    label[v]=dist;
                    for(int i =0 ;i<m;i++)
                    {
                        if(matching[i]==v)
                        {
                            #pragma omp critical
                            {
                            	label[i]= dist+1;
                            	Q_dash.push(i);
                           	i=m;
			    }
                        }
                    }
                }
            }
        }
        swap(Q,Q_dash);
        dist+= 2;     
    }
    
    t2 = clock();
    float diff ((float)t2-(float)t1);
    cout<<"Execution Time for global relabeling: "<<(diff/CLOCKS_PER_SEC)<<endl;
}

graphProcessor::graphProcessor(string filename)
{
    fstream ip(filename.c_str(), fstream::in);
    //char* s = new char [256];
    
    ip >> m ;
    ip >> n ;

    g = new int* [m];
    for(int i =0;i<m;i++)
    {
        g[i] = new int[n];
    }
    degree = new int[m];
    matching = new int[m];
    visited = new int[m+n];
    label = new int[m+n];
   
    initialize_in_parallel(matching, -1, m);
    initialize_in_parallel(visited, 0, m+n);
    initialize_in_parallel(degree, 0, m);
    
    for (int i = 0; i < m; i++)
    {
        for(int j = 0; j < n; j++)
        {
            ip >> g[i][j];
            if(g[i][j] > 0)
            {
                degree[i] ++;
            }
        }
    }
}

void graphProcessor::doSetup()
{
    degree = new int[m];
    visited = new int[m+n];
    label = new int[m+n];
    matching = new int[m];
    
    initialize_in_parallel(matching, -1, m);
    initialize_in_parallel(visited, 0, m+n);
    initialize_in_parallel(degree, 0, m);
    
    for (int i = 0; i < m; i++)
    	degree[i] = graph.findDegree(i);
    
}
void graphProcessor::doInitialMatching()
{
    queue <int> q, q_star;
    
    #pragma omp parallel for shared(q,q_star)
    for(int i = 0; i < m; i++)
    {
        if(degree[i] == 1)
        {
            #pragma omp critical	
            q.push(i);
        }
        else
        {
            #pragma omp critical
            q_star.push(i);
        }
    }
    
    int size_q = q.size();
    
    #pragma omp parallel for shared(q)
    for(int i = 0; i < size_q ; i++)
    {
	int v;
    	#pragma omp critical
        {
        	v = q.front();
        	q.pop();
        }
        
        MatchAndUpdate(v);
    }

    int size_q_star = q_star.size();
    #pragma omp parallel for shared(q_star)
    for(int i = 0; i < size_q_star ; i++)
    {
        int v;
    	#pragma omp critical
        {
        	v = q_star.front();
	        q_star.pop();
        }
        MatchAndUpdate(v);
    }

}

int graphProcessor::getSuitableVertex(int v) { 
    
    //cout<<"entereed getSuitableVertex\n";
    int minlabel=m+n; 
    queue<int> neighbours = graph.getNeighbours(v,m,n); 
    //cout<<"exited getNeighbours\n";
    int minvertex;
    if(neighbours.size()==0)
    {
       minvertex = -1; //if this vertex has no neighbour at all
       return minvertex;
    }
    
    minvertex =neighbours.front();
    
    int size = neighbours.size();
    #pragma omp parallel for
    for(int i=0;i<size;i++) 
    {     
        int u;
    	#pragma omp critical
    	{
           u = neighbours.front(); 
           neighbours.pop(); 
        }
        //cout<<"after gettin neighbour u\n";
        if(label[u]<=minlabel) 
        { 
            minlabel = label[u]; 
            minvertex = u;
        } 
    } 
    //cout<<"exited getSuitableVertex with: "<<minvertex<<"\n";
    return minvertex; 
}

int* graphProcessor::pushRelabel(int breakpoint,int RFVrow)
{   
    clock_t t1,t2;
    t1=clock();
     
    int b=breakpoint;
    doInitialMatching();
    //cout<<"doInitialMatching Finished\n";
    int colMatching[n];
    queue<int> activeQ;
    initialize_in_parallel(colMatching,0,n);
    
    globalRelabeling(); /*do not remove this - it initializes label array */
    //cout<<"gloabalRelabeling finished\n";
    
    #pragma omp parallel for
    for(int i=0;i<m;i++)
    {
	//cout << "#threads = " << omp_get_num_threads() << endl;
        if(matching[i]!=-1)
        {
            colMatching[matching[i]-m] = 1; /*indicates a matching*/
        }
    }
    
    //cout<<"after filling matching matrix in pushrelabel\n";
    
#ifdef _OPENMP
    int *flags = new int[m+n];
    initialize_in_parallel(flags,0,m+n);
    omp_lock_t locks_array[m+n]; /* should have the same size as the label array */   
    #pragma parallel for
    for(int i=0;i<m+n;i++)
       omp_init_lock(&locks_array[i]);
#endif

    /*initializing the active Queue with unmatched vertices*/
    //cout<<"initializing the active Queue with unmatched vertices\n";
    #pragma omp parallel for shared(activeQ)
    for(int j=0;j<n;j++)
    {
        if(colMatching[j]==0)
        {
            #pragma omp critical
            activeQ.push(j+m);
        }
    }
    
    int verticesProcessed=0;
    
    //cout<<"before while in pushRelabel\n";
    while(activeQ.size()>b)
    { 
        int n_size = activeQ.size();
        #pragma omp parallel
        {     
            queue<int> Q_thread;
            //cout<<"before for in push Relabel\n";
            #pragma omp for 
            for(int i=0;i<n_size;i++)
            {
                int v;
                #pragma omp critical
                {
            	    v = activeQ.front();
            	    activeQ.pop();
                }
                //cout<<"getting v in push Relabel\n";
                #pragma omp critical
                verticesProcessed++;
            
                int u = getSuitableVertex(v);
                if(u==-1)
                {
                    activeQ.pop(); //handle vertices that do not have a neighbour at all
                    continue;
                }
               // cout<"getSuitableVertex Finished\n";
                int process = 1;
#ifdef _OPENMP
	        if(!omp_test_lock(&locks_array[u]))
                {
                   #pragma omp critical
                   activeQ.push(v);
                   process = 0;
                }            	                   
#endif
	        if(process)
                {
	            if(label[u]<m+n)
                    {
                        label[v]=label[u]+1;
                        //cout<<"after getting label in PushRelabel\n";
                        if(matching[u]!=-1)
                        {
                           int w = matching[u];
                           //cout<<"after getting w in PushRelabel\n";
                           Q_thread.push(w);
                        }
                        matching[u]=v; 
                        label[u]=label[u]+2;
                        //cout<<"after updating label in PushRelabel\n";
                    }

#ifdef _OPENMP
                    omp_unset_lock(&locks_array[u]);
#endif

                } 
            
        }
        //cout<<"after for in pushRelabel\n";
        
        #pragma omp barrier
        
#ifdef _OPENMP
        int size = Q_thread.size();
        
	//cout<<"before swap in pushRelabel\n";
        #pragma omp parallel for /*nested parallel construct - check for correctness */
        for(int iter =0;iter<size;iter++)
        {
        	int val=0;
                #pragma omp critical
                {
                    if(!Q_thread.empty())
                    {
        	       val = Q_thread.front();
        	       Q_thread.pop();
        	    }
        	}
                if(omp_test_lock(&locks_array[val]))
                {
                      if(flags[val]!=1)
                      { 
                          flags[val] = 1;
                          #pragma omp critical
                          activeQ.push(val);
                      }
                      omp_unset_lock(&locks_array[val]);
                }
        }
    }
    initialize_in_parallel(flags,0,m+n);
#else
    swap(activeQ,Q_thread);
    }
#endif
//cout<<"after swap in pushRelabel\n";
    if(verticesProcessed > RFVrow)
    {
        globalRelabeling();
        verticesProcessed = 0; //restart the cycle
    }
 }
 
 t2 = clock();
 float diff ((float)t2-(float)t1);
 cout<<"Execution Time for push relable: "<<(diff/CLOCKS_PER_SEC)<<endl;
 t1 = clock();
 pushrelableSeq(RFVrow);
 t2 = clock();
 diff = ((float)t2-(float)t1);
 cout<<"Execution Time for push relable seq: "<<(diff/CLOCKS_PER_SEC)<<endl;
   
 #ifdef _OPENMP
 delete flags;
 #endif
 return matching;
}

void graphProcessor::pushrelableSeq(int RFVrow)
{
    //cout<<"entered pushRelabel seq\n";
    int colMatching[n];
    queue<int> activeQ;
    initialize_in_parallel(colMatching,0,n);
    
    for(int i=0;i<m;i++)
    {
        if(matching[i]!=-1)
        {
            colMatching[matching[i]-m] = 1; //indicates a matching
        }
    }

    //initializing the active Queue with unmatched vertices
    for(int j=0;j<n;j++)
    {
        if(colMatching[j]==0)
        {
            activeQ.push(j+m);
        }

    }
    int verticesProcessed=0;
    //cout<<"before while in push relabel seq\n";
    while(activeQ.size()>0)
    {
        queue<int> Q_thread;
        for(;!activeQ.empty();)
        {
            int v = activeQ.front();
            verticesProcessed++;
            //cout<<"before getSuitableVertex in pushRelabel seqq\n";
            int u = getSuitableVertex(v);
            if(u==-1)
            {
               activeQ.pop(); //handle vertices that do not have a neighbour 
               continue;
            }
            //cout<<"after getSuitableVertex in pushRelabel seqq\n";
            if(label[u]<m+n)
            {
                //cout<<"finding label for v in PRSeqq\n";
                label[v]=label[u]+1;
                if(matching[u]!=-1)
                {
                    //cout<<"finding w in PRSeq\n";
                    int w = matching[u];
                    Q_thread.push(w);
                }
                    //cout<<"finding matching u\n";
                    matching[u]=v; 
                    //cout<<"upating label for u\n";
                    label[u]=label[u]+2;
            }
            activeQ.pop();
        }
        swap(activeQ,Q_thread);
        if(verticesProcessed > RFVrow)
        {
            //cout<<"before gloabal relabeling\n";
            globalRelabeling();
            //cout<<"after global relabeling\n";
            verticesProcessed = 0; //restart the cycle
        }
    }
    //cout<<"exit push relabel seq\n"; 
}
int graphProcessor::getN()
{
	return n;
}
int graphProcessor::getM()
{
	return m;
}

graphProcessor::~graphProcessor()
{
    if(visited)
        delete[] visited;
    if(matching)
        delete[] matching;
    if(degree)
        delete[] degree;
}

int main()
{
    clock_t t1,t2;
    graphProcessor g;
    get_large_input(g);
    g.doSetup();
	int nthreads = atoi(getenv("OMP_NUM_THREADS"));
    t1=clock();
    int * matching = g.pushRelabel(45,2*(g.getN()));
    t2 = clock();
    float diff ((float)t2-(float)t1);
    cout<<"Execution Time with " << nthreads << " threads: "<<(diff/CLOCKS_PER_SEC/nthreads)<<endl;
}