/*matin1 matrix inversion routine used throughout the code*/
#include "cddefines.h"
#include "matin1.h"
/*matrix inversion routine used throughout the code
 * matrix inversion with accompanying solution of linear equations
 * solves the matrix equation A*X = B for the vector X
 * */
#define	DIMINDEX	100

void matin1(double avec[], 
  long int NDECLARE, 
  long int n1, 
  long int *nerror)
{
	long int i, 
	  i1, 
	  i2, 
	  i2_, 
	  i3, 
	  i_, 
	  icol, 
	  index[DIMINDEX], 
	  ipemat, 
	  ipivcl1, 
	  ipivcl2, 
	  ipivcol, 
	  jcol, 
	  lpiv=0, 
	  mn, 
	  mn_, 
	  nmin1;
	double deter, 
	  pivot, 
	  swap;

#	ifdef DEBUG_FUN
	fputs( "<+>matin1()\n", debug_fp );
#	endif

	/*     avec - 2D matrix to be inverted
	 *     NDECLARE - first dimension of avec as declared in calling program
	 *     n1 is order of matrix, so n1 <= NDECLARE
	 *
	 *     avec: a two dimensional array with n1 as column size,containing
	 *       the matrix of order n in its first n columns. the matrix of
	 *       constant vectors , b, is stored in columns n+1 through n+m of e
	 *
	 *     NERROR  - OUTPUT PARAMETER WHICH IS SET ON RETURN TO NON ZERO, IF
	 *     AT ANY ELIMINATION STEP THE CORRESPONDING COLUMN OF A CONTAINED
	 *     ONLY ZEROES. A PRINTED MESSAGE IS GIVEN ON FILE TAPE2
	 *     INDICATING THE COLUMN.
	 *
	 *     ON RETURN, avec WILL, IF NERROR = 0 CONTAIN THE INVERSE OF A IN
	 *     ITS FIRST N COLUMNS AND IF M NOT =  0 THE SOLUTION MATRIX IN
	 *     THE NEXT FOLLOWING  (N+1) TO (N+M) COLUMNS.
	 * */


	deter = 1.0f;
	if( n1 > DIMINDEX )
	{
		fprintf( ioQQQ, " matin1 called with dimIndex exceeded\n" );
		puts( "[Stop]" );
		exit(0);
	}

	ipemat = n1 + 1;
	nmin1 = n1 - 1;
	/*     the routine does its own evaluation for doubLE SUBSCRIPTING OF
	 *     array avec. */
	ipivcol = 1 - NDECLARE;
	/*     main loop to invert the matrix */
	for( mn=1; mn <= n1; mn++ )
	{
		mn_ = mn - 1;
		pivot = 0.0f;
		ipivcol += NDECLARE;
		/*     search for next pivot in column main. */
		ipivcl1 = ipivcol + mn - 1;
		ipivcl2 = ipivcol + nmin1;
		for( i=ipivcl1; i <= ipivcl2; i++ )
		{
			i_ = i - 1;
			if( fabs(avec[i_]) > fabs(pivot) )
			{
				pivot = avec[i_];
				lpiv = i;
			}
		}
		/*     is pivot different from zero? */
		if( pivot == 0 )
		{
			*nerror = -1;
			
#	ifdef DEBUG_FUN
			fputs( " <->matin1()\n", debug_fp );
#	endif
			return;
		}
		/*     get the pivot-line indicator and swap lines IF NECESSARY */
		icol = lpiv - ipivcol + 1;
		index[mn_] = icol;
		if( icol > mn )
		{
			/*     complement the determinant */
			deter = -deter;
			/*     pointer to line pivot found */
			icol -= NDECLARE;
			/*     pointer to exact pivot line */
			i3 = mn - NDECLARE;
			for( i=1; i <= ipemat; i++ )
			{
				i_ = i - 1;
				icol += NDECLARE;
				i3 += NDECLARE;
				swap = avec[i3-1];
				avec[i3-1] = avec[icol-1];
				avec[icol-1] = swap;
			}
		}
		/*     compute determinant */
		deter *= pivot;
		pivot = 1.e0/pivot;
		/*     transform pivot column */
		i3 = ipivcol + nmin1;
		for( i=ipivcol; i <= i3; i++ )
		{
			i_ = i - 1;
			avec[i_] *= -pivot;
		}
		avec[ipivcl1-1] = pivot;
		/*     pivot element transformed
		 *
		 *     now convert rest of the matrix */
		i1 = mn - NDECLARE;
		/*     pointer to pivot line elements */
		icol = 1 - NDECLARE;
		/*     general column pointer */
		for( i=1; i <= ipemat; i++ )
		{
			i_ = i - 1;
			icol += NDECLARE;
			i1 += NDECLARE;
			/*     pointers moved */
			if( i != mn )
			{
				/*     pivot column excluded */
				jcol = icol + nmin1;
				swap = avec[i1-1];
				i3 = ipivcol - 1;
				for( i2=icol; i2 <= jcol; i2++ )
				{
					i2_ = i2 - 1;
					i3 += 1;
					avec[i2_] += swap*avec[i3-1];
				}
				avec[i1-1] = swap*pivot;
			}
		}
	}

	/*     now rearrange the matrix to get right invers */
	for( i1=1; i1 <= n1; i1++ )
	{
		mn = n1 + 1 - i1;
		lpiv = index[mn-1];
		if( lpiv != mn )
		{
			icol = (lpiv - 1)*NDECLARE + 1;
			jcol = icol + nmin1;
			ipivcol = (mn - 1)*NDECLARE + 1 - icol;
			for( i2=icol; i2 <= jcol; i2++ )
			{
				i2_ = i2 - 1;
				i3 = i2 + ipivcol;
				swap = avec[i2_];
				avec[i2_] = avec[i3-1];
				avec[i3-1] = swap;
			}
		}
	}

	*nerror = 0;
	
#	ifdef DEBUG_FUN
	fputs( " <->matin1()\n", debug_fp );
#	endif
	return;

}