Advertisement
Guest User

Untitled

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