From 4d35b7d8defe6513322ba14acdac049dfd82b94a Mon Sep 17 00:00:00 2001 From: Francesco Cozzuto Date: Thu, 30 Mar 2023 01:08:14 +0200 Subject: [PATCH] changed LU to LUP and implemented matrix inversion --- src/lina.c | 195 ++++++++++++++++++++++++++++++++++++++++++++--------- src/lina.h | 8 +-- test.c | 131 ++++++++++++++++++++++++++--------- 3 files changed, 264 insertions(+), 70 deletions(-) diff --git a/src/lina.c b/src/lina.c index 33f4f88..b88f03d 100644 --- a/src/lina.c +++ b/src/lina.c @@ -679,42 +679,83 @@ void lina_conv(double *A, double *B, double *C, } } -void lina_decompLU(double *A, double *L, double *U, int n) +void lina_reallyP(int *P, double *P2, int n) +{ + memset(P2, 0, sizeof(double) * n * n); + + for (int i = 0; i < n; i++) + P2[i * n + P[i]] = 1; +} + +int lina_decompLUP(double *A, double *L, + double *U, int *P, + int n) { assert(n > 0); assert(A != L && A != U && L != U); - // TODO: Handle the case when A can not be - // decomposed. + for (int i = 0; i < n; i++) + P[i] = i; - memset(L, 0, sizeof(double) * n * n); - memset(U, 0, sizeof(double) * n * n); + int swaps = 0; + for (int i = 0; i < n; i++) { + + int v = P[i]; + double max_v = A[v * n + i]; + int max_i = i; + + for (int j = i+1; j < n; j++) { + int u = P[j]; + double abs = fabs(A[u * n + j]); + if (abs > max_v) { + max_v = abs; + max_i = j; + } + } + + if (max_i != i) { + + // Swap rows + int temp = P[i]; + P[i] = P[max_i]; + P[max_i] = temp; + + swaps++; + } + } 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]; + for (int j = 0; j < n; j++) + U[i * n + j] = A[P[i] * n + j]; - U[i * n + k] = A[i * n + k] - sum; - } + memset(L, 0, sizeof(double) * n * n); + for (int i = 0; i < n; i++) + L[i * n + i] = 1; - 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]; - } - } + for (int i = 0; i < n; i++) + for (int j = i+1; j < n; j++) { + double u = U[i * n + i]; + L[j * n + i] = U[j * n + i] / u; + for (int k = 0; k < n; k++) + U[j * n + k] -= L[j * n + i] * U[i * n + k]; } + + return swaps; +} + +static void +printSquareMatrix(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"); } /* Function: lina_det @@ -734,14 +775,20 @@ 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); + double *T = malloc(sizeof(double) * n * n * 2 + sizeof(int) * n); if (T == NULL) return false; // Do the decomposition double *L = T; - double *U = T + (n * n); - lina_decompLU(A, L, U, n); + double *U = L + (n * n); + int *P = (int*) (U + (n * n)); + + int swaps = lina_decompLUP(A, L, U, P, n); + if (swaps < 0) { + free(T); + return false; + } // Knowing that // @@ -758,9 +805,15 @@ bool lina_det(double *A, int n, double *det) // the diagonals. double prod = 1; - for (int i = 0; i < n; i++) - prod *= L[i * n + i] * U[i * n + i]; - + for (int i = 0; i < n; i++) { + double l = L[i * n + i]; + double u = U[i * n + i]; + prod *= l * u; + } + + if (swaps & 1) + prod = -prod; + if (det) *det = prod; @@ -922,7 +975,7 @@ bool lina_eig(double *M, double complex *E, int n) // // y1, y2 = (a + d)/2 +/- 1/2 sqrt{D} // - // y1 and y2 are one the conjugate of the other. Theis + // y1 and y2 are one the conjugate of the other. Their // real part is // // Re{y1, y2} = (a+d)/2 @@ -950,6 +1003,82 @@ bool lina_eig(double *M, double complex *E, int n) E[i] = A[i * n + i]; } + free(T); + return true; +} + +/* Create the n-1 by n-1 matrix D obtained by +** removing the [del_col] column and [del_row] +** frow the n by n matrix M. +*/ +static void +copyMatrixWithoutRowAndCol(double *M, double *D, int n, + int del_col, int del_row) +{ + // Copy the upper-left portion of matrix M + // that comes before the deleted column and + // row. + for (int i = 0; i < del_row; i++) + for (int j = 0; j < del_col; j++) + D[i * (n-1) + j] = M[i * n + j]; + + // Copy the lower left portion that comes + // after both the deleted column and row. + for (int i = del_row+1; i < n; i++) + for (int j = del_row+1; j < n; j++) + D[(i-1) * (n-1) + (j-1)] = M[i * n + j]; + + // Copy the bottom portion that comes after + // the deleted row but before the deleted column. + for (int i = del_row+1; i < n; i++) + for (int j = 0; j < del_col; j++) + D[(i-1) * (n-1) + j] = M[i * n + j]; + + // Copy the right portion that comes after + // the deleted column but before the deleted row. + for (int i = 0; i < del_row; i++) + for (int j = del_col+1; j < n; j++) + D[i * (n-1) + (j-1)] = M[i * n + j]; +} + +bool lina_inverse(double *M, double *D, int n) +{ + double det; + if (!lina_det(M, n, &det)) + return false; + + if (det == 0) + return false; // The matrix can't be inverted + + double *T = malloc(sizeof(double) * ((n-1) * (n-1) + n * n)); + if (T == NULL) + return false; + + double *M_t = T + (n-1) * (n-1); + lina_transpose(M, M_t, n, n); + + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) { + + copyMatrixWithoutRowAndCol(M_t, T, n, j, i); + + double det2; + if (!lina_det(T, n-1, &det2)) { + free(T); + return false; + } + + // If the determinant of M isn't zero, + // neither is this! + assert(det2 != 0); + + bool i_is_odd = i & 1; + bool j_is_odd = j & 1; + int sign = (i_is_odd == j_is_odd) ? 1 : -1; + + D[i * n + j] = sign * det2 / det; + } + free(T); return true; } \ No newline at end of file diff --git a/src/lina.h b/src/lina.h index 5ddaecf..22a8552 100644 --- a/src/lina.h +++ b/src/lina.h @@ -1,13 +1,12 @@ #include #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); -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); +bool lina_inverse(double *M, double *D, int n); void lina_conv(double *A, double *B, double *C, int Aw, int Ah, int Bw, int Bh); @@ -16,11 +15,10 @@ void lina_dot2(double *A, double *B, double *C, int m, int n, int l); bool lina_eig(double *M, double complex *E, int n); -void lina_decompLU(double *A, double *L, double *U, int n); +void lina_reallyP(int *P, double *P2, int n); +int lina_decompLUP(double *A, double *L, double *U, int *P, 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/test.c b/test.c index 90d9955..f3e214b 100644 --- a/test.c +++ b/test.c @@ -25,44 +25,111 @@ void print_vector(double complex *V, int n, FILE *stream) 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 complex E[4]; - lina_eig(A, E, 2); - print_vector(E, 2, stderr); - - double M[5*5] = { + double M[25] = { 1, 2, 3, 4, 5, 5, 1, 2, 3, 4, 4, 5, 1, 2, 3, 3, 4, 5, 1, 2, 2, 3, 4, 5, 1, }; - double complex E2[5]; - lina_eig(M, E2, 5); - print_vector(E2, 5, stderr); + + fprintf(stderr, "# --- M --- #\n"); + print_square_matrix(M, 5, stderr); + + + /* + double L[25]; + double U[25]; + int P[5]; + double P2[25]; + lina_decompLUP(M, L, U, P, 5); + lina_reallyP(P, P2, 5); + + fprintf(stderr, "# --- L --- #\n"); + print_square_matrix(L, 5, stderr); + + fprintf(stderr, "# --- U --- #\n"); + print_square_matrix(U, 5, stderr); + + fprintf(stderr, "# --- P2 --- #\n"); + print_square_matrix(P2, 5, stderr); + + double PA[25]; + lina_dot(P2, M, PA, 5, 5, 5); + fprintf(stderr, "# --- PA --- #\n"); + print_square_matrix(PA, 5, stderr); + + double LU[25]; + lina_dot(L, U, LU, 5, 5, 5); + fprintf(stderr, "# --- LU --- #\n"); + print_square_matrix(LU, 5, stderr); + + double det; + lina_det(M, 5, &det); + fprintf(stderr, "det=%2.2f\n", det); + + fprintf(stderr, "# --- eig(M) --- #\n"); + double complex E[5]; + lina_eig(M, E, 5); + print_vector(E, 5, stderr); + */ + + + double invM[25]; + lina_inverse(M, invM, 5); + + double expI[25]; + lina_dot(M, invM, expI, 5, 5, 5); + + fprintf(stderr, "# --- inv(M) --- #\n"); + print_square_matrix(invM, 5, stderr); + + fprintf(stderr, "# --- I? --- #\n"); + print_square_matrix(expI, 5, stderr); + + + /* + double M[16] = { + 1, 5, 4, 2, + 2, 1, 5, 3, + 4, 3, 2, 5, + 5, 4, 3, 1, + }; + + fprintf(stderr, "# --- M --- #\n"); + print_square_matrix(M, 4, stderr); + + double L[16]; + double U[16]; + int P[4]; + lina_decompLUP(M, L, U, P, 4); + fprintf(stderr, "# --- L,U,P --- #\n"); + print_square_matrix(L, 4, stderr); + print_square_matrix(U, 4, stderr); + + fprintf(stderr, "[ "); + for(int i = 0; i < 4; i++) + fprintf(stderr, "%d ", P[i]); + fprintf(stderr, "]\n"); + + double RP[16]; + double PM[16]; + lina_reallyP(P, RP, 4); + lina_dot(RP, M, PM, 4, 4, 4); + + double LU[16]; + lina_dot(L, U, LU, 4, 4, 4); + + fprintf(stderr, "# --- P,PM,LU --- #\n"); + print_square_matrix(RP, 4, stderr); + print_square_matrix(PM, 4, stderr); + print_square_matrix(LU, 4, stderr); + + double det; + lina_det(M, 4, &det); + + fprintf(stderr, "det(M) = %2.2f\n", det); + */ return 0; }