/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:19 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "sgei.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	ONE	1.0e0
#define	ZERO	0.0e0
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ sgei(
float *a,
long lda,
long n,
long ipvt[],
float work[])
{
#define A(I_,J_)	(*(a+(I_)*(lda)+(J_)))
	long int i, j, k, kb, kp1, l, nm1;
	float t;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	long *const Ipvt = &ipvt[0] - 1;
	float *const Work = &work[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1994-10-20 SGEI   Krogh  Changes to use M77CON
	 *>> 1987-08-18 SGEI   Lawson  Initial code.
	 *--S replaces "?": ?GEI, ?AXPY, ?SCAL, ?SWAP
	 *
	 *     This subroutine computes the inverse of the N x N matrix, A,
	 *     using the LU factorization of A given in the array A().
	 *     ------------------------------------------------------------------
	 *                        Subroutine arguments
	 *
	 *     A(,)  [inout]  An array of size at least N x N.  On entry must
	 *            contain the LU factors of an N x N  matrix, A.  It is
	 *            expected that this factorization will have been computed by
	 *            use of _GEFA, either directly or indirectly via use of
	 *            _GEFS or _GEFSC.
	 *            On return contains the N x N inverse matrix of A.
	 *
	 *     LDA  [in]  Leading dimensioning parameter for the array A(,).
	 *
	 *     N  [in]  The order of the original matrix, A.
	 *
	 *     IPVT()  [in]  An integer array of length at least N, containg a
	 *           record of the row interchanges made during factorization of
	 *           A.
	 *
	 *     WORK()  [scratch]  An array of length at least N used by this
	 *        subroutine as working space.
	 *     ------------------------------------------------------------------
	 *     ERROR CONDITION
	 *
	 *        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
	 *        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
	 *        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
	 *        setting of LDA.  The user can avoid sending a singular matrix
	 *        to this subr by testing INFO (set by _GEFS or _GEFA) or
	 *        RCOND (set by _GEFSC or _GERC) before calling this subr.
	 *        Nonsingularity is indicated by INFO .eq. 0 or RCOND .ne. ZERO.
	 *     ------------------------------------------------------------------
	 *     LINPACK. THIS VERSION DATED 08/14/78 .
	 *     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
	 *     Ref: LINPACK Users' Guide, by J. J. Dongarra, C. B. Moler,
	 *     J. R. Bunch, and G. W. Stewart, publ by Soc. for Indust. and Appl.
	 *     Math, Philadelphia, 1979.
	 *     Adapted from LINPACK for the JPL Math77 library by
	 *     C. L. Lawson, JPL, Aug 1987.
	 *     ------------------------------------------------------------------
	 *     Subprograms referenced: SAXPY, SSCAL, SSWAP
	 *     ------------------------------------------------------------------ */
	/*     ------------------------------------------------------------------
	 *     COMPUTE INVERSE(U)
	 * */
	for (k = 1; k <= n; k++)
	{
		if (A(k - 1,k - 1) != ZERO)
		{
			A(k - 1,k - 1) = ONE/A(k - 1,k - 1);
		}
		else
		{
			ermsg( "SGEI", 1, 0, "A diagonal element is zero", '.' );
			return;
		}
		t = -A(k - 1,k - 1);
		sscal( k - 1, t, &A(k - 1,0), 1 );
		kp1 = k + 1;
		for (j = kp1; j <= n; j++)
		{
			t = A(j - 1,k - 1);
			A(j - 1,k - 1) = ZERO;
			saxpy( k, t, &A(k - 1,0), 1, &A(j - 1,0), 1 );
		}
	}
 
	/*        FORM INVERSE(U)*INVERSE(L)
	 * */
	nm1 = n - 1;
	for (kb = 1; kb <= nm1; kb++)
	{
		k = n - kb;
		kp1 = k + 1;
		for (i = kp1; i <= n; i++)
		{
			Work[i] = A(k - 1,i - 1);
			A(k - 1,i - 1) = ZERO;
		}
		for (j = kp1; j <= n; j++)
		{
			t = Work[j];
			saxpy( n, t, &A(j - 1,0), 1, &A(k - 1,0), 1 );
		}
		l = Ipvt[k];
		if (l != k)
			sswap( n, &A(k - 1,0), 1, &A(l - 1,0), 1 );
	}
	return;
#undef	A
} /* end of function */