/*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 "snlsgb.h"
#include <stdlib.h>
/* PARAMETER translations */
#define AMAT 113
#define DAMAT 114
#define IN 112
#define IVNEED 3
#define L1SAV 111
#define MSAVE 115
#define NEXTIV 46
#define NEXTV 47
#define NFCALL 6
#define NFGCAL 7
#define PERM 58
#define TOOBIG 2
#define VNEED 4
/* end of PARAMETER translations */
void /*FUNCTION*/ snlsgb(
long n,
long p,
long l,
float alf[],
float b[][2],
float c[],
float y[],
void (*scalca)(long,long,long,float[],long*,float[]),
void (*scalcb)(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_)))
long int a1, da1, i, in1, iv1, k, l1, lp1, m, m0, nf;
/* 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 SNLSGB Krogh Changes to get desired C prototypes.
*>> 1994-10-20 SNLSGB Krogh Changes to use M77CON
*>> 1990-07-02 SNLSGB CLL @ JPL
*>> 1990-06-12 CLL @ JPL
*>> 1990-02-14 CLL @ JPL
*** from netlib, Fri Feb 9 13:10:09 EST 1990 ***
*--S replaces "?": ?NLSGB ,?CALCA,?CALCB,?IVSET,?RNSGB
*
* *** SOLVE SEPARABLE NONLINEAR LEAST SQUARES USING ***
* *** ANALYTICALLY COMPUTED DERIVATIVES. ***
*
* *** PARAMETER DECLARATIONS ***
* */
/* *** PURPOSE ***
*
* GIVEN A SET OF N OBSERVATIONS Y(1)....Y(N) OF A DEPENDENT VARIABLE
* T(1)...T(N), SNLSGB ATTEMPTS TO COMPUTE A LEAST SQUARES FIT
* TO A FUNCTION ETA (THE MODEL) WHICH IS A LINEAR COMBINATION
*
* L
* ETA(C,ALF,T) = SUM C * PHI(ALF,T) +PHI (ALF,T)
* J=1 J J L+1
*
* OF NONLINEAR FUNCTIONS PHI(J) DEPENDENT ON T AND ALF(1),...,ALF(P)
* (.E.G. A SUM OF EXPONENTIALS OR GAUSSIANS). THAT IS, IT DETERMINES
* NONLINEAR PARAMETERS ALF WHICH MINIMIZE
*
* 2 N 2
* NORM(RESIDUAL) = SUM (Y - ETA(C,ALF,T )) ,
* I=1 I I
*
* SUBJECT TO THE SIMPLE BOUND CONSTRAINTS B(1,I) .LE. X(I) .LE. B(2,I),
* I = 1(1)P.
*
* THE (L+1)ST TERM IS OPTIONAL.
*
* ------------------------- PARAMETER USAGE -------------------------
*
* INPUT PARAMETERS
*
* N INTEGER NUMBER OF OBSERVATIONS (MUST BE .GE. MAX(L,P)).
*
* P INTEGER NUMBER OF NONLINEAR PARAMETERS (MUST BE .GE. 1).
*
* L INTEGER NUMBER OF LINEAR PARAMETERS (MUST BE .GE. 0).
*
* ALF D.P. ARRAY P VECTOR = INITIAL ESTIMATE OF THE NONLINEAR
* PARAMETERS.
*
* SCALCA SUBROUTINE USER PROVIDED FUNCTION TO CALCULATE THE MODEL
* (I.E., TO CALCULATE PHI) -- SEE THE NOTE BELOW
* ON THE CALLING SEQUENCE FOR SCALCA.
* SCALCA MUST BE DECLARED EXTERNAL IN THE CALLING
* PROGRAM.
*
* SCALCB SUBROUTINE USER PROVIDED FUNCTION TO CALCULATE THE DERIVA-
* TIVE OF THE MODEL (I.E., OF PHI) WITH RESPECT TO
* ALF -- SEE THE NOTE BELOW ON THE CALLING
* SEQUENCE FOR SCALCB. SCALCB MUST BE DECLARED
* EXTERNAL IN THE CALLING PROGRAM.
*
* Y D.P. ARRAY VECTOR OF OBSERVATIONS.
*
* INC INTEGER ARRAY A 2 DIM. ARRAY OF DIMENSION AT LEAST (L+1,P)
* INDICATING THE POSITION OF THE NONLINEAR PARA-
* METERS IN THE MODEL. SET INC(J,K) = 1 IF ALF(K)
* APPEARS IN PHI(J). OTHERWISE SET INC(J,K) = 0.
* IF PHI((L+1)) IS NOT IN THE MODEL, SET THE L+1ST
* ROW OF INC TO ALL ZEROS. EVERY COLUMN OF INC
* MUST CONTAIN AT LEAST ONE 1.
*
* IINC INTEGER DECLARED ROW DIMENSION OF INC, WHICH MUST BE AT
* LEAST L+1.
*
* IV INTEGER ARRAY OF LENGTH AT LEAST LIV THAT CONTAINS
* VARIOUS PARAMETERS FOR THE SUBROUTINE, SUCH AS
* THE ITERATION AND FUNCTION EVALUATION LIMITS AND
* SWITCHES THAT CONTROL PRINTING. THE INPUT COM-
* PONENTS OF IV ARE DESCRIBED IN DETAIL IN THE
* PORT OPTIMIZATION DOCUMENTATION.
* IF IV(1)=0 ON INPUT, THEN DEFAULT PARAMETERS
* ARE SUPPLIED TO IV AND V. THE CALLER MAY SUPPLY
* NONDEFAULT PARAMETERS TO IV AND V BY EXECUTING A
* CALL SIVSET(1,IV,LIV,LV,V) AND THEN ASSIGNING
* NONDEFAULT VALUES TO THE APPROPRIATE COMPONENTS
* OF IV AND V BEFORE CALLING SNLSGB.
*
* LIV INTEGER LENGTH OF IV. MUST BE AT LEAST
* 115 + 4*P + L + 2*M,
* WHERE M IS THE NUMBER OF ONES IN INC.
*
* LV INTEGER LENGTH OF V. MUST BE AT LEAST
* 105 + N*(L+M+P+3) + L*(L+3)/2 + P*(2*P+21),
* WHERE M IS AS FOR LIV (SEE ABOVE).
* 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 N less than just described.
*
* V D.P. ARRAY WORK AND PARAMETER ARRAY OF LENGTH AT LEAST LV
* THAT CONTAINS SUCH INPUT COMPONENTS AS THE
* CONVERGENCE TOLERANCES. THE INPUT COMPONENTS OF
* V MAY BE SUPPLIED AS FOR IV (SEE ABOVE). NOTE
* THAT V(35) CONTAINS THE INITIAL STEP BOUND,
* WHICH, IF TOO LARGE, MAY LEAD TO OVERFLOW.
*
*
* OUTPUT PARAMETERS
*
* ALF D.P. ARRAY FINAL NONLINEAR PARAMETERS.
*
* C D.P. ARRAY L VECTOR OF LINEAR PARAMETERS -- NOTE THAT NO
* INITIAL GUESS FOR C IS REQUIRED. */
/* IV IV(1) CONTAINS A RETURN CODE DESCRIBED IN THE
* PORT OPTIMIZATION DOCUMENTATION. IF IV(1) LIES
* BETWEEN 3 AND 7, THEN THE ALGORITHM HAS
* CONVERGED (BUT IV(1) = 7 INDICATES POSSIBLE
* TROUBLE WITH THE MODEL). IV(1) = 9 OR 10 MEANS
* FUNCTION EVALUATION OR ITERATION LIMIT REACHED.
* IV(1) = 66 MEANS BAD PARAMETERS (INCLUDING A
* COLUMN OF ZEROS IN INC). NOTE THAT THE
* ALGORITHM CAN BE RESTARTED AFTER ANY RETURN WITH
* IV(1) .LT. 12 -- SEE THE PORT DOCUMENTATION.
*
* V VARIOUS ITEMS OF INTEREST, INCLUDING THE NORM OF
* THE GRADIENT(1) AND THE FUNCTION VALUE(10). SEE
* THE PORT DOCUMENTATION FOR A COMPLETE LIST.
*
*
*
* PARAMETERS FOR SCALCA(N,P,L,ALF,NF,PHI)
*
* N,L,P,ALF ARE INPUT PARAMETERS AS DESCRIBED ABOVE
*
* PHI D.P. ARRAY N*(L+1) or N*L array. SCALCA must store into
* PHI(i,j) the value of of the nonlinear function
* PHI sub j evaluated with the current values of
* ALF() and at the ith point of the data set.
* Note that if the model does not have the
* term PHI sub L+1 then one must not store to the
* (L+1)st column of PHI() as this would overwrite
* other work space.
*
* NF INTEGER CURRENT INVOCATION COUNT FOR SCALCA. IF PHI CANNOT
* BE EVALUATED AT ALF (E.G. BECAUSE AN ARGUMENT TO
* AN INTRINSIC FUNCTION IS OUT OF RANGE), THEN SCALCA
* SHOULD SIMPLY SET NF TO 0 AND RETURN. THIS
* TELLS THE ALGORITHM TO TRY A SMALLER STEP.
*
* N.B. THE DEPENDENT VARIABLE T IS NOT EXPLICITLY PASSED. IF REQUIRED,
* IT MAY BE PASSED IN NAMED COMMON.
*
*
* PARAMETERS FOR SCALCB(N,P,L,ALF,NF,DER)
*
* N,P,L,ALF,NF ARE AS FOR SCALCA
*
* DER D.P. ARRAY N*M ARRAY, WHERE M IS THE NUMBER OF ONES IN INC.
* SCALCB MUST SET DER TO THE DERIVATIVES OF THE MODEL
* WITH RESPECT TO ALF. IF THE MODEL HAS K TERMS THAT
* DEPEND ON ALF(I), THEN DER WILL HAVE K CONSECUTIVE
* COLUMNS OF DERIVATIVES WITH RESPECT TO ALF(I). THE
* COLUMNS OF DER CORRESPOND TO THE ONES IN INC WHEN
* ONE TRAVELS THROUGH INC BY COLUMNS. FOR EXAMPLE,
* IF INC HAS THE FORM...
* 1 1 0
* 0 1 0
* 1 0 0
* 0 0 1
* THEN THE FIRST TWO COLUMNS OF DER ARE FOR THE
* DERIVATIVES OF COLUMNS 1 AND 3 OF PHI WITH RESPECT
* TO ALF(1), COLUMNS 3 AND 4 OF DER ARE FOR THE
* DERIVATIVES OF COLUMNS 1 AND 2 OF PHI WITH RESPECT
* TO ALF(2), AND COLUMN 5 OF DER IS FOR THE DERIVA-
* TIVE OF COLUMN 4 OF PHI WITH RESPECT TO ALF(3).
* MORE SPECIFICALLY, DER(I,2) IS FOR THE DERIVATIVE
* OF PHI(I,3) WITH RESPECT TO ALF(1) AND DER(I,5) IS
* FOR THE DERIVATIVE OF PHI(I,4) WITH RESPECT TO
* ALF(3) (FOR I = 1,2,...,N).
* THE VALUE OF ALF PASSED TO SCALCB IS THE SAME AS
* THAT PASSED TO SCALCA THE LAST TIME IT WAS CALLED.
* (IF DER CANNOT BE EVALUATED, THEN SCALCB SHOULD SET
* NF TO 0. THIS WILL CAUSE AN ERROR RETURN.)
*
* N.B. DER IS FOR DERIVATIVES WITH RESPECT TO ALF, NOT T.
*
* ----------------------------- NOTES -------------------------------
*
* THIS PROGRAM WAS WRITTEN BY LINDA KAUFMAN AT BELL LABS, MURRAY
* HILL, N.J. IN 1977 AND EXTENSIVELY REVISED BY HER AND DAVID GAY IN
* 1980, 1981, 1983, 1984. THE WORK OF DAVID GAY WAS SUPPORTED IN PART
* BY NATIONAL SCIENCE FOUNDATION GRANT MCS-7906671.
*
* ------------------------- DECLARATIONS ----------------------------
*
*
* *** EXTERNAL SUBROUTINES ***
* */
/* SIVSET.... PROVIDES DEFAULT IV AND V VALUES.
* SRNSGB... CARRIES OUT NL2SOL ALGORITHM.
*
* *** LOCAL VARIABLES ***
* */
/* *** SUBSCRIPTS FOR IV AND V ***
* */
/* *** IV SUBSCRIPT VALUES ***
* */
/*++++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++
* */
if (Iv[1] == 0)
sivset( 1, iv, liv, lv, v );
if ((p <= 0 || l < 0) || iinc <= l)
goto L_50;
iv1 = Iv[1];
if (iv1 == 14)
goto L_90;
if (iv1 > 2 && iv1 < 12)
goto L_90;
if (iv1 == 12)
Iv[1] = 13;
if (Iv[1] != 13)
goto L_60;
if (Iv[PERM] <= MSAVE)
Iv[PERM] = MSAVE + 1;
lp1 = l + 1;
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_50;
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_50;
}
Iv[IVNEED] += 2*m;
l1 += l;
Iv[VNEED] += n*(l1 + m);
goto L_60;
L_50:
Iv[1] = 66;
L_60:
srnsgb( v, alf, b, 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[NEXTIV] = Iv[IN] + 2*m;
Iv[AMAT] = Iv[NEXTV];
Iv[DAMAT] = Iv[AMAT] + n*l1;
Iv[NEXTV] = Iv[DAMAT] + n*m;
Iv[L1SAV] = l1;
Iv[MSAVE] = m;
/* *** SET UP IN ARRAY ***
* */
in1 = Iv[IN];
for (i = 1; i <= p; i++)
{
for (k = 1; k <= lp1; k++)
{
if (INC(i - 1,k - 1) == 0)
goto L_70;
Iv[in1] = i;
Iv[in1 + 1] = k;
in1 += 2;
L_70:
;
}
}
if (iv1 == 13)
goto L_999;
L_90:
a1 = Iv[AMAT];
da1 = Iv[DAMAT];
in1 = Iv[IN];
l1 = Iv[L1SAV];
m = Iv[MSAVE];
L_100:
srnsgb( &V[a1], alf, b, c, &V[da1], (long(*)[2])&Iv[in1], iv,
l, l1, n, liv, lv, n, m, p, v, y );
switch (IARITHIF(Iv[1] - 2))
{
case -1: goto L_110;
case 0: goto L_120;
case 1: goto L_999;
}
/* *** NEW FUNCTION VALUE (R VALUE) NEEDED ***
* */
L_110:
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)) */
goto L_100;
/* *** COMPUTE DR = GRADIENT OF R COMPONENTS ***
* */
L_120:
;
(*scalcb)( n, p, l, alf, &Iv[NFGCAL], &V[da1] );
if (Iv[NFGCAL] == 0)
Iv[TOOBIG] = 1;
/* CALL SCALCB(N, P, L, ALF, IV(NFGCAL), V(DA1)) */
goto L_100;
L_999:
return;
/* *** LAST CARD OF SNLSGB FOLLOWS *** */
#undef INC
} /* end of function */