#include <stdio.h>
#include <math.h>
/* gauss22.c */
#define N 3
int main(void)
{
double A[N][N]={{1., 4., 3.}, {3., 2., 5.}, {2., 6., 6.}}, Aa[N][N]; /*簡単のため係数行列を予め指定*/
double b[N]={4., 5., 6.}, x[N], bb[N], e[N]; /*簡単のため定数ベクトルを予め指定*/
int n = N;
int i, j, k;
double akk, aik, s;
/* save original coefficients */
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++)
Aa[i][j]=A[i][j];
bb[i]=b[i];
}
/* forward operation */
for (k = 0; k < n - 1; k++) {
akk = 1 / A[k][k];
for (i = k + 1; i < n; i++) {
aik = -A[i][k] * akk;
for (j = k + 1; j <n; j++)
A[i][j] += aik * A[k][j];
b[i] += aik * b[k];
}
for (j = k + 1; j < n; j++)
A[k][j] *= akk;
b[k]*=akk;
}
/* backward operation */
x[n - 1] = b[n - 1] / A[n - 1][n - 1];
for (k = n - 2; k >= 0; k--) {
s = 0.0;
for (j = k + 1; j < n; j++)
s += A[k][j] * x[j];
x[k]=b[k]-s;
}
/* check */
for (i = 0; i < n; i++) {
s=0.0;
for (j = 0; j < n; j++)
s += Aa[i][j] * x[j];
e[i] = s - bb[i];
printf("x(%d)=%f error=%f\n",i, x[i], e[i]);
}
return 0;
}