SUBROUTINE CALCPR(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT,
     1                  PROB, IUSER, RUSER, MNPCDF)
C
C *** THIS SUBROUTINE CALCULATES A PROBABILITY FOR THE MODEL AND   ***
C *** DATA GIVEN.  FOR MULTINOMIAL PROBIT SOME ADDITIONAL STORAGE  ***
C *** CUSTOMIZATION IS REQUIRED.  THIS APPROACH CAN BE             ***
C *** USED FOR OTHER CHOICE MODELS, WITH APPROPRIATE MODIFICATIONS ***
C *** TO THE ARRAYS USED BELOW.                                    ***
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
      INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*)
      DOUBLE PRECISION X(NPAR), RCDAT(*), PROB, RUSER(*)
      EXTERNAL MNPCDF
C
C *** CALCPR PARAMETER USAGE ***
C
C IALT.... NUMBER OF CHOICES AVAILABLE IN THE CHOICE SET.
C ICDAT... VECTOR OF INTEGER DATA VALUES.
C ICH..... INTEGER INDICATING THE CHOICE.  1 <= ICH <= IALT.
C IERR.... INTEGER FOR PASSING ERROR INFORMATION.
C             IN THIS ROUTINE, IF IERR = 1 ON RETURN THEN THERE WERE
C               NO PROBLEMS.
C             IF IERR = 0 ON RETURN, THEN THE PROBABILITY COULD NOT
C               BE COMPUTED USING THE CURRENT PARAMETERS IN X.
C II...... NUMBER OF INTEGER VALUES STORED IN VECTIR ICDAT.
C IUSER... MODEL-RELATED INTEGER VALUES USED BY CALCPR.  CONTAINS
C             ARRAY POINTERS TO MANAGE DATA STORAGE, AND OTHER
C             PARAMETERS.
C MNPCDF.. SUBROUTINE WHICH CALCULATES THE CDF OF A MULTIVARIATE
C          NORMAL DISTRIBUTION.
C NPAR.... NUMBER OF PARAMETERS IN VECTOR X.
C PROB.... ON RETURN, CHOICE PROBABILITY COMPUTED USING PARAMETERS IN
C             X AND DATA IN ICDAT AND RCDAT.
C RCDAT... VECTOR OF REAL DATA VALUES.
C RUSER... MODEL-RELATED REAL VALUES USED BY CALCPR.  CAN CONTAIN
C             USEFUL PARAMETERS, AND ALSO EXTRA WORK STORAGE.
C
      EXTERNAL CALCP1
      INTEGER ISIGU, IW, NALT, NW
C
      ISIGU = IUSER(16)
      IW = IUSER(15)
      NALT = MAX(1,IUSER(6))
      NW = MAX(1, IUSER(7))
      CALL CALCP1(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT,
     1            PROB, IUSER, RUSER, NALT, RUSER(ISIGU),
     2            NW, RUSER(IW), MNPCDF)
C *** LAST LINE OF CALCPR FOLLOWS ***
      END
      SUBROUTINE CALCP1(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT,
     1                PROB, IUSER, RUSER, NALT, SIGU, NW, W, MNPCDF)
C
C *** THIS SUBROUTINE CALCULATES A PROBABILITY FOR THE MNP MODEL   ***
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
      INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*), NALT,
     1        NATTR, NW
      DOUBLE PRECISION X(NPAR), RCDAT(*), PROB, RUSER(*),
     1                 SIGU(NALT,NALT), W(NW,NALT)
      EXTERNAL MNPCDF
C
C  ***  CALCP1 PARAMETER USAGE ***
C
C IALT.... NUMBER OF CHOICES AVAILABLE IN THE CHOICE SET.
C ICDAT... VECTOR OF INTEGER DATA VALUES.
C             IN THIS SUBROUTINE, ICDAT STORES INTEGER INDEXES WHICH
C             DEFINE WHICH OF THE NOMINAL ALTERNATIVES ARE AVAILABLE
C             IN  THE CHOICE SET.  (THIS IS FOR THE CASE WHEN THERE
C             ARE NALT NOMINAL CHOICE ALTERNATIVES, BUT NOT ALL OF
C             THEM NECESSARILY APPEAR IN EVERY SUBSET.  IF ALL NALT
C             ALTERNATIVES APPEAR IN ALL SUBSETS, THEN ICSET = NALT >0
C             SHOULD BE USED WITH IDUM = 1.
C ICH..... INTEGER INDICATING THE CHOICE.  1 <= ICH <= IALT.
C ICOV.... INDICATOR FOR TYPE OF ALTERNATIVE-SPECIFIC ERRORS,
C             = 0 FOR IID ERRORS, = 1 FOR CORRELATED ERRORS.
C             IF ICSET .NE. 0, THEN THE SAME CORRELATION MATRIX IS
C             USED FOR EVERY SUBSET.  OTHERWISE, INTEGER DATA SHOULD
C             BE USED TO IDENTIFY THE ALTERNATIVES IN EACH CHOICE SET.
C             (STORED IN IUSER.)
C IDUM... INDICATOR FOR ALTERNATIVE-SPECIFIC DUMMIES,
C             = 0 FOR NO, = 1 FOR YES.  IF ICSET .NE. 0, THEN
C             THE SAME SET OF DUMMIES IS USED FOR EACH CHOICE SET.
C             OTHERWISE, INTEGER DATA SHOULD BE USED TO IDENTIFY THE
C             ALTERNATIVES IN EACH CHOICE SET (SEE NALT BELOW).
C             (STORED IN IUSER).
C IERR.... INTEGER FOR PASSING ERROR INFORMATION.
C             IN THIS ROUTINE, IF IERR = 1 ON RETURN THEN THERE WERE
C               NO PROBLEMS.
C             IF IERR = 0 ON RETURN, THEN THE PROBABILITY COULD NOT
C               BE COMPUTED USING THE CURRENT PARAMETERS IN X.
C II...... NUMBER OF INTEGER VALUES STORED IN VECTIR ICDAT.
C IUSER... MODEL-RELATED INTEGER VALUES USED BY CALCPR.  THE FIRST
C             PORTION OF IUSER CONTAINS SUCH THINGS AS ARRAY POINTERS.
C             IUSER ALSO CONTAINS STORED VALUES OF NATTR, IDUM, ETC.
C IR...... NUMBER OF REAL VALUES STORED IN VECTOR IRDAT.
C ITASTE.. INDICATOR FOR TASTE VARIATION,
C             = 0 FOR NO TASTE VARIATION, = 1 FOR UNCORRELATED TASTE
C             VARIATION, = 2 FOR CORRELATED TASTE VARIATION.
C             (STORED IN IUSER.)
C NPAR.... NUMBER OF PARAMETERS IN VECTOR X.
C PROB.... ON RETURN, CHOICE PROBABILITY COMPUTED USING PARAMETERS IN
C             X AND DATA IN ICDAT AND RCDAT.
C RCDAT... VECTOR OF REAL DATA VALUES.
C             IN THIS SUBROUTINE, THE NUMBER OF DATA VALUES SHOULD
C             BE = IALT * NATTR SO THAT THE "GENERIC" PART OF THE
C             SCALE VALUE V MAY BE COMPUTED.
C NALT.... TOTAL NUMBER OF NOMINAL CHOICE ALTERNATIVES (IF APPLICABLE).
C             IF ICSET .NE. 0 AND IDUM = 1 OR ICOV = 1 (OR BOTH), THEN
C             NALT SHOULD BE EQUAL TO ICSET.
C             OTHERWISE, NALT SHOULD BE > 0 IF EITHER IDUM OR ICOV
C             (OR BOTH) ARE > 0, AND ICDAT SHOULD BE USED TO PASS
C             INDEX INFORMATION (SEE ICDAT ABOVE).
C NATTR... NUMBER OF ATTRIBUTES (I.E., REAL DATA VARS.) PER ALTERNATIVE.
C NW...... NUMBER OF ROWS IN THE WORK-ARRAY W.
C RUSER... MODEL-RELATED REAL VALUES USED BY CALCPR.  FOR THIS MODEL,
C             IT CONTAINS A CONSTANT FOR THE COVARIANCE MATRIX SCALE,
C             AND INFORMATION USED FOR COMPUTING STEP SIZES IN FINITE-
C             DIFFERENCE CALCULATIONS.
C SIGU.... MATRIX CONTAINING THE "UNPACKED" THE FULL COVARIANCE MATRIX
C             FOR ALL NALT ALTERNATIVE-SPECIFIC ERROR TERMS.  THE
C             MATRIX IS OF DIMENSION 2 TO FACILITATE CODING.  THE
C             NORMALIZATION USED LEAVES A ROW OF ZEROS IN THE LAST
C             (NALT) ROW.  IT IS COMPUTED BEFORE THE CALL TO MINIMIZE
C             WORK WHEN CALLS ARE TO BE REPEATED.
C W....... ARRAY CONTAINING WORKSPACE FOR COVARIANCE COMPUTATIONS.
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
C
      EXTERNAL DL7VML, DV7SCP
C
      INTEGER I, IALTM1, ICOL, ICOV, ICSET, ID, IDUM, IFAULT, IIR,
     1        IOUNIT, IPCOEF, IPDUM, IPT, IPTAST, IROW, ISCALE, ISZ,
     2        ITASTE, IX, J, JP, K, KP
C
      INTEGER MAXALT, MAXAM1, LR
      PARAMETER (MAXALT=20, MAXAM1=MAXALT-1, LR=MAXAM1*(MAXAM1-1)/2)
C
      DOUBLE PRECISION SCALE, SII
      DOUBLE PRECISION V(MAXALT), SIGMA(MAXALT,MAXALT)
      DOUBLE PRECISION Z(MAXAM1), SIGZ(MAXAM1,MAXAM1), R(LR)
C
      DOUBLE PRECISION ZERO
      PARAMETER (ZERO=0.D0)
C
C  SET UP V AND SIGMA MATRIX FOR MNP SPECIFICATION.
C
C  ALTERNATIVE-SPECIFIC DUMMIES:
C
      IALTM1 = IALT - 1
      IDUM = IUSER(8)
      IF (IDUM.NE.0) THEN
         IPDUM = IUSER(12)
C        CASE 1:  ICSET = 0.
         ICSET = IUSER(5)
         IF (ICSET.EQ.0) THEN
            DO 10 I = 1, IALT
               IX = ICDAT(I)
               IF (IX.NE.NALT) THEN
                  V(I) = X(IX+IPDUM)
               ELSE
                  V(I) = ZERO
               ENDIF
 10         CONTINUE
          ELSE
C         CASE 2:  ICSET.NE.0
            V(IALT) = ZERO
            DO 20 I = 1, IALTM1
                V(I) = X(I+IPDUM)
 20         CONTINUE
          ENDIF
      ELSE
          CALL DV7SCP(IALT, V, ZERO)
      ENDIF
C
C  BETA COEFFICIENTS:
C
      NATTR = IUSER(7)
      IF (NATTR.NE.0) THEN
         IPCOEF = IUSER(11)
         ID = 0
         DO 30 I = 1, IALT
            DO 30 K = 1, NATTR
               ID = ID + 1
               V(I) = V(I) + X(IPCOEF+K)*RCDAT(ID)
 30      CONTINUE
      ENDIF

C
C  ALTERNATIVE-SPECIFIC ERRORS:
C
      ICOV = IUSER(9)
      IF (ICOV.NE.0) THEN
         ICSET = IUSER(5)
         IF (ICSET.EQ.0) THEN
            DO 40 I = 1, IALT
               IROW = ICDAT(I)
               DO 40 J = 1, I
                  ICOL = ICDAT(J)
                  IF (ICOL.LE.IROW) THEN
                     SIGMA(I,J) = SIGU(IROW,ICOL)
                  ELSE
                     SIGMA(I,J) = SIGU(ICOL,IROW)
                  ENDIF
 40         CONTINUE
         ELSE
            DO 50 I = 1, IALT
               DO 50 J = 1, I
                  SIGMA(I,J) = SIGU(I,J)
 50         CONTINUE
         ENDIF
      ELSE
         ISCALE = IUSER(18)
         SCALE = RUSER(ISCALE)
         DO 60 I = 1, IALT
            DO 60 J = 1, I
               IF (I.EQ.J) THEN
                  SIGMA(I,J) = SCALE
               ELSE
                  SIGMA(I,J) = ZERO
               ENDIF
 60      CONTINUE
       ENDIF
C
C  TASTE VARIATION:
C
      ITASTE = IUSER(10)
      IF (ITASTE.EQ.1) THEN
C        UNCORRELATED TASTE VARIATION
C        SET UP W MATRIX:
         ID = 0
         IPTAST = IUSER(14)
         DO 70 J = 1, IALT
            IPT = IPTAST
            DO 70 K = 1, NATTR
               IPT = IPT + 1
               ID = ID + 1
               W(K,J) = X(IPT) * RCDAT(ID)
 70      CONTINUE
      ENDIF
C
      IF (ITASTE.EQ.2) THEN
C        CORRELATED TASTE VARIATION
C        SET UP W MATRIX:
         ID = 1
         IPTAST = IUSER(14) + 1
         DO 80 J = 1, IALT
            CALL DL7VML(NATTR, W(1,J), X(IPTAST), RCDAT(ID))
            ID = ID + NATTR
 80      CONTINUE
      ENDIF

      IF (ITASTE.NE.0) THEN
C        TASTE VARIATION
C        ADD W(**T)W TO SIGMA:
         DO 100 I = 1, IALT
            DO 100 J = 1, I
               DO 90 K = 1, NATTR
                  SIGMA(I,J) = SIGMA(I,J) + W(K,I)*W(K,J)
 90            CONTINUE
 100     CONTINUE
      ENDIF
C
C  SYMMETRIZE SIGMA (MAY NOT BE NECESSARY???)
C
C      IF ((ICOV.NE.0).OR.(ITASTE.NE.0)) THEN
         DO 110 I = 1, IALT
            DO 110 J = 1, I
               SIGMA(J,I) = SIGMA(I,J)
 110     CONTINUE
C      ENDIF
C
C  LOWER DIMENSION VIA STANDARD TRANSFORMATION
C  (REF. PAGE 43 OF DAGANZO OR BUNCH(1991))
      ISZ = 0
      SII = SIGMA(ICH,ICH)
      DO 130 JP = 1, IALT
         IF (JP.LT.ICH) THEN
            J = JP
         ELSE
            J = JP - 1
         ENDIF
         IF (JP.NE.ICH) THEN
            Z(J) = V(JP)-V(ICH)
            DO 120 KP = 1, JP
               IF (KP.LT.ICH) THEN
                  K = KP
               ELSE
                  K = KP - 1
               ENDIF
               IF(KP.NE.ICH) THEN
                  ISZ = ISZ + 1
                  SIGZ(J,K)=SIGMA(JP,KP)-SIGMA(ICH,KP)-SIGMA(ICH,JP)+SII
               ENDIF
 120        CONTINUE
         ENDIF
 130   CONTINUE
C
      IIR = 0
      DO 150 J = 1, IALTM1
         IF (SIGZ(J,J).LE.ZERO) THEN
            IERR = 0
            RETURN
         ENDIF
         SIGZ(J,J) = SQRT(SIGZ(J,J))
         Z(J) = Z(J)/SIGZ(J,J)
         DO 140 K = 1, J-1
               IIR = IIR + 1
               R(IIR) = SIGZ(J,K)/SIGZ(J,J)/SIGZ(K,K)
 140     CONTINUE
 150  CONTINUE
C
      IERR = 1
      CALL MNPCDF(IALTM1, Z, R, PROB, IFAULT)
      IF (IFAULT.NE.0) then
         IERR = 0
         IOUNIT = IUSER(3)
         WRITE(IOUNIT,*) ' Problem evaluating mnpcdf'
      ENDIF
C *** LAST LINE OF CALCP1 FOLLOWS ***
      END
      SUBROUTINE CALCDP(NPAR, X, IERR, ICH, IALT, II, ICDAT, IR, RCDAT,
     1                  PROB0, DP, IUSER, RUSER, MNPCDF)
C
C *** THIS SUBROUTINE CALCULATES FINITE-DIFFERENCE DERIVATIVES FOR ***
C *** CHOICE PROBABILITIES.  THIS VERSION ASSUMES THAT THE CALCPR  ***
C *** BEING CALLED IS THE ONE FOR MULTINOMIAL PROBIT.  HOWEVER,    ***
C *** THE CHANGES REQUIRED FOR OTHER MODELS SHOULD BE MINOR.       ***
C *** NOTE:  THIS SUBROUTINE REQUIRES DS7GRD, AND THE ARRAYS ALPHA ***
C *** AND D SHOULD HAVE THE SAME DIMENSION AS X.                   ***
C+++++++++++++++++++++++++++  DECLARATIONS  +++++++++++++++++++++++++++
C
      INTEGER NPAR, IERR, ICH, IALT, II, ICDAT(*), IR, IUSER(*)
      DOUBLE PRECISION X(NPAR), RCDAT(*), PROB0, DP(NPAR), RUSER(*)
      EXTERNAL MNPCDF
C
      EXTERNAL CALCPR, DS7GRD, DV7SCP
      INTEGER I, ICOV, IETA0, IPCOV, IPP, IPU, IPUP, IRC, ISCALE, ISIGP,
     1        ISIGU, J, NALT, NALTM1, NFC
      DOUBLE PRECISION ETA, ETA0, PROB, SCALE, XTEMP
C
      INTEGER LX
      DOUBLE PRECISION ONE, ZERO
      PARAMETER (ZERO=0.D0, ONE=1.D0, LX=60)
C
      DOUBLE PRECISION ALPHA(LX), D(LX), WRK(6)
C
C ***  PARAMETER USAGE ***
C
C SEE CALCPR AND CALCP1
C
C *** BODY ***
C
      IERR = 1
      ICOV = IUSER(9)
      NALT = IUSER(6)
      NALTM1 = NALT - 1
      ISCALE = IUSER(18)
      SCALE = RUSER(ISCALE)
      IETA0 = IUSER(17)
      ETA0 = RUSER(IETA0)
C
      DO 10 I = 1, NPAR
          ALPHA(I) = ONE
          D(I) = ONE
 10   CONTINUE
C
      ETA = ETA0
C     ETA = ETA0/PROB
      IRC = 0
C
      PROB = PROB0
 20   CONTINUE
      CALL DS7GRD(ALPHA, D, ETA, PROB, DP, IRC,
     1             NPAR, WRK, X)
      IF (IRC.EQ.0) GO TO 40
C        IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX
         IF (ICOV.NE.0) THEN
C          SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA:
            IPCOV = IUSER(13)
            XTEMP = X(IPCOV)
            X(IPCOV) = SCALE
            ISIGP = IUSER(15)
            CALL DL7SQR(NALTM1, RUSER(ISIGP), X(IPCOV))
            X(IPCOV) = XTEMP
C          "UNPACK" FOR EASIER ACCESS IN CALCPR:
            IPP = ISIGP - 1
            ISIGU = IUSER(16)
            CALL DV7SCP(NALT*NALT, RUSER(ISIGU), ZERO)
            IPUP = ISIGU - 1
            DO 30 I = 1, NALTM1
               IPU = I + IPUP
               DO 30 J = 1, I
                  IPP = IPP + 1
                  RUSER(IPU) = RUSER(IPP)
                  IPU = IPU + NALT
 30         CONTINUE
         ENDIF
      CALL CALCPR(NPAR, X, NFC, ICH, IALT, II, ICDAT,
     1      IR, RCDAT, PROB, IUSER, RUSER, MNPCDF)
      IF (NFC.EQ.0) THEN
         IERR = 0
         RETURN
      ENDIF
      GO TO 20
 40   CONTINUE
C
C *** LAST LINE OF CALCDP FOLLOWS ***
      END
      SUBROUTINE DS7GRD (ALPHA, D, ETA0, FX, G, IRC, N, W, X)
C
C  ***  COMPUTE FINITE DIFFERENCE GRADIENT BY STWEART*S SCHEME  ***
C  ***  THIS IS SGRAD2 FROM TOMS ALGORITHM 611.
C
C     ***  PARAMETERS  ***
C
      INTEGER IRC, N
      DOUBLE PRECISION ALPHA(N), D(N), ETA0, FX, G(N), W(6), X(N)
C
C.......................................................................
C
C     ***  PURPOSE  ***
C
C        THIS SUBROUTINE USES AN EMBELLISHED FORM OF THE FINITE-DIFFER-
C     ENCE SCHEME PROPOSED BY STEWART (REF. 1) TO APPROXIMATE THE
C     GRADIENT OF THE FUNCTION F(X), WHOSE VALUES ARE SUPPLIED BY
C     REVERSE COMMUNICATION.
C
C     ***  PARAMETER DESCRIPTION  ***
C
C  ALPHA IN  (APPROXIMATE) DIAGONAL ELEMENTS OF THE HESSIAN OF F(X).
C      D IN  SCALE VECTOR SUCH THAT D(I)*X(I), I = 1,...,N, ARE IN
C             COMPARABLE UNITS.
C   ETA0 IN  ESTIMATED BOUND ON RELATIVE ERROR IN THE FUNCTION VALUE...
C             (TRUE VALUE) = (COMPUTED VALUE)*(1+E),   WHERE
C             ABS(E) .LE. ETA0.
C     FX I/O ON INPUT,  FX  MUST BE THE COMPUTED VALUE OF F(X).  ON
C             OUTPUT WITH IRC = 0, FX HAS BEEN RESTORED TO ITS ORIGINAL
C             VALUE, THE ONE IT HAD WHEN DS7GRD WAS LAST CALLED WITH
C             IRC = 0.
C      G I/O ON INPUT WITH IRC = 0, G SHOULD CONTAIN AN APPROXIMATION
C             TO THE GRADIENT OF F NEAR X, E.G., THE GRADIENT AT THE
C             PREVIOUS ITERATE.  WHEN DS7GRD RETURNS WITH IRC = 0, G IS
C             THE DESIRED FINITE-DIFFERENCE APPROXIMATION TO THE
C             GRADIENT AT X.
C    IRC I/O INPUT/RETURN CODE... BEFORE THE VERY FIRST CALL ON DS7GRD,
C             THE CALLER MUST SET IRC TO 0.  WHENEVER DS7GRD RETURNS A
C             NONZERO VALUE FOR IRC, IT HAS PERTURBED SOME COMPONENT OF
C             X... THE CALLER SHOULD EVALUATE F(X) AND CALL DS7GRD
C             AGAIN WITH FX = F(X).
C      N IN  THE NUMBER OF VARIABLES (COMPONENTS OF X) ON WHICH F
C             DEPENDS.
C      X I/O ON INPUT WITH IRC = 0, X IS THE POINT AT WHICH THE
C             GRADIENT OF F IS DESIRED.  ON OUTPUT WITH IRC NONZERO, X
C             IS THE POINT AT WHICH F SHOULD BE EVALUATED.  ON OUTPUT
C             WITH IRC = 0, X HAS BEEN RESTORED TO ITS ORIGINAL VALUE
C             (THE ONE IT HAD WHEN DS7GRD WAS LAST CALLED WITH IRC = 0)
C             AND G CONTAINS THE DESIRED GRADIENT APPROXIMATION.
C      W I/O WORK VECTOR OF LENGTH 6 IN WHICH DS7GRD SAVES CERTAIN
C             QUANTITIES WHILE THE CALLER IS EVALUATING F(X) AT A
C             PERTURBED X.
C
C     ***  APPLICATION AND USAGE RESTRICTIONS  ***
C
C        THIS ROUTINE IS INTENDED FOR USE WITH QUASI-NEWTON ROUTINES
C     FOR UNCONSTRAINED MINIMIZATION (IN WHICH CASE  ALPHA  COMES FROM
C     THE DIAGONAL OF THE QUASI-NEWTON HESSIAN APPROXIMATION).
C
C     ***  ALGORITHM NOTES  ***
C
C        THIS CODE DEPARTS FROM THE SCHEME PROPOSED BY STEWART (REF. 1)
C     IN ITS GUARDING AGAINST OVERLY LARGE OR SMALL STEP SIZES AND ITS
C     HANDLING OF SPECIAL CASES (SUCH AS ZERO COMPONENTS OF ALPHA OR G).
C
C     ***  REFERENCES  ***
C
C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION
C        METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES,
C        J. ASSOC. COMPUT. MACH. 14, PP. 72-83.
C
C     ***  HISTORY  ***
C
C     DESIGNED AND CODED BY DAVID M. GAY (SUMMER 1977/SUMMER 1980).
C
C     ***  GENERAL  ***
C
C        THIS ROUTINE WAS PREPARED IN CONNECTION WITH WORK SUPPORTED BY
C     THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS76-00324 AND
C     MCS-7906671.
C
C.......................................................................
C
C     *****  EXTERNAL FUNCTION  *****
C
      EXTERNAL DR7MDC
      DOUBLE PRECISION DR7MDC
C DR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS.
C
C     ***** LOCAL VARIABLES *****
C
      INTEGER FH, FX0, HSAVE, I, XISAVE
      DOUBLE PRECISION AAI, AFX, AFXETA, AGI, ALPHAI, AXI, AXIBAR,
     1                 DISCON, ETA, GI, H, HMIN
      DOUBLE PRECISION C2000, FOUR, HMAX0, HMIN0, H0, MACHEP, ONE, P002,
     1                 THREE, TWO, ZERO
C
      PARAMETER (C2000=2.0D+3, FOUR=4.0D+0, HMAX0=0.02D+0, HMIN0=5.0D+1,
     1     ONE=1.0D+0, P002=0.002D+0, THREE=3.0D+0,
     2     TWO=2.0D+0, ZERO=0.0D+0)
      PARAMETER (FH=3, FX0=4, HSAVE=5, XISAVE=6)
C
C---------------------------------  BODY  ------------------------------
C
      IF (IRC) 50, 10, 110
C
C     ***  FRESH START -- GET MACHINE-DEPENDENT CONSTANTS  ***
C
C     STORE MACHEP IN W(1) AND H0 IN W(2), WHERE MACHEP IS THE UNIT
C     ROUNDOFF (THE SMALLEST POSITIVE NUMBER SUCH THAT
C     1 + MACHEP .GT. 1  AND  1 - MACHEP .LT. 1),  AND  H0 IS THE
C     SQUARE-ROOT OF MACHEP.
C
 10   W(1) = DR7MDC(3)
      W(2) =  SQRT(W(1))
C
      W(FX0) = FX
C
C     ***  INCREMENT  I  AND START COMPUTING  G(I)  ***
C
 20   I =  ABS(IRC) + 1
      IF (I .GT. N) GO TO 120
         IRC = I
         AFX =  ABS(W(FX0))
         MACHEP = W(1)
         H0 = W(2)
         HMIN = HMIN0 * MACHEP
         W(XISAVE) = X(I)
         AXI =  ABS(X(I))
         AXIBAR =   MAX(AXI, ONE/D(I))
         GI = G(I)
         AGI =  ABS(GI)
         ETA =  ABS(ETA0)
         IF (AFX .GT. ZERO) ETA =   MAX(ETA, AGI*AXI*MACHEP/AFX)
         ALPHAI = ALPHA(I)
         IF (ALPHAI .EQ. ZERO) GO TO 80
         IF (GI .EQ. ZERO .OR. FX .EQ. ZERO) GO TO 90
         AFXETA = AFX*ETA
         AAI =  ABS(ALPHAI)
C
C        *** COMPUTE H = STEWART*S FORWARD-DIFFERENCE STEP SIZE.
C
         IF (GI**2 .LE. AFXETA*AAI) GO TO 30
              H = TWO* SQRT(AFXETA/AAI)
              H = H*(ONE - AAI*H/(THREE*AAI*H + FOUR*AGI))
              GO TO 40
 30      H = TWO*(AFXETA*AGI/(AAI**2))**(ONE/THREE)
         H = H*(ONE - TWO*AGI/(THREE*AAI*H + FOUR*AGI))
C
C        ***  ENSURE THAT  H  IS NOT INSIGNIFICANTLY SMALL  ***
C
 40      H =   MAX(H, HMIN*AXIBAR)
C
C        *** USE FORWARD DIFFERENCE IF BOUND ON TRUNCATION ERROR IS AT
C        *** MOST 10**-3.
C
         IF (AAI*H .LE. P002*AGI) GO TO 70
C
C        *** COMPUTE H = STEWART*S STEP FOR CENTRAL DIFFERENCE.
C
         DISCON = C2000*AFXETA
         H = DISCON/(AGI +  SQRT(GI**2 + AAI*DISCON))
C
C        ***  ENSURE THAT  H  IS NEITHER TOO SMALL NOR TOO BIG  ***
C
         H =   MAX(H, HMIN*AXIBAR)
         IF (H .GE. HMAX0*AXIBAR) H = AXIBAR * H0**(TWO/THREE)
C
C        ***  COMPUTE CENTRAL DIFFERENCE  ***
C
         IRC = -I
         GO TO 100
C
 50      H = -W(HSAVE)
         I =  ABS(IRC)
         IF (H .GT. ZERO) GO TO 60
         W(FH) = FX
         GO TO 100
C
 60      G(I) = (W(FH) - FX) / (TWO * H)
         X(I) = W(XISAVE)
         GO TO 20
C
C     ***  COMPUTE FORWARD DIFFERENCES IN VARIOUS CASES  ***
C
 70      IF (H .GE. HMAX0*AXIBAR) H = H0 * AXIBAR
         IF (ALPHAI*GI .LT. ZERO) H = -H
         GO TO 100
 80      H = AXIBAR
         GO TO 100
 90      H = H0 * AXIBAR
C
 100     X(I) = W(XISAVE) + H
         W(HSAVE) = H
         GO TO 999
C
C     ***  COMPUTE ACTUAL FORWARD DIFFERENCE  ***
C
 110     G(IRC) = (FX - W(FX0)) / W(HSAVE)
         X(IRC) = W(XISAVE)
         GO TO 20
C
C  ***  RESTORE FX AND INDICATE THAT G HAS BEEN COMPUTED  ***
C
 120  FX = W(FX0)
      IRC = 0
C
 999  RETURN
C  ***  LAST CARD OF DS7GRD FOLLOWS  ***
      END
      SUBROUTINE PCMRJ(NOBS, NPAR, X, NF, NEED, R, RP, UI, UR, UF)
      INTEGER NOBS, NPAR, NF, NEED(2), UI(*)
      DOUBLE PRECISION X(NPAR), R(NOBS), RP(NPAR,NOBS), UR(*)
      EXTERNAL UF
C
      EXTERNAL PCMRJ1
C
      INTEGER ICP, IICDAT, IICH, IIIV, IIRV, IIU, INALT, IRCDAT, IRU
C
C *** BODY ***
C
      IIU    = UI(1)
      IICH   = UI(2)
      INALT  = UI(3)
      IIIV   = UI(4)
      IIRV   = UI(5)
      IICDAT = UI(6)
C
      IRU = UI(7)
      ICP = UI(8)
C     IRW = UI(9)
      IRCDAT = UI(10)
C
      CALL PCMRJ1(NOBS, NPAR, X, NF, NEED, R, RP,
     1     UI(IIU), UI(IICH), UI(INALT), UI(IIIV), UI(IIRV), UI(IICDAT),
     2     UR(IRU), UR(ICP), UR(IRCDAT), UF)
 999  RETURN
C *** LAST LINE OF PCMRJ FOLLOWS ***
      END
      SUBROUTINE PCMRJ1(NOBS, NPAR, X, NF, NEED, R, RP,
     1     IUSER, ICHV, NALTV, IIV, IRV, ICDAT,
     2     RUSER, CPROB, RCDAT, UF)
C
C *** THIS SUBROUTINE EXPANDS THE STORAGE IN UI AND UR TO MAKE THEM ***
C *** COMPATIBLE WITH ESTIMATION OF CHOICE MODELS.                  ***
C
      INTEGER NOBS, NPAR, NF, NEED(2), IUSER(*), ICHV(NOBS),
     1        NALTV(NOBS), IIV(NOBS), IRV(NOBS), ICDAT(*)
      DOUBLE PRECISION X(NPAR), R(NOBS), RP(NPAR,NOBS), RUSER(*),
     1                 CPROB(NOBS,2), RCDAT(*)
      EXTERNAL UF
C
      EXTERNAL CALCDP, CALCPR, DL7SQR, DV7SCP
C
      INTEGER I, IALT, ICH, ICOV, IERR, II, III, IIR, IOBS, IOUNIT,
     1        IPCOV, IPP, IPU, IPUP, IR, ISCALE, ISIGP, ISIGU, J, KS,
     2        NALT, NALTM1, NFC
      DOUBLE PRECISION PROB, SCALE, XTEMP
C
      INTEGER LX
      DOUBLE PRECISION ONE, ZERO
      PARAMETER (ZERO=0.D0, ONE=1.D0, LX=60)
C
C ARRAYS:
C
C CPROB... VECTOR FOR STORING CHOICE PROBABILITIES.  CPROB(IOBS,J)
C             FOR J=1,2 STORES CHOICE PROBABILITIES FOR OBSERVATION
C             IOBS.  ONE IS THE CURRENT PROBABILITY, WHILE THE OTHER
C             ONE IS THE PROBABILITY COMPUTED AT THE PREVIOUS TRIAL
C             X.  THE CODE KEEPS TRACK OF WHICH IS WHICH USING THE
C             POINTERS STORED IN IUSER(1) AND IUSER(2).  THIS IS USED
C             IN VARIOUS WAYS TO MAKE COMPUTATION MORE EFFICIENT.
C ICHV.... VECTOR OF LENGTH NOBS.  ICHV(IOBS) CONTAINS THE INDEX OF
C             THE CHOSEN ALTERNATIVE FOR OBSERVATION IOBS.
C IIV..... VECTOR OF LENGHT NOBS.  IIV(IOBS) INDICATES THE NUMBER OF
C             INTEGER DATA VALUES STORED IN ICDAT FOR OBSERVATION IOBS.
C IRV..... VECTOR OF LENGHT NOBS.  IRV(IOBS) INDICATES THE NUMBER OF
C             REAL DATA VALUES STORED IN RCDAT FOR OBSERVATION IOBS.
C NALTV... VECTOR OF LENGHT NOBS.  NALTV(IOBS) INDICATES THE NUMBER OF
C             CHOICES AVAILABLE FOR OBSERVATION IOBS.
C
C *** BODY ***
C
      ICOV = IUSER(9)
      NALT = IUSER(6)
      NALTM1 = NALT - 1
      ISCALE = IUSER(18)
      SCALE = RUSER(ISCALE)
C
      IF (NEED(1).EQ.1) THEN
C
C *** CALCULATE RESIDUAL VECTOR ***
         KS = 1
         IF (NEED(2).EQ.IUSER(1)) KS = 2
         IUSER(KS) = NF
C
C      IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX
         IF (ICOV.NE.0) THEN
C            SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA:
             IPCOV = IUSER(13)
             XTEMP = X(IPCOV)
             X(IPCOV) = SCALE
             ISIGP = IUSER(15)
             CALL DL7SQR(NALTM1, RUSER(ISIGP), X(IPCOV))
             X(IPCOV) = XTEMP
C            "UNPACK" FOR EASIER ACCESS IN CALCPR:
             IPP = ISIGP - 1
             ISIGU = IUSER(16)
             CALL DV7SCP(NALT*NALT, RUSER(ISIGU), ZERO)
             IPUP = ISIGU - 1
             DO 10 I = 1, NALTM1
                IPU = I + IPUP
                DO 10 J = 1, I
                   IPP = IPP + 1
                   RUSER(IPU) = RUSER(IPP)
                   IPU = IPU + NALT
 10          CONTINUE
         ENDIF
         III = 1
         IIR = 1
         DO 20 IOBS = 1, NOBS
             ICH = ICHV(IOBS)
             IALT = NALTV(IOBS)
             II = IIV(IOBS)
             IR = IRV(IOBS)
             CALL CALCPR(NPAR, X, NFC, ICH, IALT, II, ICDAT(III),
     1            IR, RCDAT(IIR), PROB, IUSER, RUSER, UF)
             IF ((PROB.LE.ZERO).OR.(PROB.GT.ONE).OR.(NFC.EQ.0)) THEN
                NF = 0
                RETURN
             ENDIF
             R(IOBS) = PROB
             CPROB(IOBS,KS) = PROB
             III = III + II
             IIR = IIR + IR
 20      CONTINUE
      ELSE
C
C *** CALCULATE JACOBIAN OF RESIDUAL VECTOR ***
C
         KS = 1
         IF (IUSER(1).NE.NF) KS = 2
         IF (IUSER(KS).NE.NF) THEN
            IOUNIT = IUSER(3)
            WRITE(IOUNIT,*) ' PROBLEM WITH INITIAL ESTIMATE...'
         ENDIF
C
         III = 1
         IIR = 1
         DO 30 IOBS = 1, NOBS
             ICH = ICHV(IOBS)
             IALT = NALTV(IOBS)
             II = IIV(IOBS)
             IR = IRV(IOBS)
             PROB = CPROB(IOBS,KS)
             CALL CALCDP(NPAR, X, IERR, ICH, IALT, II, ICDAT(III),
     1             IR, RCDAT(IIR), PROB, RP(1,IOBS), IUSER, RUSER, UF)
             IF (IERR.EQ.0) THEN
                NF = 0
                RETURN
             ENDIF
             III = III + II
             IIR = IIR + IR
 30      CONTINUE
      ENDIF
 999  RETURN
C *** LAST LINE OF PCMRJ1 FOLLOWS ***
      END
      SUBROUTINE PCMRHO(NEED, F, NOBS, NF, XN, R, RD, UI, UR, W)
      INTEGER NEED(2), NOBS, NF, UI(*)
      DOUBLE PRECISION F, XN(*), R(*), RD(NOBS,*), UR(*), W(NOBS)
C
      INTEGER ICP, IOBS, IOUNIT, IRW, WEIGHT, KS
      DOUBLE PRECISION OOR, VT
C
      DOUBLE PRECISION NEGONE, ZERO
      PARAMETER (NEGONE=-1.D0, ZERO=0.D0)
C
C *** BODY ***
C
      WEIGHT = UI(14)
      IF (NEED(1).EQ.1) THEN
         VT = ZERO
         IF (WEIGHT.EQ.0) THEN
            DO 10 IOBS = 1, NOBS
                VT = VT - LOG(R(IOBS))
 10         CONTINUE
         ELSE
            IRW = UI(9)
            DO 20 IOBS = 1, NOBS
                VT = VT - UR(IRW) * LOG(R(IOBS))
                IRW = IRW + 1
 20         CONTINUE
         ENDIF
         F = VT
      ELSE
         KS = 1
         IF (UI(11).NE.NF) KS = 2
         IF (UI(10+KS).NE.NF) THEN
            IOUNIT = UI(13)
            WRITE(IOUNIT,*) ' PROBLEM WITH INITIAL POINT...'
            NF = 0
            RETURN
         ENDIF
         ICP = UI(8)
         IF (KS.EQ.2) ICP = ICP + NOBS
         IF (WEIGHT.EQ.0) THEN
            DO 30 IOBS = 1, NOBS
                OOR = NEGONE/UR(ICP)
                R(IOBS) = OOR
                W(IOBS) = R(IOBS) * OOR
                RD(IOBS,1) = W(IOBS)
                ICP = ICP + 1
 30         CONTINUE
         ELSE
            IRW = UI(9)
            DO 40 IOBS = 1, NOBS
                OOR = NEGONE/UR(ICP)
                R(IOBS) = UR(IRW) * OOR
                W(IOBS) = R(IOBS) * OOR
                RD(IOBS,1) = W(IOBS)
                ICP = ICP + 1
                IRW = IRW + 1
 40         CONTINUE
         ENDIF
      ENDIF
 999  RETURN
C *** LAST LINE OF PCMRHO FOLLOWS ***
      END
      SUBROUTINE FPRINT(NOBS, NPAR, X, NF, UI, UR, UF)
      INTEGER NOBS, NPAR, NF, UI(*)
      DOUBLE PRECISION X(NPAR), UR(*)
      EXTERNAL UF
C
      EXTERNAL FPRNT1
C
      INTEGER ICP, IICDAT, IICH, IIIV, IIRV, IIU, INALT, IRCDAT, IRU,
     1        IRW
C
C *** BODY ***
C
      IIU    = UI(1)
      IICH   = UI(2)
      INALT  = UI(3)
      IIIV   = UI(4)
      IIRV   = UI(5)
      IICDAT = UI(6)
C
      IRU = UI(7)
      ICP = UI(8)
      IRW = UI(9)
      IRCDAT = UI(10)
C
      CALL FPRNT1(NOBS, NPAR, X, NF,
     1     UI(IIU), UI(IICH), UI(INALT), UI(IIIV), UI(IIRV), UI(IICDAT),
     2     UR(IRU), UR(IRCDAT), UR(IRW), UF)
 999  RETURN
C *** LAST LINE OF FPRINT FOLLOWS ***
      END
      SUBROUTINE FPRNT1(NOBS, NPAR, X, NF,
     1     IUSER, ICHV, NALTV, IIV, IRV, ICDAT,
     2     RUSER, RCDAT, WT, UF)
C
C *** THIS SUBROUTINE EXPANDS THE STORAGE IN UI AND UR TO MAKE THEM ***
C *** COMPATIBLE WITH ESTIMATION OF CHOICE MODELS.                  ***
C *** SEE PCMRJ1 DOCUMENTATION ON ARRAYS.                           ***
C
      INTEGER NOBS, NPAR, NF, IUSER(*), ICHV(NOBS), NALTV(NOBS),
     1        IIV(NOBS), IRV(NOBS), ICDAT(*)
      DOUBLE PRECISION X(NPAR), RUSER(*), RCDAT(*), WT(NOBS)
      EXTERNAL UF
C
      EXTERNAL CALCPR, DL7SQR, DV7SCP
C
      INTEGER I, IALT, ICH, ICOV, ICSET, II, III, IIR, IOBS, IOUNIT,
     1        IPCOV, IPP, IPU, IPUP, IR, ISCALE, ISIGP, ISIGU, J, NALT,
     2        NALTM1, NFC
      DOUBLE PRECISION FPROB(20), PROB, SCALE, XTEMP
C
      INTEGER LX
      DOUBLE PRECISION ONE, ZERO
      PARAMETER (ZERO=0.D0, ONE=1.D0, LX=60)
C
C *** BODY ***
C
      ICOV = IUSER(9)
      ICSET = IUSER(5)
      IOUNIT = IUSER(3)
      NALT = IUSER(6)
      NALTM1 = NALT - 1
      ISCALE = IUSER(18)
      SCALE = RUSER(ISCALE)
C
      WRITE(IOUNIT, 10)
 10   FORMAT(//,' FINAL CHOICE SET PROBABILITIES: ',/)
C
C      IF ICOV.NE.0, SET UP AN UNPACKED SIGMA MATRIX
         IF (ICOV.NE.0) THEN
C            SQUARE THE CHOLESKY FACTOR TO GET (PACKED) SIGMA:
             IPCOV = IUSER(13)
             XTEMP = X(IPCOV)
             X(IPCOV) = SCALE
             ISIGP = IUSER(15)
             CALL DL7SQR(NALTM1, RUSER(ISIGP), X(IPCOV))
             X(IPCOV) = XTEMP
C            "UNPACK" FOR EASIER ACCESS IN CALCPR:
             IPP = ISIGP - 1
             ISIGU = IUSER(16)
             CALL DV7SCP(NALT*NALT, RUSER(ISIGU), ZERO)
             IPUP = ISIGU - 1
             DO 20 I = 1, NALTM1
                IPU = I + IPUP
                DO 20 J = 1, I
                   IPP = IPP + 1
                   RUSER(IPU) = RUSER(IPP)
                   IPU = IPU + NALT
 20          CONTINUE
         ENDIF
         III = 1
         IIR = 1
         DO 90 IOBS = 1, NOBS
             ICH = ICHV(IOBS)
             IALT = NALTV(IOBS)
             II = IIV(IOBS)
             IR = IRV(IOBS)
             DO 30 I = 1, IALT
                CALL CALCPR(NPAR, X, NFC, I, IALT, II, ICDAT(III),
     1               IR, RCDAT(IIR), PROB, IUSER, RUSER, UF)
                FPROB(I) = PROB
 30          CONTINUE
             WRITE(IOUNIT, 40) IOBS
 40          FORMAT(/,' IOBS: ',I4)
             IF (ICSET.EQ.0) WRITE(IOUNIT,50) (ICDAT(I),I=1,IALT)
 50          FORMAT('    CHOICE SET: ',20I3)
             WRITE(IOUNIT, 60) IALT, ICH, WT(IOBS)
 60          FORMAT('    NO. OF ALTS: ',I2,'    ICH: ',I2,
     1              '    WT: ',F7.3)
             WRITE(IOUNIT, 70) (FPROB(I),I=1,IALT)
 70          FORMAT('    PROBS: ',8F7.4,/,18X,8F7.4,/,18X,4F7.4)
             WRITE(IOUNIT, 80) FPROB(ICH)
 80          FORMAT('    PROB(ICH): ',F7.4)
             III = III + II
             IIR = IIR + IR
 90      CONTINUE
C
 999  RETURN
C *** LAST LINE OF FPRNT1 FOLLOWS ***
      END

.