Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <stdio.h>
- #include <stdlib.h>
- #include <mpi.h>
- #include <math.h>
- using namespace std;
- void init_vec(double *obj,int N, int is_mat);
- void matvec(double *vec, int N, int loc_rows, double *res, int MAX_ITER);
- int main(int argc, char* argv[]) {
- MPI_Init(&argc,&argv);
- int size = 0;
- int rank = 0;
- MPI_Comm_size(MPI_COMM_WORLD, &size);
- MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- for (int MAX_IT = 42; MAX_IT < 50; MAX_IT +=2) {
- for (int dim = 4992; dim <= 49920; dim += 4992){
- double t1 = MPI_Wtime();
- int N = 20000;
- N = dim;
- int *send;
- double *mat;
- double *vec = new double[N];
- double *result;
- if (N % size != 0) {
- MPI_Abort(MPI_COMM_WORLD,10);
- //printf("SIZE mod N_processors != 0\n");
- std::cout << "SIZE mod N_processors != 0" << std::endl;
- return(-1);
- }
- int mat_s;
- if (rank == 0) {
- mat_s = N*N;
- //mat = new double[mat_s];
- //init_vec(mat, mat_s);
- int elems_per_proc = mat_s / size;
- //init_vec(mat, N, 1);
- init_vec(vec, N, 0);
- result = new double[N];
- }
- matvec(vec, N, N/size, result, MAX_IT);
- //int *temp = new int[100];
- if (rank == 0) {
- //FILE * pfile;
- //pfile = fopen("results", "a");
- t1 = MPI_Wtime()- t1;
- //delete[] mat;
- delete[] vec;
- delete[] result;
- //fprintf(pfile, "NCORES: %d. MATRIX SIZE IS %d IT TOOK TIME: %lf\n", size, dim, t1);
- //fclose(pfile);
- //printf("NCORES %d MATRIX SIZE %d TIME %lf\n", size, dim, t1);
- cout << size << "," <<dim << "," << MAX_IT << "," <<t1 << endl;
- }
- }
- }
- MPI_Finalize();
- return(0);
- }
- void init_vec(double *obj,int N, int is_mat) {
- srand(time(NULL));
- if (is_mat == 0) {
- for (int k=0; k<N; ++k) {
- obj[k] = rand() % 10;
- }
- }
- else {
- for (int k = 0; k < N; ++k) {
- for (int p = k; p < N; ++p) {
- obj[k*N+p] = rand() % 10;
- obj[p*N+k] = obj[k*N+p];
- }
- }
- }
- }
- void dotproduct(double *a, double *b, double *res, int N) {
- double result = 0;
- for (int i = 0; i < N; ++i) {
- result += a[i]*b[i];
- }
- *res = result;
- }
- void constproduct(double *a, double b, double *res, int N) {
- for (int i = 0; i < N; ++i) {
- res[i] = a[i]*b;
- }
- }
- void sum(double *a, double *b, double *res, int N) {
- for (int i = 0; i < N; ++i) {
- res[i] = a[i] + b[i];
- }
- }
- void diff(double *a, double *b, double *res, int N) {
- for (int i = 0; i < N; ++i) {
- res[i] = a[i] - b[i];
- }
- }
- void ort(double *v,double *vec_storage, int curent_iter, int N) {
- double *c = (double*)calloc(N, sizeof(double));
- double *s = (double*)calloc(N, sizeof(double));
- double ress = 0;
- for (int j = 0; j<curent_iter;j++) {
- dotproduct(v, &vec_storage[N*j], &ress, N);
- constproduct(&vec_storage[N*j], ress, c, N);// v[j],scalarproduct(w,v[j],num,procnum),c,num,procnum);
- sum(s, c, s, N);
- }
- constproduct(s, -1, s, N);
- sum(v, s, v, N);
- }
- void matvec(double *vec, int N, int loc_rows, double *res, int MAX_ITER) {
- double t2 = MPI_Wtime();
- // int MAX_ITER = 20;
- double *local_mat = new double[loc_rows*N];
- double *local_res = new double[loc_rows];
- double *total_res;// = new double[N];
- double *vec_storage;
- int rank;
- double *a = new double[MAX_ITER];
- double *b = new double[MAX_ITER];
- double t3 = MPI_Wtime();
- MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- if (rank == 0) {
- vec_storage = new double[MAX_ITER*N];
- for (int k = 0; k < N; ++k) {
- vec_storage[k] = vec[k];
- }
- }
- //MPI_Scatter(&mat[0], loc_rows*N, MPI_DOUBLE, &local_mat[0], loc_rows*N, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- for (int ii = 0; ii < loc_rows; ii++) {
- for (int jj =0; jj < N; ++ jj) {
- local_mat[ii*N+jj]= (ii*jj) / sqrt(ii*ii + jj*jj);
- }
- }
- for (int iter = 1; iter < MAX_ITER-1; iter ++) {
- MPI_Bcast(&vec[0], N, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- t3 = MPI_Wtime() - t3;
- //std::cout << "DATA_TRANSFER_PROC_" << rank << ":" << t3 << std::endl;
- for (int i =0; i < loc_rows; ++i) {
- local_res[i] = 0;
- for (int j = 0; j<N;++j) {
- local_res[i] += local_mat[i*N+j] * vec[j];
- }
- }
- if (rank == 0) {
- total_res = new double[N];
- }
- MPI_Gather(&local_res[0], loc_rows, MPI_DOUBLE, &total_res[0], loc_rows, MPI_DOUBLE, 0, MPI_COMM_WORLD);
- if (rank == 0) {
- for (int i = 0; i < N; ++i) {
- vec[i] = total_res[i];
- }
- double *c = new double[N];
- dotproduct(&vec_storage[N*iter],total_res, &a[iter],N);
- constproduct(&vec_storage[N*iter], a[iter], c, N);
- diff(total_res,c,total_res,N);
- if (iter!=0) {
- constproduct(&vec_storage[N*(iter-1)],b[iter],c,N);
- diff(total_res ,c,total_res ,N);
- }
- ort(total_res,vec_storage,iter,N);
- ort(total_res,vec_storage,iter,N);
- double d;
- dotproduct(total_res,total_res,&d,N);
- b[iter+1]=sqrt(d);
- if (b[iter+1]==0) {
- //cout<<"just"<<iter<<endl;
- // k=iter;
- break;
- }
- constproduct(total_res,1/b[iter+1], &vec_storage[N*(iter+1)],N);
- }
- }
- if (rank == 0) {
- delete[] vec_storage;
- delete[] total_res;
- }
- delete[] local_res;
- delete[] local_mat;
- /*if (rank == 0){
- // dotproduct()
- for (int i =0; i < N; ++i) {
- res[i] = total_res[i];
- }
- }*/
- t2 = MPI_Wtime()- t2;
- //cout <<"INNER_PROC_" << rank <<": "<< t2 - t3 << endl;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement