explosive percolator openacc

#include<bits/stdc++.h>
#include "openacc_curand.h"
#include<curand.h>
using namespace std;

int EMPTY = 0;
int connectionNumber = 0;
//int cao = 0;

#pragma acc routine seq
int findroot(int ptr[],int i)
{
    if(ptr[i]<0) return i;

    return ptr[i]=findroot(ptr,ptr[i]);
}

/*
int findroot(int ptr[],int i)
{
    //cao = 1;
    int r,s;
    r = s = i;
    while (ptr[r]>=0)
    {
        ptr[s] = ptr[r];
        s = r;
        r = ptr[r];
    }
    return r;
}
*/

int main()
{
    int seed,vertices,m,run,filenum;

    cout<<"enter seed size: ";
    cin>>seed;
    cout<<endl;
    cout<<"enter vertice number: ";
    cin>>vertices;
    cout<<endl;
    cout<<"order number: ";
    cin>>m;
    cout<<endl;
    cout<<"enter ensemble number: ";
    cin>>run;
    cout<<endl;
    cout<<"enter filenumber: ";
    cin>>filenum;
    cout<<endl;

    int con = 0;

    for(int i=1;i<seed;i++)
    {
        con = con + i;
    }

    con = con + (vertices-seed)*m;

    int* netmap1 = new int[con+1];
    int* netmap2 = new int[con+1];

    for(int i=1;i<=con;i++)
    {
        netmap1[i] = 0;
        netmap2[i] = 0;
    }

    connectionNumber = con;

    srand(time(NULL));
    int A,B,C;
    A = vertices;
    B = run;
    C = filenum;

    string filename1;
    filename1 = "maxclus_";
    string filename2;
    filename2 = "delmx_";
    string filename3;
    filename3 = "entropy_";

    filename1 = filename1+to_string(A)+"node_"+to_string(m)+"M_"+to_string(B)+"ens"+to_string(C)+".dat";
    filename2 = filename2+to_string(vertices)+"node_"+to_string(m)+"M_"+to_string(run)+"ens"+to_string(filenum)+".dat";
    filename3 = filename3+to_string(vertices)+"node_"+to_string(m)+"M_"+to_string(run)+"ens"+to_string(filenum)+".dat";

    ofstream fout1,fout2,fout3;

    int* random = NULL;
    random = new int[connectionNumber+1];

    double* maxclus = NULL;
    maxclus = new double[connectionNumber+1];

    double* delmx = NULL;
    delmx = new double[connectionNumber+1];

    double* entropycalc = NULL;
    entropycalc = new double[connectionNumber+1];

    for(int i=0;i<=connectionNumber;i++)
    {
        maxclus[i]=0;
        delmx[i]=0;
        entropycalc[i]=0;
    }

    for(int i=0;i<=connectionNumber;i++)
    {
        random[i] = i;
    }
    int* ptr = new int[vertices+1];
    for(int i=0;i<vertices+1;i++) ptr[i]=0;


    for(int i=1;i<=25;i++)
    {
        cout<<"run : "<<i<<endl;
        vector<int> network;
        connectionNumber = 0;

        for(int i=1;i<=seed;i++)
        {
            for(int j=1;j<=seed;j++)
            {
                if(j>i)
                {

                    connectionNumber=connectionNumber + 1;
                    netmap1[connectionNumber]=i;
                    netmap2[connectionNumber]=j;
                    network.push_back(i);
                    network.push_back(j);
                }
            }
        }

        int concheck = 0;
        int ab[m];
        int countm = 0;

        for(int i = seed+1;i<=vertices; i++)
        {
            countm = 0;
            for(int k=0;k<m;k++) ab[k] = 0;
            //if(i%10000==0)random_shuffle( network.begin(),network.end() );
            //if(i%50000==0)cout<<"connection given to : "<<i<<endl;
            for(int j = 1; ;j++)
            {
                concheck = 0;
                int N1=network.size() ;
                int M1=0;
                int u = M1 + rand()/(RAND_MAX/(N1-M1+1) + 1);
                /*if(j==1)cout<<"netsize:"<<network.size()<<endl;
                if(j==1)cout<<"u:"<<u<<endl;
                if(j==1)cout<<"node:"<<network[u]<<endl;*/
                for(int n=0;n<m;n++)
                {
                    if(ab[n] == network[u]) concheck = 1;
                }

                if(concheck == 0 && network[u]!=i)
                {
                    ab[countm] = network[u];
                    countm=countm+1;



                    connectionNumber=connectionNumber+1;
                    netmap1[connectionNumber] = i;
                    netmap2[connectionNumber] = network[u];

                    network.push_back(i);
                    network.push_back(network[u]);
                }

                if(countm==m) break;

            }
        }
        random_shuffle(&random[1],&random[connectionNumber]);

        double rand_seed = time(NULL);

        #pragma acc data copy(maxclus[0:connectionNumber],delmx[0:connectionNumber],entropycalc[0:connectionNumber]), copyin(netmap1[0:connectionNumber],netmap2[0:connectionNumber],random[0:connectionNumber])
        #pragma acc enter data copyin(ptr[0:vertices+1])
        #pragma acc parallel loop independent private(random[0:connectionNumber], ptr[0:vertices+1])
        for(int rx=1;rx<=1000;rx++)
        {

            int index=0,big=0,bigtemp=0,jump=0,en1=0,en2=0;
            int node1=0,node2=0,node3=0,node4=0,nodeA=0,nodeB=0;
            int clus1=0,clus2=0,clus3=0,clus4=0;
            double entropy = log(vertices);

            curandState_t state;
            curand_init(rand_seed*rx,0,0,&state);


            for(int i=1;i<=connectionNumber;i++)
            {
                //cout<<"con number : "<<i<<endl;
                if(i!=connectionNumber)
                {

                    node1 = netmap1[random[i]];
                    node2 = netmap2[random[i]];
                    node3 = netmap1[random[i+1]];
                    node4 = netmap2[random[i+1]];

                    if(ptr[node1]==EMPTY) clus1 = 1;
                    else
                    {
                        int x = findroot(ptr,node1);
                        clus1 = -ptr[x];
                    }

                    if(ptr[node2]==EMPTY) clus2 = 1;
                    else
                    {
                        int b = findroot(ptr,node2);
                        clus2 = -ptr[b];
                    }

                    if(ptr[node3]==EMPTY) clus3 = 1;
                    else
                    {
                        int c = findroot(ptr,node3);
                        clus3 = -ptr[c];
                    }

                    if(ptr[node4]==EMPTY) clus4 = 1;
                    else
                    {
                        int d = findroot(ptr,node4);
                        clus4 = -ptr[d];
                    }

                    if(clus1*clus2<clus3*clus4)
                    {
                        nodeA = node1;
                        nodeB = node2;


                        int N1=connectionNumber;
                        int M1=i+1;
                         //int u = M1 + rand()/(RAND_MAX/(N1-M1+1) + 1);
                        int u = M1 + (curand_uniform(&state)*(N1-M1));

                        int temp = random[u];
                        random[u] = random[i+1];
                        random[i+1] = temp;
                    }
                    else
                    {
                        nodeA = node3;
                        nodeB = node4;

                        int temp = random[i+1];
                        random[i+1] = random[i];
                        random[i] = temp;
                        int N1=connectionNumber;
                        int M1=i+1;
                        //int u = M1 + rand()/(RAND_MAX/(N1-M1+1) + 1);
                        int u = M1 + (curand_uniform(&state)*(N1-M1));

                        temp = random[u];
                        random[u] = random[i+1];
                        random[i+1] = temp;
                    }
                }
                else
                {
                    nodeA = netmap1[random[i]];
                    nodeB = netmap2[random[i]];
                }



                if(ptr[nodeA]==EMPTY && ptr[nodeB]==EMPTY)
                {

                en1=1;
                en2=1;
                ptr[nodeA] = -2;
                ptr[nodeB] = nodeA;
                index = nodeA;
                entropy = (double)(entropy-(-2.0/vertices*log(1.0/vertices))+(-2.0/vertices*log(2.0/vertices)));
                if(entropy<0) entropy = 0;
                }
                else if(ptr[nodeA]==EMPTY && ptr[nodeB]!=EMPTY)
                {

                en1=1;
                int e = findroot(ptr,nodeB);
                en2 = -(ptr[e]);
                ptr[nodeA] = e;
                ptr[e] += -1;
                index = e;
                entropy = entropy-(-(double)1.0/vertices*log(1.0/(double)vertices))-(-(double)en2/vertices*log((double)en2/vertices))+(-( double)(-ptr[index])/vertices*log((-ptr[index])/(double)vertices));
                if(entropy<0) entropy = 0;


                }
                else if(ptr[nodeA]!=EMPTY && ptr[nodeB]==EMPTY)
                {

                en2 = 1;
                int f = findroot(ptr,nodeA);
                en1 = -(ptr[f]);
                ptr[nodeB] = f;
                ptr[f] += -1;
                index = f;
                entropy = entropy-(-(double)1.0/(double)vertices*log(1.0/(double)vertices))-(-(double)en1/(double)vertices*log((double)en1/vertices))+(-(double)(-ptr[index])/vertices*log((-ptr[index])/(double)vertices));
                if(entropy<0) entropy = 0;
                }
                else if(ptr[nodeA]!=EMPTY && ptr[nodeB]!=EMPTY)
                {

                    int g,h;
                    g = findroot(ptr,nodeA);
                    en1 = -(ptr[g]);
                    h = findroot(ptr,nodeB);
                    en2 = -(ptr[h]);
                    if(g!=h)
                    {
                        if(ptr[g]<ptr[h])
                        {

                            ptr[g] += ptr[h];
                            ptr[h] = g;
                            index = g;
                        }
                        else
                        {

                            ptr[h] += ptr[g];
                            ptr[g] = h;
                            index = h;
                        }
                        entropy = entropy-(-(double)en1/(double)vertices*log((double)en1/(double)vertices))-(-(double)en2/vertices*log((double)en2/(double)vertices))+(-(double)(-ptr[index])/vertices*log((double)(-ptr[index])/(double)vertices));
                    if(entropy<0) entropy = 0;
                    }
                    else
                    {

                        jump=big-bigtemp;
                        #pragma acc atomic update
                        maxclus[i] += big;
                        #pragma acc atomic update
                        delmx[i] += jump;
                        if(entropy<0) entropy = 0;
                        #pragma acc atomic update
                        entropycalc[i] += entropy;
                        bigtemp = big;


                        continue;
                    }
                }

                if(-ptr[index]>big) big = -ptr[index];

                jump = big - bigtemp;
                #pragma acc atomic update
                maxclus[i] += big;
                #pragma acc atomic update
                delmx[i] += jump;
                #pragma acc atomic update
                entropycalc[i] += entropy;
                bigtemp = big;


            }
        }
            /*for(int i=0;i<connectionNumber;i++)
            {cout<<maxclus[i]<<'\t'<<delmx[i]<<'\t'<<entropycalc[i]<<'\t'<<endl;}//<<"ptr "<<ptr[i]<<endl;}*/
    }

    fout1.open(filename1.c_str());
    fout2.open(filename2.c_str());
    fout3.open(filename3.c_str());

    connectionNumber = con;

    for(int i=1;i<=connectionNumber;i++)
    {
        fout1<<(double)i/vertices<<'\t'<<(double)maxclus[i]/vertices/run<<endl;
        fout2<<(double)i/vertices<<'\t'<<(double)delmx[i]/run<<endl;
        fout3<<(double)i/vertices<<'\t'<<(double)entropycalc[i]/run<<endl;
    }

    fout1.close();
    fout2.close();
    fout3.close();


    delete[] random;
    random = NULL;
    delete[] maxclus;
    maxclus = NULL;
    delete[] delmx;
    delmx = NULL;
    delete[] entropycalc;
    entropycalc = NULL;

    delete [] netmap1;
    netmap1 = NULL;
    delete [] netmap2;
    netmap2 = NULL;

    delete [] ptr;
    ptr = NULL;

    return 0;

}