# File rvsao/Emvel/jgauss.x
# June 14, 1994
# SPP by Doug Mink from Press et al. GAUSSJ

procedure jgauss (a,b,n,np,ierr)

# Linear equation solution by Gauss-Jordan elimination

double	a[np,np]	# Input matrix
double	b[np]		# Right hand side vectors
int	n		# Number of filled rows and columns in input matrix
int	np		# Physical dimension of input matrix
int	ierr		# 0=OK, 2=Singular matrix

int	i,j,k,l,ll,irow,icol
int	ipiv[50],indxr[50],indxc[50]
double	big,dum,pivinv

begin
	do j = 1, n {
	    ipiv[j] = 0
	    }

# Main loop of columns to be reduced
	do i = 1, n {
	    big = 0.d0
	    do j = 1, n {
		if (ipiv[j] != 1) {
		    do k = 1, n {
			if (ipiv[k] == 0) {
			    if (abs (a[j,k]) > BIG) {
				big = abs (a[j,k])
				irow = j
				icol = k
				}
			    }
			else if (ipiv[k] > 1) {
			    call eprintf ("JGAUSS: Singular matrix\n")
			    ierr = 2
			    return
			    }
			}
		    }
		}

	    ipiv[icol] = ipiv[icol] + 1

# Interchange rows to put pivot on the diagonal
	    if (irow != icol) {
		do l = 1, n {
		    dum = a[irow,l]
		    a[irow,l] = a[icol,l]
		    a[icol,l] = dum
		    }
		dum = b[irow]
		b[irow] = b[icol]
		b[icol] = dum
		}

# Divide pivot row by pivot element at (irow,icol)
	    indxr[i] = irow
	    indxc[i] = icol
	    if (a[icol,icol] == 0.d0) {
		call eprintf ("JGAUSS: Singular matrix\n")
		ierr = 2
		return
		}
	    pivinv = 1.d0 / a[icol,icol]
	    a[icol,icol] = 1.d0
	    do l = 1, n {
		a[icol,l] = a[icol,l] * pivinv
		}
	    b[icol] = b[icol] * pivinv
	    do ll = 1, n {
		if (ll != icol) {
		    dum = a[ll,icol]
		    a[ll,icol] = 0.d0
		    do l = 1, n {
			a[ll,l] = a[ll,l] - (a[icol,l] * dum)
			}
		    b[ll] = b[ll] - (b[icol] * dum)
		    }
		}
	    }

#  Unscramble column interchanges
	do l = n, 1, -1 {
	    if (indxr[l] != indxc[l]) {
		do k = 1, n {
		    dum = a[k,indxr[l]]
		    a[k,indxr[l]] = a[k,indxc[l]]
		    a[k,indxc[l]] = dum
		    }
		}
	    }
	ierr = 0
	return
end