# Untitled

a guest
Feb 1st, 2017
111
0
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 <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. }