diff options
| -rw-r--r-- | doc/software_manual.org | 40 | ||||
| -rw-r--r-- | homeworks/hw-7.org | 18 | ||||
| -rw-r--r-- | inc/lizfcm.h | 3 | ||||
| -rw-r--r-- | src/eigen.c | 19 | ||||
| -rw-r--r-- | test/eigen.t.c | 114 |
5 files changed, 128 insertions, 66 deletions
diff --git a/doc/software_manual.org b/doc/software_manual.org index 1520431..d6c7331 100644 --- a/doc/software_manual.org +++ b/doc/software_manual.org @@ -1035,22 +1035,22 @@ double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance, return lambda; } #+END_SRC -*** ~least_dominant_eigenvalue~ +*** ~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 ~tolerance~ and - ~max_iterations~ that act as stop conditions -+ Output: the least dominant eigenvalue with the lowest magnitude, approximated with the Inverse - Power Iteration Method + 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 least_dominant_eigenvalue(Matrix_double *m, Array_double *v, - double tolerance, size_t max_iterations) { +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); - double shift = 0.0; 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; @@ -1065,22 +1065,38 @@ double least_dominant_eigenvalue(Matrix_double *m, Array_double *v, 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); + Array_double *mx = m_dot_v(m, normalized_eigenvector_2); double new_lambda = - v_dot_v(mx, eigenvector_2) / v_dot_v(eigenvector_2, eigenvector_2); + 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 = eigenvector_2; + 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 + *** ~leslie_matrix~ + Author: Elizabeth Hunt + Name: ~leslie_matrix~ diff --git a/homeworks/hw-7.org b/homeworks/hw-7.org index dbdf6bb..ec8c23d 100644 --- a/homeworks/hw-7.org +++ b/homeworks/hw-7.org @@ -21,7 +21,23 @@ finds the least dominant eigenvalue on the matrix: 0 & 2 & 6 \end{bmatrix} -which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$ and should produce $\sqrt{17}$. +which has eigenvalues: $5 + \sqrt{17}, 2, 5 - \sqrt{17}$ and should thus produce $5 - \sqrt{17}$. See also the entry ~Eigen-Adjacent -> least_dominant_eigenvalue~ in the LIZFCM API documentation. +* Question Four +See ~UTEST(eigen, shifted_eigenvalue)~ in ~test/eigen.t.c~ which +finds the least dominant eigenvalue 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 thus produce $2.0$. + +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. diff --git a/inc/lizfcm.h b/inc/lizfcm.h index 295aab0..625e6bc 100644 --- a/inc/lizfcm.h +++ b/inc/lizfcm.h @@ -76,6 +76,9 @@ extern double fixed_point_secant_bisection_method(double (*f)(double), extern double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance, size_t max_iterations); +extern double shift_inverse_power_eigenvalue(Matrix_double *m, Array_double *v, + double shift, double tolerance, + size_t max_iterations); extern double least_dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance, size_t max_iterations); diff --git a/src/eigen.c b/src/eigen.c index 8fcf5c4..c4af461 100644 --- a/src/eigen.c +++ b/src/eigen.c @@ -49,12 +49,12 @@ double dominant_eigenvalue(Matrix_double *m, Array_double *v, double tolerance, return lambda; } -double least_dominant_eigenvalue(Matrix_double *m, Array_double *v, - double tolerance, size_t max_iterations) { +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); - double shift = 0.0; 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; @@ -69,17 +69,22 @@ double least_dominant_eigenvalue(Matrix_double *m, Array_double *v, 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); + Array_double *mx = m_dot_v(m, normalized_eigenvector_2); double new_lambda = - v_dot_v(mx, eigenvector_2) / v_dot_v(eigenvector_2, eigenvector_2); + 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 = eigenvector_2; + eigenvector_1 = normalized_eigenvector_2; } return lambda; } + +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 0ad0bd0..5a20d86 100644 --- a/test/eigen.t.c +++ b/test/eigen.t.c @@ -1,50 +1,7 @@ #include "lizfcm.test.h" -UTEST(eigen, leslie_matrix) { - Array_double *felicity = InitArray(double, {0.0, 1.5, 0.8}); - Array_double *survivor_ratios = InitArray(double, {0.8, 0.55}); - - Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0); - m->data[0]->data[0] = 0.0; - m->data[0]->data[1] = 1.5; - m->data[0]->data[2] = 0.8; - m->data[1]->data[0] = 0.8; - m->data[2]->data[1] = 0.55; - - Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity); - - EXPECT_TRUE(matrix_equal(leslie, m)); - - free_matrix(leslie); - free_matrix(m); - free_vector(felicity); - free_vector(survivor_ratios); -} - -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}); - Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity); - Array_double *v_guess = InitArrayWithSize(double, 3, 2.0); - double tolerance = 0.0001; - uint64_t max_iterations = 64; - - double expect_dominant_eigenvalue = 1.22005; - - double approx_dominant_eigenvalue = - dominant_eigenvalue(leslie, v_guess, tolerance, max_iterations); - - EXPECT_NEAR(expect_dominant_eigenvalue, approx_dominant_eigenvalue, - tolerance); - - free_vector(v_guess); - free_vector(survivor_ratios); - free_vector(felicity); - free_matrix(leslie); -} - -UTEST(eigen, least_dominant_eigenvalue) { - +Matrix_double *eigen_test_matrix() { + // produces a matrix that has eigenvalues [5 + sqrt{17}, 2, 5 - sqrt{17}] Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0); m->data[0]->data[0] = 2.0; m->data[0]->data[1] = 2.0; @@ -54,12 +11,17 @@ UTEST(eigen, least_dominant_eigenvalue) { m->data[1]->data[2] = 7.0; m->data[2]->data[1] = 2.0; m->data[2]->data[2] = 6.0; + return m; +} + +UTEST(eigen, least_dominant_eigenvalue) { + Matrix_double *m = eigen_test_matrix(); double expected_least_dominant_eigenvalue = 0.87689; // 5 - sqrt(17) double tolerance = 0.0001; uint64_t max_iterations = 64; - Array_double *v_guess = InitArrayWithSize(double, 3, 2.0); + Array_double *v_guess = InitArrayWithSize(double, 3, 1.0); double approx_least_dominant_eigenvalue = least_dominant_eigenvalue(m, v_guess, tolerance, max_iterations); @@ -88,3 +50,63 @@ UTEST(eigen, dominant_eigenvalue) { free_matrix(m); free_vector(v_guess); } + +UTEST(eigen, shifted_eigenvalue) { + 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 shift = (dominant_eigenvalue + least_dominant_eigenvalue) / 2.0; + + double tolerance = 0.0001; + uint64_t max_iterations = 64; + Array_double *v_guess = InitArray(double, {0.5, 1.0, 0.75}); + + double approx_middle_eigenvalue = shift_inverse_power_eigenvalue( + m, v_guess, shift, tolerance, max_iterations); + + EXPECT_NEAR(approx_middle_eigenvalue, expected_middle_eigenvalue, tolerance); +} + +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}); + Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity); + Array_double *v_guess = InitArrayWithSize(double, 3, 2.0); + double tolerance = 0.0001; + uint64_t max_iterations = 64; + + double expect_dominant_eigenvalue = 1.22005; + + double approx_dominant_eigenvalue = + dominant_eigenvalue(leslie, v_guess, tolerance, max_iterations); + + EXPECT_NEAR(expect_dominant_eigenvalue, approx_dominant_eigenvalue, + tolerance); + + free_vector(v_guess); + free_vector(survivor_ratios); + free_vector(felicity); + free_matrix(leslie); +} +UTEST(eigen, leslie_matrix) { + Array_double *felicity = InitArray(double, {0.0, 1.5, 0.8}); + Array_double *survivor_ratios = InitArray(double, {0.8, 0.55}); + + Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0); + m->data[0]->data[0] = 0.0; + m->data[0]->data[1] = 1.5; + m->data[0]->data[2] = 0.8; + m->data[1]->data[0] = 0.8; + m->data[2]->data[1] = 0.55; + + Matrix_double *leslie = leslie_matrix(survivor_ratios, felicity); + + EXPECT_TRUE(matrix_equal(leslie, m)); + + free_matrix(leslie); + free_matrix(m); + free_vector(felicity); + free_vector(survivor_ratios); +} |