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 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