Advertisement
Guest User

LanzochMPI

a guest
Dec 7th, 2016
79
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 5.31 KB | None | 0 0
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <mpi.h>
  4. #include <math.h>
  5. using namespace std;
  6. void init_vec(double *obj,int N, int is_mat);
  7. void matvec(double *vec, int N, int loc_rows, double *res, int MAX_ITER);
  8.  
  9. int main(int argc, char* argv[]) {
  10.     MPI_Init(&argc,&argv);
  11.    
  12.     int size = 0;
  13.     int rank = 0;
  14.     MPI_Comm_size(MPI_COMM_WORLD, &size);
  15.     MPI_Comm_rank(MPI_COMM_WORLD, &rank);  
  16.     for (int MAX_IT = 42; MAX_IT < 50; MAX_IT +=2) {
  17.     for (int dim = 4992; dim <= 49920; dim += 4992){
  18.             double t1 = MPI_Wtime();
  19.             int N = 20000;
  20.             N = dim;
  21.             int *send;
  22.             double *mat;
  23.             double *vec = new double[N];      
  24.             double *result;
  25.             if (N % size != 0) {
  26.                 MPI_Abort(MPI_COMM_WORLD,10);
  27.                 //printf("SIZE mod N_processors != 0\n");  
  28.                 std::cout << "SIZE mod N_processors != 0" << std::endl;
  29.                 return(-1);
  30.             }
  31.             int mat_s;
  32.             if (rank == 0) {
  33.                 mat_s = N*N;
  34.                 //mat = new double[mat_s];
  35.                 //init_vec(mat, mat_s);
  36.                 int elems_per_proc = mat_s / size;
  37.                 //init_vec(mat, N, 1);
  38.                 init_vec(vec, N, 0);
  39.                 result = new double[N];
  40.             }
  41.             matvec(vec, N, N/size, result, MAX_IT);
  42.             //int *temp = new int[100];
  43.             if (rank == 0) {
  44.                 //FILE * pfile;
  45.                 //pfile = fopen("results", "a");
  46.                 t1 = MPI_Wtime()- t1;
  47.                 //delete[] mat;
  48.                 delete[] vec;
  49.                 delete[] result;
  50.                 //fprintf(pfile, "NCORES: %d. MATRIX SIZE IS %d IT TOOK TIME: %lf\n", size, dim, t1);
  51.                 //fclose(pfile);
  52.                 //printf("NCORES %d MATRIX SIZE %d TIME %lf\n", size, dim, t1);
  53.                 cout <<  size << "," <<dim << "," << MAX_IT << "," <<t1 << endl;
  54.             }
  55.  
  56.     }
  57.     }
  58.  
  59.     MPI_Finalize();
  60.     return(0);
  61. }
  62.  
  63. void init_vec(double *obj,int N, int is_mat) {
  64.     srand(time(NULL));
  65.     if (is_mat == 0) {
  66.         for (int k=0; k<N; ++k) {
  67.             obj[k] = rand() % 10;
  68.         }
  69.     }
  70.     else {
  71.         for (int k = 0; k < N; ++k) {
  72.             for (int p = k; p < N; ++p) {
  73.                 obj[k*N+p] = rand() % 10;
  74.                 obj[p*N+k] = obj[k*N+p];
  75.             }
  76.         }
  77.     }
  78. }
  79.  
  80.  
  81. void dotproduct(double *a, double *b, double *res, int N) {
  82.     double result = 0;
  83.     for (int i = 0; i < N; ++i) {
  84.         result += a[i]*b[i];
  85.     }
  86.     *res = result;
  87. }
  88.  
  89. void constproduct(double *a, double b, double *res, int N) {
  90.     for (int i = 0; i < N; ++i) {
  91.         res[i] = a[i]*b;
  92.     }
  93. }
  94.  
  95. void sum(double *a, double *b, double *res, int N) {
  96.     for (int i = 0; i < N; ++i) {
  97.         res[i] = a[i] + b[i];
  98.     }
  99. }
  100.  
  101. void diff(double *a, double *b, double *res, int N) {
  102.     for (int i = 0; i < N; ++i) {
  103.         res[i] = a[i] - b[i];
  104.     }
  105. }
  106.  
  107. void ort(double *v,double *vec_storage, int curent_iter, int N) {
  108.     double *c = (double*)calloc(N, sizeof(double));
  109.     double *s = (double*)calloc(N, sizeof(double));
  110.     double ress = 0;
  111.     for (int j = 0; j<curent_iter;j++) {
  112.         dotproduct(v, &vec_storage[N*j], &ress, N);
  113.         constproduct(&vec_storage[N*j], ress, c, N);// v[j],scalarproduct(w,v[j],num,procnum),c,num,procnum);
  114.         sum(s, c, s, N);
  115.     }
  116.     constproduct(s, -1, s, N);
  117.     sum(v, s, v, N);
  118. }
  119.  
  120.  
  121. void matvec(double *vec, int N, int loc_rows, double *res, int MAX_ITER) {    
  122.     double t2 = MPI_Wtime();
  123. //  int MAX_ITER = 20;
  124.     double *local_mat = new double[loc_rows*N];
  125.     double *local_res = new double[loc_rows];
  126.     double *total_res;// = new double[N];
  127.     double *vec_storage;
  128.     int rank;
  129.  
  130.     double *a = new double[MAX_ITER];  
  131.     double *b = new double[MAX_ITER];
  132.  
  133.  
  134.     double t3 = MPI_Wtime();
  135.     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  136.     if (rank == 0) {
  137.         vec_storage = new double[MAX_ITER*N];
  138.         for (int k = 0; k < N; ++k) {
  139.             vec_storage[k] = vec[k];
  140.         }
  141.     }
  142.  
  143.     //MPI_Scatter(&mat[0], loc_rows*N, MPI_DOUBLE, &local_mat[0], loc_rows*N, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  144.     for (int ii = 0; ii < loc_rows; ii++) {
  145.         for (int jj =0; jj < N; ++ jj) {
  146.             local_mat[ii*N+jj]= (ii*jj) / sqrt(ii*ii + jj*jj);
  147.         }
  148.     }  
  149.  
  150.  
  151.     for (int iter = 1; iter < MAX_ITER-1; iter ++) {
  152.         MPI_Bcast(&vec[0], N, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  153.         t3 = MPI_Wtime() - t3;
  154.         //std::cout << "DATA_TRANSFER_PROC_" << rank << ":" << t3 << std::endl;
  155.         for (int i =0; i < loc_rows; ++i) {
  156.             local_res[i] = 0;
  157.             for (int j = 0; j<N;++j) {
  158.                 local_res[i] += local_mat[i*N+j] * vec[j];
  159.             }
  160.         }
  161.         if (rank == 0) {
  162.             total_res = new double[N];     
  163.         }
  164.         MPI_Gather(&local_res[0], loc_rows, MPI_DOUBLE, &total_res[0], loc_rows, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  165.         if (rank == 0) {
  166.             for (int i = 0; i < N; ++i) {
  167.                 vec[i] = total_res[i];
  168.             }
  169.             double *c = new double[N];
  170.  
  171.  
  172.             dotproduct(&vec_storage[N*iter],total_res, &a[iter],N);
  173.            
  174.             constproduct(&vec_storage[N*iter], a[iter], c, N);
  175.             diff(total_res,c,total_res,N);
  176.  
  177.             if (iter!=0) {
  178.                 constproduct(&vec_storage[N*(iter-1)],b[iter],c,N);
  179.                 diff(total_res ,c,total_res ,N);
  180.             }
  181.             ort(total_res,vec_storage,iter,N);
  182.             ort(total_res,vec_storage,iter,N);
  183.             double d;
  184.             dotproduct(total_res,total_res,&d,N);
  185.             b[iter+1]=sqrt(d);
  186.             if (b[iter+1]==0) {
  187.                 //cout<<"just"<<iter<<endl;
  188.                 // k=iter;
  189.                 break;
  190.             }
  191.             constproduct(total_res,1/b[iter+1], &vec_storage[N*(iter+1)],N);
  192.         }
  193.     }
  194.     if (rank == 0) {
  195.         delete[] vec_storage;
  196.         delete[] total_res;
  197.  
  198.     }  
  199.     delete[] local_res;
  200.     delete[] local_mat;
  201.  
  202.         /*if (rank == 0){
  203. //          dotproduct()
  204.             for (int i =0; i < N; ++i) {
  205.                 res[i] = total_res[i];
  206.             }
  207.         }*/
  208.     t2 = MPI_Wtime()- t2;
  209.     //cout <<"INNER_PROC_" << rank <<": "<< t2 - t3 << endl;
  210.  
  211. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement