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