#include <stdio.h>
#include <stdlib.h>
#define N1 3
#define N2 4
void make_square_matrix(double ***a, int n)
{
int i;
if ((*a = (double **)malloc(n * sizeof(double *))) == NULL)
exit(1);
for (i = 0; i < n; i++)
if (((*a)[i] = (double *)malloc(n * sizeof(double))) == NULL)
exit(1);
}
void free_square_matrix(double **a, int n)
{
int i;
for (i = 0; i < n; i++)
free(a[i]);
free(a);
}
void make_vector(double **a, int n)
{
if ((*a = (double *)malloc(n * sizeof(double))) == NULL)
exit(1);
}
void free_vector(double *a)
{
free(a);
}
/* LU分解 */
void LU_decomp(double **a, double **l, double **u, int n)
{
int i, j, k;
double p, s;
for (i = 0; i < n; i++) {
l[i][0] = a[i][0];
u[i][0] = 0.0;
}
p = a[0][0];
u[0][0] = 1.0;
for (j = 1; j < n; j++) {
l[0][j] = 0.0;
u[0][j] = a[0][j] / p;
}
for (i = 0; i < n; i++) {
s = a[i][i];
for (k = 0; k < i; k++)
s -= l[i][k] * u[k][i];
l[i][i] = p = s;
u[i][i] = 1.0;
for (j = i + 1; j < n; j++) {
s = a[j][i];
for (k = 0; k < i; k++)
s -= l[j][k] * u[k][i];
l[j][i] = s;
u[j][i] = 0.0;
s = a[i][j];
for (k = 0; k < i; k++)
s -= l[i][k] * u[k][j];
u[i][j] = s / p;
l[i][j] = 0.0;
}
}
}
/* LU分解を利用して一次方程式を解く */
void LU_solve(double **l, double **u, double *x, double *b, int n)
{
int i, j, m;
double s, *y;
make_vector(&y, n);
y[0] = b[0] / l[0][0];
for (i = 1; i < n; i++) {
s = b[i];
for (j = 0; j < i; j++)
s -= l[i][j] * y[j];
y[i] = s / l[i][i];
}
m = n - 1;
x[m] = y[m] / u[m][m];
for (i = m - 1; i >= 0; i--) {
s = y[i];
for (j = i + 1; j < n; j++)
s -= u[i][j] * x[j];
x[i] = s / u[i][i];
}
free_vector(y);
}
void print_matrix(double **a, int n, char *s)
{
int i, j;
printf("*** %s ***\n", s);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++)
printf("% 11.7f ", a[i][j]);
putchar('\n');
}
}
void print_vector(double *a, int n, char *s)
{
int i;
printf("*** %s ***\n", s);
for (i = 0; i < n; i++)
printf("% 11.7f\n", a[i]);
}
int main(void)
{
double **a, **l, **u, *x, *b;
int i, j;
double a1[N1][N1] = {{2.0, 3.0, -1.0}, {4.0, 4.0, -3.0}, {2.0, -3.0, 1.0}};
double b1[N1] = {5.0, 4.0, -1.0};
double a2[N2][N2] = {{1.0, 1.0, 0.0, 3.0}, {2.0, 1.0, -1.0, 1.0}, {3.0, -1.0, -1.0, 2.0}, {-1.0, 2.0, 3.0, -2.0}};
double b2[N2] = {4.0, 1.0, -2.0, 4.0};
make_square_matrix(&a, N1);
make_square_matrix(&l, N1);
make_square_matrix(&u, N1);
make_vector(&x, N1);
make_vector(&b, N1);
for (i = 0; i < N1; i++)
for (j = 0; j < N1; j++)
a[i][j] = a1[i][j];
for (i = 0; i < N1; i++)
b[i] = b1[i];
LU_decomp(a, l, u, N1);
print_matrix(a, N1, "行列 A");
print_matrix(l, N1, "行列 L");
print_matrix(u, N1, "行列 U");
print_vector(b, N1, "行列 B");
LU_solve(l, u, x, b, N1);
print_vector(x, N1, "行列 X(解)");
free_square_matrix(a, N1);
free_square_matrix(l, N1);
free_square_matrix(u, N1);
free_vector(x);
free_vector(b);
make_square_matrix(&a, N2);
make_square_matrix(&l, N2);
make_square_matrix(&u, N2);
make_vector(&x, N2);
make_vector(&b, N2);
for (i = 0; i < N2; i++)
for (j = 0; j < N2; j++)
a[i][j] = a2[i][j];
for (i = 0; i < N2; i++)
b[i] = b2[i];
LU_decomp(a, l, u, N2);
print_matrix(a, N2, "行列 A");
print_matrix(l, N2, "行列 L");
print_matrix(u, N2, "行列 U");
print_vector(b, N2, "行列 B");
LU_solve(l, u, x, b, N2);
print_vector(x, N2, "行列 X(解)");
free_square_matrix(a, N2);
free_square_matrix(l, N2);
free_square_matrix(u, N2);
free_vector(x);
free_vector(b);
return 0;
}