      SUBROUTINE   SNLSFU(N, P, L, ALF, C, Y, SCALCA, INC, IINC, IV,
     1                  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 SNLSFU Krogh  Changes to get desired C prototypes.
C>> 1994-10-20 SNLSFU Krogh  Changes to use M77CON
C>> 1990-07-02 SNLSFU CLL @ JPL
C>> 1990-06-12 CLL @ JPL
C>> 1990-02-20 CLL @ JPL
*** from netlib, Fri Feb 16 15:19:34 EST 1990 ***
c--S replaces "?": ?NLSFU,?CALCA,?IVSET,?RNSG,?V2AXY,?V7CPY,?V7SCL
C
C  ***  SOLVE SEPARABLE NONLINEAR LEAST SQUARES USING
C  ***  FINITE-DIFFERENCE DERIVATIVES.
C
C  ***  PARAMETER DECLARATIONS  ***
C
      EXTERNAL SCALCA
      INTEGER IINC, L, LIV, LV, N, P
      INTEGER INC(IINC,P), IV(LIV)
      REAL             ALF(P), C(L), V(LV), Y(N)
C
C  ***  PARAMETERS  ***
C
C      N (IN)  NUMBER OF OBSERVATIONS.
C      P (IN)  NUMBER OF NONLINEAR PARAMETERS TO BE ESTIMATED.
C      L (IN)  NUMBER OF LINEAR PARAMETERS TO BE ESTIMATED.
C    ALF (I/O) NONLINEAR PARAMETERS.
C                 INPUT = INITIAL GUESS,
C                 OUTPUT = BEST ESTIMATE FOUND.
C      C (OUT) LINEAR PARAMETERS (ESTIMATED).
C      Y (IN)  RIGHT-HAND SIDE VECTOR.
C  SCALCA (IN)  SUBROUTINE TO COMPUTE A MATRIX.
C    INC (IN)  INCIDENCE MATRIX OF DEPENDENCIES OF COLUMNS OF A ON
C                 COMPONENTS OF ALF -- INC(I,J) = 1 MEANS COLUMN I
C                 OF A DEPENDS ON ALF(J).
C   IINC (IN)  DECLARED LEAD DIMENSION OF INC.  MUST BE AT LEAST L+1.
C     IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR.
C    LIV (IN)  LENGTH OF IV.  MUST BE AT LEAST
C                 122 + 2*M + 4*P + 2*L + MAX(L+1,6*P), WHERE  M  IS
C                 THE NUMBER OF ONES IN INC.
C     LV (IN)  LENGTH OF V.  MUST BE AT LEAST
C                 105 + 2*N*(L+3) + JLEN + L*(L+3)/2 + P*(2*P + 18),
C                 WHERE  JLEN = (L+P)*(N+L+P+1),  UNLESS NEITHER A
C                 COVARIANCE MATRIX NOR REGRESSION DIAGNOSTICS ARE
C                 REQUESTED, IN WHICH CASE  JLEN = N*P.
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 4*N less than just described.
C      V (I/O) FLOATING-POINT PARAMETER AND SCRATCH VECTOR.
C
C
c -------------------------  DECLARATIONS  ----------------------------
C
C
C  ***  EXTERNAL SUBROUTINES  ***
C
      EXTERNAL SIVSET, IDSM,  SRNSG,SV2AXY,SV7CPY, SV7SCL
C
C SIVSET.... PROVIDES DEFAULT IV AND V VALUES.
C IDSM...... DETERMINES EFFICIENT ORDER FOR FINITE DIFFERENCES.
C  SRNSG... CARRIES OUT NL2SOL ALGORITHM.
C SV2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER.
C SV7CPY.... COPIES ONE VECTOR TO ANOTHER.
C SV7SCL... SCALES AND COPIES ONE VECTOR TO ANOTHER.
C
C  ***  LOCAL VARIABLES  ***
C
      LOGICAL PARTJ
      INTEGER A0, A1, AJ, ALP1, BWA1, D0, DA0, DA1, DAJ, GPTR1, GRP1,
     1        GRP2, I, I1, IN0, IN1, IN2, INI, INLEN, IPNTR1, IV1, IWA1,
     2        IWALEN, J1, JN1, JPNTR1, K, L1, LP1, M, M0, NF, NG, NGRP0,
     3        NGRP1, NGRP2, RSAVE0, RSAVE1, RSVLEN, X0I, XSAVE0, XSAVE1
      REAL             DELTA, DI, H, XI
      REAL             NEGONE, ONE, ZERO
C
C  ***  SUBSCRIPTS FOR IV AND V  ***
C
      INTEGER AMAT, COVREQ, D, DAMAT, DLTFDJ, GPTR, GRP, IN, IVNEED, J,
     1        L1SAV, MAXGRP, MODE, MSAVE, NEXTIV, NEXTV, NFCALL, NFGCAL,
     2        PERM, R, RESTOR, TOOBIG, VNEED, XSAVE
C
C  ***  IV SUBSCRIPT VALUES  ***
C
      PARAMETER (AMAT=113, COVREQ=15, D=27, DAMAT=114, DLTFDJ=43,
     1           GPTR=117, GRP=118, IN=112, IVNEED=3, J=70, L1SAV=111,
     2           MAXGRP=116, MODE=35, MSAVE=115, NEXTIV=46, NEXTV=47,
     3           NFCALL=6, NFGCAL=7, PERM=58, R=61, RESTOR=9, TOOBIG=2,
     4           VNEED=4, XSAVE=119)
      DATA NEGONE/-1.E+0/, ONE/1.E+0/, ZERO/0.E+0/
C
C++++++++++++++++++++++++++++++++  BODY  ++++++++++++++++++++++++++++++
C
      LP1 = L + 1
      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 80
      IV1 = IV(1)
      IF (IV1 .EQ. 14) GO TO 120
      IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 120
      IF (IV1 .EQ. 12) IV(1) = 13
      IF (IV(1) .NE. 13) GO TO 50
C
C  ***  FRESH START ***
C
      IF (IV(PERM) .LE. XSAVE) IV(PERM) = XSAVE + 1
C
C  ***  CHECK INC, COUNT ITS NONZEROS
C
      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 80
            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 80
 40      CONTINUE
C
C     *** NOW L1 = 1 MEANS A HAS COLUMN L+1 ***
C
C     *** COMPUTE STORAGE REQUIREMENTS ***
C
      IWALEN = max(LP1, 6*P)
      INLEN = 2 * M
      IV(IVNEED) = IV(IVNEED) + INLEN + 3*P + L + IWALEN + 3
      RSVLEN = 2 * L1 * N
      L1 = L + L1
      IV(VNEED) = IV(VNEED) + 2*N*L1 + RSVLEN + P
C
 50   CALL  SRNSG(V, ALF, C, V, IV, IV, L, 1, N, LIV, LV, N, M, P, V, Y)
      IF (IV(1) .NE. 14) GO TO 999
C
C  ***  STORAGE ALLOCATION  ***
C
      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
C
C  ***  DETERMINE HOW MANY GROUPS FOR FINITE DIFFERENCES
C  ***  (SET UP TO CALL IDSM)
C
      IN1 = IV(IN)
      JN1 = IN1 + M
      DO 70 K = 1, P
         DO 60 I = 1, LP1
            IF (INC(I,K) .EQ. 0) GO TO 60
               IV(IN1) = I
               IN1 = IN1 + 1
               IV(JN1) = K
               JN1 = JN1 + 1
 60         CONTINUE
 70      CONTINUE
      IN1 = IV(IN)
      JN1 = IN1 + M
      IWA1 = IN1 + INLEN
      NGRP1 = IWA1 + IWALEN
      BWA1 = NGRP1 + P
      IPNTR1 = BWA1 + P
      JPNTR1 = IPNTR1 + L + 2
      CALL IDSM(LP1, P, M, IV(IN1), IV(JN1), IV(NGRP1), NG, K, I,
     1         IV(IPNTR1), IV(JPNTR1), IV(IWA1), IWALEN, IV(BWA1))
      IF (I .EQ. 1) GO TO 90
         IV(1) = 69
         GO TO 50
 80   IV(1) = 66
      GO TO 50
C
C  ***  SET UP GRP AND GPTR ARRAYS FOR COMPUTING FINITE DIFFERENCES
C
C  ***  THERE ARE NG GROUPS.  GROUP I CONTAINS ALF(GRP(J)) FOR
C  ***  GPTR(I) .LE. J .LE. GPTR(I+1)-1.
C
 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
      DO 110 I = 1, NG
         IV(GPTR1) = GRP1
         GPTR1 = GPTR1 + 1
         DO 100 I1 = NGRP1, NGRP2
            IF (IV(I1) .NE. I) GO TO 100
            IV(GRP1) = I1 - NGRP0
            GRP1 = GRP1 + 1
 100        CONTINUE
 110     CONTINUE
      IV(GPTR1) = GRP1
      IF (IV1 .EQ. 13) GO TO 999
C
C  ***  INITIALIZE POINTERS  ***
C
 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) = -abs(IV(COVREQ))
C
 130  CALL  SRNSG(V(A1), ALF, C, V(DA1), IV(IN1), IV, L, L1, N, LIV, LV,
     1            N, L1, P, V, Y)
      IF (IV(1)-2) 140, 150, 999
C
C  ***  NEW FUNCTION VALUE (R VALUE) NEEDED  ***
C
 140  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
      IF (L1 .LE. L) GO TO 130
      IF (IV(RESTOR) .EQ. 2) CALL SV7CPY(N, V(RSAVE0), V(RSAVE1))
      CALL SV7CPY(N, V(RSAVE1), V(ALP1))
      GO TO 130
C
C  ***  COMPUTE DR = GRADIENT OF R COMPONENTS  ***
C
 150  IF (L1 .GT. L .AND. IV(NFGCAL) .EQ. IV(NFCALL))
     1      CALL SV7CPY(N, V(RSAVE0), V(RSAVE1))
      GPTR1 = IV(GPTR)
      DO 230 K = 1, NG
         CALL SV7CPY(P, V(XSAVE1), ALF)
         GRP1 = IV(GPTR1)
         GRP2 = IV(GPTR1+1) - 1
         GPTR1 = GPTR1 + 1
         DO 160 I1 = GRP1, GRP2
            I = IV(I1)
            XI = ALF(I)
            J1 = D0 + I
            DI = V(J1)
            IF (DI .LE. ZERO) DI = ONE
            H = DELTA * max(abs(XI), ONE/DI)
            IF (XI .LT. ZERO) H = -H
            X0I = XSAVE0 + I
            V(X0I) = XI + H
 160        CONTINUE
C%%        (*scalca)( n, p, l, &V[xsave1], &Iv[NFGCAL], &V[da1] );
         CALL SCALCA(N, P, L, V(XSAVE1), IV(NFGCAL), V(DA1))
         IF (IV(NFGCAL) .GT. 0) GO TO 170
            IV(TOOBIG) = 1
            GO TO 130
 170     JN1 = IN1
         DO 180 I = IN1, IN2
 180        IV(I) = 0
         PARTJ = IV(MODE) .LE. P
         DO 220 I1 = GRP1, GRP2
            I = IV(I1)
            DO 210 J1 = 1, L1
               IF (INC(J1,I) .EQ. 0) GO TO 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) GO TO 190
C                 *** FULL FINITE DIFFERENCE FOR COV. AND REG. DIAG. ***
                  AJ = A0 + J1*N
                  CALL SV2AXY(N, V(DAJ), NEGONE, V(AJ), V(DAJ))
                  GO TO 200
 190           IF (J1 .GT. L)
     1            CALL SV2AXY(N, V(DAJ), NEGONE, V(RSAVE0), V(DAJ))
 200           CALL SV7SCL(N, V(DAJ), H, V(DAJ))
 210           CONTINUE
 220        CONTINUE
         IF (K .GE. NG) GO TO 240
         IV(1) = -2
         CALL  SRNSG(V(A1), ALF, C, V(DA1), IV(IN1), IV, L, L1, N, LIV,
     1               LV, N, L1, P, V, Y)
         IF (-2 .NE. IV(1)) GO TO 999
 230     CONTINUE
 240  IV(1) = 2
      GO TO 130
C
 999  RETURN
C
C  ***  LAST CARD OF   SNLSFU FOLLOWS  ***
      END
