From a4907f39ff04b6667fdc68c75157e46da8331de3 Mon Sep 17 00:00:00 2001 From: Francesco Cozzuto Date: Mon, 27 Mar 2023 23:08:03 +0200 Subject: [PATCH] added QR and LU decomposition --- build.sh | 4 +- src/lina.c | 150 +++++++++++++++++++++++++++++++++++++++++++++++++++++ src/lina.h | 21 +++++--- src/qr.c | 142 ++++++++++++++++++++++++++++++++++++++++++++++++++ test.c | 58 +++++++++++++++++++++ 5 files changed, 367 insertions(+), 8 deletions(-) create mode 100644 src/qr.c create mode 100644 test.c diff --git a/build.sh b/build.sh index 74e3087..4848e0e 100755 --- a/build.sh +++ b/build.sh @@ -1,2 +1,2 @@ -gcc tests/test.c src/lina.c -o test -Wall -Wextra -g -Isrc/ -gcc tests/test_loader.c src/lina.c -o test_loader -Wall -Wextra -g -Isrc/ \ No newline at end of file +gcc tests/test.c src/lina.c -o test -Wall -Wextra -g -Isrc/ -lm +gcc tests/test_loader.c src/lina.c -o test_loader -Wall -Wextra -g -Isrc/ -lm \ No newline at end of file diff --git a/src/lina.c b/src/lina.c index 314c402..7648e92 100644 --- a/src/lina.c +++ b/src/lina.c @@ -5,6 +5,7 @@ #include #include #include +#include #include "lina.h" /* Function: lina_dot @@ -677,3 +678,152 @@ void lina_conv(double *A, double *B, double *C, C[j * Cw + i] += A[(i - Bw/2 + u) * Aw + (i - Bh/2 + v)] * B[v * Bw + u]; } } + +void lina_decompLU(double *A, double *L, double *U, int n) +{ + assert(n > 0); + assert(A != L && A != U && L != U); + + // TODO: Handle the case when A can not be + // decomposed. + + memset(L, 0, sizeof(double) * n * n); + memset(U, 0, sizeof(double) * n * n); + + /* + // Zero-out the lower half of L and the upper + // half of U. + for (int i = 0; i < n; i++) + for (int j = i+1; j < n; j++) + { + L[j * n + i] = 0; + U[i * n + j] = 0; + } + */ + + for (int i = 0; i < n; i++) + { + for (int k = i; k < n; k++) + { + int sum = 0; // L[i,j] * U[j,k] + for (int j = 0; j < i; j++) + sum += L[i * n + j] * U[j * n + k]; + + U[i * n + k] = A[i * n + k] - sum; + } + + for (int k = i; k < n; k++) + { + if (i == k) + L[i * n + i] = 1; + else + { + int sum = 0; + for (int j = 0; j < i; j++) + sum += L[k * n + j] * U[j * n + i]; + + L[k * n + i] = (A[k * n + i] - sum) / U[i * n + i]; + } + } + } +} + +bool lina_det(double *A, int n, double *det) +{ + // Allocate the space for the L,U matrices. + // I can't think of a version of this algorithm + // where a temporary buffer isn't necessary. + double *T = malloc(sizeof(double) * n * n * 2); + if (T == NULL) + return false; + + // Do the decomposition + double *L = T; + double *U = T + (n * n); + lina_decompLU(A, L, U, n); + + // Knowing that + // + // A = LU + // + // then + // + // det(A) = det(LU) = det(L)det(U) + // + // Since L and U are triangular, their + // determinant is the product of their + // diagonals, so the product of the + // determinants is the product of both + // the diagonals. + + double prod = 1; + for (int i = 0; i < n; i++) + prod *= L[i * n + i] * U[i * n + i]; + + if (det) + *det = prod; + + free(T); + return true; +} + +/* Checks that [A] is kind of upper triangular. +** +*/ +static bool isUpperTriangularEnough(double *A, int n, double eps) +{ + assert(A != NULL && n > 0 && epd > 0); + + // Check that the lower triangular portion (without + // considering the diagonal) is zero. + for (int i = 0; i < n; i++) + for (int j = 0; j < i-1; i++) + if (A[i * n + j] > eps) + return false; + + // Now check that the diagonal is also zero. Though + // since we are using the real version of the QR + // algorithm, only real eigenvalues can be found. + // Any comples eigenvalues will manifest as 2x2 blocks + // on the diagonal, so we need to allow such blocks. + // To do this, a non-zero block is allowed if it's + // not following another non-zero block. + // + // An important thing to note is that 2x2 matrices + // will always be considered upper triangular by this + // function, so the caller must manage this case. + bool flag = false; + for (int i = 0; i < n-1; i++) { + if (fabs(A[(i + 1) * n + i]) > eps) { + if (flag) + return false; + flag = true; + } else + flag = false; + } + return true; +} + +bool lina_eig(double *M, double *E, int n) +{ + double *A = malloc(sizeof(double) * n * n * 3); + if (A == NULL) + return false; + memcpy(A, M, sizeof(double) * n * n); + + double *Q = A + n * n; + double *R = Q + n * n; + + do { + for (int i = 0; i < 100; i++) { + lina_decompQR(A, Q, R, n); // A(n) = QR + lina_dot(R, Q, A, n, n, n); // A(n+1) = RQ + } + } while (!isUpperTriangularEnough(A, n, 0.1)); + + for (int i = 0; i < n; i++) + E[i] = A[i * n + i]; + + free(A); + return true; +} \ No newline at end of file diff --git a/src/lina.h b/src/lina.h index 929eaa9..d3a3aaa 100644 --- a/src/lina.h +++ b/src/lina.h @@ -1,16 +1,25 @@ +#include /* ---- Operations ---- */ -void lina_dot(double *A, double *B, double *C, int m, int n, int l); -void lina_add(double *A, double *B, double *C, int m, int n); -void lina_scale(double *A, double *B, double k, int m, int n); -void lina_transpose(double *A, double *B, int m, int n); -void lina_conv(double *A, double *B, double *C, - int Aw, int Ah, int Bw, int Bh); +void lina_dot(double *A, double *B, double *C, int m, int n, int l); +void lina_add(double *A, double *B, double *C, int m, int n); +double lina_mul(double *v, int n); +bool lina_det(double *A, int n, double *det); +void lina_scale(double *A, double *B, double k, int m, int n); +void lina_transpose(double *A, double *B, int m, int n); +void lina_conv(double *A, double *B, double *C, + int Aw, int Ah, int Bw, int Bh); void lina_dot2(double *A, double *B, double *C, int As, int Bs, int Cs, int m, int n, int l); +bool lina_eig(double *M, double *E, int n); +void lina_decompLU(double *A, double *L, double *U, int n); +void lina_decompQR(double *A, double *Q, double *R, int n); +void lina_orthoNormGramSchmidt(double *A, double *Q, int n); + + /* ---- Utilities ---- */ double *lina_loadMatrixFromStream(FILE *fp, int *width, int *height, char **error); int lina_saveMatrixToStream(FILE *fp, double *A, int width, int height, char **error); \ No newline at end of file diff --git a/src/qr.c b/src/qr.c new file mode 100644 index 0000000..2083b5d --- /dev/null +++ b/src/qr.c @@ -0,0 +1,142 @@ +#include +#include + +typedef struct { + double *items; + int size; +} square_matrix_t; + +typedef struct { + double *items; + int stride; + int size; +} vector_t; + +typedef struct { + vector_t base; + double scale; +} scaled_vector_t; + +static square_matrix_t +square_matrix_from_raw(double *M, int n) +{ + return (square_matrix_t) {.items=M, .size=n}; +} + +static vector_t +get_column_of_square_matrix(square_matrix_t M, int i) +{ + assert(i > -1 && i < M.size); + + return (vector_t) { + .items = M.items + i, + .stride = M.size, + .size = M.size + }; +} + +static void +copy_vector(vector_t V, vector_t S) +{ + assert(V.size == S.size); + for (int i = 0; i < V.size; i++) + V.items[V.stride * i] = S.items[S.stride * i]; +} + +static void +subtract_vector_inplace(vector_t V, scaled_vector_t S) +{ + assert(V.size == S.base.size); + + for (int i = 0; i < V.size; i++) + V.items[V.stride * i] -= S.scale * S.base.items[S.base.stride * i]; +} + +static void +scale_vector_inplace(vector_t V, double a) +{ + for (int i = 0; i < V.size; i++) + V.items[V.stride * i] *= a; +} + +static scaled_vector_t +scale_vector_lazily(vector_t V, double a) +{ + return (scaled_vector_t) {.base=V, .scale=a}; +} + +static double +scalar_product(vector_t V, vector_t U) +{ + assert(V.size == U.size); + + double scale = 0; + for (int i = 0; i < V.size; i++) + scale += V.items[i * V.stride] * U.items[i * U.stride]; + return scale; +} + +static double +normalize_inplace(vector_t V) +{ + // Calculate the sum of the component's squares + double sum = scalar_product(V, V); + + // Calculate the norm and scale the column + // only if the norm isn't zero. + double norm = sqrt(sum); + + if (norm != 0) + scale_vector_inplace(V, 1/norm); + + return norm; +} + +static scaled_vector_t +project(vector_t V, vector_t U) +{ + double scale_vu = scalar_product(V, U); + double scale_uu = scalar_product(U, U); + double ratio = scale_vu / scale_uu; + return scale_vector_lazily(U, ratio); +} + +/** Gram-Schmidt orthonormalization + **/ +void lina_orthoNormGramSchmidt(double *A, double *Q, int n) +{ + square_matrix_t A2 = square_matrix_from_raw(A, n); + square_matrix_t Q2 = square_matrix_from_raw(Q, n); + + for (int i = 0; i < n; i++) { + + vector_t Qi = get_column_of_square_matrix(Q2, i); + vector_t Ai = get_column_of_square_matrix(A2, i); + copy_vector(Qi, Ai); + + for (int j = 0; j < i; j++) { + vector_t Qj = get_column_of_square_matrix(Q2, j); + subtract_vector_inplace(Qi, project(Ai, Qj)); + } + + normalize_inplace(Qi); + // TODO: Handle case of zero norm + } +} + +void lina_decompQR(double *A, double *Q, double *R, int n) +{ + lina_orthoNormGramSchmidt(A, Q, n); + + square_matrix_t A2 = square_matrix_from_raw(A, n); + square_matrix_t Q2 = square_matrix_from_raw(Q, n); + + // Now calculate R by multiplying Q^t and A + for(int i = 0; i < n; i++) { // Iterate over each column i of Q.. + for(int j = 0; j < n; j++) { // ..and over each column j of A + vector_t Qi = get_column_of_square_matrix(Q2, i); + vector_t Aj = get_column_of_square_matrix(A2, j); + R[i * n + j] = scalar_product(Qi, Aj); + } + } +} \ No newline at end of file diff --git a/test.c b/test.c new file mode 100644 index 0000000..91b7e18 --- /dev/null +++ b/test.c @@ -0,0 +1,58 @@ +#include +#include "src/lina.h" + +void print_square_matrix(double *M, int n, FILE *stream) +{ + for (int i = 0; i < n; i++) + { + fprintf(stream, "| "); + for (int j = 0; j < n; j++) + { + fprintf(stderr, "%2.2f ", M[i * n + j]); + } + fprintf(stream, "|\n"); + } + fprintf(stream, "\n"); +} + +void print_vector(double *V, int n, FILE *stream) +{ + fprintf(stream, "[ "); + for (int i = 0; i < n; i++) + fprintf(stderr, "%2.2f ", V[i]); + fprintf(stream, "]\n"); +} + +int main(void) +{ + double A[4] = {1, 2, 3, 4}; + + double L[4]; + double U[4]; + double LU[4]; + lina_decompLU(A, L, U, 2); + lina_dot(L, U, LU, 2, 2, 2); + print_square_matrix(A, 2, stderr); + print_square_matrix(L, 2, stderr); + print_square_matrix(U, 2, stderr); + print_square_matrix(LU, 2, stderr); + + double det; + lina_det(A, 2, &det); + fprintf(stderr, "det=%2.2f\n", det); + + double Q[4]; + double R[4]; + double QR[4]; + lina_decompQR(A, Q, R, 2); + lina_dot(Q, R, QR, 2, 2, 2); + print_square_matrix(Q, 2, stderr); + print_square_matrix(R, 2, stderr); + print_square_matrix(QR, 2, stderr); + + double E[4]; + lina_eig(A, E, 2); + print_vector(E, 2, stderr); + + return 0; +}