diff --git a/inc/lizfcm.h b/inc/lizfcm.h index 2dc9e42..99acbf2 100644 --- a/inc/lizfcm.h +++ b/inc/lizfcm.h @@ -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 diff --git a/lib/lizfcm.a b/lib/lizfcm.a index 7a21e87..969c894 100644 Binary files a/lib/lizfcm.a and b/lib/lizfcm.a differ diff --git a/src/matrix.c b/src/matrix.c index 115e0a3..a23aad4 100644 --- a/src/matrix.c +++ b/src/matrix.c @@ -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); } } diff --git a/src/vector.c b/src/vector.c index 61692e1..e6e2544 100644 --- a/src/vector.c +++ b/src/vector.c @@ -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) diff --git a/test/main.c b/test/main.c index 798bfbf..3551af2 100644 --- a/test/main.c +++ b/test/main.c @@ -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; }