/*  @(#)inverse.c       2.1  6/26/87  */
#include <math.h>
#include <grass/libtrans.h>

#define EPSILON 1.0e-16

/* DIM_matrix is defined in "libtrans.h" */
#define N	DIM_matrix

/*
 * inverse: invert a square matrix (puts pivot elements on main diagonal).
 *          returns arg2 as the inverse of arg1.
 *
 *  This routine is based on a routine found in Andrei Rogers, "Matrix
 *  Methods in Urban and Regional Analysis", (1971), pp. 143-153.
 */
int inverse(double m[N][N])
{
    int i, j, k, l, ir = 0, ic = 0;
    int ipivot[N], itemp[N][2];
    double pivot[N], t;
    double fabs();


    if (isnull(m))
	return (-1);


    /* initialization */
    for (i = 0; i < N; i++)
	ipivot[i] = 0;

    for (i = 0; i < N; i++) {
	t = 0.0;		/* search for pivot element */

	for (j = 0; j < N; j++) {
	    if (ipivot[j] == 1)	/* found pivot */
		continue;

	    for (k = 0; k < N; k++)
		switch (ipivot[k] - 1) {
		case 0:
		    break;
		case -1:
		    if (fabs(t) < fabs(m[j][k])) {
			ir = j;
			ic = k;
			t = m[j][k];
		    }
		    break;
		case 1:
		    return (-1);
		    break;
		default:	/* shouldn't get here */
		    return (-1);
		    break;
		}
	}

	ipivot[ic] += 1;
	if (ipivot[ic] > 1) {	/* check for dependency */
	    return (-1);
	}

	/* interchange rows to put pivot element on diagonal */
	if (ir != ic)
	    for (l = 0; l < N; l++) {
		t = m[ir][l];
		m[ir][l] = m[ic][l];
		m[ic][l] = t;
	    }

	itemp[i][0] = ir;
	itemp[i][1] = ic;
	pivot[i] = m[ic][ic];

	/* check for zero pivot */
	if (fabs(pivot[i]) < EPSILON) {
	    return (-1);
	}

	/* divide pivot row by pivot element */
	m[ic][ic] = 1.0;

	for (j = 0; j < N; j++)
	    m[ic][j] /= pivot[i];

	/* reduce nonpivot rows */
	for (k = 0; k < N; k++)
	    if (k != ic) {
		t = m[k][ic];
		m[k][ic] = 0.0;

		for (l = 0; l < N; l++)
		    m[k][l] -= (m[ic][l] * t);
	    }
    }

    /* interchange columns */
    for (i = 0; i < N; i++) {
	l = N - i - 1;
	if (itemp[l][0] == itemp[l][1])
	    continue;

	ir = itemp[l][0];
	ic = itemp[l][1];

	for (k = 0; k < N; k++) {
	    t = m[k][ir];
	    m[k][ir] = m[k][ic];
	    m[k][ic] = t;
	}
    }

    return 1;
}




#define ZERO 1.0e-8

/*
 * isnull: returns 1 if matrix is null, else 0.
 */

int isnull(double a[N][N])
{
    register int i, j;
    double fabs();


    for (i = 0; i < N; i++)
	for (j = 0; j < N; j++)
	    if ((fabs(a[i][j]) - ZERO) > ZERO)
		return 0;

    return 1;
}