Initial commit

1 files changed, 110 insertions, 0 deletions

diff --git a/math/gsurfit/gs_chomat.gx b/math/gsurfit/gs_chomat.gx
new file mode 100644
index 00000000..023b3c12
--- /dev/null
+++ b/math/gsurfit/gs_chomat.gx
@@ -0,0 +1,110 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+include <math/gsurfit.h>
+$if (datatype == r)
+include "gsurfitdef.h"
+$else
+include "dgsurfitdef.h"
+$endif
+
+# GSCHOFAC -- Routine to calculate the Cholesky factorization of a
+# symmetric, positive semi-definite banded matrix. This routines was
+# adapted from the bchfac.f routine described in "A Practical Guide
+# to Splines", Carl de Boor (1978).
+
+procedure $tgschofac (matrix, nbands, nrows, matfac, ier)
+
+PIXEL matrix[nbands, nrows]	# data matrix
+int	nbands			# number of bands
+int	nrows			# number of rows
+PIXEL matfac[nbands, nrows]	# Cholesky factorization
+int	ier			# error code
+
+int	i, n, j, imax, jmax
+PIXEL	ratio
+
+begin
+	if (nrows == 1) {
+	    if (matrix[1,1] > 0.)
+	        matfac[1,1] = 1. / matrix[1,1]
+	    return
+	}
+
+	# copy matrix into matfac
+	do n = 1, nrows {
+	    do j = 1, nbands
+		matfac[j,n] = matrix[j,n]
+	}
+
+	do n = 1, nrows {
+
+	    # test to see if matrix is singular
+	    if (((matfac[1,n]+matrix[1,n])-matrix[1,n]) <= 1000/MAX_PIXEL) {
+		do j = 1, nbands
+		    matfac[j,n] = 0.
+		ier = SINGULAR
+		next
+	    }
+
+	    matfac[1,n] = 1. / matfac[1,n]
+	    imax = min (nbands - 1, nrows - n)
+	    if (imax < 1)
+		next
+
+	    jmax = imax
+	    do i = 1, imax {
+		ratio = matfac[i+1,n] * matfac[1,n]
+		do j = 1, jmax
+		    matfac[j,n+i] = matfac[j,n+i] - matfac[j+i,n] * ratio
+		jmax = jmax - 1
+		matfac[i+1,n] = ratio
+	    }
+	}
+end
+
+
+# GSCHOSLV -- Solve the matrix whose Cholesky factorization was calculated in
+# GSCHOFAC for the coefficients. This routine was adapted from bchslv.f
+# described in "A Practical Guide to Splines", by Carl de Boor (1978).
+
+procedure $tgschoslv (matfac, nbands, nrows, vector, coeff)
+
+PIXEL matfac[nbands,nrows] 		# Cholesky factorization
+int	nbands				# number of bands
+int	nrows				# number of rows
+PIXEL vector[nrows]			# right side of matrix equation
+PIXEL coeff[nrows]			# coefficients
+
+int	i, n, j, jmax, nbndm1
+
+begin
+	if (nrows == 1) {
+	    coeff[1] = vector[1] * matfac[1,1]
+	    return
+	}
+
+	# copy vector to coefficients
+	do i = 1, nrows
+	    coeff[i] = vector[i]
+
+	# forward substitution
+	nbndm1 = nbands - 1
+	do n = 1, nrows {
+	    jmax = min (nbndm1, nrows - n)
+	    if (jmax >= 1) {
+	        do j = 1, jmax
+		    coeff[j+n] = coeff[j+n] - matfac[j+1,n] * coeff[n]
+	    }
+	}
+
+	# back substitution
+	for (n = nrows; n >= 1; n = n - 1) {
+	    coeff[n] = coeff[n] * matfac[1,n]
+	    jmax = min (nbndm1, nrows - n)
+	    if (jmax >= 1) {
+		do j = 1, jmax
+		    coeff[n] = coeff[n] - matfac[j+1,n] * coeff[j+n]
+	    }
+	}
+end