simponic/lizfcm
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

lu factorization

Browse Source
This commit is contained in:
Elizabeth Hunt 2023-10-11 15:56:20 -06:00
parent faaef032d8
commit 78bd8c4a95
Signed by: simponic
GPG Key ID: 52B3774857EB24B1
5 changed files with 153 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
inc/lizfcm.h
View File

@ -12,22 +12,26 @@ extern double forward_derivative_at(double (*f)(double), double a, double h);
extern double backward_derivative_at(double (*f)(double), double a, double h);
extern double sum_v(Array_double *v);
extern Array_double *scale_v(Array_double *v1, double m);
extern Array_double *add_v(Array_double *v1, Array_double *v2);
extern Array_double *minus_v(Array_double *v1, Array_double *v2);
extern double dot_v(Array_double *v1, Array_double *v2);
extern double l2_norm(Array_double *v);
extern double l1_norm(Array_double *v);
extern double linf_norm(Array_double *v);
extern double l2_distance(Array_double *v1, Array_double *v2);
extern double l1_distance(Array_double *v1, Array_double *v2);
extern double linf_distance(Array_double *v1, Array_double *v2);
extern Array_double *copy_vector(Array_double *v1);
extern void free_vector(Array_double *v);
extern void format_vector_into(Array_double *v, char *s);
extern void format_matrix_into(Matrix_double *m, char *s);
extern void put_identity_diagonal(Matrix_double *m);
extern Line *least_squares_lin_reg(Array_double *x, Array_double *y);
extern void put_identity_diagonal(Matrix_double *m);
extern Matrix_double *copy_matrix(Matrix_double *m);
extern void free_matrix(Matrix_double *m);
extern Matrix_double **put_lu_decomp(Matrix_double *m);
extern void format_matrix_into(Matrix_double *m, char *s);
#endif // LIZFCM_H

BIN
lib/lizfcm.a
View File

Binary file not shown.

62
src/matrix.c
View File

@ -9,14 +9,74 @@ void put_identity_diagonal(Matrix_double *m) {
m->data[y]->data[y] = 1.0;
}
Matrix_double *copy_matrix(Matrix_double *m) {
Matrix_double *copy = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
for (size_t y = 0; y < copy->rows; y++) {
free_vector(copy->data[y]);
copy->data[y] = copy_vector(m->data[y]);
}
return copy;
}
Matrix_double **put_lu_decomp(Matrix_double *m) {
assert(m->cols == m->rows);
Matrix_double *u = copy_matrix(m);
Matrix_double *l = InitMatrixWithSize(double, m->rows, m->cols, 0.0);
put_identity_diagonal(l);
Matrix_double **u_l = malloc(sizeof(Matrix_double *) * 2);
for (size_t y = 0; y < m->rows; y++) {
if (u->data[y]->data[y] == 0) {
printf("ERROR: a pivot is zero in given matrix\n");
exit(-1);
}
}
if (u && l) {
for (size_t x = 0; x < m->cols; x++) {
for (size_t y = x + 1; y < m->rows; y++) {
double denom = u->data[x]->data[x];
if (denom == 0) {
printf("ERROR: non-factorable matrix\n");
exit(-1);
}
double factor = -(u->data[y]->data[x] / denom);
Array_double *scaled = scale_v(u->data[x], factor);
Array_double *added = add_v(scaled, u->data[y]);
free_vector(scaled);
free_vector(u->data[y]);
u->data[y] = added;
l->data[y]->data[x] = -factor;
}
}
}
u_l[0] = u;
u_l[1] = l;
return u_l;
}
void free_matrix(Matrix_double *m) {
for (size_t y = 0; y < m->rows; ++y)
free_vector(m->data[y]);
free(m);
}
void format_matrix_into(Matrix_double *m, char *s) {
sprintf(s, "");
if (m->rows == 0)
sprintf(s, "empty");
for (size_t y = 0; y < m->rows; ++y) {
char row_s[256];
char *row_s = malloc(sizeof(char) * 256);
format_vector_into(m->data[y], row_s);
sprintf(s, "%s %s \n", s, row_s);
free(row_s);
}
}

23
src/vector.c
View File

@ -41,12 +41,19 @@ double sum_v(Array_double *v) {
return sum;
}
Array_double *scale_v(Array_double *v, double m) {
Array_double *copy = copy_vector(v);
for (size_t i = 0; i < v->size; i++)
copy->data[i] *= m;
return copy;
}
Array_double *add_v(Array_double *v1, Array_double *v2) {
assert(v1->size == v2->size);
Array_double *sum = InitArrayWithSize(double, v1->size, 0);
Array_double *sum = copy_vector(v1);
for (size_t i = 0; i < v1->size; i++)
sum->data[i] = v1->data[i] + v2->data[i];
sum->data[i] += v2->data[i];
return sum;
}
@ -80,6 +87,18 @@ double linf_distance(Array_double *v1, Array_double *v2) {
return dist;
}
Array_double *copy_vector(Array_double *v) {
Array_double *copy = InitArrayWithSize(double, v->size, 0.0);
for (size_t i = 0; i < copy->size; ++i)
copy->data[i] = v->data[i];
return copy;
}
void free_vector(Array_double *v) {
free(v->data);
free(v);
}
void format_vector_into(Array_double *v, char *s) {
sprintf(s, "");
if (v->size == 0)

106
test/main.c
View File

@ -5,56 +5,72 @@
double f(double x) { return (x - 1) / (x + 1); }
int main() {
printf("smaceps(): %.10e\n", smaceps());
printf("dmaceps(): %.10e\n", dmaceps());
Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3});
Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5});
char v_s[256];
char w_s[256];
format_vector_into(v, v_s);
format_vector_into(w, w_s);
printf("v: %s\n", v_s);
printf("w: %s\n", w_s);
printf("l1_norm(v): %f\n", l1_norm(v));
printf("l2_norm(v): %f\n", l2_norm(v));
printf("linf_norm(v): %f\n", linf_norm(v));
printf("l1_dist(v, w): %f\n", l1_distance(v, w));
printf("l2_dist(v, w): %f\n", l2_distance(v, w));
printf("linf_dist(v, w): %f\n", linf_distance(v, w));
double h = 0.001;
printf("approx f'(1) w/ c.d.: %f\n", central_derivative_at(&f, 1, h));
printf("approx f'(1) w/ fw.d.: %f\n", forward_derivative_at(&f, 1, h));
printf("approx f'(1) w/ bw.d.: %f\n", backward_derivative_at(&f, 1, h));
v = InitArray(double, {1, 2, 3, 4, 5});
w = InitArray(double, {2, 3, 4, 5, 6});
format_vector_into(v, v_s);
format_vector_into(w, w_s);
printf("v: %s\n", v_s);
printf("w: %s\n", w_s);
Line *line = least_squares_lin_reg(v, w);
printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
v = InitArray(double, {1, 2, 3, 4, 5, 6, 7});
w = InitArray(double, {0.5, 3, 2, 3.5, 5, 6, 7.5});
format_vector_into(v, v_s);
format_vector_into(w, w_s);
printf("v: %s\n", v_s);
printf("w: %s\n", w_s);
line = least_squares_lin_reg(v, w);
printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
// printf("smaceps(): %.10e\n", smaceps());
// printf("dmaceps(): %.10e\n", dmaceps());
//
// Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3});
// Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5});
//
// char v_s[256];
// char w_s[256];
// format_vector_into(v, v_s);
// format_vector_into(w, w_s);
//
// printf("v: %s\n", v_s);
// printf("w: %s\n", w_s);
// printf("l1_norm(v): %f\n", l1_norm(v));
// printf("l2_norm(v): %f\n", l2_norm(v));
// printf("linf_norm(v): %f\n", linf_norm(v));
//
// printf("l1_dist(v, w): %f\n", l1_distance(v, w));
// printf("l2_dist(v, w): %f\n", l2_distance(v, w));
// printf("linf_dist(v, w): %f\n", linf_distance(v, w));
//
// double h = 0.001;
// printf("approx f'(1) w/ c.d.: %f\n", central_derivative_at(&f, 1, h));
// printf("approx f'(1) w/ fw.d.: %f\n", forward_derivative_at(&f, 1, h));
// printf("approx f'(1) w/ bw.d.: %f\n", backward_derivative_at(&f, 1, h));
//
// v = InitArray(double, {1, 2, 3, 4, 5});
// w = InitArray(double, {2, 3, 4, 5, 6});
// format_vector_into(v, v_s);
// format_vector_into(w, w_s);
// printf("v: %s\n", v_s);
// printf("w: %s\n", w_s);
//
// Line *line = least_squares_lin_reg(v, w);
// printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
//
// v = InitArray(double, {1, 2, 3, 4, 5, 6, 7});
// w = InitArray(double, {0.5, 3, 2, 3.5, 5, 6, 7.5});
//
// format_vector_into(v, v_s);
// format_vector_into(w, w_s);
// printf("v: %s\n", v_s);
// printf("w: %s\n", w_s);
// line = least_squares_lin_reg(v, w);
// printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
char m_s[256];
Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
put_identity_diagonal(m);
for (int i = 0; i < 8; i++) {
for (int j = 0; j < 8; j++) {
m->data[i]->data[j] = (i + 1.0) + j * 3 + (rand() % 12);
}
}
format_matrix_into(m, m_s);
printf("%s", m_s);
printf("%s\n", m_s);
Matrix_double **u_l = put_lu_decomp(m);
Matrix_double *u = u_l[0];
Matrix_double *l = u_l[1];
format_matrix_into(u, m_s);
printf("%s\n", m_s);
format_matrix_into(l, m_s);
printf("%s\n", m_s);
return 0;
}