1 files changed, 99 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/stxtools/ludcmd.x b/pkg/utilities/nttools/stxtools/ludcmd.x
new file mode 100644
index 00000000..708a9df9
--- /dev/null
+++ b/pkg/utilities/nttools/stxtools/ludcmd.x
@@ -0,0 +1,99 @@
+define	TINY	1.d-20
+
+# ludcmd -- lower-upper decomposition
+# Double-precision version of ludcmp from Numerical Recipes.
+# This differs from the Numerical Recipes version in the following ways:
+# (1) the calling sequence also includes an ISTAT parameter, (2) memory
+# is allocated instead of using the fixed array VV, and (3) double
+# precision is used.
+# This routine decomposes a matrix (in-place) into lower and upper
+# triangular portions.  Use lubksd to obtain a solution to A * X = B
+# or to compute the inverse of the matrix A.
+# If the matrix is singular, ISTAT is set to one.
+#
+# Phil Hodge, 28-Oct-1988  Subroutine copied from Numerical Recipes.
+# Phil Hodge, 10-Sep-1992  Convert to double precision and rename from ludcmp.
+
+procedure ludcmd (a, n, np, indx, d, istat)
+
+double	a[np,np]	# io: input a, output decomposed a
+int	n		# i: logical size of a is n x n
+int	np		# i: space allocated for a
+int	indx[n]		# o: index to be used by xt_lubksb
+double	d		# o: +1 or -1
+int	istat		# o: OK if no problem; 1 if matrix is singular
+#--
+pointer sp
+pointer vv		# scratch space
+double	aamax
+double	sum
+double	dum
+int	i, j, k
+int	imax
+
+begin
+	istat = OK				# initial value
+
+	call smark (sp)
+	call salloc (vv, n, TY_DOUBLE)
+
+	d = 1.d0
+	do i = 1, n {
+	    aamax = 0.d0
+	    do j = 1, n
+		if (abs(a[i,j]) > aamax)
+		    aamax = abs(a[i,j])
+	    if (aamax == 0.d0) {
+		istat = 1
+		return
+	    }
+	    Memd[vv+i-1] = 1.d0 / aamax
+	}
+	do j = 1, n {
+	    if (j > 1) {
+		do i = 1, j-1 {
+		    sum = a[i,j]
+		    if (i > 1) {
+			do k = 1, i-1
+			    sum = sum - a[i,k] * a[k,j]
+			a[i,j] = sum
+		    }
+		}
+	    }
+	    aamax = 0.d0
+	    do i = j, n {
+		sum = a[i,j]
+		if (j > 1) {
+		    do k = 1, j-1
+			sum = sum - a[i,k] * a[k,j]
+		    a[i,j] = sum
+		}
+		dum = Memd[vv+i-1] * abs[sum]
+		if (dum >= aamax) {
+		    imax = i
+		    aamax = dum
+		}
+	    }
+	    if (j != imax) {
+		do k = 1, n {
+		    dum = a[imax,k]
+		    a[imax,k] = a[j,k]
+		    a[j,k] = dum
+		}
+		d = -d
+		Memd[vv+imax-1] = Memd[vv+j-1]
+	    }
+	    indx[j] = imax
+	    if (j != n) {
+		if (a[j,j] == 0.d0)
+		    a[j,j] = TINY
+		dum = 1.d0 / a[j,j]
+		do i = j+1, n
+		    a[i,j] = a[i,j] * dum
+	    }
+	}
+	if (a[n,n] == 0.d0)
+	    a[n,n] = TINY
+
+	call sfree (sp)
+end