conjuguate gradient cleaned and commented

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>

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;
 
+    // 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,7 +90,7 @@ 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);
 
@@ -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);
 }