q5 hw7

6 files changed, 124 insertions, 3 deletions

diff --git a/doc/software_manual.org b/doc/software_manual.org
index d6c7331..e12032d 100644
--- a/doc/software_manual.org
+++ b/doc/software_manual.org
@@ -1096,7 +1096,44 @@ double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
   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~
diff --git a/homeworks/hw-7.org b/homeworks/hw-7.org
index ec8c23d..e18d1ee 100644
--- a/homeworks/hw-7.org
+++ b/homeworks/hw-7.org
@@ -41,3 +41,21 @@ With the initial guess: $[0.5, 1.0, 0.75]$.
 
 See also the entry ~Eigen-Adjacent -> shift_inverse_power_eigenvalue~ in the LIZFCM API
 documentation.
+* Question Five
+See ~UTEST(eigen, partition_find_eigenvalues)~ in ~test/eigen.t.c~ which
+finds the eigenvalues in a partition of 10 on the matrix:
+
+\begin{bmatrix}
+2 & 2 & 4 \\
+1 & 4 & 7 \\
+0 & 2 & 6 
+\end{bmatrix}
+
+which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$, and should produce all three from
+the partitions when given the guesses $[0.5, 1.0, 0.75]$ from the questions above.
+
+See also the entry ~Eigen-Adjacent -> partition_find_eigenvalues~ in the LIZFCM API
+documentation.
+
+* Question Six
+
diff --git a/inc/lizfcm.h b/inc/lizfcm.h
index 625e6bc..1bb5322 100644
--- a/inc/lizfcm.h
+++ b/inc/lizfcm.h
@@ -82,6 +82,10 @@ extern double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v,
 extern double least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
                                         double tolerance,
                                         size_t max_iterations);
+extern Array_double *partition_find_eigenvalues(Matrix_double *m,
+                                                Matrix_double *guesses,
+                                                double tolerance,
+                                                size_t max_iterations);
 extern Matrix_double *leslie_matrix(Array_double *age_class_surivor_ratio,
                                     Array_double *age_class_offspring);
 #endif // LIZFCM_H
diff --git a/src/eigen.c b/src/eigen.c
index c4af461..92ef88c 100644
--- a/src/eigen.c
+++ b/src/eigen.c
@@ -84,6 +84,32 @@ double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v,
   return lambda;
 }
 
+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;
+}
+
 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);
diff --git a/test/eigen.t.c b/test/eigen.t.c
index 5a20d86..dc01aa7 100644
--- a/test/eigen.t.c
+++ b/test/eigen.t.c
@@ -1,4 +1,5 @@
 #include "lizfcm.test.h"
+#include <math.h>
 
 Matrix_double *eigen_test_matrix() {
   // produces a matrix that has eigenvalues [5 + sqrt{17}, 2, 5 - sqrt{17}]
@@ -69,6 +70,40 @@ UTEST(eigen, shifted_eigenvalue) {
   EXPECT_NEAR(approx_middle_eigenvalue, expected_middle_eigenvalue, tolerance);
 }
 
+UTEST(eigen, partition_find_eigenvalues) {
+  Matrix_double *m = eigen_test_matrix();
+
+  double least_dominant_eigenvalue = 0.87689; // 5 - sqrt{17}
+  double dominant_eigenvalue = 9.12311;       // 5 + sqrt{17}
+  double expected_middle_eigenvalue = 2.0;
+  double expected_eigenvalues[3] = {least_dominant_eigenvalue,
+                                    expected_middle_eigenvalue,
+                                    dominant_eigenvalue};
+
+  size_t partitions = 10;
+  Matrix_double *guesses = InitMatrixWithSize(double, partitions, 3, 0.0);
+  for (size_t y = 0; y < guesses->rows; y++) {
+    free_vector(guesses->data[y]);
+    guesses->data[y] = InitArray(double, {0.5, 1.0, 0.75});
+  }
+
+  double tolerance = 0.0001;
+  uint64_t max_iterations = 64;
+
+  int eigenvalues_found[3] = {false, false, false};
+  Array_double *partition_eigenvalues =
+      partition_find_eigenvalues(m, guesses, tolerance, max_iterations);
+
+  for (size_t i = 0; i < partition_eigenvalues->size; i++)
+    for (size_t eigenvalue_i = 0; eigenvalue_i < 3; eigenvalue_i++)
+      if (fabs(partition_eigenvalues->data[i] - expected_eigenvalues[i]) <=
+          tolerance)
+        eigenvalues_found[eigenvalue_i] = true;
+
+  for (size_t eigenvalue_i = 0; eigenvalue_i < 3; eigenvalue_i++)
+    EXPECT_TRUE(eigenvalues_found[eigenvalue_i]);
+}
+
 UTEST(eigen, leslie_matrix_dominant_eigenvalue) {
   Array_double *felicity = InitArray(double, {0.0, 1.5, 0.8});
   Array_double *survivor_ratios = InitArray(double, {0.8, 0.55});
diff --git a/test/lizfcm.test.h b/test/lizfcm.test.h
index 2374b83..9819d46 100644
--- a/test/lizfcm.test.h
+++ b/test/lizfcm.test.h
@@ -1,7 +1,8 @@
-#include "lizfcm.h"
-#include "utest.h"
-
 #ifndef LIZFCM_TEST_H
 #define LIZFCM_TEST_H
 
+#include "lizfcm.h"
+
+#include "utest.h"
+
 #endif // LIZFCM_TEST_H