diff --git a/grad_conj/algebre.c b/grad_conj/algebre.c index 37371db236bb0de28468cb1f5d725d3f0a329077..8aa2588d1e584b2ed31b853fb3bbedca05bf1da5 100644 --- a/grad_conj/algebre.c +++ b/grad_conj/algebre.c @@ -21,6 +21,21 @@ void prod_mat_vec(Vector * vector_res, Matrix * matrix, Vector * vector){ } } } +// parallelisable +void prod_mat_vec_p(Vector * vector_res, Matrix * matrix, Vector * vector, int beg, int nbLigne){ + if(vector->t!=0) + exit(TRANSPOSED_VECTOR_USED); + if(matrix->width!=vector->size) + exit(UNCOMPATIBLE_ALGEBRIC_MULTIPLICATION); + if(nbLigne>vector_res->size) + exit(CONTAINER_RES_DO_NOT_FIT); + for(int i=0; i!=nbLigne; i++){ + vector_res->c[i]=0; + for(int j=0; j!=matrix->width; j++){ + vector_res->c[i]+=matrix->c[beg+i][j]*vector->c[j]; + } + } +} /* VECTOR_t x MATRIX */ @@ -47,6 +62,14 @@ void sum_vec_alpha_vec(Vector * vector_res, Vector * vector1, float alpha, Vecto for(int i=0;i!=vector1->size;i++) vector_res->c[i]=vector1->c[i]+alpha*vector2->c[i]; } +// parallelisable +void sum_vec_alpha_vec_p(Vector * vector_res, Vector * vector1, float alpha, Vector * vector2, int beg, int nbLigne){ + if(vector_res->size<nbLigne) + exit(CONTAINER_RES_DO_NOT_FIT); + for(int i=0;i!=nbLigne;i++){ + vector_res->c[i]=vector1->c[beg+i]+alpha*vector2->c[beg+i]; + } +} /* VECTOR_t x MATRIX x VECTOR */ @@ -69,13 +92,9 @@ float norm_vec_2(Vector * vector){ return res; } -/* VECTOR_t x VECTOR */ +/* VECTOR_t . VECTOR */ float prod_vec_t_vec(Vector * vector_t, Vector * vector){ - if(vector_t->t!=1) - exit(NON_TRANSPOSED_VECTOR_USED); - if(vector->t!=0) - exit(TRANSPOSED_VECTOR_USED); if(vector->size!=vector_t->size) exit(UNCOMPATIBLE_ALGEBRIC_MULTIPLICATION); float res=0; @@ -84,6 +103,16 @@ float prod_vec_t_vec(Vector * vector_t, Vector * vector){ return res; } +// parallelisable +float prod_vec_t_vec_p(Vector * vector_t, Vector * vector, int beg, int nbLigne){ + if(vector->size<beg+nbLigne) + exit(UNCOMPATIBLE_ALGEBRIC_MULTIPLICATION); + float res=0; + for(int i=beg;i!=beg+nbLigne;i++) + res+=vector->c[i]*vector_t->c[i]; + return res; +} + /* VECTOR x Scal without creating new object */ void prod_vec_scal(Vector * vector_res, Vector * vector, float scal){ diff --git a/grad_conj/algebre.h b/grad_conj/algebre.h index b1a61f8bec64460e4876257878e4749f4ee87abd..8ee4f9d6d5260bbf249d1b6f32974753827bd8b5 100644 --- a/grad_conj/algebre.h +++ b/grad_conj/algebre.h @@ -4,11 +4,15 @@ void prod_mat_vec(Vector * vector_res, Matrix * matrix, Vector * vector); void prod_vec_mat(Vector * vector_res, Vector * vector, Matrix * matrix); +void prod_mat_vec_p(Vector * vector_res, Matrix * matrix, Vector * vector, int beg, int nbLigne); + void sum_vec_alpha_vec(Vector * vector_res, Vector * vector1, float alpha, Vector * vector2); +void sum_vec_alpha_vec_p(Vector * vector_res, Vector * vector1, float alpha, Vector * vector2, int beg, int nbLigne); float prod_vec_t_mat_vec(Vector * vector, Matrix * matrix); -float prod_vec_t_vec(Vector * vector_t, Vector * vector); float norm_vec_2(Vector * vector); +float prod_vec_t_vec(Vector * vector_t, Vector * vector); +float prod_vec_t_vec_p(Vector * vector_t, Vector * vector, int beg, int nbLigne); void prod_vec_scal(Vector * vector_res, Vector * vector, float scal); diff --git a/grad_conj/bin/Debug/grad_conj.exe b/grad_conj/bin/Debug/grad_conj.exe index 9bf973f8c18faa718cf3d1a71920c3af1a2e7010..d331d272654698c2e7412c4f1c4f81e709f52501 100644 Binary files a/grad_conj/bin/Debug/grad_conj.exe and b/grad_conj/bin/Debug/grad_conj.exe differ diff --git a/grad_conj/grad_conj.depend b/grad_conj/grad_conj.depend index cce5b3e4564f2b643d4ae8c0206e138cb846a664..3a982e49c642fd6dd1fba3ed14da3a4fd3b414fa 100644 --- a/grad_conj/grad_conj.depend +++ b/grad_conj/grad_conj.depend @@ -1,15 +1,15 @@ # depslib dependency file v1.0 -1489617942 source:b:\mes documents\progra\calcint\grad_conj\mathobject.c +1489764264 source:b:\mes documents\progra\calcint\grad_conj\mathobject.c <stdlib.h> <stdio.h> "mathObject.h" "customError.h" -1489612623 b:\mes documents\progra\calcint\grad_conj\mathobject.h +1489763407 b:\mes documents\progra\calcint\grad_conj\mathobject.h 1489596942 b:\mes documents\progra\calcint\grad_conj\customerror.h -1489619150 source:b:\mes documents\progra\calcint\grad_conj\main.c +1489766951 source:b:\mes documents\progra\calcint\grad_conj\main.c <stdio.h> <stdlib.h> <time.h> @@ -21,21 +21,23 @@ "algebre.h" "iterative_methods.h" -1489617314 source:b:\mes documents\progra\calcint\grad_conj\algebre.c +1489742437 source:b:\mes documents\progra\calcint\grad_conj\algebre.c <stdlib.h> <stdio.h> "mathObject.h" "customError.h" "algebre.h" -1489617252 b:\mes documents\progra\calcint\grad_conj\algebre.h +1489742440 b:\mes documents\progra\calcint\grad_conj\algebre.h -1489619088 source:b:\mes documents\progra\calcint\grad_conj\iterative_methods.c +1489762392 source:b:\mes documents\progra\calcint\grad_conj\iterative_methods.c <stdlib.h> <stdio.h> + <math.h> + <mpi.h> "mathObject.h" "customError.h" "algebre.h" -1489612728 b:\mes documents\progra\calcint\grad_conj\iterative_methods.h +1489681364 b:\mes documents\progra\calcint\grad_conj\iterative_methods.h diff --git a/grad_conj/iterative_methods.c b/grad_conj/iterative_methods.c index a2a62ae302f404157c09766df6cc108652c361bf..d020f8f7917532968fd480c3bf000eb6930d5ca3 100644 --- a/grad_conj/iterative_methods.c +++ b/grad_conj/iterative_methods.c @@ -1,5 +1,7 @@ #include <stdlib.h> #include <stdio.h> +#include <math.h> +#include <mpi.h> #include "mathObject.h" #include "customError.h" @@ -7,6 +9,21 @@ #define MAX_ITERATION 1000000 +void grad_conj_parallel(Vector * b, Matrix * A, Vector * x, float epsilon, int nbProc, int rank, MPI_Status * Status){ + + + //printf("nb_line_for_proc : \n"); + //printVector(nb_line_for_proc); + + + //printf("first_line_for_proc : \n"); + //printVector(first_line_for_proc); + + + + +} + void grad_conj(Vector * b, Matrix * A, Vector * x, float epsilon){ int size_vec=b->size; @@ -30,9 +47,16 @@ void grad_conj(Vector * b, Matrix * A, Vector * x, float epsilon){ norm_r=norm_vec_2(r); while(norm_r>epsilon && k < MAX_ITERATION){ alpha=norm_r/prod_vec_t_mat_vec(p,A); + printf("alpha : %f\n",alpha); sum_vec_alpha_vec(x,x,alpha,p); + printf("x : \n"); + printVector(x); prod_mat_vec(Ax,A,p); + printf("Ax : \n"); + printVector(Ax); sum_vec_alpha_vec(r_next,r,-1.0*alpha,Ax); + printf("r_next : \n"); + printVector(r_next); norm_r_next=norm_vec_2(r_next); if(norm_r_next>epsilon){ beta=norm_r_next/norm_r; @@ -42,10 +66,13 @@ void grad_conj(Vector * b, Matrix * A, Vector * x, float epsilon){ } norm_r=norm_r_next; moyenne_r+=(norm_r-moyenne_r)/k; + printf("p : \n"); + printVector(p); + /* if(min_r>norm_r){ printf("k : %d, ||r|| = %10.5f, m = %f\n",k,norm_r,moyenne_r); min_r=norm_r; - } + }*/ } if(k>=MAX_ITERATION) printf("Solution has not reached the asked precision : %f in %d step.",epsilon,MAX_ITERATION); diff --git a/grad_conj/iterative_methods.c.save-failed b/grad_conj/iterative_methods.c.save-failed new file mode 100644 index 0000000000000000000000000000000000000000..c84b3b94571c323f09d9b0c4567de4bb7ded4ce5 --- /dev/null +++ b/grad_conj/iterative_methods.c.save-failed @@ -0,0 +1,72 @@ +#include <stdlib.h> +#include <stdio.h> +#include <math.h> + +#include "mathObject.h" +#include "customError.h" +#include "algebre.h" + +#define MAX_ITERATION 1000000 + +void grad_conj_parallelisable(Vector * b, Matrix * A, Vector * x, float epsilon, int nbProc, int k){ + int nb_ss_matrix_tot = (int)pow(2.,(float)k); + Vector * nb_ss_matrix_par_pb = allocateVector(nbProc); + + for(int i=0;i!=nb_ss_matrix_par_pb->size;i++){ + nb_ss_matrix_par_pb[i]=nb_ss_matrix_tot/nbProc; + if(i<nb_ss_matrix_tot%nb_ss_matrix_par_pb) + h_prob+=1; + } + printVector(nb_ss_matrix_par_pb) +} + +void grad_conj(Vector * b, Matrix * A, Vector * x, float epsilon){ + int size_vec=b->size; + + Vector * r = allocateVector(size_vec); + Vector * r_next = allocateVector(size_vec); + Vector * p = allocateVector(size_vec); + Vector * Ax = allocateVector(size_vec); + + int k=0; + float alpha; + float norm_r; + float norm_r_next; + float beta; + float min_r=1000000.0; + float moyenne_r=0; + + prod_mat_vec(Ax,A,x); + sum_vec_alpha_vec(r,b,-1.0,Ax); + copyVector(r,p); + + norm_r=norm_vec_2(r); + while(norm_r>epsilon && k < MAX_ITERATION){ + alpha=norm_r/prod_vec_t_mat_vec(p,A); + sum_vec_alpha_vec(x,x,alpha,p); + prod_mat_vec(Ax,A,p); + sum_vec_alpha_vec(r_next,r,-1.0*alpha,Ax); + norm_r_next=norm_vec_2(r_next); + if(norm_r_next>epsilon){ + beta=norm_r_next/norm_r; + copyVector(r_next,r); + sum_vec_alpha_vec(p,r,beta,p); + k+=1; + } + norm_r=norm_r_next; + moyenne_r+=(norm_r-moyenne_r)/k; + if(min_r>norm_r){ + printf("k : %d, ||r|| = %10.5f, m = %f\n",k,norm_r,moyenne_r); + min_r=norm_r; + } + } + if(k>=MAX_ITERATION) + printf("Solution has not reached the asked precision : %f in %d step.",epsilon,MAX_ITERATION); + else + printf("Solution found with %d steps.\n",k); + + freeVector(r); + freeVector(r_next); + freeVector(p); + freeVector(Ax); +} diff --git a/grad_conj/iterative_methods.h b/grad_conj/iterative_methods.h index 3f3723366b19929159c955c4d717af898c5aaa9c..4925c6a88a095be77fcc973f1178685b654dc404 100644 --- a/grad_conj/iterative_methods.h +++ b/grad_conj/iterative_methods.h @@ -2,5 +2,6 @@ #define ITERATIVE_METHODS_H_INCLUDED void grad_conj(Vector * b, Matrix * A, Vector * x, float epsilon); +void grad_conj_parallel(Vector * b, Matrix * A, Vector * x, float epsilon, int nbProc, int rank, MPI_Status * Status); #endif // ITERATIVE_METHODS_H_INCLUDED diff --git a/grad_conj/main.c b/grad_conj/main.c index d291b55070c3c2908c84a734f3250f19810284eb..eaf535384a32c8b2157996a145da0f0a53805d75 100644 --- a/grad_conj/main.c +++ b/grad_conj/main.c @@ -10,38 +10,165 @@ #include "algebre.h" #include "iterative_methods.h" +/* +Commands : +cd 'B:\Mes Documents\progra\calcInt\' +C:\MPICH2\bin\mpiexec.exe -localonly 10 +*/ + int main (int argc, char *argv[]) { int rank, size; MPI_Status status; - int test=0; + int sizeMatrix=3; + float epsilon = 0.0001; + + float norm_r_2_in_use,norm_r_2_reduced; + float norm_r_next_2_in_use,norm_r_next_2_reduced; + + MPI_Init (&argc, &argv); /* starts MPI */ + MPI_Comm_rank(MPI_COMM_WORLD, &rank); /* get current process id */ + MPI_Comm_size(MPI_COMM_WORLD, &size); /* get number of processes */ srand(time(NULL)); - testMatrix(); - testVector(); + /* definition of the problem */ + int nbProc=size; + int nb_ligne_tot = sizeMatrix; + int * nb_line_for_proc = malloc(sizeof(int)*nbProc); + int * first_line_for_proc = malloc(sizeof(int)*nbProc); + + /* initialisation of the domains */ + int nb_ligne_par_proc_default=nb_ligne_tot/nbProc; + for(int i=0;i!=nbProc;i++) + nb_line_for_proc[i]=nb_ligne_par_proc_default; + for(int i=0;i!=nb_ligne_tot-nb_ligne_par_proc_default*nbProc;i++){ + nb_line_for_proc[i]+=1; + } + first_line_for_proc[0]=0; + for(int i=1;i!=nbProc;i++){ + first_line_for_proc[i]=first_line_for_proc[i-1]+nb_line_for_proc[i-1]; + } - testAlgebric(); + int nb_line_proc_in_use = nb_line_for_proc[rank]; - Vector * b = allocateVector(10000); - Matrix * A = allocateMatrix(10000,10000); - Vector * x = allocateVector(10000); + Vector * b_tot = allocateVector(sizeMatrix); + Vector * b = allocateVector(nb_line_proc_in_use); + Matrix * A = allocateMatrix(nb_line_proc_in_use,sizeMatrix); + Vector * x_tot = allocateVector(sizeMatrix); + fillVector(x_tot,0); + Matrix * A_tot; + if(rank==0){ + A_tot = allocateMatrix(sizeMatrix,sizeMatrix); + initProbVector(b_tot,1.0); + initProbMatrixSymetric(A_tot,1.0); + for(int i=1;i!=A_tot->width;i++){ + A_tot->c[i][i]+=0;//A_tot->width; + } + //printMatrix(A_tot); + for(int p=1;p!=nbProc;p++){ + for(int i=first_line_for_proc[p]; + i!=first_line_for_proc[p]+nb_line_for_proc[p]; + i++){ + MPI_Send(A_tot->c[i],sizeMatrix,MPI_FLOAT,p,0,MPI_COMM_WORLD); + } + } + for(int i=0;i!=nb_line_proc_in_use;i++){ + A->c[i]=A_tot->c[i]; + } + }else + for(int i=0;i!=nb_line_proc_in_use;i++) + MPI_Recv(A->c[i],sizeMatrix,MPI_FLOAT,0,0,MPI_COMM_WORLD,&status); + MPI_Scatterv(b_tot->c,nb_line_for_proc,first_line_for_proc,MPI_FLOAT,b->c,nb_line_proc_in_use,MPI_FLOAT,0,MPI_COMM_WORLD); - initProbVector(b,5.0); + /** resolution of the conjugate gradient **/ + + /* initialisation - allocation */ + Vector * r = allocateVector(nb_line_proc_in_use); + Vector * p = allocateVector(nb_line_proc_in_use); + Vector * p_tot = allocateVector(sizeMatrix); + Vector * Ax = allocateVector(nb_line_proc_in_use); + Vector * x = allocateVector(nb_line_proc_in_use); fillVector(x,0); - initProbMatrixSymetric(A,5.0); - for(int i=0;i!=A->width;i++){ - A->c[i][i]+=A->width/15; + + int k=0; + float alpha; + float beta; + float pAp_in_use; + float pAp_reduced; + int done=0; + + double time_proc_in_use=MPI_Wtime(); + double time_reduced=0; + + prod_mat_vec_p(Ax,A,x_tot,0,nb_line_proc_in_use); + sum_vec_alpha_vec_p(r,b,-1.0,Ax,0,nb_line_proc_in_use); + copyVector(r,p); + norm_r_2_in_use = norm_vec_2(r); + MPI_Allreduce(&norm_r_2_in_use,&norm_r_2_reduced,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD); + const int MAX_ITERATION_MAIN=100; + while(!done && k<MAX_ITERATION_MAIN){ + MPI_Allgatherv(p->c,nb_line_proc_in_use,MPI_FLOAT,p_tot->c,nb_line_for_proc,first_line_for_proc,MPI_FLOAT,MPI_COMM_WORLD); + prod_mat_vec_p(Ax,A,p_tot,0,nb_line_proc_in_use); + /*printf("Ax :\n"); + printVector(Ax); + printf("p :\n"); + printVector(p);*/ + pAp_in_use=prod_vec_t_vec_p(p,Ax,0,nb_line_proc_in_use); + //printf("pAp : %f\n",pAp_in_use); + MPI_Allreduce(&pAp_in_use,&pAp_reduced,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD); + //printf("pAp_reduc : %f\n",pAp_reduced); + alpha=norm_r_2_reduced/pAp_reduced; + //printf("pAp_reduc : %f\n",pAp_reduced); + sum_vec_alpha_vec_p(x,x,alpha,p,0,nb_line_proc_in_use); + /*printf("x :\n"); + printVector(x);*/ + sum_vec_alpha_vec_p(r,r,-1.0*alpha,Ax,0,nb_line_proc_in_use); + /*printf("r :\n"); + printVector(r);*/ + norm_r_next_2_in_use = norm_vec_2(r); + MPI_Allreduce(&norm_r_next_2_in_use,&norm_r_next_2_reduced,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD); + if(norm_r_next_2_reduced<epsilon) + done=1; + else{ + beta=norm_r_next_2_reduced/norm_r_2_reduced; + sum_vec_alpha_vec_p(p,r,beta,p,0,nb_line_proc_in_use); + norm_r_2_reduced=norm_r_next_2_reduced; + k++; + } + /*if(rank==0){ + printf("p\n"); + printVector(p_tot); + printf("%f\n",norm_r_2_reduced); + }*/ } - /* - printf("\n\n########\nproblem : b-Ax=0 \n"); - printf("A : \n"); - printMatrix(A); - printf("b : \n"); - printVector(b);*/ - - printf("calculation : \n"); - grad_conj(b,A,x,0.000001); - - printf("grad conj result : \n"); - //printVector(x); + if(k==MAX_ITERATION_MAIN) + printf("failed to converge\n"); + MPI_Allgatherv(x->c,nb_line_proc_in_use,MPI_FLOAT,x_tot->c,nb_line_for_proc,first_line_for_proc,MPI_FLOAT,MPI_COMM_WORLD); + time_proc_in_use=MPI_Wtime()-time_proc_in_use; + MPI_Reduce(&time_proc_in_use,&time_reduced,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD); + printf("time for proc %d : %f\n",rank,time_proc_in_use); + if(rank==0){ + printSystem(A_tot,b_tot); + printf("grad conj result : in %d step and %f s\n",k,time_reduced); + printf("x resultat : \n"); + printVector(x_tot); + printf("b = Ax : \n"); + prod_mat_vec(b_tot,A_tot,x_tot); + printVector(b_tot); + } + + if(0){ + fillVector(x_tot,0.); + printf("Iterative :\n"); + grad_conj(b_tot,A_tot,x_tot,0.01); + printf("b : \n"); + printVector(b_tot); + printf("x : \n"); + printVector(x_tot); + } + + MPI_Finalize(); + + if(rank==0) + freeMatrix(A_tot); } diff --git a/grad_conj/mathObject.c b/grad_conj/mathObject.c index 5ff139ae01b5dd17856487f0c09cca03071c5489..16f2b52b8687e5deda20b40a5370e26813996b77 100644 --- a/grad_conj/mathObject.c +++ b/grad_conj/mathObject.c @@ -59,7 +59,7 @@ void freeMatrix(Matrix * matrix){ void printMatrix(Matrix * matrix){ for(int i=0; i<matrix->height; i++){ for(int j=0; j<matrix->width; j++){ - printf("%5.2f, ",matrix->c[i][j]); + printf("%f, ",matrix->c[i][j]); } printf("\n"); } @@ -101,6 +101,7 @@ Vector * allocateVector(int size){ return vector; } void initProbVector(Vector * vector, int max_value){ + rand(); for(int i=0; i < vector->size; i++){ vector->c[i]=((float)rand())/32768.0*(float)max_value; } @@ -125,7 +126,7 @@ void freeVector(Vector * vector){ void printVector(Vector * vector){ for(int i=0; i<vector->size; i++){ - printf("%5.2f \n",vector->c[i]); + printf("%f \n",vector->c[i]); } if(vector->t) printf("is transposed\n"); @@ -147,3 +148,17 @@ void testVector(){ freeVector(vector); printf("Ok\nThis is a random (from 0 to 10) 8 sized generated vector, it works well !\n"); } + +void printSystem(Matrix * A,Vector * b){ + for(int i=0;i!=A->height;i++){ + for(int j=0;j!=A->width-1;j++) + printf("x%d*%f + ",(j+1),A->c[i][j]); + printf("x%d*%f = %f,\n",A->width,A->c[i][A->width-1],b->c[i]); + } +} + +/* OTHER FUNCTIONALITY */ + +int cInt(float a){ + return (int)(a+0.00001); +} diff --git a/grad_conj/mathObject.h b/grad_conj/mathObject.h index 2799838e9cf881a2ce9f48624cec2a230f982e68..2a5b47e218b217dda24672a03baae916f947100c 100644 --- a/grad_conj/mathObject.h +++ b/grad_conj/mathObject.h @@ -32,4 +32,7 @@ void freeVector(Vector *); void printVector(Vector *); void testVector(); +int cInt(float); +void printSystem(Matrix * A,Vector * b); + #endif // MATHOBJECT_H_INCLUDED diff --git a/grad_conj/obj/Debug/algebre.o b/grad_conj/obj/Debug/algebre.o index a15ad74ba1e829337ddc7b94d2ba14f5501ddca8..2d112836bdbd1b587e1ee47f88b4a33b19932dc5 100644 Binary files a/grad_conj/obj/Debug/algebre.o and b/grad_conj/obj/Debug/algebre.o differ diff --git a/grad_conj/obj/Debug/iterative_methods.o b/grad_conj/obj/Debug/iterative_methods.o index 519f3eb08eb188be76b7fd5afa2ba6cef9faed4b..3425d61448d3975fda686a3bbf30bd29767d233a 100644 Binary files a/grad_conj/obj/Debug/iterative_methods.o and b/grad_conj/obj/Debug/iterative_methods.o differ diff --git a/grad_conj/obj/Debug/main.o b/grad_conj/obj/Debug/main.o index c5611a3ed116db893bfc856ffcd9d3680ac00b64..0afdbdb4f64a05b197a92d55d5aec4ff85645881 100644 Binary files a/grad_conj/obj/Debug/main.o and b/grad_conj/obj/Debug/main.o differ diff --git a/grad_conj/obj/Debug/mathObject.o b/grad_conj/obj/Debug/mathObject.o index 2b419c5fb1c3ebcf4acee0206d84d84d17ddf09d..6eaedbdc6b4f32cbde9bd2183c7d52a66f6eabda 100644 Binary files a/grad_conj/obj/Debug/mathObject.o and b/grad_conj/obj/Debug/mathObject.o differ