Skip to content
Snippets Groups Projects
Commit 7d140d87 authored by Nicolas Fley's avatar Nicolas Fley
Browse files

compiling and converging parallel conjuguate gradient done

parent 7cdce459
Branches
No related tags found
No related merge requests found
......@@ -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){
......
......@@ -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);
......
No preview for this file type
# 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
#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);
......
#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);
}
......@@ -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
......@@ -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];
}
int nb_line_proc_in_use = nb_line_for_proc[rank];
testAlgebric();
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);
Vector * b = allocateVector(10000);
Matrix * A = allocateMatrix(10000,10000);
Vector * x = allocateVector(10000);
/** resolution of the conjugate gradient **/
initProbVector(b,5.0);
/* 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++;
}
/*
printf("\n\n########\nproblem : b-Ax=0 \n");
printf("A : \n");
printMatrix(A);
/*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");
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);*/
printVector(b_tot);
printf("x : \n");
printVector(x_tot);
}
printf("calculation : \n");
grad_conj(b,A,x,0.000001);
MPI_Finalize();
printf("grad conj result : \n");
//printVector(x);
if(rank==0)
freeMatrix(A_tot);
}
......@@ -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);
}
......@@ -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
No preview for this file type
No preview for this file type
No preview for this file type
No preview for this file type
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment