◆ N2g1_r()

Sub N2g1_r	(	M As	Long,
		N As	Long,
		X() As	Double,
		Info As	Long,
		YY() As	Double,
		YYp() As	Double,
		IRev As	Long,
		Optional NFcall As	Long,
		Optional NFjcall As	Long,
		Optional Niter As	Long
	)

Nonlinear least squares approximation (adaptive algorithm) (reverse communication version)

Purpose: This routine minimizes the sum of the squares of M nonlinear functions in N variables by the adaptive algorithm which combines and augments a Gauss-Newton, Levenberg-Marquardt and other techniques for better convergence.
minimize the sum of fi(x1, x2, ..., xn)^2 (sum for i = 1 to M)

The user must provide the function and Jacobian values in accordance with IRev.

N2g1_r is equivalent to using N2g_r with setting default parameters. Covariance matrix and regression diagnostic vector are not computed.

Parameters

[in]	M	Number of data. (M > 0)
[in]	N	Number of parameters. (0 < N <= M)
[in,out]	X()	Array X(LX - 1) (LX >= N) [in] X must contain an initial estimate of the solution vector. [out] IRev = 0: Solution vector, i.e. best point so far found. IRev = 1, 2: The abscissa where the function value or derivatives shoule be evaluated.
[out]	Info	= 0: Successful exit. = -1: The argument M had an illegal value. (M < N) = -2: The argument N had an illegal value. (N < 1) = -3: The argument X() is invalid. = -5: The argument YY() is invalid. = -6: The argument YYp() is invalid. = 7: Singular convergence. (Hessian near the current iterate appears to be singular) = 8: False convergence. (Iterate appears to be converging to a noncritical point. Tolerances may be too small) = 9: Function evaluation limit reached. = 10: Iteration limit reached. = 63: F(X) cannot be computed at the initial X. = 65: The gradient could not be computed at X.
[in]	YY()	Array YYp(LYY - 1) (LYY >= M) When returned with IRev = 1, the function value at X() should be given in YY() in the next call.
[in]	YYp()	Array YYp(LYYp1 - 1, LYYp2 - 1) (LYYp1 >= M, LYYp2 >= N) When returned with IRev = 2, Jacobian matrix at X() should be given in YYp() in the next call.
[in,out]	IRev	Control variable for reverse communication. [in] Before first call, IRev should be initialized to zero. On succeeding calls, IRev should not be altered. [out] If IRev is not zero, set the required values to the specified variables or display the intermediate result as follows and call the routine again. = 0: Computation finished. See return code in Info. = 1: User should set the function values at X() in YY(). Do not alter any variables other than YY(). = 2: User should set the derivatives at X() in YYp(). Do not alter any variables other than YYp().
[out]	NFcall	(Optional) Number of function evaluations with IRev = 1.
[out]	NFjcall	(Optional) Number of Jacobian evaluations with IRev = 2.
[out]	Niter	(Optional) Number of iterations.

Reference: netlib/port

Example Program: Approximate the following data with model function f(x) = c1*(1 - exp(-c2*x)). Determine two parameters c1 and c2 by the nonlinear least squares method.
f(x) x

10.07 77.6

29.61 239.9

50.76 434.8

81.78 760.0

The initial estimates of the solution are c1 = 500 and c2 = 0.0001.
Sub FN2g(M As Long, N As Long, X() As Double, Nf As Long, Fvec() As Double)

Dim Xdata(3) As Double, Ydata(3) As Double, I As Long

Ydata(0) = 10.07: Xdata(0) = 77.6

Ydata(1) = 29.61: Xdata(1) = 239.9

Ydata(2) = 50.76: Xdata(2) = 434.8

Ydata(3) = 81.78: Xdata(3) = 760

For I = 0 To M - 1

Fvec(I) = Ydata(I) - X(0) * (1 - Exp(-Xdata(I) * X(1)))

Next

End Sub

Sub JN2g(M As Long, N As Long, X() As Double, Nf As Long, Fjac() As Double)

Dim Xdata(3) As Double, Ydata(3) As Double, I As Long

Ydata(0) = 10.07: Xdata(0) = 77.6

Ydata(1) = 29.61: Xdata(1) = 239.9

Ydata(2) = 50.76: Xdata(2) = 434.8

Ydata(3) = 81.78: Xdata(3) = 760

For I = 0 To M - 1

Fjac(I, 0) = Exp(-Xdata(I) * X(1)) - 1

Fjac(I, 1) = -Xdata(I) * X(0) * Exp(-X(1) * Xdata(I))

Next

End Sub

Sub Ex_N2g1_r()

Const M = 4, N = 2

Dim X(N - 1) As Double, Info As Long

Dim YY(M - 1) As Double, YYp(M - 1, N - 1) As Double, IRev As Long

X(0) = 500: X(1) = 0.0001

IRev = 0

Do

Call N2g1_r(M, N, X(), Info, YY(), YYp(), IRev)

If IRev = 1 Then

Call FN2g(M, N, X(), 0, YY())

ElseIf IRev = 2 Then

Call JN2g(M, N, X(), 0, YYp())

End If

Loop While IRev <> 0

Debug.Print "C1, C2 =", X(0), X(1)

Debug.Print "Info =", Info

End Sub

Example Results: C1, C2 = 241.084896112856 5.44942234058364E-04

Info = 0