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

convergence of iterative gradient conjugate

parent b220072c
Branches
No related tags found
No related merge requests found
......@@ -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);
}
#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
No preview for this file type
......@@ -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
# 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
......@@ -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);
}
#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
......@@ -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);
}
......@@ -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);
......
......@@ -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();
......
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