C
C           SAMPLE DRIVER PROGRAM FOR NONLINEAR LEAST SQUARES ROUTINE
C
C
C            SOLUTION TO SAMPLE PROBLEM (TO 3 SIGNIFICANT DIGITS):
C               LINEAR PARAMETERS  .375, 1.94, -1.46
C               NONLINEAR PARAMETERS  .0129, .0221
C
      DOUBLE PRECISION Y(50), T(50, 1), ALF(12), BETA(8), A(50, 13),
     X W(50), DSQRT, ETA
      INTEGER P, OUTPUT
      EXTERNAL ADA
      DATA INPUT /5/, OUTPUT /6/
C
C              READ AND LIST THE DATA
C
      NMAX = 50
      IPRINT = 1
      READ (INPUT, 100) N, L, NL, P, IV
      WRITE (OUTPUT, 101) N, L, NL, P, IV
      DO 12 I = 1, N
         W(I) = 1.0
   12    READ (INPUT, 105) (T(I,J), J = 1, IV), Y(I)
      DO 15 I = 1, N
   15    WRITE (OUTPUT, 213) I, (T(I,J), J = 1, IV), Y(I)
C
      READ (INPUT, 102) (ALF(K), K = 1, NL)
      WRITE (OUTPUT, 103) (ALF(K), K = 1, NL)
      WRITE (OUTPUT, 104)
C
      CALL VARPRO (L, NL, N, NMAX, L+P+2, IV, T, Y, W, ADA, A, IPRINT,
     X ALF, BETA, IERR)
C
C              PRINT LOWER TRIANGLE OF COVARIANCE MATRIX
C
      IF (IERR .LE. -4) GO TO 80
      LPNL = L + NL
      WRITE (OUTPUT, 106)
      DO 65 I = 1, LPNL
         WRITE (OUTPUT, 107) (A(I, J), J = 1, I)
   65    A(I, I) = DSQRT(A(I, I))
C
C              PRINT STANDARD DEVIATIONS OF PARAMETER ESTIMATES
C
      WRITE (OUTPUT, 108)
      WRITE (OUTPUT, 107) (A(J, J), J = 1, LPNL)
C
C              CALCULATE AND PRINT LOWER TRIANG. OF CORRELATION MATRIX
C
      IF (LPNL .EQ. 1) GO TO 72
      DO 68 I = 2, LPNL
         IM1 = I - 1
         DO 68 J = 1, IM1
   68       A(I,J) = A(I,J)/ (A(I,I)*A(J,J))
C
   72 WRITE (OUTPUT, 109)
      DO 74 I = 1, LPNL
         A(I,I) = 1.
   74    WRITE (OUTPUT, 107) (A(I, J), J = 1, I)
C
C              PRINT RESIDUALS
C
      WRITE (OUTPUT, 212)
      DO 75 I = 1, N
         ETA = Y(I) - A(I, LPNL+1)/DSQRT(W(I))
   75    WRITE (OUTPUT, 213) I, W(I), (T(I,J), J = 1, IV), Y(I),
     X   ETA, A(I, LPNL+1)
C
   80 STOP
C
  100 FORMAT (5I5)
  101 FORMAT (36H1   NON-LINEAR LEAST SQUARES PROBLEM //
     1 25H NUMBER OF OBSERVATIONS =, I5, 32H   NUMBER OF LINEAR PARAMETE
     2RS = , I4  // 33H NUMBER OF NONLINEAR PARAMETERS =, I4,
     3 45H NUMBER OF NONVANISHING PARTIAL DERIVATIVES =, I4 //
     4 34H NUMBER OF INDEPENDENT VARIABLES =, I4 //
     5 33H    I        T(I)            Y(I) //)
  102 FORMAT (4E20.7)
  103 FORMAT (30H0 INITIAL NONLINEAR PARAMETERS //(4E20.7))
  104 FORMAT (1H0, 50(1H*))
  105 FORMAT (5E15.7)
  106 FORMAT (/ 24H       COVARIANCE MATRIX  /)
  107 FORMAT (8E15.7)
  108 FORMAT (47H0    STANDARD DEVIATIONS OF PARAMETER ESTIMATES /)
  109 FORMAT (/ 25H       CORRELATION MATRIX  /)
  212 FORMAT (/89H   I         W(I)            T(I)            Y(I)
     X  PREDICTED Y     WEIGHTED RESIDUAL //)
  213 FORMAT (I5, 7E16.7)
      END
C
C
      SUBROUTINE ADA (LP1, NL, N, NMAX, LPP2, IV, A, INC, T, ALF, ISEL)
C
C        EVALUATE THE FUNCTIONS PHI AND THEIR PARTIAL DERIVATIVES
C        D PHI(J)/D ALF(K), AT THE SAMPLE POINTS T(I).
C        ISEL = 1  MEANS COMPUTE BOTH FUNCTIONS AND DERIVATIVES AND
C                  INITIALIZE INC AND CONSTANT PHI'S
C             = 2  MEANS COMPUTE ONLY THE NONCONSTANT FUNCTIONS PHI
C             = 3  MEANS COMPUTE ONLY THE DERIVATIVES.
C
C        THIS PARTICULAR ROUTINE IS FOR FITTING TWO EXPONENTIALS AND
C        ONE CONSTANT TERM:
C
C        ETA(C, ALF; T) = C  + C  * EXP(-ALF *T) + C * EXP(-ALF * T)
C                          1    2           1       3          2
C
C        WHERE C IS THE VECTOR OF LINEAR PARAMETERS TO BE DETERMINED,
C        AND ALF IS THE VECTOR OF NONLINEAR PARAMETERS TO BE DETERMINED.
C        HERE N = 33, L = 3, NL = P = 2, NCFUN = 1, PHI(1) = 1, PHI(2) =
C        EXP(-ALF(1)*T), AND PHI(3) = EXP(-ALF(2)*T).  THIS EXAMPLE AND
C        TEST DATA ARE TAKEN FROM REFERENCE 3 (SEE VARPRO).
C
C     ..................................................................
C
      DOUBLE PRECISION A(NMAX, LPP2), ALF(NL), T(NMAX), DEXP
      INTEGER INC(12, LP1)
C
      GO TO (10, 20, 30), ISEL
C
C              ISEL = 1, IN THIS CASE THE INCIDENCE MATRIX INC IS:
C
C               0  1  0  0     <- ALF(1)
C               0  0  1  0     <- ALF(2)
C
   10 INC(1, 2) = 1
      INC(2, 3) = 1
C
C              ISEL = 1, SET CONSTANT FUNCTION PHI(1) (REMOVE IF NO
C              CONSTANT FUNCTIONS)
C
      DO 12 I = 1, N
   12    A(I,1) = 1.0
C
C              ISEL = 1 OR 2, COMPUTE NONCONSTANT FUNCTIONS PHI
C
   20 DO 22 I = 1, N
         A(I,2) = DEXP(-ALF(1)*T(I))
   22    A(I,3) = DEXP(-ALF(2)*T(I))
C                           COLUMN L+1 = 4 IS LEFT FOR WORKSPACE.
      IF (ISEL .EQ. 2) GO TO 35
C
C              ISEL = 1 OR 3, COMPUTE DERIVATIVES DPHI(I) / D ALF(J)
C
   30 DO 32 I = 1, N
         A(I,5) = -T(I)*DEXP(-ALF(1)*T(I))
   32    A(I,6) = -T(I)*DEXP(-ALF(2)*T(I))
   35 RETURN
      END