compute dominant eigenvalue

5 files changed, 83 insertions, 0 deletions

diff --git a/inc/lizfcm.h b/inc/lizfcm.h
index 2e12a50..85a29bb 100644
--- a/inc/lizfcm.h
+++ b/inc/lizfcm.h
@@ -40,6 +40,7 @@ extern Matrix_double *gaussian_elimination(Matrix_double *m);
 extern Array_double *solve_matrix_gaussian(Matrix_double *m, Array_double *b);
 extern Array_double *m_dot_v(Matrix_double *m, Array_double *v);
 extern Matrix_double *m_dot_m(Matrix_double *a, Matrix_double *b);
+extern Matrix_double *transpose(Matrix_double *m);
 extern Array_double *col_v(Matrix_double *m, size_t x);
 extern Matrix_double *copy_matrix(Matrix_double *m);
 extern Matrix_double *add_column(Matrix_double *m, Array_double *col);
@@ -72,4 +73,7 @@ extern double fixed_point_secant_bisection_method(double (*f)(double),
                                                   double x_0, double x_1,
                                                   double tolerance,
                                                   size_t max_iterations);
+
+extern double dominant_eigenvalue(Matrix_double *m, Array_double *v,
+                                  double tolerance, size_t max_iterations);
 #endif // LIZFCM_H
diff --git a/src/eigen.c b/src/eigen.c
new file mode 100644
index 0000000..aaacedc
--- /dev/null
+++ b/src/eigen.c
@@ -0,0 +1,31 @@
+#include "lizfcm.h"
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+double 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 error = tolerance;
+  size_t iter = max_iterations;
+  double lambda = 0.0;
+  Array_double *eigenvector_1 = copy_vector(v);
+
+  while (error >= tolerance && (--iter) > 0) {
+    Array_double *eigenvector_2 = m_dot_v(m, eigenvector_1);
+
+    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/src/matrix.c b/src/matrix.c
index 0891734..04c5adc 100644
--- a/src/matrix.c
+++ b/src/matrix.c
@@ -42,6 +42,16 @@ Matrix_double *m_dot_m(Matrix_double *a, Matrix_double *b) {
   return prod;
 }
 
+Matrix_double *transpose(Matrix_double *m) {
+  Matrix_double *transposed = InitMatrixWithSize(double, m->cols, m->rows, 0.0);
+
+  for (size_t x = 0; x < m->rows; x++)
+    for (size_t y = 0; y < m->cols; y++)
+      transposed->data[y]->data[x] = m->data[x]->data[y];
+
+  return transposed;
+}
+
 Matrix_double *put_identity_diagonal(Matrix_double *m) {
   assert(m->rows == m->cols);
   Matrix_double *copy = copy_matrix(m);
diff --git a/test/eigen.t.c b/test/eigen.t.c
new file mode 100644
index 0000000..985a304
--- /dev/null
+++ b/test/eigen.t.c
@@ -0,0 +1,23 @@
+#include "lizfcm.test.h"
+
+UTEST(eigen, dominant_eigenvalue) {
+  Matrix_double *m = InitMatrixWithSize(double, 2, 2, 0.0);
+  m->data[0]->data[0] = 2.0;
+  m->data[0]->data[1] = -12.0;
+  m->data[1]->data[0] = 1.0;
+  m->data[1]->data[1] = -5.0;
+
+  Array_double *v_guess = InitArrayWithSize(double, 2, 1.0);
+  double tolerance = 0.0001;
+  uint64_t max_iterations = 64;
+
+  double expect_dominant_eigenvalue = -2.0;
+
+  double approx_dominant_eigenvalue =
+      dominant_eigenvalue(m, v_guess, tolerance, max_iterations);
+
+  EXPECT_NEAR(expect_dominant_eigenvalue, approx_dominant_eigenvalue,
+              tolerance);
+  free_matrix(m);
+  free_vector(v_guess);
+}
diff --git a/test/matrix.t.c b/test/matrix.t.c
index 1c72b85..597d6e0 100644
--- a/test/matrix.t.c
+++ b/test/matrix.t.c
@@ -230,3 +230,18 @@ UTEST(matrix, m_dot_m) {
   free_matrix(b);
   free_matrix(prod);
 }
+
+UTEST(matrix, transpose) {
+  Matrix_double *a = InitMatrixWithSize(double, 1, 3, 12.0);
+  a->data[0]->data[1] = 13.0;
+  Matrix_double *b = InitMatrixWithSize(double, 3, 1, 12.0);
+  b->data[1]->data[0] = 13.0;
+
+  Matrix_double *a_t = transpose(a);
+
+  EXPECT_TRUE(matrix_equal(a_t, b));
+
+  free_matrix(a_t);
+  free_matrix(a);
+  free_matrix(b);
+}