/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:47 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include <math.h>
#include "fcrt.h"
#include "snlsfu.h"
#include <stdlib.h>
		/* PARAMETER translations */
#define	AMAT	113
#define	COVREQ	15
#define	D	27
#define	DAMAT	114
#define	DLTFDJ	43
#define	GPTR	117
#define	GRP	118
#define	IN	112
#define	IVNEED	3
#define	L1SAV	111
#define	MAXGRP	116
#define	MODE	35
#define	MSAVE	115
#define	NEXTIV	46
#define	NEXTV	47
#define	NFCALL	6
#define	NFGCAL	7
#define	PERM	58
#define	RESTOR	9
#define	TOOBIG	2
#define	VNEED	4
#define	XSAVE	119
		/* end of PARAMETER translations */
 
void /*FUNCTION*/ snlsfu(
long n,
long p,
long l,
float alf[],
float c[],
float y[],
void (*scalca)(long,long,long,float[],long*,float[]),
long *inc,
long iinc,
long iv[],
long liv,
long lv,
float v[])
{
#define INC(I_,J_)	(*(inc+(I_)*(iinc)+(J_)))
	LOGICAL32 partj;
	long int a0, a1, aj, alp1, bwa1, d0, da0, da1, daj, gptr1, grp1,
	 grp2, i, i1, in0, in1, in2, ini, inlen, ipntr1, iv1, iwa1, iwalen,
	 j1, jn1, jpntr1, k, l1, lp1, m, m0, nf, ng, ngrp0, ngrp1, ngrp2,
	 rsave0, rsave1, rsvlen, x0i, xsave0, xsave1;
	float delta, di, h, xi;
	static float negone = -1.e0;
	static float one = 1.e0;
	static float zero = 0.e0;
		/* OFFSET Vectors w/subscript range: 1 to dimension */
	float *const Alf = &alf[0] - 1;
	float *const C = &c[0] - 1;
	long *const Iv = &iv[0] - 1;
	float *const V = &v[0] - 1;
	float *const Y = &y[0] - 1;
		/* end of OFFSET VECTORS */
 
	/* Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
	 * ALL RIGHTS RESERVED.
	 * Based on Government Sponsored Research NAS7-03001.
	 *>> 1996-04-27 SNLSFU Krogh  Changes to get desired C prototypes.
	 *>> 1994-10-20 SNLSFU Krogh  Changes to use M77CON
	 *>> 1990-07-02 SNLSFU CLL @ JPL
	 *>> 1990-06-12 CLL @ JPL
	 *>> 1990-02-20 CLL @ JPL
	 *** from netlib, Fri Feb 16 15:19:34 EST 1990 ***
	 *--S replaces "?": ?NLSFU,?CALCA,?IVSET,?RNSG,?V2AXY,?V7CPY,?V7SCL
	 *
	 *  ***  SOLVE SEPARABLE NONLINEAR LEAST SQUARES USING
	 *  ***  FINITE-DIFFERENCE DERIVATIVES.
	 *
	 *  ***  PARAMETER DECLARATIONS  ***
	 * */
 
	/*  ***  PARAMETERS  ***
	 *
	 *      N (IN)  NUMBER OF OBSERVATIONS.
	 *      P (IN)  NUMBER OF NONLINEAR PARAMETERS TO BE ESTIMATED.
	 *      L (IN)  NUMBER OF LINEAR PARAMETERS TO BE ESTIMATED.
	 *    ALF (I/O) NONLINEAR PARAMETERS.
	 *                 INPUT = INITIAL GUESS,
	 *                 OUTPUT = BEST ESTIMATE FOUND.
	 *      C (OUT) LINEAR PARAMETERS (ESTIMATED).
	 *      Y (IN)  RIGHT-HAND SIDE VECTOR.
	 *  SCALCA (IN)  SUBROUTINE TO COMPUTE A MATRIX.
	 *    INC (IN)  INCIDENCE MATRIX OF DEPENDENCIES OF COLUMNS OF A ON
	 *                 COMPONENTS OF ALF -- INC(I,J) = 1 MEANS COLUMN I
	 *                 OF A DEPENDS ON ALF(J).
	 *   IINC (IN)  DECLARED LEAD DIMENSION OF INC.  MUST BE AT LEAST L+1.
	 *     IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR.
	 *    LIV (IN)  LENGTH OF IV.  MUST BE AT LEAST
	 *                 122 + 2*M + 4*P + 2*L + MAX(L+1,6*P), WHERE  M  IS
	 *                 THE NUMBER OF ONES IN INC.
	 *     LV (IN)  LENGTH OF V.  MUST BE AT LEAST
	 *                 105 + 2*N*(L+3) + JLEN + L*(L+3)/2 + P*(2*P + 18),
	 *                 WHERE  JLEN = (L+P)*(N+L+P+1),  UNLESS NEITHER A
	 *                 COVARIANCE MATRIX NOR REGRESSION DIAGNOSTICS ARE
	 *                 REQUESTED, IN WHICH CASE  JLEN = N*P.
	 *                 If row L+1 of INC() contains only zeros, meaning
	 *                 the term PHI sub (L+1) is absent from the model,
	 *                 then LV can be 4*N less than just described.
	 *      V (I/O) FLOATING-POINT PARAMETER AND SCRATCH VECTOR.
	 *
	 *
	 * -------------------------  DECLARATIONS  ----------------------------
	 *
	 *
	 *  ***  EXTERNAL SUBROUTINES  ***
	 * */
 
	/* SIVSET.... PROVIDES DEFAULT IV AND V VALUES.
	 * IDSM...... DETERMINES EFFICIENT ORDER FOR FINITE DIFFERENCES.
	 *  SRNSG... CARRIES OUT NL2SOL ALGORITHM.
	 * SV2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER.
	 * SV7CPY.... COPIES ONE VECTOR TO ANOTHER.
	 * SV7SCL... SCALES AND COPIES ONE VECTOR TO ANOTHER.
	 *
	 *  ***  LOCAL VARIABLES  ***
	 * */
 
	/*  ***  SUBSCRIPTS FOR IV AND V  ***
	 * */
 
	/*  ***  IV SUBSCRIPT VALUES  ***
	 * */
 
	/*++++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++
	 * */
	lp1 = l + 1;
	if (Iv[1] == 0)
		sivset( 1, iv, liv, lv, v );
	if ((p <= 0 || l < 0) || iinc <= l)
		goto L_80;
	iv1 = Iv[1];
	if (iv1 == 14)
		goto L_120;
	if (iv1 > 2 && iv1 < 12)
		goto L_120;
	if (iv1 == 12)
		Iv[1] = 13;
	if (Iv[1] != 13)
		goto L_50;
 
	/*  ***  FRESH START ***
	 * */
	if (Iv[PERM] <= XSAVE)
		Iv[PERM] = XSAVE + 1;
 
	/*  ***  CHECK INC, COUNT ITS NONZEROS
	 * */
	l1 = 0;
	m = 0;
	for (i = 1; i <= p; i++)
	{
		m0 = m;
		if (l == 0)
			goto L_20;
		for (k = 1; k <= l; k++)
		{
			if (INC(i - 1,k - 1) < 0 || INC(i - 1,k - 1) > 1)
				goto L_80;
			if (INC(i - 1,k - 1) == 1)
				m += 1;
		}
L_20:
		if (INC(i - 1,lp1 - 1) != 1)
			goto L_30;
		m += 1;
		l1 = 1;
L_30:
		if ((m == m0 || INC(i - 1,lp1 - 1) < 0) || INC(i - 1,lp1 - 1) >
		 1)
			goto L_80;
	}
 
	/*     *** NOW L1 = 1 MEANS A HAS COLUMN L+1 ***
	 *
	 *     *** COMPUTE STORAGE REQUIREMENTS ***
	 * */
	iwalen = max( lp1, 6*p );
	inlen = 2*m;
	Iv[IVNEED] += inlen + 3*p + l + iwalen + 3;
	rsvlen = 2*l1*n;
	l1 += l;
	Iv[VNEED] += 2*n*l1 + rsvlen + p;
 
L_50:
	srnsg( v, alf, c, v, (long(*)[2])iv, iv, l, 1, n, liv, lv, n,
	 m, p, v, y );
	if (Iv[1] != 14)
		goto L_999;
 
	/*  ***  STORAGE ALLOCATION  ***
	 * */
	Iv[IN] = Iv[NEXTIV];
	Iv[AMAT] = Iv[NEXTV];
	Iv[DAMAT] = Iv[AMAT] + n*l1;
	Iv[XSAVE] = Iv[DAMAT] + n*l1;
	Iv[NEXTV] = Iv[XSAVE] + p + rsvlen;
	Iv[L1SAV] = l1;
	Iv[MSAVE] = m;
 
	/*  ***  DETERMINE HOW MANY GROUPS FOR FINITE DIFFERENCES
	 *  ***  (SET UP TO CALL IDSM)
	 * */
	in1 = Iv[IN];
	jn1 = in1 + m;
	for (k = 1; k <= p; k++)
	{
		for (i = 1; i <= lp1; i++)
		{
			if (INC(k - 1,i - 1) == 0)
				goto L_60;
			Iv[in1] = i;
			in1 += 1;
			Iv[jn1] = k;
			jn1 += 1;
L_60:
			;
		}
	}
	in1 = Iv[IN];
	jn1 = in1 + m;
	iwa1 = in1 + inlen;
	ngrp1 = iwa1 + iwalen;
	bwa1 = ngrp1 + p;
	ipntr1 = bwa1 + p;
	jpntr1 = ipntr1 + l + 2;
	idsm( lp1, p, m, &Iv[in1], &Iv[jn1], &Iv[ngrp1], &ng, &k, &i,
	 &Iv[ipntr1], &Iv[jpntr1], &Iv[iwa1], iwalen, &Iv[bwa1] );
	if (i == 1)
		goto L_90;
	Iv[1] = 69;
	goto L_50;
L_80:
	Iv[1] = 66;
	goto L_50;
 
	/*  ***  SET UP GRP AND GPTR ARRAYS FOR COMPUTING FINITE DIFFERENCES
	 *
	 *  ***  THERE ARE NG GROUPS.  GROUP I CONTAINS ALF(GRP(J)) FOR
	 *  ***  GPTR(I) .LE. J .LE. GPTR(I+1)-1.
	 * */
L_90:
	Iv[MAXGRP] = ng;
	Iv[GPTR] = in1 + 2*l1;
	gptr1 = Iv[GPTR];
	Iv[GRP] = gptr1 + ng + 1;
	Iv[NEXTIV] = Iv[GRP] + p;
	grp1 = Iv[GRP];
	ngrp0 = ngrp1 - 1;
	ngrp2 = ngrp0 + p;
	for (i = 1; i <= ng; i++)
	{
		Iv[gptr1] = grp1;
		gptr1 += 1;
		for (i1 = ngrp1; i1 <= ngrp2; i1++)
		{
			if (Iv[i1] != i)
				goto L_100;
			Iv[grp1] = i1 - ngrp0;
			grp1 += 1;
L_100:
			;
		}
	}
	Iv[gptr1] = grp1;
	if (iv1 == 13)
		goto L_999;
 
	/*  ***  INITIALIZE POINTERS  ***
	 * */
L_120:
	a1 = Iv[AMAT];
	a0 = a1 - n;
	da1 = Iv[DAMAT];
	da0 = da1 - n;
	in1 = Iv[IN];
	in0 = in1 - 2;
	l1 = Iv[L1SAV];
	in2 = in1 + 2*l1 - 1;
	d0 = Iv[D] - 1;
	ng = Iv[MAXGRP];
	xsave1 = Iv[XSAVE];
	xsave0 = xsave1 - 1;
	rsave1 = xsave1 + p;
	rsave0 = rsave1 + n;
	alp1 = a1 + l*n;
	delta = V[DLTFDJ];
	Iv[COVREQ] = -labs( Iv[COVREQ] );
 
L_130:
	srnsg( &V[a1], alf, c, &V[da1], (long(*)[2])&Iv[in1], iv, l, l1,
	 n, liv, lv, n, l1, p, v, y );
	switch (IARITHIF(Iv[1] - 2))
	{
		case -1: goto L_140;
		case  0: goto L_150;
		case  1: goto L_999;
	}
 
	/*  ***  NEW FUNCTION VALUE (R VALUE) NEEDED  ***
	 * */
L_140:
	nf = Iv[NFCALL];
     (*scalca)( n, p, l, alf, &nf, &V[a1] );
	if (nf <= 0)
		Iv[TOOBIG] = 1;
	/*      CALL SCALCA(N, P, L, ALF, NF, V(A1)) */
	if (l1 <= l)
		goto L_130;
	if (Iv[RESTOR] == 2)
		sv7cpy( n, &V[rsave0], &V[rsave1] );
	sv7cpy( n, &V[rsave1], &V[alp1] );
	goto L_130;
 
	/*  ***  COMPUTE DR = GRADIENT OF R COMPONENTS  ***
	 * */
L_150:
	if (l1 > l && Iv[NFGCAL] == Iv[NFCALL])
		sv7cpy( n, &V[rsave0], &V[rsave1] );
	gptr1 = Iv[GPTR];
	for (k = 1; k <= ng; k++)
	{
		sv7cpy( p, &V[xsave1], alf );
		grp1 = Iv[gptr1];
		grp2 = Iv[gptr1 + 1] - 1;
		gptr1 += 1;
		for (i1 = grp1; i1 <= grp2; i1++)
		{
			i = Iv[i1];
			xi = Alf[i];
			j1 = d0 + i;
			di = V[j1];
			if (di <= zero)
				di = one;
			h = delta*fmaxf( fabsf( xi ), one/di );
			if (xi < zero)
				h = -h;
			x0i = xsave0 + i;
			V[x0i] = xi + h;
		}
        (*scalca)( n, p, l, &V[xsave1], &Iv[NFGCAL], &V[da1] );
		if (Iv[NFGCAL] > 0)
			goto L_170;
		/*         CALL SCALCA(N, P, L, V(XSAVE1), IV(NFGCAL), V(DA1)) */
		Iv[TOOBIG] = 1;
		goto L_130;
L_170:
		jn1 = in1;
		for (i = in1; i <= in2; i++)
		{
			Iv[i] = 0;
		}
		partj = Iv[MODE] <= p;
		for (i1 = grp1; i1 <= grp2; i1++)
		{
			i = Iv[i1];
			for (j1 = 1; j1 <= l1; j1++)
			{
				if (INC(i - 1,j1 - 1) == 0)
					goto L_210;
				ini = in0 + 2*j1;
				Iv[ini] = i;
				Iv[ini + 1] = j1;
				x0i = xsave0 + i;
				h = one/(V[x0i] - Alf[i]);
				daj = da0 + j1*n;
				if (partj)
					goto L_190;
				/*                 *** FULL FINITE DIFFERENCE FOR COV. AND REG. DIAG. *** */
				aj = a0 + j1*n;
				sv2axy( n, &V[daj], negone, &V[aj], &V[daj] );
				goto L_200;
L_190:
				if (j1 > l)
					sv2axy( n, &V[daj], negone, &V[rsave0], &V[daj] );
L_200:
				sv7scl( n, &V[daj], h, &V[daj] );
L_210:
				;
			}
		}
		if (k >= ng)
			goto L_240;
		Iv[1] = -2;
		srnsg( &V[a1], alf, c, &V[da1], (long(*)[2])&Iv[in1], iv,
		 l, l1, n, liv, lv, n, l1, p, v, y );
		if (-2 != Iv[1])
			goto L_999;
	}
L_240:
	Iv[1] = 2;
	goto L_130;
 
L_999:
	return;
 
	/*  ***  LAST CARD OF   SNLSFU FOLLOWS  *** */
#undef	INC
} /* end of function */