/********************************/
/*	gauss3  valdes 8 30 82	*/
/*				*/
/*	Gaussian elimination	*/
/*	matrix inversion	*/
/*	Written by J.F. Jarvis	*/
/********************************/

gauss3 (n, ep, a, x)		/* 	n = number of rows ( = columns) ep =
				   error value (typically = .00001) a = pointer
				   to doubleing point matrix to be inverted x =
				   pointer to inverted matrix to be returned
				   note---original matrix a is distroyed in
				   thhe inversion. */
int     n;
double  ep, *a, *x;
{
    int     i, j, k, kp, col, row, rc, rp1, cp1, n2, rpc, op, rpn;
    double  z, t, ratio, absf (), *fp;

    n2 = n * n;
    col = 0;
    cp1 = 1;

    for (fp = x; fp < x + n2; fp++)
	*fp = 0.0;
    for (fp = x; fp < x + n2; fp += (n + 1))
	*fp = 1.0;

    rpn = 0;
    for (row = 0; row < n2; row += n) {
	z = 0.0;
	rpn += n;
	kp = 0;
	rc = (row + col);
	for (k = rc; k < n2; k += n)
	    if ((z - (t = absf (*(a + k)))) < 0) {
		z = t;
		kp = (k - col);
	    }

	if ((row - kp) < 0) {
	    op = kp + col;
	    for (j = rc; j < rpn; j++) {
		z = *(a + j);
		*(a + j) = *(a + op);
		*(a + op) = z;
		op++;
	    }
	    for (j = row; j < rpn; j++) {
		z = *(x + j);
		*(x + j) = *(x + kp);
		*(x + kp) = z;
		kp++;
	    }
	}

	if ((absf (*(a + row + col)) - ep) <= 0)
	    return (0);

	if ((col - n + 1) == 0)
	    break;

	for (k = rpn; k < n2; k += n) {
	    if ((t = (*(a + k + col))) == 0.0)
		continue;
	    ratio = t / (*(a + row + col));
	    for (j = cp1; j < n; j++)
		*(a + k + j) -= (ratio * (*(a + row + j)));
	    for (j = 0; j < n; j++)
		*(x + k + j) -= (ratio * (*(x + row + j)));
	}

	col++;
	cp1++;
    }

    op = n - 1;
    for (row = n2 - n; row >= 0; row -= n) {
	ratio = *(a + row + op);
	for (col = 0; col < n; col++) {
	    z = 0.0;
	    rc = row + col;
	    cp1 = op + 1;
	    if ((op + 1 - n) != 0) {
		for (k = row + n; k < n2; k += n) {
		    z += (*(a + row + cp1)) * (*(x + k + col));
		    cp1++;
		}
	    }

	    *(x + rc) = (*(x + rc) - z) / ratio;

	}
	op--;
    }

    return (1);
}

/*Absolute value*/
double
        absf (a)
double  a;
{
    double  b;

    b = a < 0.? -a : a;
    return (b);
}