#* HISTORY *
#* B.Simon	09-Nov-92	converted to spp from asurv 1.2

# MATINV -- Matrix inverse routine

# This program is adopted from Philip R. Bevington, "Data Reduction 
# and Error Analysis for the Physical Sciences', 1969 McGraw Hill (NY:NY), 
# program b-2 p. 302. Several minor modifications were done by T. Isobe.

procedure matinv (array, nvar, det)

double	array[nvar,nvar]	# u: symmetric matrix, inverse on output
int	nvar			# i: size of matrix
double	det			# o: matrix determinant (0 indicates singular)
#--
double	amax, save
int	i, j, k, l
pointer	sp, ik, jk

begin
	# Allocate dynamic memory for temporary arrays

	call smark (sp)
	call salloc (ik, nvar, TY_INT)
	call salloc (jk, nvar, TY_INT)

	call aclri (Memi[ik], nvar)
	call aclri (Memi[jk], nvar)

	det = 1.0
	do k= 1, nvar {
	    
	    # Find largest element array(i,j) in rest of matrix
	    
	    amax = 0.0       
	    repeat {
		do i= k, nvar {
		    do j = k, nvar {
			if (abs (amax) - abs (array[i,j]) > 0.0) 
			    break
			amax = array[i,j]
			Memi[ik+k-1] = i
			Memi[jk+k-1] = j
		    }
		}
		
		# Interchange rows and columns to put amax in array[k,k]
		
		if (amax == 0.0) {
		    det = 0.0
		    return
		}
		
		i = Memi[ik+k-1]
		if (i - k < 0)
		    next
		
		if ((i - k) > 0) {
		    do j = 1, nvar {
			save = array[k,j]
			array[k,j] = array[i,j]
			array[i,j] = - save
		    }
		}
		
		
		j = Memi[jk+k-1]
		if (j - k < 0)
		    next
		
		if ((j - k) > 0) {
		    do i = 1, nvar {
			save = array[i,k]
			array[i,k] = array[i,j]
			array[i,j] = - save
		    }
		}
		
		break
	    }

	    # Accumulate elements of inverse matrix

	    do i = 1, nvar {
		if (i != k)
		    array[i,k] = - array[i,k] / amax
	    }

	    do i = 1, nvar {
		do j = 1, nvar {
		    if (i != k && j != k)
			array[i,j] = array[i,j] + array[i,k] * array[k,j]
		}
	    }

	    do j = 1, nvar {
		if(j != k) 
		    array[k,j] = array[k,j] / amax
	    }

	    array[k,k] = 1.0 / amax
	    det = det * amax
	}

	# Restore ordering of matrix

	do l = 1, nvar {
	    k = nvar - l + 1
	    j = Memi[ik+k-1]

	    if (j > k) {
		do i = 1, nvar {
		    save = array[i,k]
		    array[i,k] = - array[i,j]
		    array[i,j] = save
		}
	    }

	    i = Memi[jk+k-1]
	    if (i > k) {
		do j = 1, nvar {
		    save = array[k,j]
		    array[k,j] = - array[i,j]
		    array[i,j] = save
		}
	    }
	}

	call sfree (sp)
	return
end