/*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 */