#include
#define MAX_N 10
int ReadData(double A[][MAX_N], double b[], int *N)
{
int i, j;
if (scanf("%d", N) == EOF)
return 0;
for (i=0; i<*N; i++)
for (j=0; j<*N; j++)
scanf("%lf", &A[i][j]);
for (i=0; i<*N; i++)
scanf("%lf", &b[i]);
return 1;
}
void LRDecomposition(double L[][MAX_N], double R[][MAX_N], double A[][MAX_N], int N)
{
int i, j, k, m;
/* 初始化L */
for (i=0; i for (j=i; j {
if (j == i)
L[i][j] = 1;
else
L[i][j] = 0;
}
/* 初始化R */
for (i=0; i for (j=0; j R[i][j] = 0;
/* 计算第1步 */
for (j=0; j R[0][j] = A[0][j];
for (i=1; i L[i][0] = A[i][0]/R[0][0];
/* 计算第2至N步 */
for (k=1; k {
/* 计算R */
for (j=k; j {
R[k][j] = A[k][j];
for (m=0; m R[k][j] -= L[k][m]*R[m][j];
}
/* 计算L */
for (i=k+1; i {
L[i][k] = A[i][k];
for (m=0; m L[i][k] -= L[i][m]*R[m][k];
L[i][k] /= R[k][k];
}
}
}
void SolveEquation(double x[], double y[], double L[][MAX_N], double R[][MAX_N], double b[], int N)
{
int i, j;
/* 计算y */
y[0] = b[0];
for (i=1; i {
y[i] = b[i];
for (j=0; j y[i] -= L[i][j]*y[j];
}
/* 计算x */
x[N-1] = y[N-1]/R[N-1][N-1];
for (i=N-2; i>=0; i--)
{
x[i] = y[i];
for (j=i+1; j x[i] -= R[i][j]*x[j];
x[i] /= R[i][i];
}
}
void DisplayResult(double L[][MAX_N], double R[][MAX_N], double x[], double y[], int N)
{
int i, j;
printf("L = \n");
for (i=0; i {
for (j=0; j printf("%10.6f", L[i][j]);
printf("\n");
}
printf("R = \n");
for (i=0; i {
for (j=0; j printf("%10.6f", R[i][j]);
printf("\n");
}
printf("y = \n");
for (i=0; i printf("%10.6f", y[i]);
printf("\nx = \n");
for (i=0; i printf("%10.6f", x[i]);
printf("\n\n");
}
int main()
{
int N, T = 0;
double A[MAX_N][MAX_N];
double L[MAX_N][MAX_N];
double R[MAX_N][MAX_N];
double b[MAX_N];
double y[MAX_N];
double x[MAX_N];
freopen("in.txt", "r", stdin);
freopen("out.txt", "w", stdout);
while (ReadData(A, b, &N))
{
LRDecomposition(L, R, A, N);
SolveEquation(x, y, L, R, b, N);
printf("Case # %d\n", ++T);
DisplayResult(L, R, x, y, N);
}
return 0;
}
-------------------------------------->
