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 +#include +#include +#include + +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); +}