compute dominant eigenvalue

This commit is contained in:
Elizabeth Hunt 2023-11-15 14:16:15 -07:00
parent 3f1f18b149
commit 1a6b952736
Signed by: simponic
GPG Key ID: 52B3774857EB24B1
5 changed files with 83 additions and 0 deletions

View File

@ -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 *solve_matrix_gaussian(Matrix_double *m, Array_double *b);
extern Array_double *m_dot_v(Matrix_double *m, Array_double *v); 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 *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 Array_double *col_v(Matrix_double *m, size_t x);
extern Matrix_double *copy_matrix(Matrix_double *m); extern Matrix_double *copy_matrix(Matrix_double *m);
extern Matrix_double *add_column(Matrix_double *m, Array_double *col); 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 x_0, double x_1,
double tolerance, double tolerance,
size_t max_iterations); size_t max_iterations);
extern double dominant_eigenvalue(Matrix_double *m, Array_double *v,
double tolerance, size_t max_iterations);
#endif // LIZFCM_H #endif // LIZFCM_H

31
src/eigen.c Normal file
View File

@ -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;
}

View File

@ -42,6 +42,16 @@ Matrix_double *m_dot_m(Matrix_double *a, Matrix_double *b) {
return prod; 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) { Matrix_double *put_identity_diagonal(Matrix_double *m) {
assert(m->rows == m->cols); assert(m->rows == m->cols);
Matrix_double *copy = copy_matrix(m); Matrix_double *copy = copy_matrix(m);

23
test/eigen.t.c Normal file
View File

@ -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);
}

View File

@ -230,3 +230,18 @@ UTEST(matrix, m_dot_m) {
free_matrix(b); free_matrix(b);
free_matrix(prod); 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);
}