Skip to content
Snippets Groups Projects
Commit e484dc44 authored by Lucas Christesen Ahler's avatar Lucas Christesen Ahler
Browse files

Debuggin linear equations

parent 11da16af
No related branches found
No related tags found
No related merge requests found
CFLAGS = -Wall -O1 -std=gnu11 $(shell gsl-config --cflags) CFLAGS = -Wall -O3 -std=gnu11 $(shell gsl-config --cflags)
LDLIBS = -lm $(shell gsl-config --libs) LDLIBS = -lm $(shell gsl-config --libs)
default: out.txt time_graph.png default: out.txt time_graph.png
...@@ -6,6 +6,8 @@ default: out.txt time_graph.png ...@@ -6,6 +6,8 @@ default: out.txt time_graph.png
out.txt time_test.txt: main out.txt time_test.txt: main
./$< > $@ ./$< > $@
main: main.c functions.c
time_graph.png: time_test.txt Makefile time_graph.png: time_test.txt Makefile
echo '\ echo '\
set terminal png;\ set terminal png;\
......
#include "functions.h"
void print_vector(char s[], gsl_vector* v) {
printf("%s\n",s);
for(int i=0;i<v->size;i++) {
printf("%15g",gsl_vector_get(v,i));
printf("\n");
}
}
void print_matrix(gsl_matrix* M) {
for(int i = 0; i < M->size1; i++) {
for(int j = 0; j < M->size2; j++) {
printf("%15g\t", gsl_matrix_get(M, i, j));
}
printf("\n");
}
}
void GS_decomp(gsl_matrix* A, gsl_matrix* R) {
int n = A -> size1;
int m = A -> size2;
for(int i=0;i<m;i++) {
gsl_vector_view ai = gsl_matrix_column(A,i);
double Rii = gsl_blas_dnrm2(&(ai.vector));
gsl_matrix_set(R,i,i,Rii);
gsl_vector_scale(&(ai.vector),1/Rii);
for(int j=i+1;j<m;j++) {
gsl_vector_view aj = gsl_matrix_column(A,j);
double Rij;
gsl_blas_ddot(&(ai.vector),&(aj.vector),&Rij);
gsl_matrix_set(R,i,j,Rij);
gsl_vector* qiRij = gsl_vector_alloc(n);
gsl_vector_memcpy(qiRij,&(ai.vector));
gsl_vector_scale(qiRij,Rij);
gsl_vector_sub(&(aj.vector),qiRij);
}
}
}
void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x) {
assert(Q->size1==Q->size2);
int m = Q -> size2;
gsl_blas_dgemv(CblasTrans,1,Q,b,0,x);
for(int i=m-1;i>=0;i--) {
double xi = gsl_vector_get(x,i);
for(int j=i+1;j<m;j++) {
xi += -gsl_matrix_get(R,i,j)*gsl_vector_get(x,j);
}
xi /= gsl_matrix_get(R,i,i);
gsl_vector_set(x,i,xi);
}
}
void GS_inverse(gsl_matrix* Q, gsl_matrix* R, gsl_matrix* B) {
for(int j=0;j<Q->size2;j++) {
gsl_vector* ej = gsl_vector_alloc(Q->size1);
gsl_vector_set_basis(ej,j);
gsl_vector_view Bj = gsl_matrix_column(B,j);
GS_solve(Q,R,ej,&(Bj.vector));
gsl_vector_free(ej);
}
}
\ No newline at end of file
#ifndef FUNCTIONS_H
#define FUNCTIONS_H
#include<math.h>
#include<time.h>
#include<stdio.h>
#include<assert.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_linalg.h>
#define RND (double)rand()/RAND_MAX
void print_vector(char s[], gsl_vector* v);
void print_matrix(gsl_matrix* M);
void GS_decomp(gsl_matrix* A, gsl_matrix* R);
void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x);
void GS_inverse(gsl_matrix* Q, gsl_matrix* R, gsl_matrix* B);
#endif
\ No newline at end of file
#include<math.h> #include "functions.h"
#include<time.h>
#include<stdio.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_linalg.h>
#define RND (double)rand()/RAND_MAX
void print_vector(char s[], gsl_vector* v) {
printf("%s\n",s);
for(int i=0;i<v->size;i++) {
printf("%15g",gsl_vector_get(v,i));
printf("\n");
}
}
void print_matrix(gsl_matrix* M) {
for(int i = 0; i < M->size1; i++) {
for(int j = 0; j < M->size2; j++) {
printf("%15g\t", gsl_matrix_get(M, i, j));
}
printf("\n");
}
}
void GS_decomp(gsl_matrix* A, gsl_matrix* R) {
int n = A -> size1;
int m = A -> size2;
for(int i=0;i<m;i++) {
gsl_vector_view ai = gsl_matrix_column(A,i);
double Rii = gsl_blas_dnrm2(&(ai.vector));
gsl_matrix_set(R,i,i,Rii);
gsl_vector_scale(&(ai.vector),1/Rii);
for(int j=i+1;j<m;j++) {
gsl_vector_view aj = gsl_matrix_column(A,j);
double Rij;
gsl_blas_ddot(&(ai.vector),&(aj.vector),&Rij);
gsl_matrix_set(R,i,j,Rij);
gsl_vector* qiRij = gsl_vector_alloc(n);
gsl_vector_memcpy(qiRij,&(ai.vector));
gsl_vector_scale(qiRij,Rij);
gsl_vector_sub(&(aj.vector),qiRij);
}
}
}
void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x) {
int m = Q -> size2;
gsl_blas_dgemv(CblasTrans,1,Q,b,0,x);
for(int i=m-1;i>=0;i--) {
double xi = gsl_vector_get(x,i);
for(int j=i+1;j<m;j++) {
xi += -gsl_matrix_get(R,i,j)*gsl_vector_get(x,j);
}
xi /= gsl_matrix_get(R,i,i);
gsl_vector_set(x,i,xi);
}
}
void GS_inverse(gsl_matrix* Q, gsl_matrix* R, gsl_matrix* B) {
for(int j=0;j<Q->size2;j++) {
gsl_vector* ej = gsl_vector_alloc(Q->size1);
gsl_vector_set_basis(ej,j);
gsl_vector_view Bj = gsl_matrix_column(B,j);
GS_solve(Q,R,ej,&(Bj.vector));
gsl_vector_free(ej);
}
}
int main() { int main() {
int n =3; int n =3;
...@@ -139,42 +71,35 @@ int main() { ...@@ -139,42 +71,35 @@ int main() {
print_matrix(AB); print_matrix(AB);
FILE* time_test_out_stream = fopen("time_test.txt","w"); FILE* time_test_out_stream = fopen("time_test.txt","w");
for(int N=2;N<100;N++) { for(int N=300;N<500;N++) {
int iterations = 200;
gsl_matrix* My_time_test_A = gsl_matrix_alloc(N,N); gsl_matrix* My_time_test_A = gsl_matrix_alloc(N,N);
gsl_matrix* Gsl_time_test_A = gsl_matrix_alloc(N,N);
gsl_matrix* My_time_test_R = gsl_matrix_alloc(N,N); gsl_matrix* My_time_test_R = gsl_matrix_alloc(N,N);
gsl_vector* Gsl_time_test_v = gsl_vector_alloc(N); gsl_vector* Gsl_time_test_v = gsl_vector_alloc(N);
for(int k=0; k<My_time_test_A->size1;k++) {
for(int j=0; j<My_time_test_A->size2;j++) {
double kj=RND;
gsl_matrix_set(My_time_test_A,k,j,kj);
}
}
gsl_matrix_memcpy(Gsl_time_test_A,My_time_test_A);
clock_t my_time1 = clock(); clock_t my_time1 = clock();
for(int i=0;i<iterations;i++) { GS_decomp(My_time_test_A,My_time_test_R);
for(int k=0; k<My_time_test_A->size1;k++) {
for(int j=0; j<My_time_test_A->size2;j++) {
double kj=RND;
gsl_matrix_set(My_time_test_A,k,j,kj);
}
}
GS_decomp(My_time_test_A,My_time_test_R);
}
clock_t my_time2 = clock(); clock_t my_time2 = clock();
double my_time = (double)(my_time2-my_time1); double my_time = (double)(my_time2-my_time1);
clock_t gsl_time1 = clock(); clock_t gsl_time1 = clock();
for(int i=0;i<iterations;i++) { gsl_linalg_QR_decomp(Gsl_time_test_A,Gsl_time_test_v);
for(int k=0; k<My_time_test_A->size1;k++) {
for(int j=0; j<My_time_test_A->size2;j++) {
double kj=RND;
gsl_matrix_set(My_time_test_A,k,j,kj);
}
}
gsl_linalg_QR_decomp(My_time_test_A,Gsl_time_test_v);
}
clock_t gsl_time2 = clock(); clock_t gsl_time2 = clock();
double gsl_time = (double)(gsl_time2-gsl_time1); double gsl_time = (double)(gsl_time2-gsl_time1);
my_time /= iterations; my_time /= CLOCKS_PER_SEC; my_time /= CLOCKS_PER_SEC;
gsl_time /= iterations; gsl_time /= CLOCKS_PER_SEC; gsl_time /= CLOCKS_PER_SEC;
fprintf(time_test_out_stream,"%i %g %g\n",N,my_time/iterations,gsl_time/iterations); fprintf(time_test_out_stream,"%i %g %g\n",N,my_time,gsl_time);
gsl_matrix_free(My_time_test_A); gsl_matrix_free(My_time_test_A);
gsl_matrix_free(Gsl_time_test_A);
gsl_matrix_free(My_time_test_R); gsl_matrix_free(My_time_test_R);
gsl_vector_free(Gsl_time_test_v); gsl_vector_free(Gsl_time_test_v);
} }
......
...@@ -55,5 +55,4 @@ A*B: ...@@ -55,5 +55,4 @@ A*B:
-1.73472e-16 -5.41234e-16 1 -1.73472e-16 -5.41234e-16 1
B*A: B*A:
1 -3.33067e-16 1.38778e-17 1 -3.33067e-16 1.38778e-17
3.88578e-16 1 -6.76542e-17 3.88578e-16
1.7486e-15 -8.60423e-16 1 \ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment