summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--inc/lizfcm.h4
-rw-r--r--lib/lizfcm.abin14128 -> 15472 bytes
-rw-r--r--notes/Oct-13.org8
-rw-r--r--src/matrix.c30
-rw-r--r--test/main.c100
5 files changed, 102 insertions, 40 deletions
diff --git a/inc/lizfcm.h b/inc/lizfcm.h
index 99acbf2..40046d2 100644
--- a/inc/lizfcm.h
+++ b/inc/lizfcm.h
@@ -30,7 +30,9 @@ 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 Matrix_double **lu_decomp(Matrix_double *m);
+extern Array_double *bsubst(Matrix_double *u, Array_double *b);
+extern Array_double *fsubst(Matrix_double *l, Array_double *b);
extern void format_matrix_into(Matrix_double *m, char *s);
diff --git a/lib/lizfcm.a b/lib/lizfcm.a
index 78ca79e..7ea59ed 100644
--- a/lib/lizfcm.a
+++ b/lib/lizfcm.a
Binary files differ
diff --git a/notes/Oct-13.org b/notes/Oct-13.org
new file mode 100644
index 0000000..853a6d6
--- /dev/null
+++ b/notes/Oct-13.org
@@ -0,0 +1,8 @@
+* Root Finding
+Finx x \in R such that f(x) = 0
+
+If g(x) is a function and we want to find x such that g(x) is
+extremal then we find x \in R such that g'(x) = 0
+
+
+Hanging Cable Problem y(x) = c_1 cosh(\frac{x- c_2}{c_2})
diff --git a/src/matrix.c b/src/matrix.c
index 953d50f..438a468 100644
--- a/src/matrix.c
+++ b/src/matrix.c
@@ -1,7 +1,7 @@
#include "lizfcm.h"
#include <assert.h>
-#include <string.h>
#include <stdio.h>
+#include <string.h>
void put_identity_diagonal(Matrix_double *m) {
assert(m->rows == m->cols);
@@ -19,7 +19,7 @@ Matrix_double *copy_matrix(Matrix_double *m) {
return copy;
}
-Matrix_double **put_lu_decomp(Matrix_double *m) {
+Matrix_double **lu_decomp(Matrix_double *m) {
assert(m->cols == m->rows);
Matrix_double *u = copy_matrix(m);
@@ -63,6 +63,32 @@ Matrix_double **put_lu_decomp(Matrix_double *m) {
return u_l;
}
+Array_double *bsubst(Matrix_double *u, Array_double *b) {
+ assert(u->rows == b->size && u->cols == u->rows);
+
+ Array_double *x = copy_vector(b);
+ for (int64_t row = b->size - 1; row >= 0; row--) {
+ for (size_t col = b->size - 1; col > row; col--)
+ x->data[row] -= x->data[col] * u->data[row]->data[col];
+ x->data[row] /= u->data[row]->data[row];
+ }
+ return x;
+}
+
+Array_double *fsubst(Matrix_double *l, Array_double *b) {
+ assert(l->rows == b->size && l->cols == l->rows);
+
+ Array_double *x = copy_vector(b);
+
+ for (size_t row = 0; row < b->size; row++) {
+ for (size_t col = 0; col < row; col++)
+ x->data[row] -= x->data[col] * l->data[row]->data[col];
+ x->data[row] /= l->data[row]->data[row];
+ }
+
+ return x;
+}
+
void free_matrix(Matrix_double *m) {
for (size_t y = 0; y < m->rows; ++y)
free_vector(m->data[y]);
diff --git a/test/main.c b/test/main.c
index 75e60ee..2364b3d 100644
--- a/test/main.c
+++ b/test/main.c
@@ -10,19 +10,18 @@ int main() {
printf("dmaceps(): %.10e\n", dmaceps());
printf("========\n");
+ printf("Norm, Distance\n");
Array_double *v = InitArray(double, {3, 1, -4, 1, 5, -9, 3});
- char v_s[256];
- strcpy(v_s, "");
- format_vector_into(v, v_s);
+ char s[2048];
+ strcpy(s, "");
+ format_vector_into(v, s);
+ printf("v: %s", s);
Array_double *w = InitArray(double, {-2, 7, 1, -8, -2, 8, 5});
- char w_s[256];
- strcpy(w_s, "");
- format_vector_into(w, w_s);
+ strcpy(s, "");
+ format_vector_into(w, s);
+ printf("w: %s", s);
- printf("Norm, Distance\n");
- printf("v: %s", v_s);
- printf("w: %s", 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));
@@ -38,52 +37,79 @@ int main() {
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));
printf("========\n");
+ printf("Least Squares\n");
v = InitArray(double, {1, 2, 3, 4, 5});
- strcpy(v_s, "");
- format_vector_into(v, v_s);
+ strcpy(s, "");
+ format_vector_into(v, s);
+ printf("v: %s", s);
w = InitArray(double, {2, 3, 4, 5, 6});
- strcpy(w_s, "");
- format_vector_into(w, w_s);
+ strcpy(s, "");
+ format_vector_into(w, s);
+ printf("w: %s", s);
+
Line *line = least_squares_lin_reg(v, w);
- printf("Least Squares\n");
- printf("v: %s", v_s);
- printf("w: %s", w_s);
+
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});
- strcpy(v_s, "");
- format_vector_into(v, v_s);
+ strcpy(s, "");
+ format_vector_into(v, s);
+ printf("v: %s", s);
w = InitArray(double, {0.5, 3, 2, 3.5, 5, 6, 7.5});
- strcpy(w_s, "");
- format_vector_into(w, w_s);
- printf("v: %s", v_s);
- printf("w: %s", w_s);
+ strcpy(s, "");
+ format_vector_into(w, s);
+ printf("w: %s", s);
+
line = least_squares_lin_reg(v, w);
printf("least_squares_lin_reg(v, w): (%f)x + %f\n", line->m, line->a);
printf("========\n");
printf("LU Decomp\n");
- char m_s[2048];
- Matrix_double *m = InitMatrixWithSize(double, 8, 8, 0.0);
- 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);
+ Matrix_double *m = InitMatrixWithSize(double, 10, 10, 0.0);
+ for (int i = 0; i < 10; i++) {
+ for (int j = 0; j < 10; j++)
+ m->data[i]->data[j] = (100 - rand() % 200);
}
- format_matrix_into(m, m_s);
- printf("m = %s", m_s);
+ strcpy(s, "");
+ format_matrix_into(m, s);
+ printf("m = %s", s);
- Matrix_double **u_l = put_lu_decomp(m);
+ Array_double *b = InitArrayWithSize(double, 10, 100.0);
+ Matrix_double **u_l = lu_decomp(m);
Matrix_double *u = u_l[0];
Matrix_double *l = u_l[1];
- strcpy(m_s, "");
- format_matrix_into(u, m_s);
- printf("u = %s", m_s);
- strcpy(m_s, "");
- format_matrix_into(l, m_s);
- printf("l = %s", m_s);
+ strcpy(s, "");
+ format_matrix_into(u, s);
+ printf("u = %s", s);
+ strcpy(s, "");
+ format_matrix_into(l, s);
+ printf("l = %s", s);
+ strcpy(s, "");
+ format_vector_into(b, s);
+ printf("b = %s", s);
printf("========\n");
- printf("Back Substitution\n");
+ printf("Backward -> Forward Substitution\n");
+
+ Array_double *b_fsub = fsubst(l, b);
+ strcpy(s, "");
+ format_vector_into(b_fsub, s);
+ printf("x: %s\n", s);
+
+ Array_double *x_bsub = bsubst(u, b_fsub);
+ strcpy(s, "");
+ format_vector_into(x_bsub, s);
+ printf("x: %s", s);
+
+ free_vector(b_fsub);
+
+ printf("Verifications (each should be approximately 100)\n");
+ for (size_t row = 0; row < m->rows; row++) {
+ double curr = 0;
+ for (size_t col = 0; col < m->cols; col++)
+ curr += m->data[row]->data[col] * x_bsub->data[col];
+ printf("Substitution for row %zu = %f\n", row, curr);
+ }
return 0;
}
generated by cgit v1.2.3-70-g09d2 (git 2.46.0) at 2025-09-21 19:24:17 +0000