add least dominant eigenvalue

3 files changed, 62 insertions, 0 deletions

diff --git a/inc/lizfcm.h b/inc/lizfcm.h
index 4cdba8d..295aab0 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 least_dominant_eigenvalue(Matrix_double *m, Array_double *v,
+                                        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 36ccc92..8fcf5c4 100644
--- a/src/eigen.c
+++ b/src/eigen.c
@@ -48,3 +48,38 @@ 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) {
+  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;
+
+  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);
+    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;
+}
diff --git a/test/eigen.t.c b/test/eigen.t.c
index f271bf2..0ad0bd0 100644
--- a/test/eigen.t.c
+++ b/test/eigen.t.c
@@ -43,6 +43,30 @@ UTEST(eigen, leslie_matrix_dominant_eigenvalue) {
   free_matrix(leslie);
 }
 
+UTEST(eigen, least_dominant_eigenvalue) {
+
+  Matrix_double *m = InitMatrixWithSize(double, 3, 3, 0.0);
+  m->data[0]->data[0] = 2.0;
+  m->data[0]->data[1] = 2.0;
+  m->data[0]->data[2] = 4.0;
+  m->data[1]->data[0] = 1.0;
+  m->data[1]->data[1] = 4.0;
+  m->data[1]->data[2] = 7.0;
+  m->data[2]->data[1] = 2.0;
+  m->data[2]->data[2] = 6.0;
+
+  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);
+  double approx_least_dominant_eigenvalue =
+      least_dominant_eigenvalue(m, v_guess, tolerance, max_iterations);
+
+  EXPECT_NEAR(expected_least_dominant_eigenvalue,
+              approx_least_dominant_eigenvalue, tolerance);
+}
+
 UTEST(eigen, dominant_eigenvalue) {
   Matrix_double *m = InitMatrixWithSize(double, 2, 2, 0.0);
   m->data[0]->data[0] = 2.0;