/*	matinv.cp */

/*
	a procedure to invert a matrix
	translated from FORTRAN code by
	Robert M. Delaney
	5/2/91
*/

/*
	array indices run from 0 to n-1 where n is the order of the square matrix
*/

/*
	we need to dynamically declare the matrix array and pass it to this procedure as per
	
	array = (double **)malloc(n*sizeof(double *));
	for (i=0;i<n;i++)
		array[i] = (double *)malloc(long)(n*sizeof(double));
*/

#include <stdlib.h>
#include "matinv.h"

static double dabs(double);

double matinv(double **array, long n)
{
	long	i, j, k, L, *ik, *jk;
	double	amax, save, zero, one, determinant;
	
	zero = 0;
	one = 1;
	
	ik = (long *)malloc(n*sizeof(long));
	jk = (long *)malloc(n*sizeof(long));
	
	determinant = one;
	
	for(k=0; k<n; k++)
	{
		/* find largest element in abs value in rest of matrix */
		amax = zero;
		
		for(i=k; i<n; i++)
		{
			for(j=k; j<n; j++)
			{
				if(dabs(array[i][j]) > dabs(amax))
				{
					amax = array[i][j];
					ik[k] = i;
					jk[k] = j;
				}
			}
		}
		
		if(dabs(amax) == zero)
		{
			determinant = 0;
			return determinant;
		}
		
		/* interchange rows and columns to put amax in array(k,k) */
		i = ik[k];
		if(i>k)
		{
			for(j=0; j<n; j++)
			{
				save = array[k][j];
				array[k][j] = array[i][j];
				array[i][j] = -save;
			}
		}
		
		j = jk[k];
		if(j>k)
		{
			for(i=0; i<n; i++)
			{
				save = array[i][k];
				array[i][k] = array[i][j];
				array[i][j] = -save;
			}
		}
		
		/* accumulate elements of inverse matrix */
		for(i=0; i<n; i++)
		{
			if(i!=k)
				array[i][k] = -array[i][k]/amax;
		}
		
		for(i=0; i<n; i++)
		{
			for(j=0; j<n; j++)
			{
				if(i!=k && j!=k)
					array[i][j] = array[i][j] + array[i][k]*array[k][j];
			}
		}
		
		for(j=0; j<n; j++)
		{
			if(j!=k)
				array[k][j] = array[k][j]/amax;
		}
		
		array[k][k] = one/amax;
		determinant = determinant*amax;
	}
		
	/* restore ordering of the matrix */
	for(L=0; L<n; L++)
	{
		k = n - L - 1;
		j = ik[k];
		
		if(j>k)
		{
			for(i=0; i<n; i++)
			{
				save = array[i][k];
				array[i][k] = -array[i][j];
				array[i][j] = save;
			}
		}
		
		i = jk[k];
		
		if(i>k)
		{
			for(j=0; j<n; j++)
			{
				save = array[k][j];
				array[k][j] = -array[i][j];
				array[i][j] = save;
			}
		}
	}
	
	free(ik);
	free(jk);
	
	return determinant;
	
} /* matinv */


static double dabs(double x)
{
	if (x<0)
		x = -x;
	return(x);
	
} /* dabs */
