diff --git a/grad_conj/algebre.c b/grad_conj/algebre.c
index 6a1f08d7797182b5bf59d7a8158762508d7857eb..37371db236bb0de28468cb1f5d725d3f0a329077 100644
--- a/grad_conj/algebre.c
+++ b/grad_conj/algebre.c
@@ -5,22 +5,100 @@
#include "customError.h"
#include "algebre.h"
-Vector * prod_mat_vec(Matrix * matrix, Vector * vector){
+/* MATRIX x VECTOR */
+
+void prod_mat_vec(Vector * vector_res, Matrix * matrix, Vector * vector){
+ if(vector->t!=0)
+ exit(TRANSPOSED_VECTOR_USED);
if(matrix->width!=vector->size)
exit(UNCOMPATIBLE_ALGEBRIC_MULTIPLICATION);
- Vector * vector_res = allocateVector(matrix->height);
+ if(vector_res->size!=matrix->height)
+ exit(CONTAINER_RES_DO_NOT_FIT);
for(int i=0; i!=matrix->height; i++){
vector_res->c[i]=0;
for(int j=0; j!=matrix->width; j++){
vector_res->c[i]+=matrix->c[i][j]*vector->c[j];
}
}
- return vector_res;
+}
+
+/* VECTOR_t x MATRIX */
+
+void prod_vec_mat(Vector * vector_res, Vector * vector, Matrix * matrix){
+ if(vector->t!=1)
+ exit(NON_TRANSPOSED_VECTOR_USED);
+ if(matrix->height!=vector->size)
+ exit(UNCOMPATIBLE_ALGEBRIC_MULTIPLICATION);
+ if(vector_res->size!=matrix->width)
+ exit(CONTAINER_RES_DO_NOT_FIT);
+ vector_res->t=1;
+ for(int i=0; i!=matrix->width; i++){
+ vector_res->c[i]=0;
+ for(int j=0; j!=matrix->height; j++){
+ vector_res->c[i]+=matrix->c[j][i]*vector->c[j];
+ }
+ }
+}
+
+/* VECTOR + alpha x VECTOR */
+void sum_vec_alpha_vec(Vector * vector_res, Vector * vector1, float alpha, Vector * vector2){
+ if(vector_res->size!=vector1->size || vector1->size != vector2->size)
+ exit(UNCOMPATIBLE_ALGEBRIC_ADDITION);
+ for(int i=0;i!=vector1->size;i++)
+ vector_res->c[i]=vector1->c[i]+alpha*vector2->c[i];
+}
+
+/* VECTOR_t x MATRIX x VECTOR */
+
+float prod_vec_t_mat_vec(Vector * vector, Matrix * matrix){
+ vector->t=1;
+ Vector * subVec = allocateVector(matrix->width);
+ prod_vec_mat(subVec,vector,matrix);
+ vector->t=0;
+ float res = prod_vec_t_vec(subVec,vector);
+ freeVector(subVec);
+ return res;
+}
+
+/* NORM */
+
+float norm_vec_2(Vector * vector){
+ float res=0;
+ for(int i=0;i!=vector->size;i++)
+ res+=vector->c[i]*vector->c[i];
+ return res;
+}
+
+/* VECTOR_t x 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;
+ for(int i=0;i!=vector->size;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){
+ if(vector_res->size!=vector->size)
+ exit(CONTAINER_RES_DO_NOT_FIT);
+ for(int i=0; i!=vector->size; i++){
+ vector_res->c[i]=vector->c[i]*scal;
+ }
}
void testAlgebric(){
Vector * vector = allocateVector(8);
Vector * vector2 = allocateVector(8);
+ Vector * vector_res_5 = allocateVector(5);
+ Vector * vector_res_8 = allocateVector(8);
Matrix * matrix = allocateMatrix(5,8);
Matrix * matrix2 = allocateMatrix(8,8);
initProbVector(vector,10);
@@ -28,14 +106,31 @@ void testAlgebric(){
initProbVector(vector2,10);
initProbMatrix(matrix2,10);
- Vector * vector_res_prod_mat_vec = prod_mat_vec(matrix,vector);
+ prod_mat_vec(vector_res_5,matrix,vector);
+ printf("prod_mat_vec : \n");
+ printVector(vector_res_5);
+
+ vector->t=1;
+ prod_vec_mat(vector_res_8,vector,matrix2);
+ printf("prod_vec_mat : \n");
+ printVector(vector_res_8);
+
+ printf("prod_vec_t_vec : %f\n",prod_vec_t_vec(vector,vector2));
+ printf("prod_vec_t_mat_vec : %f\n",prod_vec_t_mat_vec(vector,matrix2));
+ printf("norm_vec_2 : %f\n",norm_vec_2(vector));
- printVector(vector_res_prod_mat_vec);
+ printf("sum_vec_alpha_vec : \n");
+ sum_vec_alpha_vec(vector,vector,-1,vector);
+ printVector(vector);
- freeVector(vector_res_prod_mat_vec);
+ printf("prod_vec_scal_u : \n");
+ prod_vec_scal(vector_res_8,vector_res_8,10.0);
+ printVector(vector_res_8);
freeVector(vector);
freeVector(vector2);
+ freeVector(vector_res_5);
+ freeVector(vector_res_8);
freeMatrix(matrix);
freeMatrix(matrix2);
}
diff --git a/grad_conj/algebre.h b/grad_conj/algebre.h
index 2482e958b08f501d8f52cedd8f16686dae648702..b1a61f8bec64460e4876257878e4749f4ee87abd 100644
--- a/grad_conj/algebre.h
+++ b/grad_conj/algebre.h
@@ -1,6 +1,17 @@
#ifndef ALGEBRE_H_INCLUDED
#define ALGEBRE_H_INCLUDED
-Vector * prod_mat_vec(Matrix * matrix, Vector * vector);
+void prod_mat_vec(Vector * vector_res, Matrix * matrix, Vector * vector);
+void prod_vec_mat(Vector * vector_res, Vector * vector, Matrix * matrix);
+
+void sum_vec_alpha_vec(Vector * vector_res, Vector * vector1, float alpha, Vector * vector2);
+
+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);
+
+void prod_vec_scal(Vector * vector_res, Vector * vector, float scal);
+
+void testAlgebric();
#endif // ALGEBRE_H_INCLUDED
diff --git a/grad_conj/bin/Debug/grad_conj.exe b/grad_conj/bin/Debug/grad_conj.exe
index 7be3e4cba3ea6c3c1c4ea50e36694bc24fde84b3..9bf973f8c18faa718cf3d1a71920c3af1a2e7010 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/customError.h b/grad_conj/customError.h
index 5047c438f77f22c9696cf559929191bb83d0cc08..825e95d036fa00e66e0e0df62e3a00e74c58c06d 100644
--- a/grad_conj/customError.h
+++ b/grad_conj/customError.h
@@ -4,5 +4,9 @@
#define ERROR_NO_MEMORY_LEFT 1000
#define UNCOMPATIBLE_ALGEBRIC_MULTIPLICATION 1001
#define NON_SQUARE_MATRIX 1002
+#define NON_TRANSPOSED_VECTOR_USED 1003
+#define TRANSPOSED_VECTOR_USED 1004
+#define CONTAINER_RES_DO_NOT_FIT 1005
+#define UNCOMPATIBLE_ALGEBRIC_ADDITION 1006
#endif // CUSTOMERROR_H_INCLUDED
diff --git a/grad_conj/grad_conj.depend b/grad_conj/grad_conj.depend
index 9f7b888c3a56dc53897d045e80aace4cbfd088d7..cce5b3e4564f2b643d4ae8c0206e138cb846a664 100644
--- a/grad_conj/grad_conj.depend
+++ b/grad_conj/grad_conj.depend
@@ -1,40 +1,41 @@
# depslib dependency file v1.0
-1489421217 source:b:\mes documents\progra\calcint\grad_conj\mathobject.c
+1489617942 source:b:\mes documents\progra\calcint\grad_conj\mathobject.c
<stdlib.h>
<stdio.h>
"mathObject.h"
"customError.h"
-1489420631 b:\mes documents\progra\calcint\grad_conj\mathobject.h
+1489612623 b:\mes documents\progra\calcint\grad_conj\mathobject.h
-1489421052 b:\mes documents\progra\calcint\grad_conj\customerror.h
+1489596942 b:\mes documents\progra\calcint\grad_conj\customerror.h
-1489420337 source:b:\mes documents\progra\calcint\grad_conj\main.c
+1489619150 source:b:\mes documents\progra\calcint\grad_conj\main.c
<stdio.h>
<stdlib.h>
<time.h>
<mpi.h>
<time.h>
+ <math.h>
"customError.h"
"mathObject.h"
"algebre.h"
"iterative_methods.h"
-1489419800 source:b:\mes documents\progra\calcint\grad_conj\algebre.c
+1489617314 source:b:\mes documents\progra\calcint\grad_conj\algebre.c
<stdlib.h>
<stdio.h>
"mathObject.h"
"customError.h"
"algebre.h"
-1489419451 b:\mes documents\progra\calcint\grad_conj\algebre.h
+1489617252 b:\mes documents\progra\calcint\grad_conj\algebre.h
-1489420232 source:b:\mes documents\progra\calcint\grad_conj\iterative_methods.c
+1489619088 source:b:\mes documents\progra\calcint\grad_conj\iterative_methods.c
<stdlib.h>
<stdio.h>
"mathObject.h"
"customError.h"
"algebre.h"
-1489420308 b:\mes documents\progra\calcint\grad_conj\iterative_methods.h
+1489612728 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 9666c0cf7f85aef71677fe907aca03f400051f37..a2a62ae302f404157c09766df6cc108652c361bf 100644
--- a/grad_conj/iterative_methods.c
+++ b/grad_conj/iterative_methods.c
@@ -5,3 +5,55 @@
#include "customError.h"
#include "algebre.h"
+#define MAX_ITERATION 1000000
+
+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 4e6cc232e6a0fbe28e9191ab70d7c336bbe311f9..3f3723366b19929159c955c4d717af898c5aaa9c 100644
--- a/grad_conj/iterative_methods.h
+++ b/grad_conj/iterative_methods.h
@@ -1,6 +1,6 @@
#ifndef ITERATIVE_METHODS_H_INCLUDED
#define ITERATIVE_METHODS_H_INCLUDED
-
+void grad_conj(Vector * b, Matrix * A, Vector * x, float epsilon);
#endif // ITERATIVE_METHODS_H_INCLUDED
diff --git a/grad_conj/main.c b/grad_conj/main.c
index badc9601dfa0f6dbf07e6b4575e7081cb016c8b3..d291b55070c3c2908c84a734f3250f19810284eb 100644
--- a/grad_conj/main.c
+++ b/grad_conj/main.c
@@ -3,6 +3,7 @@
#include <time.h>
#include <mpi.h>
#include <time.h>
+#include <math.h>
#include "customError.h"
#include "mathObject.h"
@@ -20,4 +21,27 @@ int main (int argc, char *argv[]) {
testVector();
testAlgebric();
+
+ Vector * b = allocateVector(10000);
+ Matrix * A = allocateMatrix(10000,10000);
+ Vector * x = allocateVector(10000);
+
+ initProbVector(b,5.0);
+ fillVector(x,0);
+ initProbMatrixSymetric(A,5.0);
+ for(int i=0;i!=A->width;i++){
+ A->c[i][i]+=A->width/15;
+ }
+ /*
+ 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);
}
diff --git a/grad_conj/mathObject.c b/grad_conj/mathObject.c
index 246086c36bf963e42ca1759750e6959790a8e98b..5ff139ae01b5dd17856487f0c09cca03071c5489 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.1f, ",matrix->c[i][j]);
+ printf("%5.2f, ",matrix->c[i][j]);
}
printf("\n");
}
@@ -94,6 +94,7 @@ void testMatrix(){
Vector * allocateVector(int size){
Vector * vector = malloc(sizeof(Vector));
vector->size = size;
+ vector->t = 0;
vector->c = (float *)malloc(size * sizeof(float));
if(vector->c == NULL)
exit(ERROR_NO_MEMORY_LEFT);
@@ -104,6 +105,17 @@ void initProbVector(Vector * vector, int max_value){
vector->c[i]=((float)rand())/32768.0*(float)max_value;
}
}
+void fillVector(Vector * vector, float val){
+ for(int i = 0;i!=vector->size;i++)
+ vector->c[i]=val;
+}
+void copyVector(Vector * vec_original, Vector * vec_copy){
+ if(vec_original->size!=vec_copy->size)
+ exit(CONTAINER_RES_DO_NOT_FIT);
+ for(int i=0;i!=vec_original->size;i++){
+ vec_copy->c[i]=vec_original->c[i];
+ }
+}
void freeVector(Vector * vector){
free(vector->c);
free(vector);
@@ -113,8 +125,10 @@ void freeVector(Vector * vector){
void printVector(Vector * vector){
for(int i=0; i<vector->size; i++){
- printf("%5.1f \n",vector->c[i]);
+ printf("%5.2f \n",vector->c[i]);
}
+ if(vector->t)
+ printf("is transposed\n");
}
/* TEST BASIC FUNCTIONALITY */
@@ -126,6 +140,8 @@ void testVector(){
printf("OK\nInitialisation : ");
initProbVector(vector,10);
printf("OK\n");
+ copyVector(vector,vector);
+ printf("Copy : OK\n");
printVector(vector);
printf("Printing : OK\nFree : ");
freeVector(vector);
diff --git a/grad_conj/mathObject.h b/grad_conj/mathObject.h
index 7a17847d2b3ea19de55a0af992bb8367c6511fc0..2799838e9cf881a2ce9f48624cec2a230f982e68 100644
--- a/grad_conj/mathObject.h
+++ b/grad_conj/mathObject.h
@@ -5,6 +5,7 @@ typedef struct vector{
int size;
float * c;
+ int t;
}Vector;
@@ -25,6 +26,8 @@ void testMatrix();
Vector * allocateVector(int);
void initProbVector(Vector *, int);
+void fillVector(Vector * vector, float val);
+void copyVector(Vector *, Vector *);
void freeVector(Vector *);
void printVector(Vector *);
void testVector();
diff --git a/grad_conj/obj/Debug/algebre.o b/grad_conj/obj/Debug/algebre.o
index d7a488fe512475e6c801609a47c3cc1e17b198fb..a15ad74ba1e829337ddc7b94d2ba14f5501ddca8 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 797a45098251cbfc9931813348776a3bae0bd086..519f3eb08eb188be76b7fd5afa2ba6cef9faed4b 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 4e88e17f0dd7d1fed42d8b72bd526874f303500f..c5611a3ed116db893bfc856ffcd9d3680ac00b64 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 65e0dc38b65d1f77b1afa3f54ffce01b019a3a9a..2b419c5fb1c3ebcf4acee0206d84d84d17ddf09d 100644
Binary files a/grad_conj/obj/Debug/mathObject.o and b/grad_conj/obj/Debug/mathObject.o differ