From 6bf4b90c90f15f4ab60833bddf5b5756d1a6b1f6 Mon Sep 17 00:00:00 2001 From: Elizabeth Alexander Hunt Date: Thu, 2 Jul 2026 11:55:17 -0700 Subject: Init --- Homework/math4610/doc/software_manual.org | 1376 +++++++++++++++++++++++++ Homework/math4610/doc/software_manual.pdf | Bin 0 -> 227425 bytes Homework/math4610/doc/software_manual.tex | 1583 +++++++++++++++++++++++++++++ 3 files changed, 2959 insertions(+) create mode 100644 Homework/math4610/doc/software_manual.org create mode 100644 Homework/math4610/doc/software_manual.pdf create mode 100644 Homework/math4610/doc/software_manual.tex (limited to 'Homework/math4610/doc') diff --git a/Homework/math4610/doc/software_manual.org b/Homework/math4610/doc/software_manual.org new file mode 100644 index 0000000..b41ec03 --- /dev/null +++ b/Homework/math4610/doc/software_manual.org @@ -0,0 +1,1376 @@ +#+TITLE: LIZFCM Software Manual (v0.6) +#+AUTHOR: Elizabeth Hunt +#+LATEX_HEADER: \notindent \notag \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry} +#+LATEX: \setlength\parindent{0pt} +#+STARTUP: entitiespretty fold inlineimages + +* Design +The LIZFCM static library (at [[https://github.com/Simponic/math-4610]]) is a successor to my +attempt at writing codes for the Fundamentals of Computational Mathematics course in Common +Lisp, but the effort required to meet the requirement of creating a static library became +too difficult to integrate outside of the ~ASDF~ solution that Common Lisp already brings +to the table. + +All of the work established in ~deprecated-cl~ has been painstakingly translated into +the C programming language. I have a couple tenets for its design: + ++ Implementations of routines should all be done immutably in respect to arguments. ++ Functional programming is good (it's... rough in C though). ++ Routines are separated into "modules" that follow a form of separation of concerns + in files, and not individual files per function. + +* Compilation +A provided ~Makefile~ is added for convencience. It has been tested on an ~arm~-based M1 machine running +MacOS as well as ~x86~ Arch Linux. + +1. ~cd~ into the root of the repo +2. ~make~ + +Then, as of homework 5, the testing routines are provided in ~test~ and utilize the +~utest~ "micro"library. They compile to a binary in ~./dist/lizfcm.test~. + +Execution of the Makefile will perform compilation of individual routines. + +But, in the requirement of manual intervention (should the little alien workers +inside the computer fail to do their job), one can use the following command to +produce an object file: + +\begin{verbatim} + gcc -Iinc/ -lm -Wall -c src/.c -o build/.o +\end{verbatim} + +Which is then bundled into a static library in ~lib/lizfcm.a~ and can be linked +in the standard method. + +* The LIZFCM API +** Simple Routines +*** ~smaceps~ ++ Author: Elizabeth Hunt ++ Name: ~smaceps~ ++ Location: ~src/maceps.c~ ++ Input: none ++ Output: a ~float~ returning the specific "Machine Epsilon" of a machine on a + single precision floating point number at which it becomes "indistinguishable". + +#+BEGIN_SRC c +float smaceps() { + float one = 1.0; + float machine_epsilon = 1.0; + float one_approx = one + machine_epsilon; + + while (fabsf(one_approx - one) > 0) { + machine_epsilon /= 2; + one_approx = one + machine_epsilon; + } + + return machine_epsilon; +} +#+END_SRC + +*** ~dmaceps~ ++ Author: Elizabeth Hunt ++ Name: ~dmaceps~ ++ Location: ~src/maceps.c~ ++ Input: none ++ Output: a ~double~ returning the specific "Machine Epsilon" of a machine on a + double precision floating point number at which it becomes "indistinguishable". + +#+BEGIN_SRC c +double dmaceps() { + double one = 1.0; + double machine_epsilon = 1.0; + double one_approx = one + machine_epsilon; + + while (fabs(one_approx - one) > 0) { + machine_epsilon /= 2; + one_approx = one + machine_epsilon; + } + + return machine_epsilon; +} +#+END_SRC + +** Derivative Routines +*** ~central_derivative_at~ ++ Author: Elizabeth Hunt ++ Name: ~central_derivative_at~ ++ Location: ~src/approx_derivative.c~ ++ Input: + - ~f~ is a pointer to a one-ary function that takes a double as input and produces + a double as output + - ~a~ is the domain value at which we approximate ~f'~ + - ~h~ is the step size ++ Output: a ~double~ of the approximate value of ~f'(a)~ via the central difference + method. + +#+BEGIN_SRC c +double central_derivative_at(double (*f)(double), double a, double h) { + assert(h > 0); + + double x2 = a + h; + double x1 = a - h; + + double y2 = f(x2); + double y1 = f(x1); + + return (y2 - y1) / (x2 - x1); +} +#+END_SRC + +*** ~forward_derivative_at~ ++ Author: Elizabeth Hunt ++ Name: ~forward_derivative_at~ ++ Location: ~src/approx_derivative.c~ ++ Input: + - ~f~ is a pointer to a one-ary function that takes a double as input and produces + a double as output + - ~a~ is the domain value at which we approximate ~f'~ + - ~h~ is the step size ++ Output: a ~double~ of the approximate value of ~f'(a)~ via the forward difference + method. + +#+BEGIN_SRC c +double forward_derivative_at(double (*f)(double), double a, double h) { + assert(h > 0); + + double x2 = a + h; + double x1 = a; + + double y2 = f(x2); + double y1 = f(x1); + + return (y2 - y1) / (x2 - x1); +} +#+END_SRC + +*** ~backward_derivative_at~ ++ Author: Elizabeth Hunt ++ Name: ~backward_derivative_at~ ++ Location: ~src/approx_derivative.c~ ++ Input: + - ~f~ is a pointer to a one-ary function that takes a double as input and produces + a double as output + - ~a~ is the domain value at which we approximate ~f'~ + - ~h~ is the step size ++ Output: a ~double~ of the approximate value of ~f'(a)~ via the backward difference + method. + +#+BEGIN_SRC c +double backward_derivative_at(double (*f)(double), double a, double h) { + assert(h > 0); + + double x2 = a; + double x1 = a - h; + + double y2 = f(x2); + double y1 = f(x1); + + return (y2 - y1) / (x2 - x1); +} +#+END_SRC + +** Vector Routines +*** Vector Arithmetic: ~add_v, minus_v~ ++ Author: Elizabeth Hunt ++ Name(s): ~add_v~, ~minus_v~ ++ Location: ~src/vector.c~ ++ Input: two pointers to locations in memory wherein ~Array_double~'s lie ++ Output: a pointer to a new ~Array_double~ as the result of addition or subtraction + of the two input ~Array_double~'s + +#+BEGIN_SRC c +Array_double *add_v(Array_double *v1, Array_double *v2) { + assert(v1->size == v2->size); + + Array_double *sum = copy_vector(v1); + for (size_t i = 0; i < v1->size; i++) + sum->data[i] += v2->data[i]; + return sum; +} + +Array_double *minus_v(Array_double *v1, Array_double *v2) { + assert(v1->size == v2->size); + + Array_double *sub = InitArrayWithSize(double, v1->size, 0); + for (size_t i = 0; i < v1->size; i++) + sub->data[i] = v1->data[i] - v2->data[i]; + return sub; +} +#+END_SRC + +*** Norms: ~l1_norm~, ~l2_norm~, ~linf_norm~ ++ Author: Elizabeth Hunt ++ Name(s): ~l1_norm~, ~l2_norm~, ~linf_norm~ ++ Location: ~src/vector.c~ ++ Input: a pointer to a location in memory wherein an ~Array_double~ lies ++ Output: a ~double~ representing the value of the norm the function applies + +#+BEGIN_SRC c +double l1_norm(Array_double *v) { + double sum = 0; + for (size_t i = 0; i < v->size; ++i) + sum += fabs(v->data[i]); + return sum; +} + +double l2_norm(Array_double *v) { + double norm = 0; + for (size_t i = 0; i < v->size; ++i) + norm += v->data[i] * v->data[i]; + return sqrt(norm); +} + +double linf_norm(Array_double *v) { + assert(v->size > 0); + double max = v->data[0]; + for (size_t i = 0; i < v->size; ++i) + max = c_max(v->data[i], max); + return max; +} +#+END_SRC + +*** ~vector_distance~ ++ Author: Elizabeth Hunt ++ Name: ~vector_distance~ ++ Location: ~src/vector.c~ ++ Input: two pointers to locations in memory wherein ~Array_double~'s lie, and a pointer to a + one-ary function ~norm~ taking as input a pointer to an ~Array_double~ and returning a double + representing the norm of that ~Array_double~ + +#+BEGIN_SRC c +double vector_distance(Array_double *v1, Array_double *v2, + double (*norm)(Array_double *)) { + Array_double *minus = minus_v(v1, v2); + double dist = (*norm)(minus); + free(minus); + return dist; +} +#+END_SRC + +*** Distances: ~l1_distance~, ~l2_distance~, ~linf_distance~ ++ Author: Elizabeth Hunt ++ Name(s): ~l1_distance~, ~l2_distance~, ~linf_distance~ ++ Location: ~src/vector.c~ ++ Input: two pointers to locations in memory wherein ~Array_double~'s lie, and the distance + via the corresponding ~l1~, ~l2~, or ~linf~ norms ++ Output: A ~double~ representing the distance between the two ~Array_doubles~'s by the given + norm. + +#+BEGIN_SRC c +double l1_distance(Array_double *v1, Array_double *v2) { + return vector_distance(v1, v2, &l1_norm); +} + +double l2_distance(Array_double *v1, Array_double *v2) { + return vector_distance(v1, v2, &l2_norm); +} + +double linf_distance(Array_double *v1, Array_double *v2) { + return vector_distance(v1, v2, &linf_norm); +} +#+END_SRC + +*** ~sum_v~ ++ Author: Elizabeth Hunt ++ Name: ~sum_v~ ++ Location: ~src/vector.c~ ++ Input: a pointer to an ~Array_double~ ++ Output: a ~double~ representing the sum of all the elements of an ~Array_double~ + +#+BEGIN_SRC c +double sum_v(Array_double *v) { + double sum = 0; + for (size_t i = 0; i < v->size; i++) + sum += v->data[i]; + return sum; +} +#+END_SRC + +*** ~scale_v~ ++ Author: Elizabeth Hunt ++ Name: ~scale_v~ ++ Location: ~src/vector.c~ ++ Input: a pointer to an ~Array_double~ and a scalar ~double~ to scale the vector ++ Output: a pointer to a new ~Array_double~ of the scaled input ~Array_double~ + +#+BEGIN_SRC c +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; +} +#+END_SRC + +*** ~free_vector~ ++ Author: Elizabeth Hunt ++ Name: ~free_vector~ ++ Location: ~src/vector.c~ ++ Input: a pointer to an ~Array_double~ ++ Output: nothing. ++ Side effect: free the memory of the reserved ~Array_double~ on the heap + +#+BEGIN_SRC c +void free_vector(Array_double *v) { + free(v->data); + free(v); +} +#+END_SRC + +*** ~add_element~ ++ Author: Elizabeth Hunt ++ Name: ~add_element~ ++ Location: ~src/vector.c~ ++ Input: a pointer to an ~Array_double~ ++ Output: a new ~Array_double~ with element ~x~ appended. + +#+BEGIN_SRC c +Array_double *add_element(Array_double *v, double x) { + Array_double *pushed = InitArrayWithSize(double, v->size + 1, 0.0); + for (size_t i = 0; i < v->size; ++i) + pushed->data[i] = v->data[i]; + pushed->data[v->size] = x; + return pushed; +} +#+END_SRC + +*** ~slice_element~ ++ Author: Elizabeth Hunt ++ Name: ~slice_element~ ++ Location: ~src/vector.c~ ++ Input: a pointer to an ~Array_double~ ++ Output: a new ~Array_double~ with element ~x~ sliced. + +#+BEGIN_SRC c +Array_double *slice_element(Array_double *v, size_t x) { + Array_double *sliced = InitArrayWithSize(double, v->size - 1, 0.0); + for (size_t i = 0; i < v->size - 1; ++i) + sliced->data[i] = i >= x ? v->data[i + 1] : v->data[i]; + return sliced; +} +#+END_SRC + +*** ~copy_vector~ ++ Author: Elizabeth Hunt ++ Name: ~copy_vector~ ++ Location: ~src/vector.c~ ++ Input: a pointer to an ~Array_double~ ++ Output: a pointer to a new ~Array_double~ whose ~data~ and ~size~ are copied from the input + ~Array_double~ + +#+BEGIN_SRC c +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; +} +#+END_SRC + +*** ~format_vector_into~ ++ Author: Elizabeth Hunt ++ Name: ~format_vector_into~ ++ Location: ~src/vector.c~ ++ Input: a pointer to an ~Array_double~ and a pointer to a c-string ~s~ to "print" the vector out + into ++ Output: nothing. ++ Side effect: overwritten memory into ~s~ + +#+BEGIN_SRC c +void format_vector_into(Array_double *v, char *s) { + if (v->size == 0) { + strcat(s, "empty"); + return; + } + + for (size_t i = 0; i < v->size; ++i) { + char num[64]; + strcpy(num, ""); + + sprintf(num, "%f,", v->data[i]); + strcat(s, num); + } + strcat(s, "\n"); +} +#+END_SRC + +** Matrix Routines +*** ~lu_decomp~ ++ Author: Elizabeth Hunt ++ Name: ~lu_decomp~ ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ $m$ to decompose into a lower triangular and upper triangular + matrix $L$, $U$, respectively such that $LU = m$. ++ Output: a pointer to the location in memory in which two ~Matrix_double~'s reside: the first + representing $L$, the second, $U$. ++ Errors: Fails assertions when encountering a matrix that cannot be + decomposed + +#+BEGIN_SRC c +Matrix_double **lu_decomp(Matrix_double *m) { + assert(m->cols == m->rows); + + Matrix_double *u = copy_matrix(m); + Matrix_double *l_empt = InitMatrixWithSize(double, m->rows, m->cols, 0.0); + Matrix_double *l = put_identity_diagonal(l_empt); + free_matrix(l_empt); + + 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"); + assert(false); + } + } + + 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"); + assert(false); + } + + 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; +} +#+END_SRC +*** ~bsubst~ ++ Author: Elizabeth Hunt ++ Name: ~bsubst~ ++ Location: ~src/matrix.c~ ++ Input: a pointer to an upper-triangular ~Matrix_double~ $u$ and a ~Array_double~ + $b$ ++ Output: a pointer to a new ~Array_double~ whose entries are given by performing + back substitution + +#+BEGIN_SRC c +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; +} +#+END_SRC +*** ~fsubst~ ++ Author: Elizabeth Hunt ++ Name: ~fsubst~ ++ Location: ~src/matrix.c~ ++ Input: a pointer to a lower-triangular ~Matrix_double~ $l$ and a ~Array_double~ + $b$ ++ Output: a pointer to a new ~Array_double~ whose entries are given by performing + forward substitution + +#+BEGIN_SRC c +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; +} +#+END_SRC + +*** ~solve_matrix_lu_bsubst~ ++ Author: Elizabeth Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ $m$ and a pointer to an ~Array_double~ $b$ ++ Output: $x$ such that $mx = b$ if such a solution exists (else it's non LU-factorable as discussed + above) + +Here we make use of forward substitution to first solve $Ly = b$ given $L$ as the $L$ factor in +~lu_decomp~. Then we use back substitution to solve $Ux = y$ for $x$ similarly given $U$. + +Then, $LUx = b$, thus $x$ is a solution. + +#+BEGIN_SRC c +Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) { + assert(b->size == m->rows); + assert(m->rows == m->cols); + + Array_double *x = copy_vector(b); + Matrix_double **u_l = lu_decomp(m); + Matrix_double *u = u_l[0]; + Matrix_double *l = u_l[1]; + + Array_double *b_fsub = fsubst(l, b); + x = bsubst(u, b_fsub); + free_vector(b_fsub); + + free_matrix(u); + free_matrix(l); + free(u_l); + + return x; +} +#+END_SRC + +*** ~gaussian_elimination~ ++ Author: Elizabeth Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ $m$ ++ Output: a pointer to a copy of $m$ in reduced echelon form + +This works by finding the row with a maximum value in the column $k$. Then, it uses that as a pivot, and +applying reduction to all other rows. The general idea is available at [[https://en.wikipedia.org/wiki/Gaussian_elimination]]. + +#+BEGIN_SRC c +Matrix_double *gaussian_elimination(Matrix_double *m) { + uint64_t h = 0, k = 0; + + Matrix_double *m_cp = copy_matrix(m); + + while (h < m_cp->rows && k < m_cp->cols) { + uint64_t max_row = h; + double max_val = 0.0; + + for (uint64_t row = h; row < m_cp->rows; row++) { + double val = fabs(m_cp->data[row]->data[k]); + if (val > max_val) { + max_val = val; + max_row = row; + } + } + + if (max_val == 0.0) { + k++; + continue; + } + + if (max_row != h) { + Array_double *swp = m_cp->data[max_row]; + m_cp->data[max_row] = m_cp->data[h]; + m_cp->data[h] = swp; + } + + for (uint64_t row = h + 1; row < m_cp->rows; row++) { + double factor = m_cp->data[row]->data[k] / m_cp->data[h]->data[k]; + m_cp->data[row]->data[k] = 0.0; + + for (uint64_t col = k + 1; col < m_cp->cols; col++) { + m_cp->data[row]->data[col] -= m_cp->data[h]->data[col] * factor; + } + } + + h++; + k++; + } + + return m_cp; +} +#+END_SRC + +*** ~solve_matrix_gaussian~ ++ Author: Elizabeth Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ $m$ and a target ~Array_double~ $b$ ++ Output: a pointer to a vector $x$ being the solution to the equation $mx = b$ + +We first perform ~gaussian_elimination~ after augmenting $m$ and $b$. Then, as $m$ is in reduced echelon form, it's an upper +triangular matrix, so we can perform back substitution to compute $x$. + +#+BEGIN_SRC c +Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) { + assert(b->size == m->rows); + assert(m->rows == m->cols); + + Matrix_double *m_augment_b = add_column(m, b); + Matrix_double *eliminated = gaussian_elimination(m_augment_b); + + Array_double *b_gauss = col_v(eliminated, m->cols); + Matrix_double *u = slice_column(eliminated, m->rows); + + Array_double *solution = bsubst(u, b_gauss); + + free_matrix(m_augment_b); + free_matrix(eliminated); + free_matrix(u); + free_vector(b_gauss); + + return solution; +} +#+END_SRC + + +*** ~m_dot_v~ ++ Author: Elizabeth Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ $m$ and ~Array_double~ $v$ ++ Output: the dot product $mv$ as an ~Array_double~ + +#+BEGIN_SRC c +Array_double *m_dot_v(Matrix_double *m, Array_double *v) { + assert(v->size == m->cols); + + Array_double *product = copy_vector(v); + + for (size_t row = 0; row < v->size; ++row) + product->data[row] = v_dot_v(m->data[row], v); + + return product; +} +#+END_SRC + +*** ~put_identity_diagonal~ ++ Author: Elizabeth Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ ++ Output: a pointer to a copy to ~Matrix_double~ whose diagonal is full of 1's + +#+BEGIN_SRC c +Matrix_double *put_identity_diagonal(Matrix_double *m) { + assert(m->rows == m->cols); + Matrix_double *copy = copy_matrix(m); + for (size_t y = 0; y < m->rows; ++y) + copy->data[y]->data[y] = 1.0; + return copy; +} +#+END_SRC + +*** ~slice_column~ ++ Author: Elizabeth Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ ++ Output: a pointer to a copy of the given ~Matrix_double~ with column at ~x~ sliced + +#+BEGIN_SRC c +Matrix_double *slice_column(Matrix_double *m, size_t x) { + Matrix_double *sliced = copy_matrix(m); + + for (size_t row = 0; row < m->rows; row++) { + Array_double *old_row = sliced->data[row]; + sliced->data[row] = slice_element(old_row, x); + free_vector(old_row); + } + sliced->cols--; + + return sliced; +} +#+END_SRC + +*** ~add_column~ ++ Author: Elizabet Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ and a new vector representing the appended column ~x~ ++ Output: a pointer to a copy of the given ~Matrix_double~ with a new column ~x~ + +#+BEGIN_SRC c +Matrix_double *add_column(Matrix_double *m, Array_double *v) { + Matrix_double *pushed = copy_matrix(m); + + for (size_t row = 0; row < m->rows; row++) { + Array_double *old_row = pushed->data[row]; + pushed->data[row] = add_element(old_row, v->data[row]); + free_vector(old_row); + } + + pushed->cols++; + return pushed; +} +#+END_SRC + +*** ~copy_matrix~ ++ Author: Elizabeth Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ ++ Output: a pointer to a copy of the given ~Matrix_double~ + +#+BEGIN_SRC c +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; +} +#+END_SRC + +*** ~free_matrix~ ++ Author: Elizabeth Hunt ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ ++ Output: none. ++ Side Effects: frees memory reserved by a given ~Matrix_double~ and its member + ~Array_double~ vectors describing its rows. + +#+BEGIN_SRC c +void free_matrix(Matrix_double *m) { + for (size_t y = 0; y < m->rows; ++y) + free_vector(m->data[y]); + free(m); +} +#+END_SRC + +*** ~format_matrix_into~ ++ Author: Elizabeth Hunt ++ Name: ~format_matrix_into~ ++ Location: ~src/matrix.c~ ++ Input: a pointer to a ~Matrix_double~ and a pointer to a c-string ~s~ to "print" the vector out + into ++ Output: nothing. ++ Side effect: overwritten memory into ~s~ + +#+BEGIN_SRC c +void format_matrix_into(Matrix_double *m, char *s) { + if (m->rows == 0) + strcpy(s, "empty"); + + for (size_t y = 0; y < m->rows; ++y) { + char row_s[5192]; + strcpy(row_s, ""); + + format_vector_into(m->data[y], row_s); + strcat(s, row_s); + } + strcat(s, "\n"); +} +#+END_SRC +** Root Finding Methods +*** ~find_ivt_range~ ++ Author: Elizabeth Hunt ++ Name: ~find_ivt_range~ ++ Location: ~src/roots.c~ ++ Input: a pointer to a oneary function taking a double and producing a double, the beginning point + in $R$ to search for a range, a ~delta~ step that is taken, and a ~max_steps~ number of maximum + iterations to perform. ++ Output: a pair of ~double~'s in an ~Array_double~ representing a closed closed interval ~[beginning, end]~ + +#+BEGIN_SRC c +// f is well defined at start_x + delta*n for all n on the integer range [0, +// max_iterations] +Array_double *find_ivt_range(double (*f)(double), double start_x, double delta, + size_t max_iterations) { + double a = start_x; + + while (f(a) * f(a + delta) >= 0 && max_iterations > 0) { + max_iterations--; + a += delta; + } + + double end = a + delta; + double begin = a - delta; + + if (max_iterations == 0 && f(begin) * f(end) >= 0) + return NULL; + return InitArray(double, {begin, end}); +} +#+END_SRC +*** ~bisect_find_root~ ++ Author: Elizabeth Hunt ++ Name(s): ~bisect_find_root~ ++ Input: a one-ary function taking a double and producing a double, a closed interval represented + by ~a~ and ~b~: ~[a, b]~, a ~tolerance~ at which we return the estimated root once $b-a < \text{tolerance}$, and a + ~max_iterations~ to break us out of a loop if we can never reach the ~tolerance~. ++ Output: a vector of size of 3, ~double~'s representing first the range ~[a,b]~ and then the midpoint, + ~c~ of the range. ++ Description: recursively uses binary search to split the interval until we reach ~tolerance~. We + also assume the function ~f~ is continuous on ~[a, b]~. + +#+BEGIN_SRC c +// f is continuous on [a, b] +Array_double *bisect_find_root(double (*f)(double), double a, double b, + double tolerance, size_t max_iterations) { + assert(a <= b); + // guarantee there's a root somewhere between a and b by IVT + assert(f(a) * f(b) < 0); + + double c = (1.0 / 2) * (a + b); + if (b - a < tolerance || max_iterations == 0) + return InitArray(double, {a, b, c}); + + if (f(a) * f(c) < 0) + return bisect_find_root(f, a, c, tolerance, max_iterations - 1); + return bisect_find_root(f, c, b, tolerance, max_iterations - 1); +} +#+END_SRC +*** ~bisect_find_root_with_error_assumption~ ++ Author: Elizabeth Hunt ++ Name: ~bisect_find_root_with_error_assumption~ ++ Input: a one-ary function taking a double and producing a double, a closed interval represented + by ~a~ and ~b~: ~[a, b]~, and a ~tolerance~ equivalent to the above definition in ~bisect_find_root~ ++ Output: a ~double~ representing the estimated root ++ Description: using the bisection method we know that $e_k \le (\frac{1}{2})^k (b_0 - a_0)$. So we can + calculate $k$ at the worst possible case (that the error is exactly the tolerance) to be + $\frac{log(tolerance) - log(b_0 - a_0)}{log(\frac{1}{2})}$. We pass this value into the ~max_iterations~ + of ~bisect_find_root~ as above. +#+BEGIN_SRC c +double bisect_find_root_with_error_assumption(double (*f)(double), double a, + double b, double tolerance) { + assert(a <= b); + + uint64_t max_iterations = + (uint64_t)ceil((log(tolerance) - log(b - a)) / log(1 / 2.0)); + + Array_double *a_b_root = bisect_find_root(f, a, b, tolerance, max_iterations); + double root = a_b_root->data[2]; + free_vector(a_b_root); + + return root; +} +#+END_SRC + +*** ~fixed_point_iteration_method~ ++ Author: Elizabeth Hunt ++ Name: ~fixed_point_iteration_method~ ++ Location: ~src/roots.c~ ++ Input: a pointer to a oneary function $f$ taking a double and producing a double of which we are + trying to find a root, a guess $x_0$, and a function $g$ of the same signature of $f$ at which we + "step" our guesses according to the fixed point iteration method: $x_k = g(x_{k-1})$. Additionally, a + ~max_iterations~ representing the maximum number of "steps" to take before arriving at our + approximation and a ~tolerance~ to return our root if it becomes within [0 - tolerance, 0 + tolerance]. ++ Assumptions: $g(x)$ must be a function such that at the point $x^*$ (the found root) the derivative + $|g'(x^*)| \lt 1$ ++ Output: a double representing the found approximate root $\approx x^*$. + +#+BEGIN_SRC c +double fixed_point_iteration_method(double (*f)(double), double (*g)(double), + double x_0, double tolerance, + size_t max_iterations) { + if (max_iterations <= 0) + return x_0; + + double root = g(x_0); + if (tolerance >= fabs(f(root))) + return root; + + return fixed_point_iteration_method(f, g, root, tolerance, + max_iterations - 1); +} +#+END_SRC + +*** ~fixed_point_newton_method~ ++ Author: Elizabeth Hunt ++ Name: ~fixed_point_newton_method~ ++ Location: ~src/roots.c~ ++ Input: a pointer to a oneary function $f$ taking a double and producing a double of which we are + trying to find a root and another pointer to a function fprime of the same signature, a guess $x_0$, + and a ~max_iterations~ and ~tolerance~ as defined in the above method are required inputs. ++ Description: continually computes elements in the sequence $x_n = x_{n-1} - \frac{f(x_{n-1})}{f'p(x_{n-1})}$ ++ Output: a double representing the found approximate root $\approx x^*$ recursively applied to the sequence + given +#+BEGIN_SRC c +double fixed_point_newton_method(double (*f)(double), double (*fprime)(double), + double x_0, double tolerance, + size_t max_iterations) { + if (max_iterations <= 0) + return x_0; + + double root = x_0 - f(x_0) / fprime(x_0); + if (tolerance >= fabs(f(root))) + return root; + + return fixed_point_newton_method(f, fprime, root, tolerance, + max_iterations - 1); +} +#+END_SRC + +*** ~fixed_point_secant_method~ ++ Author: Elizabeth Hunt ++ Name: ~fixed_point_secant_method~ ++ Location: ~src/roots.c~ ++ Input: a pointer to a oneary function $f$ taking a double and producing a double of which we are + trying to find a root, a guess $x_0$ and $x_1$ in which a root lies between $[x_0, x_1]$; applying the + sequence $x_n = x_{n-1} - f(x_{n-1}) \frac{x_{n-1} - x_{n-2}}{f(x_{n-1}) - f(x_{n-2})}$. + Additionally, a ~max_iterations~ and ~tolerance~ as defined in the above method are required + inputs. ++ Output: a double representing the found approximate root $\approx x^*$ recursively applied to the sequence. +#+BEGIN_SRC c +double fixed_point_secant_method(double (*f)(double), double x_0, double x_1, + double tolerance, size_t max_iterations) { + if (max_iterations == 0) + return x_1; + + double root = x_1 - f(x_1) * ((x_1 - x_0) / (f(x_1) - f(x_0))); + + if (tolerance >= fabs(f(root))) + return root; + + return fixed_point_secant_method(f, x_1, root, tolerance, max_iterations - 1); +} +#+END_SRC +*** ~fixed_point_secant_bisection_method~ ++ Author: Elizabeth Hunt ++ Name: ~fixed_point_secant_method~ ++ Location: ~src/roots.c~ ++ Input: a pointer to a oneary function $f$ taking a double and producing a double of which we are + trying to find a root, a guess $x_0$, and a $x_1$ of which we define our first interval $[x_0, x_1]$. + Then, we perform a single iteration of the ~fixed_point_secant_method~ on this interval; if it + produces a root outside, we refresh the interval and root respectively with the given + ~bisect_find_root~ method. Additionally, a ~max_iterations~ and ~tolerance~ as defined in the above method are required + inputs. ++ Output: a double representing the found approximate root $\approx x^*$ continually applied with the + constraints defined. + +#+BEGIN_SRC c +double fixed_point_secant_bisection_method(double (*f)(double), double x_0, + double x_1, double tolerance, + size_t max_iterations) { + double begin = x_0; + double end = x_1; + double root = x_0; + + while (tolerance < fabs(f(root)) && max_iterations > 0) { + max_iterations--; + + double secant_root = fixed_point_secant_method(f, begin, end, tolerance, 1); + + if (secant_root < begin || secant_root > end) { + Array_double *range_root = bisect_find_root(f, begin, end, tolerance, 1); + + begin = range_root->data[0]; + end = range_root->data[1]; + root = range_root->data[2]; + + free_vector(range_root); + continue; + } + + root = secant_root; + + if (f(root) * f(begin) < 0) + end = secant_root; // the root exists in [begin, secant_root] + else + begin = secant_root; + } + + return root; +} +#+END_SRC + +** Linear Routines +*** ~least_squares_lin_reg~ ++ Author: Elizabeth Hunt ++ Name: ~least_squares_lin_reg~ ++ Location: ~src/lin.c~ ++ Input: two pointers to ~Array_double~'s whose entries correspond two ordered + pairs in R^2 ++ Output: a linear model best representing the ordered pairs via least squares + regression + +#+BEGIN_SRC c +Line *least_squares_lin_reg(Array_double *x, Array_double *y) { + assert(x->size == y->size); + + uint64_t n = x->size; + double sum_x = sum_v(x); + double sum_y = sum_v(y); + double sum_xy = v_dot_v(x, y); + double sum_xx = v_dot_v(x, x); + double denom = ((n * sum_xx) - (sum_x * sum_x)); + + Line *line = malloc(sizeof(Line)); + line->m = ((sum_xy * n) - (sum_x * sum_y)) / denom; + line->a = ((sum_y * sum_xx) - (sum_x * sum_xy)) / denom; + + return line; +} +#+END_SRC + +** Eigen-Adjacent +*** ~dominant_eigenvalue~ ++ Author: Elizabeth Hunt ++ Name: ~dominant_eigenvalue~ ++ Location: ~src/eigen.c~ ++ Input: a pointer to an invertible matrix ~m~, an initial eigenvector guess ~v~ (that is non + zero or orthogonal to an eigenvector with the dominant eigenvalue), a ~tolerance~ and + ~max_iterations~ that act as stop conditions ++ Output: the dominant eigenvalue with the highest magnitude, approximated with the Power + Iteration Method + +#+BEGIN_SRC c +double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance, + size_t max_iterations) { + assert(m->rows == m->cols); + assert(m->rows == v->size); + + double error = tolerance; + size_t iter = max_iterations; + double lambda = 0.0; + Array_double *eigenvector_1 = copy_vector(v); + + while (error >= tolerance && (--iter) > 0) { + Array_double *eigenvector_2 = m_dot_v(m, eigenvector_1); + Array_double *normalized_eigenvector_2 = + scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2)); + free_vector(eigenvector_2); + eigenvector_2 = normalized_eigenvector_2; + + Array_double *mx = m_dot_v(m, eigenvector_2); + double new_lambda = + v_dot_v(mx, eigenvector_2) / v_dot_v(eigenvector_2, eigenvector_2); + + error = fabs(new_lambda - lambda); + lambda = new_lambda; + free_vector(eigenvector_1); + eigenvector_1 = eigenvector_2; + } + + return lambda; +} +#+END_SRC +*** ~shift_inverse_power_eigenvalue~ ++ Author: Elizabeth Hunt ++ Name: ~least_dominant_eigenvalue~ ++ Location: ~src/eigen.c~ ++ Input: a pointer to an invertible matrix ~m~, an initial eigenvector guess ~v~ (that is non + zero or orthogonal to an eigenvector with the dominant eigenvalue), a ~shift~ to act as the + shifted \delta, and ~tolerance~ and ~max_iterations~ that act as stop conditions. ++ Output: the eigenvalue closest to ~shift~ with the lowest magnitude closest to 0, approximated + with the Inverse Power Iteration Method +#+BEGIN_SRC c +double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v, + double shift, double tolerance, + size_t max_iterations) { + assert(m->rows == m->cols); + assert(m->rows == v->size); + + Matrix_double *m_c = copy_matrix(m); + for (size_t y = 0; y < m_c->rows; ++y) + m_c->data[y]->data[y] = m_c->data[y]->data[y] - shift; + + double error = tolerance; + size_t iter = max_iterations; + double lambda = shift; + Array_double *eigenvector_1 = copy_vector(v); + + while (error >= tolerance && (--iter) > 0) { + Array_double *eigenvector_2 = solve_matrix_lu_bsubst(m_c, eigenvector_1); + Array_double *normalized_eigenvector_2 = + scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2)); + free_vector(eigenvector_2); + + Array_double *mx = m_dot_v(m, normalized_eigenvector_2); + double new_lambda = + v_dot_v(mx, normalized_eigenvector_2) / + v_dot_v(normalized_eigenvector_2, normalized_eigenvector_2); + + error = fabs(new_lambda - lambda); + lambda = new_lambda; + free_vector(eigenvector_1); + eigenvector_1 = normalized_eigenvector_2; + } + + return lambda; +} +#+END_SRC + +*** ~least_dominant_eigenvalue~ ++ Author: Elizabeth Hunt ++ Name: ~least_dominant_eigenvalue~ ++ Location: ~src/eigen.c~ ++ Input: a pointer to an invertible matrix ~m~, an initial eigenvector guess ~v~ (that is non + zero or orthogonal to an eigenvector with the dominant eigenvalue), a ~tolerance~ and + ~max_iterations~ that act as stop conditions. ++ Output: the least dominant eigenvalue with the lowest magnitude closest to 0, approximated + with the Inverse Power Iteration Method. +#+BEGIN_SRC c +double least_dominant_eigenvalue(Matrix_double *m, Array_double *v, + double tolerance, size_t max_iterations) { + return shift_inverse_power_eigenvalue(m, v, 0.0, tolerance, max_iterations); +} +#+END_SRC +*** ~partition_find_eigenvalues~ ++ Author: Elizabeth Hunt ++ Name: ~partition_find_eigenvalues~ ++ Location: ~src/eigen.c~ ++ Input: a pointer to an invertible matrix ~m~, a matrix whose rows correspond to initial + eigenvector guesses at each "partition" which is computed from a uniform distribution + between the number of rows this "guess matrix" has and the distance between the least + dominant eigenvalue and the most dominant. Additionally, a ~max_iterations~ and a ~tolerance~ + that act as stop conditions. ++ Output: a vector of ~doubles~ corresponding to the "nearest" eigenvalue at the midpoint of + each partition, via the given guess of that partition. +#+BEGIN_SRC c +Array_double *partition_find_eigenvalues(Matrix_double *m, + Matrix_double *guesses, + double tolerance, + size_t max_iterations) { + assert(guesses->rows >= + 2); // we need at least, the most and least dominant eigenvalues + + double end = dominant_eigenvalue(m, guesses->data[guesses->rows - 1], + tolerance, max_iterations); + double begin = + least_dominant_eigenvalue(m, guesses->data[0], tolerance, max_iterations); + + double delta = (end - begin) / guesses->rows; + Array_double *eigenvalues = InitArrayWithSize(double, guesses->rows, 0.0); + for (size_t i = 0; i < guesses->rows; i++) { + double box_midpoint = ((delta * i) + (delta * (i + 1))) / 2; + + double nearest_eigenvalue = shift_inverse_power_eigenvalue( + m, guesses->data[i], box_midpoint, tolerance, max_iterations); + + eigenvalues->data[i] = nearest_eigenvalue; + } + + return eigenvalues; +} +#+END_SRC +*** ~leslie_matrix~ ++ Author: Elizabeth Hunt ++ Name: ~leslie_matrix~ ++ Location: ~src/eigen.c~ ++ Input: two pointers to ~Array_double~'s representing the ratio of individuals in an age class + $x$ getting to the next age class $x+1$ and the number of offspring that individuals in an age + class create in age class 0. ++ Output: the leslie matrix generated from the input vectors. + +#+BEGIN_SRC c +Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio, + Array_double *age_class_offspring) { + assert(age_class_surivor_ratio->size + 1 == age_class_offspring->size); + + Matrix_double *leslie = InitMatrixWithSize(double, age_class_offspring->size, + age_class_offspring->size, 0.0); + + free_vector(leslie->data[0]); + leslie->data[0] = age_class_offspring; + + for (size_t i = 0; i < age_class_surivor_ratio->size; i++) + leslie->data[i + 1]->data[i] = age_class_surivor_ratio->data[i]; + return leslie; +} +#+END_SRC +** Jacobi / Gauss-Siedel +*** ~jacobi_solve~ ++ Author: Elizabeth Hunt ++ Name: ~jacobi_solve~ ++ Location: ~src/matrix.c~ ++ Input: a pointer to a diagonally dominant square matrix $m$, a vector representing + the value $b$ in $mx = b$, a double representing the maximum distance between + the solutions produced by iteration $i$ and $i+1$ (by L2 norm a.k.a cartesian + distance), and a ~max_iterations~ which we force stop. ++ Output: the converged-upon solution $x$ to $mx = b$ +#+BEGIN_SRC c +Array_double *jacobi_solve(Matrix_double *m, Array_double *b, + double l2_convergence_tolerance, + size_t max_iterations) { + assert(m->rows == m->cols); + assert(b->size == m->cols); + size_t iter = max_iterations; + + Array_double *x_k = InitArrayWithSize(double, b->size, 0.0); + Array_double *x_k_1 = + InitArrayWithSize(double, b->size, rand_from(0.1, 10.0)); + + while ((--iter) > 0 && l2_distance(x_k_1, x_k) > l2_convergence_tolerance) { + for (size_t i = 0; i < m->rows; i++) { + double delta = 0.0; + for (size_t j = 0; j < m->cols; j++) { + if (i == j) + continue; + delta += m->data[i]->data[j] * x_k->data[j]; + } + x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i]; + } + + Array_double *tmp = x_k; + x_k = x_k_1; + x_k_1 = tmp; + } + + free_vector(x_k); + return x_k_1; +} +#+END_SRC + +*** ~gauss_siedel_solve~ ++ Author: Elizabeth Hunt ++ Name: ~gauss_siedel_solve~ ++ Location: ~src/matrix.c~ ++ Input: a pointer to a [[https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method][diagonally dominant or symmetric and positive definite]] + square matrix $m$, a vector representing + the value $b$ in $mx = b$, a double representing the maximum distance between + the solutions produced by iteration $i$ and $i+1$ (by L2 norm a.k.a cartesian + distance), and a ~max_iterations~ which we force stop. ++ Output: the converged-upon solution $x$ to $mx = b$ ++ Description: we use almost the exact same method as ~jacobi_solve~ but modify + only one array in accordance to the Gauss-Siedel method, but which is necessarily + copied before due to the convergence check. +#+BEGIN_SRC c + +Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b, + double l2_convergence_tolerance, + size_t max_iterations) { + assert(m->rows == m->cols); + assert(b->size == m->cols); + size_t iter = max_iterations; + + Array_double *x_k = InitArrayWithSize(double, b->size, 0.0); + Array_double *x_k_1 = + InitArrayWithSize(double, b->size, rand_from(0.1, 10.0)); + + while ((--iter) > 0) { + for (size_t i = 0; i < x_k->size; i++) + x_k->data[i] = x_k_1->data[i]; + + for (size_t i = 0; i < m->rows; i++) { + double delta = 0.0; + for (size_t j = 0; j < m->cols; j++) { + if (i == j) + continue; + delta += m->data[i]->data[j] * x_k_1->data[j]; + } + x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i]; + } + + if (l2_distance(x_k_1, x_k) <= l2_convergence_tolerance) + break; + } + + free_vector(x_k); + return x_k_1; +} +#+END_SRC + +** Appendix / Miscellaneous +*** Random ++ Author: Elizabeth Hunt ++ Name: ~rand_from~ ++ Location: ~src/rand.c~ ++ Input: a pair of doubles, min and max to generate a random number min + \le x \le max ++ Output: a random double in the constraints shown + +#+BEGIN_SRC c +double rand_from(double min, double max) { + return min + (rand() / (RAND_MAX / (max - min))); +} +#+END_SRC +*** Data Types +**** ~Line~ ++ Author: Elizabeth Hunt ++ Location: ~inc/types.h~ + +#+BEGIN_SRC c +typedef struct Line { + double m; + double a; +} Line; +#+END_SRC +**** The ~Array_~ and ~Matrix_~ ++ Author: Elizabeth Hunt ++ Location: ~inc/types.h~ + +We define two Pre processor Macros ~DEFINE_ARRAY~ and ~DEFINE_MATRIX~ that take +as input a type, and construct a struct definition for the given type for +convenient access to the vector or matrices dimensions. + +Such that ~DEFINE_ARRAY(int)~ would expand to: + +#+BEGIN_SRC c + typedef struct { + int* data; + size_t size; + } Array_int +#+END_SRC + +And ~DEFINE_MATRIX(int)~ would expand a to ~Matrix_int~; containing a pointer to +a collection of pointers of ~Array_int~'s and its dimensions. + +#+BEGIN_SRC c + typedef struct { + Array_int **data; + size_t cols; + size_t rows; + } Matrix_int +#+END_SRC + +*** Macros +**** ~c_max~ and ~c_min~ ++ Author: Elizabeth Hunt ++ Location: ~inc/macros.h~ ++ Input: two structures that define an order measure ++ Output: either the larger or smaller of the two depending on the measure + +#+BEGIN_SRC c +#define c_max(x, y) (((x) >= (y)) ? (x) : (y)) +#define c_min(x, y) (((x) <= (y)) ? (x) : (y)) +#+END_SRC + +**** ~InitArray~ ++ Author: Elizabeth Hunt ++ Location: ~inc/macros.h~ ++ Input: a type and array of values to initialze an array with such type ++ Output: a new ~Array_type~ with the size of the given array and its data + +#+BEGIN_SRC c +#define InitArray(TYPE, ...) \ + ({ \ + TYPE temp[] = __VA_ARGS__; \ + Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ + arr->size = sizeof(temp) / sizeof(temp[0]); \ + arr->data = malloc(arr->size * sizeof(TYPE)); \ + memcpy(arr->data, temp, arr->size * sizeof(TYPE)); \ + arr; \ + }) +#+END_SRC + +**** ~InitArrayWithSize~ ++ Author: Elizabeth Hunt ++ Location: ~inc/macros.h~ ++ Input: a type, a size, and initial value ++ Output: a new ~Array_type~ with the given size filled with the initial value + +#+BEGIN_SRC c +#define InitArrayWithSize(TYPE, SIZE, INIT_VALUE) \ + ({ \ + Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ + arr->size = SIZE; \ + arr->data = malloc(arr->size * sizeof(TYPE)); \ + for (size_t i = 0; i < arr->size; i++) \ + arr->data[i] = INIT_VALUE; \ + arr; \ + }) +#+END_SRC + +**** ~InitMatrixWithSize~ ++ Author: Elizabeth Hunt ++ Location: ~inc/macros.h~ ++ Input: a type, number of rows, columns, and initial value ++ Output: a new ~Matrix_type~ of size ~rows x columns~ filled with the initial + value + +#+BEGIN_SRC c +#define InitMatrixWithSize(TYPE, ROWS, COLS, INIT_VALUE) \ + ({ \ + Matrix_##TYPE *matrix = malloc(sizeof(Matrix_##TYPE)); \ + matrix->rows = ROWS; \ + matrix->cols = COLS; \ + matrix->data = malloc(matrix->rows * sizeof(Array_##TYPE *)); \ + for (size_t y = 0; y < matrix->rows; y++) \ + matrix->data[y] = InitArrayWithSize(TYPE, COLS, INIT_VALUE); \ + matrix; \ + }) +#+END_SRC + diff --git a/Homework/math4610/doc/software_manual.pdf b/Homework/math4610/doc/software_manual.pdf new file mode 100644 index 0000000..1772e7c Binary files /dev/null and b/Homework/math4610/doc/software_manual.pdf differ diff --git a/Homework/math4610/doc/software_manual.tex b/Homework/math4610/doc/software_manual.tex new file mode 100644 index 0000000..dac099b --- /dev/null +++ b/Homework/math4610/doc/software_manual.tex @@ -0,0 +1,1583 @@ +% Created 2023-12-11 Mon 19:22 +% Intended LaTeX compiler: pdflatex +\documentclass[11pt]{article} +\usepackage[utf8]{inputenc} +\usepackage[T1]{fontenc} +\usepackage{graphicx} +\usepackage{longtable} +\usepackage{wrapfig} +\usepackage{rotating} +\usepackage[normalem]{ulem} +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{capt-of} +\usepackage{hyperref} +\notindent \notag \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry} +\author{Elizabeth Hunt} +\date{\today} +\title{LIZFCM Software Manual (v0.6)} +\hypersetup{ + pdfauthor={Elizabeth Hunt}, + pdftitle={LIZFCM Software Manual (v0.6)}, + pdfkeywords={}, + pdfsubject={}, + pdfcreator={Emacs 28.2 (Org mode 9.7-pre)}, + pdflang={English}} +\begin{document} + +\maketitle +\tableofcontents + +\setlength\parindent{0pt} + +\section{Design} +\label{sec:org63deaf6} +The LIZFCM static library (at \url{https://github.com/Simponic/math-4610}) is a successor to my +attempt at writing codes for the Fundamentals of Computational Mathematics course in Common +Lisp, but the effort required to meet the requirement of creating a static library became +too difficult to integrate outside of the \texttt{ASDF} solution that Common Lisp already brings +to the table. + +All of the work established in \texttt{deprecated-cl} has been painstakingly translated into +the C programming language. I have a couple tenets for its design: + +\begin{itemize} +\item Implementations of routines should all be done immutably in respect to arguments. +\item Functional programming is good (it's\ldots{} rough in C though). +\item Routines are separated into "modules" that follow a form of separation of concerns +in files, and not individual files per function. +\end{itemize} + +\section{Compilation} +\label{sec:org7291327} +A provided \texttt{Makefile} is added for convencience. It has been tested on an \texttt{arm}-based M1 machine running +MacOS as well as \texttt{x86} Arch Linux. + +\begin{enumerate} +\item \texttt{cd} into the root of the repo +\item \texttt{make} +\end{enumerate} + +Then, as of homework 5, the testing routines are provided in \texttt{test} and utilize the +\texttt{utest} "micro"library. They compile to a binary in \texttt{./dist/lizfcm.test}. + +Execution of the Makefile will perform compilation of individual routines. + +But, in the requirement of manual intervention (should the little alien workers +inside the computer fail to do their job), one can use the following command to +produce an object file: + +\begin{verbatim} + gcc -Iinc/ -lm -Wall -c src/.c -o build/.o +\end{verbatim} + +Which is then bundled into a static library in \texttt{lib/lizfcm.a} and can be linked +in the standard method. + +\section{The LIZFCM API} +\label{sec:org1ebe7fa} +\subsection{Simple Routines} +\label{sec:orgff18c6b} +\subsubsection{\texttt{smaceps}} +\label{sec:org443df5e} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{smaceps} +\item Location: \texttt{src/maceps.c} +\item Input: none +\item Output: a \texttt{float} returning the specific "Machine Epsilon" of a machine on a +single precision floating point number at which it becomes "indistinguishable". +\end{itemize} + +\begin{verbatim} +float smaceps() { + float one = 1.0; + float machine_epsilon = 1.0; + float one_approx = one + machine_epsilon; + + while (fabsf(one_approx - one) > 0) { + machine_epsilon /= 2; + one_approx = one + machine_epsilon; + } + + return machine_epsilon; +} +\end{verbatim} + +\subsubsection{\texttt{dmaceps}} +\label{sec:org5121603} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{dmaceps} +\item Location: \texttt{src/maceps.c} +\item Input: none +\item Output: a \texttt{double} returning the specific "Machine Epsilon" of a machine on a +double precision floating point number at which it becomes "indistinguishable". +\end{itemize} + +\begin{verbatim} +double dmaceps() { + double one = 1.0; + double machine_epsilon = 1.0; + double one_approx = one + machine_epsilon; + + while (fabs(one_approx - one) > 0) { + machine_epsilon /= 2; + one_approx = one + machine_epsilon; + } + + return machine_epsilon; +} +\end{verbatim} + +\subsection{Derivative Routines} +\label{sec:org6fd324c} +\subsubsection{\texttt{central\_derivative\_at}} +\label{sec:orge9f0821} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{central\_derivative\_at} +\item Location: \texttt{src/approx\_derivative.c} +\item Input: +\begin{itemize} +\item \texttt{f} is a pointer to a one-ary function that takes a double as input and produces +a double as output +\item \texttt{a} is the domain value at which we approximate \texttt{f'} +\item \texttt{h} is the step size +\end{itemize} +\item Output: a \texttt{double} of the approximate value of \texttt{f'(a)} via the central difference +method. +\end{itemize} + +\begin{verbatim} +double central_derivative_at(double (*f)(double), double a, double h) { + assert(h > 0); + + double x2 = a + h; + double x1 = a - h; + + double y2 = f(x2); + double y1 = f(x1); + + return (y2 - y1) / (x2 - x1); +} +\end{verbatim} + +\subsubsection{\texttt{forward\_derivative\_at}} +\label{sec:org8720f28} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{forward\_derivative\_at} +\item Location: \texttt{src/approx\_derivative.c} +\item Input: +\begin{itemize} +\item \texttt{f} is a pointer to a one-ary function that takes a double as input and produces +a double as output +\item \texttt{a} is the domain value at which we approximate \texttt{f'} +\item \texttt{h} is the step size +\end{itemize} +\item Output: a \texttt{double} of the approximate value of \texttt{f'(a)} via the forward difference +method. +\end{itemize} + +\begin{verbatim} +double forward_derivative_at(double (*f)(double), double a, double h) { + assert(h > 0); + + double x2 = a + h; + double x1 = a; + + double y2 = f(x2); + double y1 = f(x1); + + return (y2 - y1) / (x2 - x1); +} +\end{verbatim} + +\subsubsection{\texttt{backward\_derivative\_at}} +\label{sec:org1589b19} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{backward\_derivative\_at} +\item Location: \texttt{src/approx\_derivative.c} +\item Input: +\begin{itemize} +\item \texttt{f} is a pointer to a one-ary function that takes a double as input and produces +a double as output +\item \texttt{a} is the domain value at which we approximate \texttt{f'} +\item \texttt{h} is the step size +\end{itemize} +\item Output: a \texttt{double} of the approximate value of \texttt{f'(a)} via the backward difference +method. +\end{itemize} + +\begin{verbatim} +double backward_derivative_at(double (*f)(double), double a, double h) { + assert(h > 0); + + double x2 = a; + double x1 = a - h; + + double y2 = f(x2); + double y1 = f(x1); + + return (y2 - y1) / (x2 - x1); +} +\end{verbatim} + +\subsection{Vector Routines} +\label{sec:org493841e} +\subsubsection{Vector Arithmetic: \texttt{add\_v, minus\_v}} +\label{sec:org3912c29} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name(s): \texttt{add\_v}, \texttt{minus\_v} +\item Location: \texttt{src/vector.c} +\item Input: two pointers to locations in memory wherein \texttt{Array\_double}'s lie +\item Output: a pointer to a new \texttt{Array\_double} as the result of addition or subtraction +of the two input \texttt{Array\_double}'s +\end{itemize} + +\begin{verbatim} +Array_double *add_v(Array_double *v1, Array_double *v2) { + assert(v1->size == v2->size); + + Array_double *sum = copy_vector(v1); + for (size_t i = 0; i < v1->size; i++) + sum->data[i] += v2->data[i]; + return sum; +} + +Array_double *minus_v(Array_double *v1, Array_double *v2) { + assert(v1->size == v2->size); + + Array_double *sub = InitArrayWithSize(double, v1->size, 0); + for (size_t i = 0; i < v1->size; i++) + sub->data[i] = v1->data[i] - v2->data[i]; + return sub; +} +\end{verbatim} + +\subsubsection{Norms: \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm}} +\label{sec:orged74cfb} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name(s): \texttt{l1\_norm}, \texttt{l2\_norm}, \texttt{linf\_norm} +\item Location: \texttt{src/vector.c} +\item Input: a pointer to a location in memory wherein an \texttt{Array\_double} lies +\item Output: a \texttt{double} representing the value of the norm the function applies +\end{itemize} + +\begin{verbatim} +double l1_norm(Array_double *v) { + double sum = 0; + for (size_t i = 0; i < v->size; ++i) + sum += fabs(v->data[i]); + return sum; +} + +double l2_norm(Array_double *v) { + double norm = 0; + for (size_t i = 0; i < v->size; ++i) + norm += v->data[i] * v->data[i]; + return sqrt(norm); +} + +double linf_norm(Array_double *v) { + assert(v->size > 0); + double max = v->data[0]; + for (size_t i = 0; i < v->size; ++i) + max = c_max(v->data[i], max); + return max; +} +\end{verbatim} + +\subsubsection{\texttt{vector\_distance}} +\label{sec:org20a5773} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{vector\_distance} +\item Location: \texttt{src/vector.c} +\item Input: two pointers to locations in memory wherein \texttt{Array\_double}'s lie, and a pointer to a +one-ary function \texttt{norm} taking as input a pointer to an \texttt{Array\_double} and returning a double +representing the norm of that \texttt{Array\_double} +\end{itemize} + +\begin{verbatim} +double vector_distance(Array_double *v1, Array_double *v2, + double (*norm)(Array_double *)) { + Array_double *minus = minus_v(v1, v2); + double dist = (*norm)(minus); + free(minus); + return dist; +} +\end{verbatim} + +\subsubsection{Distances: \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance}} +\label{sec:orgac16178} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name(s): \texttt{l1\_distance}, \texttt{l2\_distance}, \texttt{linf\_distance} +\item Location: \texttt{src/vector.c} +\item Input: two pointers to locations in memory wherein \texttt{Array\_double}'s lie, and the distance +via the corresponding \texttt{l1}, \texttt{l2}, or \texttt{linf} norms +\item Output: A \texttt{double} representing the distance between the two \texttt{Array\_doubles}'s by the given +norm. +\end{itemize} + +\begin{verbatim} +double l1_distance(Array_double *v1, Array_double *v2) { + return vector_distance(v1, v2, &l1_norm); +} + +double l2_distance(Array_double *v1, Array_double *v2) { + return vector_distance(v1, v2, &l2_norm); +} + +double linf_distance(Array_double *v1, Array_double *v2) { + return vector_distance(v1, v2, &linf_norm); +} +\end{verbatim} + +\subsubsection{\texttt{sum\_v}} +\label{sec:org876aafa} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{sum\_v} +\item Location: \texttt{src/vector.c} +\item Input: a pointer to an \texttt{Array\_double} +\item Output: a \texttt{double} representing the sum of all the elements of an \texttt{Array\_double} +\end{itemize} + +\begin{verbatim} +double sum_v(Array_double *v) { + double sum = 0; + for (size_t i = 0; i < v->size; i++) + sum += v->data[i]; + return sum; +} +\end{verbatim} + +\subsubsection{\texttt{scale\_v}} +\label{sec:orgf1d236c} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{scale\_v} +\item Location: \texttt{src/vector.c} +\item Input: a pointer to an \texttt{Array\_double} and a scalar \texttt{double} to scale the vector +\item Output: a pointer to a new \texttt{Array\_double} of the scaled input \texttt{Array\_double} +\end{itemize} + +\begin{verbatim} +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; +} +\end{verbatim} + +\subsubsection{\texttt{free\_vector}} +\label{sec:org2ca163d} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{free\_vector} +\item Location: \texttt{src/vector.c} +\item Input: a pointer to an \texttt{Array\_double} +\item Output: nothing. +\item Side effect: free the memory of the reserved \texttt{Array\_double} on the heap +\end{itemize} + +\begin{verbatim} +void free_vector(Array_double *v) { + free(v->data); + free(v); +} +\end{verbatim} + +\subsubsection{\texttt{add\_element}} +\label{sec:org7a99233} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{add\_element} +\item Location: \texttt{src/vector.c} +\item Input: a pointer to an \texttt{Array\_double} +\item Output: a new \texttt{Array\_double} with element \texttt{x} appended. +\end{itemize} + +\begin{verbatim} +Array_double *add_element(Array_double *v, double x) { + Array_double *pushed = InitArrayWithSize(double, v->size + 1, 0.0); + for (size_t i = 0; i < v->size; ++i) + pushed->data[i] = v->data[i]; + pushed->data[v->size] = x; + return pushed; +} +\end{verbatim} + +\subsubsection{\texttt{slice\_element}} +\label{sec:org6c07c99} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{slice\_element} +\item Location: \texttt{src/vector.c} +\item Input: a pointer to an \texttt{Array\_double} +\item Output: a new \texttt{Array\_double} with element \texttt{x} sliced. +\end{itemize} + +\begin{verbatim} +Array_double *slice_element(Array_double *v, size_t x) { + Array_double *sliced = InitArrayWithSize(double, v->size - 1, 0.0); + for (size_t i = 0; i < v->size - 1; ++i) + sliced->data[i] = i >= x ? v->data[i + 1] : v->data[i]; + return sliced; +} +\end{verbatim} + +\subsubsection{\texttt{copy\_vector}} +\label{sec:org81f7cc1} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{copy\_vector} +\item Location: \texttt{src/vector.c} +\item Input: a pointer to an \texttt{Array\_double} +\item Output: a pointer to a new \texttt{Array\_double} whose \texttt{data} and \texttt{size} are copied from the input +\texttt{Array\_double} +\end{itemize} + +\begin{verbatim} +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; +} +\end{verbatim} + +\subsubsection{\texttt{format\_vector\_into}} +\label{sec:orgd168171} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{format\_vector\_into} +\item Location: \texttt{src/vector.c} +\item Input: a pointer to an \texttt{Array\_double} and a pointer to a c-string \texttt{s} to "print" the vector out +into +\item Output: nothing. +\item Side effect: overwritten memory into \texttt{s} +\end{itemize} + +\begin{verbatim} +void format_vector_into(Array_double *v, char *s) { + if (v->size == 0) { + strcat(s, "empty"); + return; + } + + for (size_t i = 0; i < v->size; ++i) { + char num[64]; + strcpy(num, ""); + + sprintf(num, "%f,", v->data[i]); + strcat(s, num); + } + strcat(s, "\n"); +} +\end{verbatim} + +\subsection{Matrix Routines} +\label{sec:org5c45c12} +\subsubsection{\texttt{lu\_decomp}} +\label{sec:orgf1e0ac3} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{lu\_decomp} +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} \(m\) to decompose into a lower triangular and upper triangular +matrix \(L\), \(U\), respectively such that \(LU = m\). +\item Output: a pointer to the location in memory in which two \texttt{Matrix\_double}'s reside: the first +representing \(L\), the second, \(U\). +\item Errors: Fails assertions when encountering a matrix that cannot be +decomposed +\end{itemize} + +\begin{verbatim} +Matrix_double **lu_decomp(Matrix_double *m) { + assert(m->cols == m->rows); + + Matrix_double *u = copy_matrix(m); + Matrix_double *l_empt = InitMatrixWithSize(double, m->rows, m->cols, 0.0); + Matrix_double *l = put_identity_diagonal(l_empt); + free_matrix(l_empt); + + 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"); + assert(false); + } + } + + 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"); + assert(false); + } + + 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; +} +\end{verbatim} +\subsubsection{\texttt{bsubst}} +\label{sec:orgec7e4b5} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{bsubst} +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to an upper-triangular \texttt{Matrix\_double} \(u\) and a \texttt{Array\_double} +\(b\) +\item Output: a pointer to a new \texttt{Array\_double} whose entries are given by performing +back substitution +\end{itemize} + +\begin{verbatim} +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; +} +\end{verbatim} +\subsubsection{\texttt{fsubst}} +\label{sec:org72ff2ed} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{fsubst} +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a lower-triangular \texttt{Matrix\_double} \(l\) and a \texttt{Array\_double} +\(b\) +\item Output: a pointer to a new \texttt{Array\_double} whose entries are given by performing +forward substitution +\end{itemize} + +\begin{verbatim} +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; +} +\end{verbatim} + +\subsubsection{\texttt{solve\_matrix\_lu\_bsubst}} +\label{sec:orga735557} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} \(m\) and a pointer to an \texttt{Array\_double} \(b\) +\item Output: \(x\) such that \(mx = b\) if such a solution exists (else it's non LU-factorable as discussed +above) +\end{itemize} + +Here we make use of forward substitution to first solve \(Ly = b\) given \(L\) as the \(L\) factor in +\texttt{lu\_decomp}. Then we use back substitution to solve \(Ux = y\) for \(x\) similarly given \(U\). + +Then, \(LUx = b\), thus \(x\) is a solution. + +\begin{verbatim} +Array_double *solve_matrix_lu_bsubst(Matrix_double *m, Array_double *b) { + assert(b->size == m->rows); + assert(m->rows == m->cols); + + Array_double *x = copy_vector(b); + Matrix_double **u_l = lu_decomp(m); + Matrix_double *u = u_l[0]; + Matrix_double *l = u_l[1]; + + Array_double *b_fsub = fsubst(l, b); + x = bsubst(u, b_fsub); + free_vector(b_fsub); + + free_matrix(u); + free_matrix(l); + free(u_l); + + return x; +} +\end{verbatim} + +\subsubsection{\texttt{gaussian\_elimination}} +\label{sec:org71d2519} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} \(m\) +\item Output: a pointer to a copy of \(m\) in reduced echelon form +\end{itemize} + +This works by finding the row with a maximum value in the column \(k\). Then, it uses that as a pivot, and +applying reduction to all other rows. The general idea is available at \url{https://en.wikipedia.org/wiki/Gaussian\_elimination}. + +\begin{verbatim} +Matrix_double *gaussian_elimination(Matrix_double *m) { + uint64_t h = 0, k = 0; + + Matrix_double *m_cp = copy_matrix(m); + + while (h < m_cp->rows && k < m_cp->cols) { + uint64_t max_row = h; + double max_val = 0.0; + + for (uint64_t row = h; row < m_cp->rows; row++) { + double val = fabs(m_cp->data[row]->data[k]); + if (val > max_val) { + max_val = val; + max_row = row; + } + } + + if (max_val == 0.0) { + k++; + continue; + } + + if (max_row != h) { + Array_double *swp = m_cp->data[max_row]; + m_cp->data[max_row] = m_cp->data[h]; + m_cp->data[h] = swp; + } + + for (uint64_t row = h + 1; row < m_cp->rows; row++) { + double factor = m_cp->data[row]->data[k] / m_cp->data[h]->data[k]; + m_cp->data[row]->data[k] = 0.0; + + for (uint64_t col = k + 1; col < m_cp->cols; col++) { + m_cp->data[row]->data[col] -= m_cp->data[h]->data[col] * factor; + } + } + + h++; + k++; + } + + return m_cp; +} +\end{verbatim} + +\subsubsection{\texttt{solve\_matrix\_gaussian}} +\label{sec:org230915f} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} \(m\) and a target \texttt{Array\_double} \(b\) +\item Output: a pointer to a vector \(x\) being the solution to the equation \(mx = b\) +\end{itemize} + +We first perform \texttt{gaussian\_elimination} after augmenting \(m\) and \(b\). Then, as \(m\) is in reduced echelon form, it's an upper +triangular matrix, so we can perform back substitution to compute \(x\). + +\begin{verbatim} +Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) { + assert(b->size == m->rows); + assert(m->rows == m->cols); + + Matrix_double *m_augment_b = add_column(m, b); + Matrix_double *eliminated = gaussian_elimination(m_augment_b); + + Array_double *b_gauss = col_v(eliminated, m->cols); + Matrix_double *u = slice_column(eliminated, m->rows); + + Array_double *solution = bsubst(u, b_gauss); + + free_matrix(m_augment_b); + free_matrix(eliminated); + free_matrix(u); + free_vector(b_gauss); + + return solution; +} +\end{verbatim} + + +\subsubsection{\texttt{m\_dot\_v}} +\label{sec:org83c8351} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} \(m\) and \texttt{Array\_double} \(v\) +\item Output: the dot product \(mv\) as an \texttt{Array\_double} +\end{itemize} + +\begin{verbatim} +Array_double *m_dot_v(Matrix_double *m, Array_double *v) { + assert(v->size == m->cols); + + Array_double *product = copy_vector(v); + + for (size_t row = 0; row < v->size; ++row) + product->data[row] = v_dot_v(m->data[row], v); + + return product; +} +\end{verbatim} + +\subsubsection{\texttt{put\_identity\_diagonal}} +\label{sec:orge3fcb3e} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} +\item Output: a pointer to a copy to \texttt{Matrix\_double} whose diagonal is full of 1's +\end{itemize} + +\begin{verbatim} +Matrix_double *put_identity_diagonal(Matrix_double *m) { + assert(m->rows == m->cols); + Matrix_double *copy = copy_matrix(m); + for (size_t y = 0; y < m->rows; ++y) + copy->data[y]->data[y] = 1.0; + return copy; +} +\end{verbatim} + +\subsubsection{\texttt{slice\_column}} +\label{sec:org95e39ba} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} +\item Output: a pointer to a copy of the given \texttt{Matrix\_double} with column at \texttt{x} sliced +\end{itemize} + +\begin{verbatim} +Matrix_double *slice_column(Matrix_double *m, size_t x) { + Matrix_double *sliced = copy_matrix(m); + + for (size_t row = 0; row < m->rows; row++) { + Array_double *old_row = sliced->data[row]; + sliced->data[row] = slice_element(old_row, x); + free_vector(old_row); + } + sliced->cols--; + + return sliced; +} +\end{verbatim} + +\subsubsection{\texttt{add\_column}} +\label{sec:org9a2ad93} +\begin{itemize} +\item Author: Elizabet Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} and a new vector representing the appended column \texttt{x} +\item Output: a pointer to a copy of the given \texttt{Matrix\_double} with a new column \texttt{x} +\end{itemize} + +\begin{verbatim} +Matrix_double *add_column(Matrix_double *m, Array_double *v) { + Matrix_double *pushed = copy_matrix(m); + + for (size_t row = 0; row < m->rows; row++) { + Array_double *old_row = pushed->data[row]; + pushed->data[row] = add_element(old_row, v->data[row]); + free_vector(old_row); + } + + pushed->cols++; + return pushed; +} +\end{verbatim} + +\subsubsection{\texttt{copy\_matrix}} +\label{sec:org63765c0} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} +\item Output: a pointer to a copy of the given \texttt{Matrix\_double} +\end{itemize} + +\begin{verbatim} +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; +} +\end{verbatim} + +\subsubsection{\texttt{free\_matrix}} +\label{sec:orgc337967} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} +\item Output: none. +\item Side Effects: frees memory reserved by a given \texttt{Matrix\_double} and its member +\texttt{Array\_double} vectors describing its rows. +\end{itemize} + +\begin{verbatim} +void free_matrix(Matrix_double *m) { + for (size_t y = 0; y < m->rows; ++y) + free_vector(m->data[y]); + free(m); +} +\end{verbatim} + +\subsubsection{\texttt{format\_matrix\_into}} +\label{sec:org6b188b4} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{format\_matrix\_into} +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \texttt{Matrix\_double} and a pointer to a c-string \texttt{s} to "print" the vector out +into +\item Output: nothing. +\item Side effect: overwritten memory into \texttt{s} +\end{itemize} + +\begin{verbatim} +void format_matrix_into(Matrix_double *m, char *s) { + if (m->rows == 0) + strcpy(s, "empty"); + + for (size_t y = 0; y < m->rows; ++y) { + char row_s[5192]; + strcpy(row_s, ""); + + format_vector_into(m->data[y], row_s); + strcat(s, row_s); + } + strcat(s, "\n"); +} +\end{verbatim} +\subsection{Root Finding Methods} +\label{sec:org352ccdf} +\subsubsection{\texttt{find\_ivt\_range}} +\label{sec:orgb9a0d16} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{find\_ivt\_range} +\item Location: \texttt{src/roots.c} +\item Input: a pointer to a oneary function taking a double and producing a double, the beginning point +in \(R\) to search for a range, a \texttt{delta} step that is taken, and a \texttt{max\_steps} number of maximum +iterations to perform. +\item Output: a pair of \texttt{double}'s in an \texttt{Array\_double} representing a closed closed interval \texttt{[beginning, end]} +\end{itemize} + +\begin{verbatim} +// f is well defined at start_x + delta*n for all n on the integer range [0, +// max_iterations] +Array_double *find_ivt_range(double (*f)(double), double start_x, double delta, + size_t max_iterations) { + double a = start_x; + + while (f(a) * f(a + delta) >= 0 && max_iterations > 0) { + max_iterations--; + a += delta; + } + + double end = a + delta; + double begin = a - delta; + + if (max_iterations == 0 && f(begin) * f(end) >= 0) + return NULL; + return InitArray(double, {begin, end}); +} +\end{verbatim} +\subsubsection{\texttt{bisect\_find\_root}} +\label{sec:org25382b3} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name(s): \texttt{bisect\_find\_root} +\item Input: a one-ary function taking a double and producing a double, a closed interval represented +by \texttt{a} and \texttt{b}: \texttt{[a, b]}, a \texttt{tolerance} at which we return the estimated root once \(b-a < \text{tolerance}\), and a +\texttt{max\_iterations} to break us out of a loop if we can never reach the \texttt{tolerance}. +\item Output: a vector of size of 3, \texttt{double}'s representing first the range \texttt{[a,b]} and then the midpoint, +\texttt{c} of the range. +\item Description: recursively uses binary search to split the interval until we reach \texttt{tolerance}. We +also assume the function \texttt{f} is continuous on \texttt{[a, b]}. +\end{itemize} + +\begin{verbatim} +// f is continuous on [a, b] +Array_double *bisect_find_root(double (*f)(double), double a, double b, + double tolerance, size_t max_iterations) { + assert(a <= b); + // guarantee there's a root somewhere between a and b by IVT + assert(f(a) * f(b) < 0); + + double c = (1.0 / 2) * (a + b); + if (b - a < tolerance || max_iterations == 0) + return InitArray(double, {a, b, c}); + + if (f(a) * f(c) < 0) + return bisect_find_root(f, a, c, tolerance, max_iterations - 1); + return bisect_find_root(f, c, b, tolerance, max_iterations - 1); +} +\end{verbatim} +\subsubsection{\texttt{bisect\_find\_root\_with\_error\_assumption}} +\label{sec:org4b9cb72} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{bisect\_find\_root\_with\_error\_assumption} +\item Input: a one-ary function taking a double and producing a double, a closed interval represented +by \texttt{a} and \texttt{b}: \texttt{[a, b]}, and a \texttt{tolerance} equivalent to the above definition in \texttt{bisect\_find\_root} +\item Output: a \texttt{double} representing the estimated root +\item Description: using the bisection method we know that \(e_k \le (\frac{1}{2})^k (b_0 - a_0)\). So we can +calculate \(k\) at the worst possible case (that the error is exactly the tolerance) to be +\(\frac{log(tolerance) - log(b_0 - a_0)}{log(\frac{1}{2})}\). We pass this value into the \texttt{max\_iterations} +of \texttt{bisect\_find\_root} as above. +\end{itemize} +\begin{verbatim} +double bisect_find_root_with_error_assumption(double (*f)(double), double a, + double b, double tolerance) { + assert(a <= b); + + uint64_t max_iterations = + (uint64_t)ceil((log(tolerance) - log(b - a)) / log(1 / 2.0)); + + Array_double *a_b_root = bisect_find_root(f, a, b, tolerance, max_iterations); + double root = a_b_root->data[2]; + free_vector(a_b_root); + + return root; +} +\end{verbatim} + +\subsubsection{\texttt{fixed\_point\_iteration\_method}} +\label{sec:org4cee2bd} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{fixed\_point\_iteration\_method} +\item Location: \texttt{src/roots.c} +\item Input: a pointer to a oneary function \(f\) taking a double and producing a double of which we are +trying to find a root, a guess \(x_0\), and a function \(g\) of the same signature of \(f\) at which we +"step" our guesses according to the fixed point iteration method: \(x_k = g(x_{k-1})\). Additionally, a +\texttt{max\_iterations} representing the maximum number of "steps" to take before arriving at our +approximation and a \texttt{tolerance} to return our root if it becomes within [0 - tolerance, 0 + tolerance]. +\item Assumptions: \(g(x)\) must be a function such that at the point \(x^*\) (the found root) the derivative +\(|g'(x^*)| \lt 1\) +\item Output: a double representing the found approximate root \(\approx x^*\). +\end{itemize} + +\begin{verbatim} +double fixed_point_iteration_method(double (*f)(double), double (*g)(double), + double x_0, double tolerance, + size_t max_iterations) { + if (max_iterations <= 0) + return x_0; + + double root = g(x_0); + if (tolerance >= fabs(f(root))) + return root; + + return fixed_point_iteration_method(f, g, root, tolerance, + max_iterations - 1); +} +\end{verbatim} + +\subsubsection{\texttt{fixed\_point\_newton\_method}} +\label{sec:org93e3999} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{fixed\_point\_newton\_method} +\item Location: \texttt{src/roots.c} +\item Input: a pointer to a oneary function \(f\) taking a double and producing a double of which we are +trying to find a root and another pointer to a function fprime of the same signature, a guess \(x_0\), +and a \texttt{max\_iterations} and \texttt{tolerance} as defined in the above method are required inputs. +\item Description: continually computes elements in the sequence \(x_n = x_{n-1} - \frac{f(x_{n-1})}{f'p(x_{n-1})}\) +\item Output: a double representing the found approximate root \(\approx x^*\) recursively applied to the sequence +given +\end{itemize} +\begin{verbatim} +double fixed_point_newton_method(double (*f)(double), double (*fprime)(double), + double x_0, double tolerance, + size_t max_iterations) { + if (max_iterations <= 0) + return x_0; + + double root = x_0 - f(x_0) / fprime(x_0); + if (tolerance >= fabs(f(root))) + return root; + + return fixed_point_newton_method(f, fprime, root, tolerance, + max_iterations - 1); +} +\end{verbatim} + +\subsubsection{\texttt{fixed\_point\_secant\_method}} +\label{sec:orgf3f0711} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{fixed\_point\_secant\_method} +\item Location: \texttt{src/roots.c} +\item Input: a pointer to a oneary function \(f\) taking a double and producing a double of which we are +trying to find a root, a guess \(x_0\) and \(x_1\) in which a root lies between \([x_0, x_1]\); applying the +sequence \(x_n = x_{n-1} - f(x_{n-1}) \frac{x_{n-1} - x_{n-2}}{f(x_{n-1}) - f(x_{n-2})}\). +Additionally, a \texttt{max\_iterations} and \texttt{tolerance} as defined in the above method are required +inputs. +\item Output: a double representing the found approximate root \(\approx x^*\) recursively applied to the sequence. +\end{itemize} +\begin{verbatim} +double fixed_point_secant_method(double (*f)(double), double x_0, double x_1, + double tolerance, size_t max_iterations) { + if (max_iterations == 0) + return x_1; + + double root = x_1 - f(x_1) * ((x_1 - x_0) / (f(x_1) - f(x_0))); + + if (tolerance >= fabs(f(root))) + return root; + + return fixed_point_secant_method(f, x_1, root, tolerance, max_iterations - 1); +} +\end{verbatim} +\subsubsection{\texttt{fixed\_point\_secant\_bisection\_method}} +\label{sec:orgeaef048} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{fixed\_point\_secant\_method} +\item Location: \texttt{src/roots.c} +\item Input: a pointer to a oneary function \(f\) taking a double and producing a double of which we are +trying to find a root, a guess \(x_0\), and a \(x_1\) of which we define our first interval \([x_0, x_1]\). +Then, we perform a single iteration of the \texttt{fixed\_point\_secant\_method} on this interval; if it +produces a root outside, we refresh the interval and root respectively with the given +\texttt{bisect\_find\_root} method. Additionally, a \texttt{max\_iterations} and \texttt{tolerance} as defined in the above method are required +inputs. +\item Output: a double representing the found approximate root \(\approx x^*\) continually applied with the +constraints defined. +\end{itemize} + +\begin{verbatim} +double fixed_point_secant_bisection_method(double (*f)(double), double x_0, + double x_1, double tolerance, + size_t max_iterations) { + double begin = x_0; + double end = x_1; + double root = x_0; + + while (tolerance < fabs(f(root)) && max_iterations > 0) { + max_iterations--; + + double secant_root = fixed_point_secant_method(f, begin, end, tolerance, 1); + + if (secant_root < begin || secant_root > end) { + Array_double *range_root = bisect_find_root(f, begin, end, tolerance, 1); + + begin = range_root->data[0]; + end = range_root->data[1]; + root = range_root->data[2]; + + free_vector(range_root); + continue; + } + + root = secant_root; + + if (f(root) * f(begin) < 0) + end = secant_root; // the root exists in [begin, secant_root] + else + begin = secant_root; + } + + return root; +} +\end{verbatim} + +\subsection{Linear Routines} +\label{sec:orge3b6d97} +\subsubsection{\texttt{least\_squares\_lin\_reg}} +\label{sec:orgcc90c4a} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{least\_squares\_lin\_reg} +\item Location: \texttt{src/lin.c} +\item Input: two pointers to \texttt{Array\_double}'s whose entries correspond two ordered +pairs in R\textsuperscript{2} +\item Output: a linear model best representing the ordered pairs via least squares +regression +\end{itemize} + +\begin{verbatim} +Line *least_squares_lin_reg(Array_double *x, Array_double *y) { + assert(x->size == y->size); + + uint64_t n = x->size; + double sum_x = sum_v(x); + double sum_y = sum_v(y); + double sum_xy = v_dot_v(x, y); + double sum_xx = v_dot_v(x, x); + double denom = ((n * sum_xx) - (sum_x * sum_x)); + + Line *line = malloc(sizeof(Line)); + line->m = ((sum_xy * n) - (sum_x * sum_y)) / denom; + line->a = ((sum_y * sum_xx) - (sum_x * sum_xy)) / denom; + + return line; +} +\end{verbatim} + +\subsection{Eigen-Adjacent} +\label{sec:orga3c637f} +\subsubsection{\texttt{dominant\_eigenvalue}} +\label{sec:org0306c8a} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{dominant\_eigenvalue} +\item Location: \texttt{src/eigen.c} +\item Input: a pointer to an invertible matrix \texttt{m}, an initial eigenvector guess \texttt{v} (that is non +zero or orthogonal to an eigenvector with the dominant eigenvalue), a \texttt{tolerance} and +\texttt{max\_iterations} that act as stop conditions +\item Output: the dominant eigenvalue with the highest magnitude, approximated with the Power +Iteration Method +\end{itemize} + +\begin{verbatim} +double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance, + size_t max_iterations) { + assert(m->rows == m->cols); + assert(m->rows == v->size); + + double error = tolerance; + size_t iter = max_iterations; + double lambda = 0.0; + Array_double *eigenvector_1 = copy_vector(v); + + while (error >= tolerance && (--iter) > 0) { + Array_double *eigenvector_2 = m_dot_v(m, eigenvector_1); + Array_double *normalized_eigenvector_2 = + scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2)); + free_vector(eigenvector_2); + eigenvector_2 = normalized_eigenvector_2; + + Array_double *mx = m_dot_v(m, eigenvector_2); + double new_lambda = + v_dot_v(mx, eigenvector_2) / v_dot_v(eigenvector_2, eigenvector_2); + + error = fabs(new_lambda - lambda); + lambda = new_lambda; + free_vector(eigenvector_1); + eigenvector_1 = eigenvector_2; + } + + return lambda; +} +\end{verbatim} +\subsubsection{\texttt{shift\_inverse\_power\_eigenvalue}} +\label{sec:orgc29637a} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{least\_dominant\_eigenvalue} +\item Location: \texttt{src/eigen.c} +\item Input: a pointer to an invertible matrix \texttt{m}, an initial eigenvector guess \texttt{v} (that is non +zero or orthogonal to an eigenvector with the dominant eigenvalue), a \texttt{shift} to act as the +shifted \(\delta\), and \texttt{tolerance} and \texttt{max\_iterations} that act as stop conditions. +\item Output: the eigenvalue closest to \texttt{shift} with the lowest magnitude closest to 0, approximated +with the Inverse Power Iteration Method +\end{itemize} +\begin{verbatim} +double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v, + double shift, double tolerance, + size_t max_iterations) { + assert(m->rows == m->cols); + assert(m->rows == v->size); + + Matrix_double *m_c = copy_matrix(m); + for (size_t y = 0; y < m_c->rows; ++y) + m_c->data[y]->data[y] = m_c->data[y]->data[y] - shift; + + double error = tolerance; + size_t iter = max_iterations; + double lambda = shift; + Array_double *eigenvector_1 = copy_vector(v); + + while (error >= tolerance && (--iter) > 0) { + Array_double *eigenvector_2 = solve_matrix_lu_bsubst(m_c, eigenvector_1); + Array_double *normalized_eigenvector_2 = + scale_v(eigenvector_2, 1.0 / linf_norm(eigenvector_2)); + free_vector(eigenvector_2); + + Array_double *mx = m_dot_v(m, normalized_eigenvector_2); + double new_lambda = + v_dot_v(mx, normalized_eigenvector_2) / + v_dot_v(normalized_eigenvector_2, normalized_eigenvector_2); + + error = fabs(new_lambda - lambda); + lambda = new_lambda; + free_vector(eigenvector_1); + eigenvector_1 = normalized_eigenvector_2; + } + + return lambda; +} +\end{verbatim} + +\subsubsection{\texttt{least\_dominant\_eigenvalue}} +\label{sec:org5df73a2} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{least\_dominant\_eigenvalue} +\item Location: \texttt{src/eigen.c} +\item Input: a pointer to an invertible matrix \texttt{m}, an initial eigenvector guess \texttt{v} (that is non +zero or orthogonal to an eigenvector with the dominant eigenvalue), a \texttt{tolerance} and +\texttt{max\_iterations} that act as stop conditions. +\item Output: the least dominant eigenvalue with the lowest magnitude closest to 0, approximated +with the Inverse Power Iteration Method. +\end{itemize} +\begin{verbatim} +double least_dominant_eigenvalue(Matrix_double *m, Array_double *v, + double tolerance, size_t max_iterations) { + return shift_inverse_power_eigenvalue(m, v, 0.0, tolerance, max_iterations); +} +\end{verbatim} +\subsubsection{\texttt{partition\_find\_eigenvalues}} +\label{sec:org3dde7af} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{partition\_find\_eigenvalues} +\item Location: \texttt{src/eigen.c} +\item Input: a pointer to an invertible matrix \texttt{m}, a matrix whose rows correspond to initial +eigenvector guesses at each "partition" which is computed from a uniform distribution +between the number of rows this "guess matrix" has and the distance between the least +dominant eigenvalue and the most dominant. Additionally, a \texttt{max\_iterations} and a \texttt{tolerance} +that act as stop conditions. +\item Output: a vector of \texttt{doubles} corresponding to the "nearest" eigenvalue at the midpoint of +each partition, via the given guess of that partition. +\end{itemize} +\begin{verbatim} +Array_double *partition_find_eigenvalues(Matrix_double *m, + Matrix_double *guesses, + double tolerance, + size_t max_iterations) { + assert(guesses->rows >= + 2); // we need at least, the most and least dominant eigenvalues + + double end = dominant_eigenvalue(m, guesses->data[guesses->rows - 1], + tolerance, max_iterations); + double begin = + least_dominant_eigenvalue(m, guesses->data[0], tolerance, max_iterations); + + double delta = (end - begin) / guesses->rows; + Array_double *eigenvalues = InitArrayWithSize(double, guesses->rows, 0.0); + for (size_t i = 0; i < guesses->rows; i++) { + double box_midpoint = ((delta * i) + (delta * (i + 1))) / 2; + + double nearest_eigenvalue = shift_inverse_power_eigenvalue( + m, guesses->data[i], box_midpoint, tolerance, max_iterations); + + eigenvalues->data[i] = nearest_eigenvalue; + } + + return eigenvalues; +} +\end{verbatim} +\subsubsection{\texttt{leslie\_matrix}} +\label{sec:orgca10ed3} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{leslie\_matrix} +\item Location: \texttt{src/eigen.c} +\item Input: two pointers to \texttt{Array\_double}'s representing the ratio of individuals in an age class +\(x\) getting to the next age class \(x+1\) and the number of offspring that individuals in an age +class create in age class 0. +\item Output: the leslie matrix generated from the input vectors. +\end{itemize} + +\begin{verbatim} +Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio, + Array_double *age_class_offspring) { + assert(age_class_surivor_ratio->size + 1 == age_class_offspring->size); + + Matrix_double *leslie = InitMatrixWithSize(double, age_class_offspring->size, + age_class_offspring->size, 0.0); + + free_vector(leslie->data[0]); + leslie->data[0] = age_class_offspring; + + for (size_t i = 0; i < age_class_surivor_ratio->size; i++) + leslie->data[i + 1]->data[i] = age_class_surivor_ratio->data[i]; + return leslie; +} +\end{verbatim} +\subsection{Jacobi / Gauss-Siedel} +\label{sec:org91c563c} +\subsubsection{\texttt{jacobi\_solve}} +\label{sec:org2cd6098} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{jacobi\_solve} +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a diagonally dominant square matrix \(m\), a vector representing +the value \(b\) in \(mx = b\), a double representing the maximum distance between +the solutions produced by iteration \(i\) and \(i+1\) (by L2 norm a.k.a cartesian +distance), and a \texttt{max\_iterations} which we force stop. +\item Output: the converged-upon solution \(x\) to \(mx = b\) +\end{itemize} +\begin{verbatim} +Array_double *jacobi_solve(Matrix_double *m, Array_double *b, + double l2_convergence_tolerance, + size_t max_iterations) { + assert(m->rows == m->cols); + assert(b->size == m->cols); + size_t iter = max_iterations; + + Array_double *x_k = InitArrayWithSize(double, b->size, 0.0); + Array_double *x_k_1 = + InitArrayWithSize(double, b->size, rand_from(0.1, 10.0)); + + while ((--iter) > 0 && l2_distance(x_k_1, x_k) > l2_convergence_tolerance) { + for (size_t i = 0; i < m->rows; i++) { + double delta = 0.0; + for (size_t j = 0; j < m->cols; j++) { + if (i == j) + continue; + delta += m->data[i]->data[j] * x_k->data[j]; + } + x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i]; + } + + Array_double *tmp = x_k; + x_k = x_k_1; + x_k_1 = tmp; + } + + free_vector(x_k); + return x_k_1; +} +\end{verbatim} + +\subsubsection{\texttt{gauss\_siedel\_solve}} +\label{sec:org6633923} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{gauss\_siedel\_solve} +\item Location: \texttt{src/matrix.c} +\item Input: a pointer to a \href{https://en.wikipedia.org/wiki/Gauss\%E2\%80\%93Seidel\_method}{diagonally dominant or symmetric and positive definite} +square matrix \(m\), a vector representing +the value \(b\) in \(mx = b\), a double representing the maximum distance between +the solutions produced by iteration \(i\) and \(i+1\) (by L2 norm a.k.a cartesian +distance), and a \texttt{max\_iterations} which we force stop. +\item Output: the converged-upon solution \(x\) to \(mx = b\) +\item Description: we use almost the exact same method as \texttt{jacobi\_solve} but modify +only one array in accordance to the Gauss-Siedel method, but which is necessarily +copied before due to the convergence check. +\end{itemize} +\begin{verbatim} + +Array_double *gauss_siedel_solve(Matrix_double *m, Array_double *b, + double l2_convergence_tolerance, + size_t max_iterations) { + assert(m->rows == m->cols); + assert(b->size == m->cols); + size_t iter = max_iterations; + + Array_double *x_k = InitArrayWithSize(double, b->size, 0.0); + Array_double *x_k_1 = + InitArrayWithSize(double, b->size, rand_from(0.1, 10.0)); + + while ((--iter) > 0) { + for (size_t i = 0; i < x_k->size; i++) + x_k->data[i] = x_k_1->data[i]; + + for (size_t i = 0; i < m->rows; i++) { + double delta = 0.0; + for (size_t j = 0; j < m->cols; j++) { + if (i == j) + continue; + delta += m->data[i]->data[j] * x_k_1->data[j]; + } + x_k_1->data[i] = (b->data[i] - delta) / m->data[i]->data[i]; + } + + if (l2_distance(x_k_1, x_k) <= l2_convergence_tolerance) + break; + } + + free_vector(x_k); + return x_k_1; +} +\end{verbatim} + +\subsection{Appendix / Miscellaneous} +\label{sec:orga72494e} +\subsubsection{Random} +\label{sec:org4940c39} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Name: \texttt{rand\_from} +\item Location: \texttt{src/rand.c} +\item Input: a pair of doubles, min and max to generate a random number min +\(\le\) x \(\le\) max +\item Output: a random double in the constraints shown +\end{itemize} + +\begin{verbatim} +double rand_from(double min, double max) { + return min + (rand() / (RAND_MAX / (max - min))); +} +\end{verbatim} +\subsubsection{Data Types} +\label{sec:org8d3f6e1} +\begin{enumerate} +\item \texttt{Line} +\label{sec:orgc0df901} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{inc/types.h} +\end{itemize} + +\begin{verbatim} +typedef struct Line { + double m; + double a; +} Line; +\end{verbatim} +\item The \texttt{Array\_} and \texttt{Matrix\_} +\label{sec:org435e816} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{inc/types.h} +\end{itemize} + +We define two Pre processor Macros \texttt{DEFINE\_ARRAY} and \texttt{DEFINE\_MATRIX} that take +as input a type, and construct a struct definition for the given type for +convenient access to the vector or matrices dimensions. + +Such that \texttt{DEFINE\_ARRAY(int)} would expand to: + +\begin{verbatim} +typedef struct { + int* data; + size_t size; +} Array_int +\end{verbatim} + +And \texttt{DEFINE\_MATRIX(int)} would expand a to \texttt{Matrix\_int}; containing a pointer to +a collection of pointers of \texttt{Array\_int}'s and its dimensions. + +\begin{verbatim} +typedef struct { + Array_int **data; + size_t cols; + size_t rows; +} Matrix_int +\end{verbatim} +\end{enumerate} + +\subsubsection{Macros} +\label{sec:orga2161be} +\begin{enumerate} +\item \texttt{c\_max} and \texttt{c\_min} +\label{sec:org16ca9c3} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{inc/macros.h} +\item Input: two structures that define an order measure +\item Output: either the larger or smaller of the two depending on the measure +\end{itemize} + +\begin{verbatim} +#define c_max(x, y) (((x) >= (y)) ? (x) : (y)) +#define c_min(x, y) (((x) <= (y)) ? (x) : (y)) +\end{verbatim} + +\item \texttt{InitArray} +\label{sec:orgcaff993} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{inc/macros.h} +\item Input: a type and array of values to initialze an array with such type +\item Output: a new \texttt{Array\_type} with the size of the given array and its data +\end{itemize} + +\begin{verbatim} +#define InitArray(TYPE, ...) \ + ({ \ + TYPE temp[] = __VA_ARGS__; \ + Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ + arr->size = sizeof(temp) / sizeof(temp[0]); \ + arr->data = malloc(arr->size * sizeof(TYPE)); \ + memcpy(arr->data, temp, arr->size * sizeof(TYPE)); \ + arr; \ + }) +\end{verbatim} + +\item \texttt{InitArrayWithSize} +\label{sec:orga925ddb} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{inc/macros.h} +\item Input: a type, a size, and initial value +\item Output: a new \texttt{Array\_type} with the given size filled with the initial value +\end{itemize} + +\begin{verbatim} +#define InitArrayWithSize(TYPE, SIZE, INIT_VALUE) \ + ({ \ + Array_##TYPE *arr = malloc(sizeof(Array_##TYPE)); \ + arr->size = SIZE; \ + arr->data = malloc(arr->size * sizeof(TYPE)); \ + for (size_t i = 0; i < arr->size; i++) \ + arr->data[i] = INIT_VALUE; \ + arr; \ + }) +\end{verbatim} + +\item \texttt{InitMatrixWithSize} +\label{sec:orgf90d7c8} +\begin{itemize} +\item Author: Elizabeth Hunt +\item Location: \texttt{inc/macros.h} +\item Input: a type, number of rows, columns, and initial value +\item Output: a new \texttt{Matrix\_type} of size \texttt{rows x columns} filled with the initial +value +\end{itemize} + +\begin{verbatim} +#define InitMatrixWithSize(TYPE, ROWS, COLS, INIT_VALUE) \ + ({ \ + Matrix_##TYPE *matrix = malloc(sizeof(Matrix_##TYPE)); \ + matrix->rows = ROWS; \ + matrix->cols = COLS; \ + matrix->data = malloc(matrix->rows * sizeof(Array_##TYPE *)); \ + for (size_t y = 0; y < matrix->rows; y++) \ + matrix->data[y] = InitArrayWithSize(TYPE, COLS, INIT_VALUE); \ + matrix; \ + }) +\end{verbatim} +\end{enumerate} +\end{document} \ No newline at end of file -- cgit v1.3