daily pastebin goal
24%
SHARE
TWEET

Алгоритм Ланцоша с выборочной ортогонализацией

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