c

/*

    Skeleton file for 433-677 Project 1

    aharwood, 2007

    Student logins:

*/

#include <omp.h>

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

char *probname[]={"prefix sum",

     "matrix multiplication",

     "Gaussian elimination",

     "Sieve of Eratosthenes",

     "maximum contiguous sum"};

/* function to dynamically allocate a matrix of size n*n floats */

float **matrix(int n){

    float **M;

    int i;

    M=(float **)malloc(sizeof(float *)*n);

    if(M==NULL){

        printf("failed to allocate matrix rows\n");

        exit(1);

    }

    for(i=0;i<n;i++){

        M[i]=(float *)malloc(sizeof(float)*n);

        if(M[i]==NULL){

            printf("failed to allocate matrix columns\n");

            exit(1);

        }

    }

    return M;

}

int main(int argc, char** argv){

    int t;  // problem type

    int n;  // size of the problem

    int P;    // number of threads

    int *ia; // array for ints

    float *fa; // array for floats

    float **A,**B,**C; // matrices

    int sqrtn; // the square root of n

    int numprimes; // the number of primes <= n

    float m;    // for Guassian elimination

    float max;    // for maximum contiguous sum

    float sum;    // for maximum contiguous sum

    int start,end; // for maximum contiguous sum

    int i,j,k;  // dummy index variables

    int psd;

    int psk;

    if(argc != 4) {

        printf("usage: proj1 problem_type num_threads problem_size\n");

        printf("problem_type: 1 - prefix sum\n");

        printf("problem_type: 2 - matrix multiplication\n");

        printf("problem_type: 3 - Gaussian elimination\n");

        printf("problem_type: 4 - Sieve of Eratosthenes\n");

        printf("problem_type: 5 - maximum contiguous sum\n");

        exit(1);

    }    

    sscanf(argv[1],"%i",&t);

    sscanf(argv[2],"%i",&P);

    sscanf(argv[3],"%i",&n);        

    printf("number of processors: %d\n", omp_get_num_procs());

    omp_set_num_threads(P);

        printf("max threads: %d\n", omp_get_max_threads());

    printf("problem: %s\n",probname[t-1]);

    printf("array size %i\n",n);

    sqrtn=(int)ceil(sqrt(n));

    if(P==1){

        /*************************/

        /* Sequential algorithms */

        /*************************/

        switch(t){

        case 1:

            /* Initialize the array */

            fa=(float *)malloc(sizeof(float)*n);

            if(fa==NULL){

                printf("failed to allocate array\n");

                exit(1);

            }

            for(i=0;i<n;i++) fa[i]=i;    

            /* Compute the prefix sum */

            for(i=1;i<n;i++) fa[i]=fa[i]+fa[i-1];

            // uncomment to check

            printf("last value %f\n",fa[n-1]);

            break;

        case 2:

            /* Initialize the arrays */

            A=matrix(n);

            B=matrix(n);

            C=matrix(n);

            for(i=0;i<n;i++)

                for(j=0;j<n;j++){

                    A[i][j]=i+j;    // some dummy values

                    B[i][j]=i-j;    // some dummy values

                }

            /* Compute the matrix multiplication */

            for(i=0;i<n;i++)

                for(j=0;j<n;j++){

                      C[i][j]=0;

                      for(k=0;k<n;k++)

                        C[i][j]+=A[i][k]*B[k][j];

                }

            // uncomment to print answer for small arrays

            //for(i=0;i<n;i++)

            //    for(j=0;j<n;j++)

            //        printf("C[%d][%d]=%f\n",i,j,C[i][j]);

            break;

        case 3:

            /* Initialize the arrays */

            A=matrix(n);

            fa=(float *)malloc(sizeof(float)*n);

            if(fa==NULL){

                printf("failed to allocate array\n");

                exit(1);

            }

            for(i=0;i<n;i++){

                fa[i]=2.0*(1.0+i);    // dummy value

                for(j=0;j<n;j++)

                    A[i][j]=10.0+(i+2.0)/(j+1.0); // dummy value

            }

            /* Do Gaussian elimination */

            for(i=0;i<n-1;i++)

                  for(j=i+1;j<n;j++){

                    m=A[j][i]/A[i][i];

                    for(k=i;k<n;k++)

                          A[j][k]=A[j][k]-A[i][k]*m;

                    fa[j]=fa[j]-fa[i]*m;

                  }

            // uncomment to print resulting array for small arrays

            //for(i=0;i<n;i++)

                //for(j=0;j<n;j++)

                    //printf("A[%d][%d]=%f\n",i,j,A[i][j]);

            break;

        case 4:

            /* Initialize the array */

            ia=(int *)malloc(sizeof(int)*n);

            if(ia==NULL){

                printf("failed to allocate array\n");

                exit(1);

            }

            for(i=0;i<n;i++) ia[i]=i;

            /* Use the Sieve of Eratosthenes to find prime numbers */

            for(i=2;i<=sqrtn;i++){

                if(ia[i]>0)

                    for(j=i*i;j<n;j+=i)

                        ia[j]=0;

            }        

            //Uncomment to print the primes, for small sizes of n only

            //for(i=0;i<n;i++) if(ia[i]>0) printf("%i is a prime\n",ia[i]);

            /* Compact the array */

            numprimes=0;

            for(j=0;j<n;j++)

                if(ia[j]>0){

                    ia[numprimes]=ia[j];

                    numprimes++;

                }

            //Uncomment to print the number of primes

            //printf("%d primes less than or equal to %d\n", numprimes, n);

            break;

        case 5:

            /* Initialize the array */

            fa=(float *)malloc(sizeof(float)*n);

            if(fa==NULL){

                printf("failed to allocate array\n");

                exit(1);

            }

            srand48(100); // dummy seed value for random numbers

            for(i=0;i<n;i++) fa[i]=drand48()-0.5;

            /* Find maximum contiguous sum */

            sum=0.0;

            max=0.0;

            start=0;

            end=0;    

            i=0;

            for(j=0;j<n;j++){

                sum+=fa[j];

                if(sum>max){

                    max=sum;

                    start=i;

                    end=j;

                } else if (sum<0.0){

                    i=j+1;

                    sum=0;

                }

            }

            //uncomment to print answer    

            //printf("maximum contiguous sum = %f\n",max);

            break;        

        }

    } else {

        /***********************/

        /* Parallel algorithms */

        /***********************/

        switch(t){

        case 1:    /* Put parallel prefix sum algorithm here */

            /* Initialize the array */

            omp_set_num_threads(P);

            //printf("Parallel Prefix Sum\n");

            fa=(float *)malloc(sizeof(float)*n);

            if(fa==NULL){

                printf("failed to allocate array\n");

                exit(1);

            }

            /* Some initializations */

            for (i=0; i <n; i++)

                fa[i] = i;

            for(psd=0;psd<=(int)log2(n)-1;psd++)

            {

                int var1=(int)pow(2,psd+1);

                #pragma omp parallel for

                //for(psk=0;psk<n-1;psk+=(int)pow(2,psd+1))

                for(psk=0;psk<n-1;psk+=var1)

                {

                    fa[psk+var1-1]=fa[psk+(int)pow(2,psd)-1]+fa[psk+var1-1];

                    //printf("fa[%d]=%f\n",psk+(int)pow(2,psd+1)-1,fa[psk+(int)pow(2,psd+1)-1]);

                }

            }

            fa[n-1]=0.0;

            for(psd=(int)log2(n);psd>=0;psd--)

            {

                int var2=(int)pow(2,psd+1);

                int var3=(int)pow(2,psd);

                #pragma omp parallel for

                //for(psk=0;psk<n-1;psk+=(int)pow(2,psd+1))

                for(psk=0;psk<n-1;psk+=var2)

                {

                    float t = fa[psk+var3-1];

                    fa[psk+var3-1]=fa[psk+var2-1];

                    fa[psk+var2-1]= t + fa[psk+var2-1];

                    //printf("fa[%d]=%f\n",psk+(int)pow(2,psd)-1,fa[psk+(int)pow(2,psd+1)-1]);

                }

            }

            printf("last value %f\n",fa[n-1]);

            //for(int jj=0;jj<8;jj++)

            //{

                //printf("fa[%d]=%f\n",jj,fa[jj]);

            //}

            break;

        case 2: /* Put parallel matrix multplication algorithm here */

            break;

        case 3: /* Put parallel Gaussian elimination algorithm here */

            /* Initialize the arrays */

            omp_set_num_threads(P);

            A=matrix(n);

            fa=(float *)malloc(sizeof(float)*n);

            if(fa==NULL){

                printf("failed to allocate array\n");

                exit(1);

            }

            for(i=0;i<n;i++){

                fa[i]=2.0*(1.0+i);    // dummy value

                for(j=0;j<n;j++)

                    A[i][j]=10.0+(i+2.0)/(j+1.0); // dummy value

            }

            /* Do Gaussian elimination */

            for(i=0;i<n-1;i++)

                #pragma omp parallel for private(j,k,m)

                  for(j=i+1;j<n;j++){

                    m=A[j][i]/A[i][i];

                    //#pragma omp parallel for

                    for(k=i;k<n;k++)

                          A[j][k]=A[j][k]-A[i][k]*m;

                    fa[j]=fa[j]-fa[i]*m;

                  }

            // uncomment to print resulting array for small arrays

            //for(i=0;i<n;i++)

                //for(j=0;j<n;j++)

                    //printf("A[%d][%d]=%f\n",i,j,A[i][j]);

            break;

        case 4: /* Put parallel Sieve of Eratosthenes algoritm here */

            break;

        case 5: /* Put parallel maximum contiguous sum algorithm here */

            break;

        }

    }

    exit(0);

}

v