Initial commit

1 files changed, 59 insertions, 0 deletions

diff --git a/math/nlfit/nliter.gx b/math/nlfit/nliter.gx
new file mode 100644
index 00000000..7857f852
--- /dev/null
+++ b/math/nlfit/nliter.gx
@@ -0,0 +1,59 @@
+include <mach.h>
+include	<math/nlfit.h>
+$if (datatype == r)
+include "nlfitdefr.h"
+$else
+include "nlfitdefd.h"
+$endif
+
+# NLITER - Routine to compute one iteration of the fitting process
+
+procedure nliter$t (nl, x, z, w, npts, nvars, ier)
+
+pointer	nl		# pointer to nl fitting structure
+PIXEL	x[ARB]		# independent variables (npts * nvars)
+PIXEL	z[ARB]		# function values (npts)
+PIXEL	w[ARB]		# weights (npts)
+int	npts		# number of points
+int	nvars		# number of independent variables
+int	ier		# error code
+
+int	i, index
+PIXEL	nlacpts$t(), nlresid$t()
+
+begin
+	# Do the initial accumulation.
+   	NL_OLDSQ(nl) = nlacpts$t (nl, x, z, w, npts, nvars)
+
+	# Set up temporary parameter array.
+	call amov$t (PARAM(NL_PARAM(nl)), TRY(NL_TRY(nl)), NL_NPARAMS(nl))
+
+	repeat {
+
+	    # Solve the matrix.
+	    call nlsolve$t (nl, ier)
+	    if (ier == NO_DEG_FREEDOM)
+	        return
+
+	    # Increment the parameters.
+	    do i = 1, NL_NFPARAMS(nl) {
+	        index = PLIST(NL_PLIST(nl)+i-1)
+	        TRY(NL_TRY(nl)+index-1) = TRY(NL_TRY(nl)+index-1) +
+		    DPARAM(NL_DPARAM(nl)+i-1)
+	    }
+
+	    # Reaccumulate the residuals and increment lambda.
+	    NL_SUMSQ(nl) = nlresid$t (nl, x, z, w, npts, nvars)
+	    #if (NL_OLDSQ(nl) > (NL_SUMSQ(nl) + NL_TOL(nl) * NL_SUMSQ(nl))) {
+	    if (NL_OLDSQ(nl) >= NL_SUMSQ(nl)) {
+		call amov$t (TRY(NL_TRY(nl)), PARAM(NL_PARAM(nl)),
+		    NL_NPARAMS(nl))
+		NL_LAMBDA(nl) = PIXEL (0.10) * NL_LAMBDA(nl)
+		break
+	    } else
+	        NL_LAMBDA(nl) = min (PIXEL(LAMBDAMAX),
+		    PIXEL (10.0) * NL_LAMBDA(nl))
+
+	} until (NL_LAMBDA(nl) <= PIXEL(0.0) ||
+	    NL_LAMBDA(nl) >= PIXEL(LAMBDAMAX))
+end