changed LU to LUP and implemented matrix inversion
This commit is contained in:
Francesco Cozzuto
+158
-29
@@ -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;
|
||||
}
|
||||
|
||||
for (int k = i; k < n; k++)
|
||||
{
|
||||
if (i == k)
|
||||
memset(L, 0, sizeof(double) * n * n);
|
||||
for (int i = 0; i < n; i++)
|
||||
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,8 +805,14 @@ 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
|
||||
@@ -953,3 +1006,79 @@ bool lina_eig(double *M, double complex *E, int n)
|
||||
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;
|
||||
}
|
||||
+3
-5
@@ -1,13 +1,12 @@
|
||||
#include <complex.h>
|
||||
#include <stdbool.h>
|
||||
|
||||
/* ---- 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);
|
||||
@@ -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;
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user