Advertisement
bhavyagoyal

Untitled

Nov 16th, 2015
392
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 18.89 KB | None | 0 0
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <stdlib.h>
  4. #include <sys/time.h>
  5. #include <time.h>
  6. #include <math.h>
  7. #include "evaluate.h"
  8. #include <omp.h>
  9.  
  10. void initialize(double** A, double** B, double** C, int* n);
  11. void* tune_matrix_multiply();
  12. void matrix_multiply(double* A, double* B, double* C, int n, void* tunning_params);
  13. void getUserID(char userID[20],char userName[30]);
  14.  
  15.  
  16. int main(int argc, char *argv[]){
  17.     srand(time(NULL));
  18.     void* tunning_params = tune_matrix_multiply();
  19.  
  20.     int n;
  21.     double *A, *B, *C, tolerance = pow(10,-7);
  22.     char userID[20];
  23.     char userName[30];
  24.     getUserID(userID,userName);
  25.  
  26.     if (debug_run(argc, argv,userID,userName)) {  /* This is debug run for students to test their code. Do not edit this. */
  27.         initialize(&A, &B, &C, &n);
  28.         double etOpenMPst, etOpenMPsp;
  29.         etOpenMPst = omp_get_wtime();
  30.         matrix_multiply(A, B, C, n, tunning_params); /* Multiply B and C to store the results into A <- B * C */
  31.         etOpenMPsp = omp_get_wtime();
  32.         double total_flops = myCalcFlops(etOpenMPst, etOpenMPsp, n);
  33.         printf("Flops : %lf\n",total_flops);
  34.  
  35.         int check = check_Correctness(A, B, C, n, tolerance);
  36.         if(check) printf("Correct : true \n");
  37.         else printf("Correct : false \n");
  38.         free(A);
  39.         free(B);
  40.         free(C);
  41.        
  42.     } else{
  43.         /* This is an evaluation run
  44.          * and it will be used for evaluation of student answers
  45.          * so dont wory about the code in the else block.
  46.         */
  47.         int sizes[4];
  48.  
  49.         find_matrix_sizes(argc, argv, sizes);
  50.        
  51.         double limit = 100.0, flops[4]= {0.0,0.0,0.0,0.0};
  52.         int s, correctness[4] = {0,0,0,0};
  53.         for(s=0;s<4;s++){
  54.             int n0=sizes[s];
  55.             //allocate memory of A, B, C
  56.             A = (double *) malloc(n0*n0*sizeof(double));
  57.             B = (double *) malloc(n0*n0*sizeof(double));
  58.             C = (double *) malloc(n0*n0*sizeof(double));
  59.             int j, i, k;
  60.  
  61.             for(j=0;j<5;j++){
  62.                 // Randomly initializing elements B & C from 0-limit
  63.                 /* B is a Lower Triangular Matrix and C is a n/16 banded Matrix */
  64.                 for(i=0;i<n0;i++){
  65.                     for(k=0;k<n0;k++){
  66.                         *(A + n0*i + k) = 0.0;
  67.                        
  68.                         if (k>i) *(B + n0*i + k) = 0.0;
  69.                         else *(B + n0*i + k) = (double)rand()/(double)(RAND_MAX/limit);
  70.                        
  71.                         if (k>i+n0/16-1 || i>k+n0/16-1) *(C + n0*i + k) = 0.0;
  72.                         else *(C + n0*i + k) = (double)rand()/(double)(RAND_MAX/limit);
  73.                     }
  74.                 }
  75.                
  76.                 if(j==0){
  77.                     double etOpenMPst, etOpenMPsp;
  78.                     etOpenMPst = omp_get_wtime();
  79.                     matrix_multiply(A, B, C, n0, tunning_params);
  80.                     etOpenMPsp = omp_get_wtime();
  81.                 }
  82.                 else{
  83.                     double etOpenMPst, etOpenMPsp;
  84.                     etOpenMPst = omp_get_wtime();
  85.                     matrix_multiply(A, B, C, n0, tunning_params);
  86.                     etOpenMPsp = omp_get_wtime();
  87.  
  88.                     flops[s] += myCalcFlops(etOpenMPst,etOpenMPsp, n0);
  89.                    
  90.                     int check = check_Correctness(A, B, C, n0, tolerance);
  91.                     if(check==1){
  92.                         correctness[s]++;
  93.                     }
  94.                 }
  95.             }
  96.             flops[s] = flops[s]/4.0;
  97.             free(A);
  98.             free(B);
  99.             free(C);
  100.         }
  101.         double overall_flop_score = (flops[0] + 2 * flops[1] + 3 * flops[2] + 4 * flops[3])/10;
  102.         int q,overall_correctness = 1;
  103.         for(q=0;q<4;q++){
  104.             if(correctness[q] != 4){
  105.                 overall_correctness = 0;
  106.             }
  107.         }
  108.             printf("UserName:%s\tuserID:%s\n",userName,userID);
  109.         printf("Size\tcorrectness\tflops\n");
  110.         printf("small \t\t%d\t%lf\n",correctness[0],flops[0]);
  111.         printf("medium \t\t%d\t%lf\n",correctness[1],flops[1]);
  112.         printf("large \t\t%d\t%lf\n",correctness[2],flops[2]);
  113.         printf("extralarge \t%d\t%lf\n",correctness[3],flops[3]);
  114.         printf("\n");
  115.         printf("Score = %lf\n",overall_correctness * overall_flop_score);
  116.        
  117.         FILE* outputfile=fopen("output.csv","w");
  118.         fprintf(outputfile, "\"%s\",\"%s\",%d,%lf,%d,%lf,%d,%lf,%d,%lf,%lf",userID,userName,correctness[0],flops[0],correctness[1],flops[1],correctness[2],flops[2],correctness[3],flops[3],overall_correctness * overall_flop_score);
  119.         fclose(outputfile);
  120.         FILE* outputfile2=fopen("output2.csv","w");
  121.         fprintf(outputfile, "\"%s\",%d,%lf,%d,%lf,%d,%lf,%d,%lf,%lf",userName,correctness[0],flops[0],correctness[1],flops[1],correctness[2],flops[2],correctness[3],flops[3],overall_correctness * overall_flop_score);
  122.         fclose(outputfile2);
  123.     }
  124.     return 0;
  125. }
  126. //Add only code in the functions
  127. void getUserID(char userID[20],char userName[30]){
  128.     //Enter your moodle userID here i.e. your entry number
  129.     strcpy(userID,"cs1120202");
  130.     //Enter your name here
  131.     strcpy(userName,"Aashish Goyal");
  132. }
  133.  
  134. void initialize(double** A, double** B, double** C, int* n){
  135.     /*
  136.     Initialize matrix A , B & C where each of them is a square matrix of size n*n
  137.     Try 4 sizes of matrix:
  138.     a) small - should fit in the L1 cache
  139.     b) medium - should fit in the L2 cache
  140.     b) large - should fit in the L3 cache
  141.     b) extra-large - should fit in the memory
  142.     */
  143.     *n=1024;
  144.     *A = (double *) malloc(*n**n*sizeof(double));
  145.     *B = (double *) malloc(*n**n*sizeof(double));
  146.     *C = (double *) malloc(*n**n*sizeof(double));
  147.     int i, k;
  148.     for(i=0;i<*n;i++){
  149.         for(k=0;k<*n;k++){
  150.             *(*A + *n*i + k) = 0.0;
  151.                        
  152.             if(k>i) *(*B + *n*i + k) = 0.0;
  153.             else *(*B + *n*i + k) = (double)rand()/(double)(RAND_MAX/100.0);
  154.                        
  155.             if(k>i+*n/16-1 || i>k+*n/16-1) *(*C + *n*i + k) = 0.0;
  156.             else *(*C + *n*i + k) = (double)rand()/(double)(RAND_MAX/100.0);
  157.         }
  158.                    
  159.     }
  160. }
  161.  
  162. void* tune_matrix_multiply(){
  163.     /* Write code to test various methods to optimize the matrix multiplication
  164.        and store these parameters in a data structure(for ex, array of values ). This
  165.        data structure can be returned by this function and will be passed to the
  166.        matrix multiply function so that it can efficiently multiply the matrices.
  167.     */
  168.     printf("Number of procs %d \n",omp_get_num_procs() );
  169.     // omp_set_num_threads(omp_get_num_procs() - 2);
  170.    
  171.    
  172.        
  173.    
  174.         int n = 1024;
  175.    
  176.         int i,j,k;
  177.         int* mapper = (int*)malloc(sizeof(int) *(n*(n+1)));
  178.  
  179.  
  180.         #pragma omp parallel for private(i,j,k) shared(n, mapper) default(none) collapse(1)
  181.         for(i = 0; i < n; i++){
  182.             for(j = 0, k = (i*(i+1))/2 ; j < i+1; j++){
  183.                 mapper[(k+j)*2] = i;
  184.                 mapper[(k+j)*2 + 1] = j;
  185.             }
  186.         }
  187.  
  188.         mapper[0] = omp_get_num_procs() + 1;
  189.        
  190.         // return (void*)mapper;
  191.         return mapper;
  192.    
  193.    
  194. }
  195.  
  196.  
  197. #define CR(x,y) ((x*n)+y)
  198. #define CR2(x,y) ((x<<d)+y)
  199.  
  200. #define MIN(x,y) (((x)<(y)) ? (x):(y) )
  201. #define MAX(x,y) (((x)>(y)) ? (x):(y) )
  202.  
  203.  
  204.  
  205. void matrix_multiply(double* A, double* B, double* C, int n, void* tunning_params){
  206.     // Multiply matrix B & C and store the results in A
  207.  
  208.         // CORRECTNESS IS ENSURED AS PARALLELISM ONLY ON i
  209.  
  210.  
  211.             // printf("number of threads %d \n", omp_get_num_threads());
  212.  
  213.         // on gcl computer able to get 23 Gflops with this code
  214.  
  215.         // register int i,j,k;
  216.         // #pragma omp parallel private(i,j,k) shared(A,B,C,n) default(none)
  217.         // {
  218.         //     #pragma omp for schedule(auto)
  219.         //     for(i = 0; i < n; i++){
  220.         //         for(k = 0; k < n; k++){
  221.         //             for(j = 0; j < n; j++){
  222.         //                 A[CR(i,j)] += B[CR(i,k)] * C[CR(k,j)];
  223.         //             }
  224.         //         }
  225.         //     }
  226.         // }
  227.  
  228.         // int i,j,k;
  229.         // double *D, *E, *F;
  230.    
  231.         // printf("%d\n" , n);
  232.        
  233.         // #pragma omp parallel private(i,j,k,D,E,F) shared(A,B,C,n) default(none)
  234.         // {
  235.         //     // #pragma omp for schedule(auto)
  236.         //     #pragma omp parallel for schedule(static) private(i,j,k,D,E,F) shared(A,B,C,n) default(none) collapse(2)
  237.         //     for(i = 0; i < n; i++){
  238.         //         for(k = 0, E = (B + (i*n)), F = C; k < n; k++, E++){
  239.         //             for(j = 0, D = (A + (i*n)); j < n; j++,D++,F++){
  240.         //                     *D += (*E) * (*F);
  241.         //             }
  242.         //         }
  243.         //     }
  244.         // }
  245.  
  246.    
  247.         int *tune = tunning_params;
  248.  
  249.  
  250.         // current best
  251.         int  i,j,k;
  252.         int m = n;
  253.        
  254.  
  255.        
  256.         // fresh code
  257.        
  258.         switch (n){
  259.             case 1024:{
  260.                 // #pragma omp parallel for private(i,j,k) shared(A,B,C,n) default(none) schedule(guided)
  261.    
  262.                 // #pragma omp parallel for schedule(dynamic) num_threads(tune[0])
  263.                 // for(i = n-1; i >= 0; i--){
  264.                 //     F = A + (i << d);
  265.                 //     G = B + (i << d);
  266.                 //     for(k = 0; k < i + 1; k++){
  267.                 //         // double reg = B[CR2(i,k)];
  268.                 //         double reg = *(G + k);
  269.                 //         const int val1 = MAX(0, k-n_16+1);
  270.                 //         const int val2 = MIN(n, k+n_16);
  271.                 //         for(j = val1, E = F + val1, H = C + CR2(k,val1); j < val2; j++, E++, H++){
  272.                 //             // *E += reg * C[CR2(k,j)];
  273.                 //             *E += reg * (*H);
  274.                 //         }
  275.                 //     }
  276.                 // }
  277.    
  278.                 int n_16 = 1024/16;
  279.                 const int n = 1024;
  280.                 #pragma omp parallel for schedule(dynamic) num_threads(omp_get_num_procs()+1)
  281.                 for(i = n-1; i >= 0; i--){
  282.                     for( k = 0; k <= (i); k++){
  283.                         const int val1 = MAX(0, k-n_16+1);
  284.                         const int val2 = MIN(n, k+n_16);
  285.                         // for( j = 0; j < (val2-val1); j++){
  286.                         for( j = val1; j < (val2); j++){
  287.                             A[(i*n)+(j)] = A[(i*n)+(j)] + B[(i*n)+k]*C[(k*n)+(j)];
  288.                         }
  289.                     }
  290.                 }
  291.                
  292.                
  293.                 return;            
  294.             }
  295.             case 256:{
  296.                 // const int n = 256;
  297.                 int d = 8;
  298.                 // const int p = 256/16;
  299.                 #pragma omp parallel for collapse(1) num_threads(9) schedule(dynamic)
  300.                 for(i = 0; i < n; i++){
  301.                     for(k = 0; k < i+1; k++){
  302.                         for(j = MAX(0, k-n/16+1); j < MIN(n, k+n/16); j++){
  303.                             A[CR(i,j)] += B[CR(i,k)] * C[CR(k,j)];
  304.                         }
  305.                     }
  306.                 }
  307.                 return;
  308.             }
  309.             case 64:{
  310.                 int l , q , p;
  311.                 for(i = 0 ; i < n ; i++){
  312.                     for(j = 0 , l = MIN(i+2,n) ,q = MIN(i+3,n) , p = n*i ; j < i + 1 ; j++){
  313.                         A[p+j] += B[p + j] * C[n*j + j];
  314.                         if(i && j < i)  A[p + j] += B[p + j + 1] * C[n*(j+1) + j];
  315.                         if(i-1 && j < i-1)  A[p + j] += B[p + j + 2] * C[n*(j+2) + j];
  316.                         if(i-2 && j < i-2)  A[p + j] += B[p + j + 3] * C[n*(j+3) + j];
  317.                         if(j + 1 && j + 1 < l)  A[p + j + 1] += B[p + j]*C[n*j + j + 1];
  318.                         if(j + 2 && j + 2 < q)  A[p + j + 2] += B[p + j]*C[n*j + j + 2];
  319.                         if(j + 3 && j + 3 < MIN(i+4,n))  A[p + j + 3] += B[p + j]*C[n*j + j + 3];
  320.                     }
  321.                 }
  322.                 return;
  323.             }
  324.             case 32:{
  325.                 for(i = 0 ; i < n ; i++){
  326.                     int l, p ;
  327.    
  328.                     int d = 5;
  329.    
  330.                     for(j = 0 , l = MIN(i+2,n) , p = i<<d ; j < i + 1 ; j++){
  331.                         A[p+j] += B[p + j] * C[n*j + j];
  332.                         if(i && j < i)  A[p + j] += B[p + j + 1] * C[((j+1)<<d) + j];
  333.                         if(j + 1 && j + 1 < l)  A[p + j + 1] += B[p + j]*C[(j<<d) + j + 1];
  334.                     }
  335.                 }
  336.                 return;
  337.             }
  338.             default:{
  339.                 const int p = n/16;
  340.                 #pragma omp parallel for private(i,j,k) shared(A,B,C,n) default(none) collapse(1) schedule(dynamic) num_threads(tune[0])
  341.                 for(i = 0; i < n; i++){
  342.                     for(k = 0; k < i+1; k++){
  343.                         for(j = MAX(0, k-p+1); j < MIN(n, k+p); j++){
  344.                             A[CR(i,j)] += B[CR(i,k)] * C[CR(k,j)];
  345.                         }
  346.                     }
  347.                 }
  348.                 return;
  349.             }
  350.    
  351.         }
  352.  
  353.         // if (n > 300)
  354.         // {
  355.         //     const int n  = m;
  356.         //     const int n_16 = n/16;
  357.         //     // #pragma omp parallel for private(i,j,k) shared(A,B,C,n) default(none) schedule(guided)
  358.         //     #pragma omp parallel for schedule(dynamic) num_threads(tune[0])
  359.         //     for(i = n-1; i >= 0 ; i--){
  360.         //         for(k = 0; k < i+1; k++){
  361.         //             double register reg = B[CR(i,k)];
  362.         //             const int val1 = MAX(0, k-n_16+1);
  363.         //             const int val2 = MIN(n, k+n_16);
  364.         //             for(j = val1; j < val2; j++){
  365.         //                 A[CR(i,j)] += reg * C[CR(k,j)];
  366.         //             }
  367.         //         }
  368.         //     }
  369.         // }
  370.         // else if (n > 100){
  371.         //     const int p = n/16;
  372.         //     #pragma omp parallel for private(i,j,k) shared(A,B,C,n) default(none) collapse(1) schedule(dynamic) num_threads(tune[0])
  373.         //     for(i = 0; i < n; i++){
  374.         //         for(k = 0; k < i+1; k++){
  375.         //             for(j = MAX(0, k-p+1); j < MIN(n, k+p); j++){
  376.         //                 A[CR(i,j)] += B[CR(i,k)] * C[CR(k,j)];
  377.         //             }
  378.         //         }
  379.         //     }
  380.         // }
  381.         // else {
  382.  
  383.         //     // #pragma omp parallel for private(i,j,k) shared(A,B,C,n) default(none) collapse(1) schedule(dynamic)
  384.         //     // for(i = 0; i < n; i++){
  385.         //     //     for(k = 0; k < i+1; k++){
  386.         //     //         for(j = MAX(0, k-n/16+1); j < MIN(n, k+n/16); j++){
  387.         //     //             A[CR(i,j)] += B[CR(i,k)] * C[CR(k,j)];
  388.         //     //         }
  389.         //     //     }
  390.         //     // }
  391.  
  392.         //     if(n == 32){
  393.         //             for(i = 0 ; i < n ; i++){
  394.         //                 int l, p ;
  395.         //                 int d = 0;
  396.         //                 int m = 1;
  397.         //                 while(m < n){
  398.         //                     d += 1;
  399.         //                     m *= 2;
  400.         //                 }
  401.                        
  402.         //                 for(j = 0 , l = MIN(i+2,n) , p = i<<d ; j < i + 1 ; j++){
  403.         //                     A[p+j] += B[p + j] * C[n*j + j];
  404.         //                     if(i && j < i)  A[p + j] += B[p + j + 1] * C[((j+1)<<d) + j];
  405.         //                     if(j + 1 && j + 1 < l)  A[p + j + 1] += B[p + j]*C[(j<<d) + j + 1];
  406.         //                 }
  407.         //             }
  408.         //             return;
  409.         //         }
  410.              
  411.         //     else if(n == 64){
  412.         //             int l , q , p;
  413.         //             for(i = 0 ; i < n ; i++){
  414.         //                 for(j = 0 , l = MIN(i+2,n) ,q = MIN(i+3,n) , p = n*i ; j < i + 1 ; j++){
  415.         //                     A[p+j] += B[p + j] * C[n*j + j];
  416.         //                     if(i && j < i)  A[p + j] += B[p + j + 1] * C[n*(j+1) + j];
  417.         //                     if(i-1 && j < i-1)  A[p + j] += B[p + j + 2] * C[n*(j+2) + j];
  418.         //                     if(i-2 && j < i-2)  A[p + j] += B[p + j + 3] * C[n*(j+3) + j];
  419.         //                     if(j + 1 && j + 1 < l)  A[p + j + 1] += B[p + j]*C[n*j + j + 1];
  420.         //                     if(j + 2 && j + 2 < q)  A[p + j + 2] += B[p + j]*C[n*j + j + 2];
  421.         //                     if(j + 3 && j + 3 < MIN(i+4,n))  A[p + j + 3] += B[p + j]*C[n*j + j + 3];
  422.         //                 }
  423.         //             }
  424.         //             return;
  425.         //         }
  426.         //     else{
  427.         //         #pragma omp parallel for private(i,j,k) shared(A,B,C,n) default(none) collapse(1) schedule(dynamic)
  428.         //         for(i = 0; i < n; i++){
  429.         //             for(k = 0; k < i+1; k++){
  430.         //                 for(j = MAX(0, k-n/16+1); j < MIN(n, k+n/16); j++){
  431.         //                     A[CR(i,j)] += B[CR(i,k)] * C[CR(k,j)];
  432.         //                 }
  433.         //             }
  434.         //         }
  435.         //     }
  436.  
  437.         // }
  438.  
  439.  
  440.         // next
  441.         // int i,j,k;
  442.         // // int mapper[(n*(n+1))];
  443.  
  444.  
  445.         // // #pragma omp parallel for private(i,j,k) shared(n, mapper) default(none) collapse(1)
  446.         // // for(i = 0; i < n; i++){
  447.         // //     for(j = 0, k = (i*(i+1))/2 ; j < i+1; j++){
  448.         // //         mapper[(k+j)*2] = i;
  449.         // //         mapper[(k+j)*2 + 1] = j;
  450.         // //     }
  451.         // // }
  452.  
  453.         // int *mapper  = (int*) tunning_params;
  454.        
  455.         // // for(i = 0; i < 10; i++){
  456.         // //     printf("%d ", mapper[2*i + 1]);
  457.         // // }printf("\n");
  458.  
  459.         // int l;
  460.         // int m;
  461.         // #pragma omp parallel for private(i,j,k,l) shared(A,B,C,n,mapper) default(none) collapse(1) schedule()
  462.         // for(l = 0; l < n*(n+1)/2; l++){
  463.         //     i = mapper[2*(l)];
  464.         //     k = mapper[2*(l) + 1];
  465.            
  466.         //     for(j = MAX(0, k-n/16+1); j < MIN(n, k+n/16); j++){
  467.         //         #pragma omp atomic
  468.         //         A[CR(i,j)] += B[CR(i,k)] * C[CR(k,j)];
  469.         //     }
  470.         // }
  471.  
  472.        
  473.        
  474.  
  475.  
  476.         // int i,j,k;
  477.  
  478.         // for(i = 0; i < n; i++){
  479.         //     for(j = MAX(0, i - n/16 + 1); j < i; j++){
  480.         //         double tmp = C[CR(i,j)];
  481.         //         C[CR(i,j)] = C[CR(j,i)];
  482.         //         C[CR(j,i)] = tmp;
  483.         //     }
  484.         // }
  485.  
  486.         // {
  487.         //     #pragma omp parallel for private(i,j,k) shared(A,B,C,n) default(none) collapse(2) schedule(dynamic)
  488.         //     for(i = 0; i < n; i++){
  489.         //         for( j = 0; j < (n); j++ ){
  490.         //             double register val = 0;
  491.         //             for( k = MAX(0, j-n/16+1); k < MIN(i+1, j+n/16); k++){
  492.         //                 val += B[CR(i,k)] * C[CR(j,k)];
  493.         //             }
  494.         //             A[CR(i,j)] = val;
  495.         //         }
  496.         //     }
  497.         // }
  498.  
  499.         // for(i = 0; i < n; i++){
  500.         //     for(j = MAX(0, i - n/16 + 1); j < i; j++){
  501.         //         double tmp = C[CR(i,j)];
  502.         //         C[CR(i,j)] = C[CR(j,i)];
  503.         //         C[CR(j,i)] = tmp;
  504.         //     }
  505.         // }
  506.            
  507.         // int i,j,k;
  508.  
  509.         // int m = n;
  510.         // {
  511.         //     const int n = m;
  512.         //     double *D, *E, *F;
  513.    
  514.         //     #pragma omp parallel for default(none) private(i, j, k, D, E, F) shared(A, B, C)
  515.         //     for(i = 0; i < n ; i++){
  516.         //             for(k = 0, E = (B + (i*n)); k < i + 1; k++, E++){
  517.         //                     for (j = MAX(0, k-n/16+1), D = (A + (i * n) + j), F = (C + (k * n) +j) ; j < MIN(n,k+n/1); j++, D++, F++){
  518.         //                             *D += (*E) * (*F);
  519.         //                     }
  520.         //             }
  521.         //     }
  522.            
  523.         // }
  524.  
  525.  
  526. }
  527.  
  528. // #pragma GCC pop_options
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement