/*
 * gauss1.c --- Gauss の消去法で連立1次方程式を解く
 */

#include "matutil.h"

gauss1(int n, matrix a, vector x)
{
    int i, j, k;
    double q, s;

    /* 前進消去 */
    for (k = 1; k <= n-1; k++)
        for (i = k + 1; i <= n; i++) {
            q = a[i][k] / a[k][k];
#ifdef  ORIGINAL
            for (j = 1; j <= n; j++)
#else
            for (j = k + 1; j <= n; j++)
#endif
                a[i][j] -= q * a[k][j];
            x[i] -= q * x[k];
        }

    /* 後退代入 */
    x[n] /= a[n][n];
    for (k = n - 1; k >= 1; k--) {
        s = x[k];
        for (j = k + 1; j <= n; j++)
            s -= a[k][j] * x[j];
        x[k] = s / a[k][k];
    }
}