add shifted eigenvalue procedure and unit test

1 files changed, 28 insertions, 12 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~