hw8 checkpoint

7 files changed, 113 insertions, 4 deletions

diff --git a/homeworks/hw-8.org b/homeworks/hw-8.org
new file mode 100644
index 0000000..42d9dac
--- /dev/null
+++ b/homeworks/hw-8.org
@@ -0,0 +1,17 @@
+#+TITLE: Homework 7
+#+AUTHOR: Elizabeth Hunt
+#+LATEX_HEADER: \notindent \notag  \usepackage{amsmath} \usepackage[a4paper,margin=1in,portrait]{geometry}
+#+LATEX: \setlength\parindent{0pt}
+#+OPTIONS: toc:nil
+
+TODO: Update LIZFCM org file with jacobi solve, format_matrix_into, rand
+
+* Question One
+See ~UTEST(jacobi, solve_jacobi)~ in ~test/jacobi.t.c~ and the entry
+~Jacobi -> solve_jacobi~ in the LIZFCM API documentation.
+* Question Two
+A problem arises when using the Jacobi method to solve for the previous population
+distribution, $n_k$, from $Ln_{k} = n_{k+1}$, because a Leslie matrix is not diagonally
+dominant and will cause a division by zero. Likewise, we cannot factor it into $L$
+and $U$ terms and apply back substitution because pivot points are zero.
+* Question Three
diff --git a/inc/lizfcm.h b/inc/lizfcm.h
index 1bb5322..70ebd2f 100644
--- a/inc/lizfcm.h
+++ b/inc/lizfcm.h
@@ -88,4 +88,10 @@ extern Array_double *partition_find_eigenvalues(Matrix_double *m,
                                                 size_t max_iterations);
 extern Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio,
                                     Array_double *age_class_offspring);
+
+extern double rand_from(double min, double max);
+
+extern Array_double *jacobi_solve(Matrix_double *m, Array_double *b,
+                                  double tolerance, size_t max_iterations);
+
 #endif // LIZFCM_H
diff --git a/src/matrix.c b/src/matrix.c
index 04c5adc..5f36d12 100644
--- a/src/matrix.c
+++ b/src/matrix.c
@@ -2,9 +2,9 @@
 #include <assert.h>
 #include <math.h>
 #include <stdio.h>
-#include <string.h>
+n #include<string.h>
 
-Array_double *m_dot_v(Matrix_double *m, Array_double *v) {
+    Array_double *m_dot_v(Matrix_double *m, Array_double *v) {
   assert(v->size == m->cols);
 
   Array_double *product = copy_vector(v);
@@ -222,6 +222,35 @@ Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b) {
   return solution;
 }
 
+Array_double *jacobi_solve(Matrix_double *m, Array_double *b,
+                           double l2_convergence_tolerance,
+                           size_t max_iterations) {
+  size_t iter = max_iterations;
+
+  Array_double *x_k = InitArrayWithSize(double, b->size, rand_from(0.1, 10.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 < x_k->size; i++) {
+      double delta = 0.0;
+      for (size_t j = 0; j < x_k->size; 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;
+}
+
 Matrix_double *slice_column(Matrix_double *m, size_t x) {
   Matrix_double *sliced = copy_matrix(m);
 
@@ -259,7 +288,7 @@ void format_matrix_into(Matrix_double *m, char *s) {
     strcpy(s, "empty");
 
   for (size_t y = 0; y < m->rows; ++y) {
-    char row_s[256];
+    char row_s[5192];
     strcpy(row_s, "");
 
     format_vector_into(m->data[y], row_s);
diff --git a/src/rand.c b/src/rand.c
new file mode 100644
index 0000000..ff375ed
--- /dev/null
+++ b/src/rand.c
@@ -0,0 +1,7 @@
+#include "lizfcm.h"
+
+double rand_from(double min, double max) {
+  double range = (max - min);
+  double div = RAND_MAX / range;
+  return min + (rand() / div);
+}
diff --git a/test/jacobi.t.c b/test/jacobi.t.c
new file mode 100644
index 0000000..dc13d6e
--- /dev/null
+++ b/test/jacobi.t.c
@@ -0,0 +1,33 @@
+#include "lizfcm.test.h"
+#include <math.h>
+
+Matrix_double *generate_ddm(size_t n) {
+  Matrix_double *m = InitMatrixWithSize(double, n, n, rand_from(0.0, 1.0));
+
+  for (size_t y = 0; y < m->rows; y++) {
+    m->data[y]->data[y] += sum_v(m->data[y]);
+  }
+
+  return m;
+}
+
+UTEST(jacobi, jacobi_solve) {
+  Matrix_double *m = generate_ddm(2);
+
+  Array_double *b_1 = InitArrayWithSize(double, m->rows, 1.0);
+  Array_double *b = m_dot_v(m, b_1);
+
+  double tolerance = 0.001;
+  size_t max_iter = 400;
+  Array_double *solution = jacobi_solve(m, b, tolerance, max_iter);
+
+  for (size_t y = 0; y < m->rows; y++) {
+    double dot = v_dot_v(m->data[y], solution);
+    EXPECT_NEAR(b->data[y], dot, 0.1);
+  }
+
+  free_matrix(m);
+  free_vector(b_1);
+  free_vector(b);
+  free_vector(solution);
+}
diff --git a/test/main.c b/test/main.c
index 5d16fe7..6a92d4a 100644
--- a/test/main.c
+++ b/test/main.c
@@ -1,5 +1,12 @@
 #include "lizfcm.test.h"
+#include <stdlib.h>
+#include <time.h>
 
 UTEST(basic, unit_tests) { ASSERT_TRUE(1); }
 
-UTEST_MAIN();
+UTEST_STATE();
+int main(int argc, const char *const argv[]) {
+  srand(time(NULL));
+
+  return utest_main(argc, argv);
+}
diff --git a/test/rand.t.c b/test/rand.t.c
new file mode 100644
index 0000000..9ed2e9c
--- /dev/null
+++ b/test/rand.t.c
@@ -0,0 +1,10 @@
+#include "lizfcm.test.h"
+
+UTEST(rand, rand_from) {
+  double min = -2.0;
+  double max = 5.0;
+  for (size_t i = 0; i < 1000; i++) {
+    double r = rand_from(min, max);
+    ASSERT_TRUE(min <= r && r <= max);
+  }
+}