diff --git a/grad_conj/bin/Debug/grad_conj.exe b/grad_conj/bin/Debug/grad_conj.exe
index d331d272654698c2e7412c4f1c4f81e709f52501..1a0a7040d85ca72ac14c6265ffaeb0943bd13721 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 3a982e49c642fd6dd1fba3ed14da3a4fd3b414fa..07345927b86b95938db746aa1a23f4c5c483b857 100644
--- a/grad_conj/grad_conj.depend
+++ b/grad_conj/grad_conj.depend
@@ -9,7 +9,7 @@
1489596942 b:\mes documents\progra\calcint\grad_conj\customerror.h
-1489766951 source:b:\mes documents\progra\calcint\grad_conj\main.c
+1489771819 source:b:\mes documents\progra\calcint\grad_conj\main.c
<stdio.h>
<stdlib.h>
<time.h>
diff --git a/grad_conj/main.c b/grad_conj/main.c
index eaf535384a32c8b2157996a145da0f0a53805d75..d48a216d9c5054886c10c5535c41f41ecf2e3417 100644
--- a/grad_conj/main.c
+++ b/grad_conj/main.c
@@ -17,15 +17,25 @@ C:\MPICH2\bin\mpiexec.exe -localonly 10
*/
int main (int argc, char *argv[]) {
+ // MPI variable
int rank, size;
MPI_Status status;
- int sizeMatrix=3;
+ // parameters of the problem
+ 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;
+ // get user asked size of the problem
+ if(argc == 2){
+ int bufEntU=atoi(argv[1]);
+ if(bufEntU>0){
+ sizeMatrix=bufEntU;
+ if(rank==0)
+ printf("You have defined a %d sized problem.\n",sizeMatrix);
+ }
+ }
+ // initialisation of the MPI context
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 */
@@ -37,7 +47,8 @@ int main (int argc, char *argv[]) {
int * nb_line_for_proc = malloc(sizeof(int)*nbProc);
int * first_line_for_proc = malloc(sizeof(int)*nbProc);
- /* initialisation of the domains */
+ /* initialisation of the domains according to his size
+ and the number of processus */
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;
@@ -51,24 +62,27 @@ int main (int argc, char *argv[]) {
int nb_line_proc_in_use = nb_line_for_proc[rank];
+ /** Creating problem and scattering it to each processus **/
+
+ /* System's problem variables */
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){
+
+ if(rank==0){ // only the root create the problem
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;
+ A_tot->c[i][i]+=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++){
+ i++){ // Scattering of the Matrix
MPI_Send(A_tot->c[i],sizeMatrix,MPI_FLOAT,p,0,MPI_COMM_WORLD);
}
}
@@ -76,9 +90,9 @@ int main (int argc, char *argv[]) {
A->c[i]=A_tot->c[i];
}
}else
- for(int i=0;i!=nb_line_proc_in_use;i++)
+ for(int i=0;i!=nb_line_proc_in_use;i++) // receiving scattered matrix from the root
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);
+ 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);
/** resolution of the conjugate gradient **/
@@ -90,85 +104,95 @@ int main (int argc, char *argv[]) {
Vector * x = allocateVector(nb_line_proc_in_use);
fillVector(x,0);
+ /* intern problem's variables */
int k=0;
float alpha;
float beta;
float pAp_in_use;
float pAp_reduced;
+ float norm_r_2_in_use,norm_r_2_reduced;
+ float norm_r_next_2_in_use,norm_r_next_2_reduced;
int done=0;
double time_proc_in_use=MPI_Wtime();
double time_reduced=0;
+ //r0=b-Ax0
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);
+ //p0=r0
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;
+ k=0;
while(!done && k<MAX_ITERATION_MAIN){
+ // get all p to be able to do A.p
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);
+ /** alpha = <Rk,Rk>/<Pk,A.Pk> **/
+ /* <Pk,A.Pk> */
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);
+ // <r,r> already done
alpha=norm_r_2_reduced/pAp_reduced;
- //printf("pAp_reduc : %f\n",pAp_reduced);
+ //Xk+1=Xk+alpha*Pk
sum_vec_alpha_vec_p(x,x,alpha,p,0,nb_line_proc_in_use);
- /*printf("x :\n");
- printVector(x);*/
+ //Rk+1=Rk-alpha*A.Pk
sum_vec_alpha_vec_p(r,r,-1.0*alpha,Ax,0,nb_line_proc_in_use);
- /*printf("r :\n");
- printVector(r);*/
+ // get <Rk+1,Rk+1> to check epsilon end condition
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 the norm less than epsilon, it's done
if(norm_r_next_2_reduced<epsilon)
done=1;
else{
+ // Beta = <Rk+1,Rk+1>/<Rk,Rk> every thing is already calculated
beta=norm_r_next_2_reduced/norm_r_2_reduced;
+ // Pk+1 = Rk+1 + Beta*Pk
sum_vec_alpha_vec_p(p,r,beta,p,0,nb_line_proc_in_use);
+ // Rk = Rk+1 : avoiding some useless calcul
norm_r_2_reduced=norm_r_next_2_reduced;
+ // k = k + 1
k++;
}
- /*if(rank==0){
- printf("p\n");
- printVector(p_tot);
- printf("%f\n",norm_r_2_reduced);
- }*/
}
if(k==MAX_ITERATION_MAIN)
printf("failed to converge\n");
+ // Reconstruct full x solution
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);
+ // get effective execution time of the process
time_proc_in_use=MPI_Wtime()-time_proc_in_use;
+ // reduce time's execution of each process
MPI_Reduce(&time_proc_in_use,&time_reduced,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
+ /** Print some results information **/
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");
+ printf("x result : \n");
printVector(x_tot);
- printf("b = Ax : \n");
+ printf("recaculating b = Ax gives : \n");
prod_mat_vec(b_tot,A_tot,x_tot);
printVector(b_tot);
}
+ // we may want to free some stuff isn'it ?
- 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);
+ if(rank==0){
+ freeMatrix(A_tot);
+ freeVector(b_tot);
+ freeVector(x_tot);
}
-
+ freeMatrix(A);
+ freeVector(b_tot);
+ freeVector(b);
+ freeVector(x_tot);
+ freeVector(x);
+ freeVector(p_tot);
+ freeVector(p);
+ freeVector(r);
+ freeVector(Ax);
+
+ // let's finalize
MPI_Finalize();
-
- if(rank==0)
- freeMatrix(A_tot);
}
diff --git a/grad_conj/obj/Debug/main.o b/grad_conj/obj/Debug/main.o
index 0afdbdb4f64a05b197a92d55d5aec4ff85645881..eb05314bae3f5ef1a4492076ab97cde93f096c30 100644
Binary files a/grad_conj/obj/Debug/main.o and b/grad_conj/obj/Debug/main.o differ