!
! *********************************************************************************
!  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!  &&&&& gfortran -o housym3 houysym3.f95 -lcblas -lclapack &&&&&
!        export PATH=$PATH:/Users/keithtpoole/usr/local/bin/
!Macintosh HD/Users/usr/local/bin/gfortran
!
!        gfortran -o dw-brent dw-brent.f95 -lcblas -lclapack
!        ./dw-brent
!
!  EPSILON SCROLL BAR -- Command-X toggle-scroll-bar
!  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!
!  DW-NOMINATE-HANDS-STANDALONE.F95----DYNAMIC-WEIGHTED NOMINATE -- THIS PROGRAM 
!                     IMPLEMENTS THE DW-NOMINATE MODEL USING A GREATLY
!                     SIMPLIFIED STRUCTURE ON THE COMBINED HOUSE AND SENATE
!
!                        24 JUNE 2014              
!               REVISED: 19 JULY 2015
!                 
! *********************************************************************************
!
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION M2ND,MACHEP
      INTEGER ICONG(54001),ID1(54001),ISTATE(54001),IDIST(54001),&
              IPARTY(54001),LVOTE(3800),KWHERE(99999,127),&
              NUMCONG(127),NUMCONGRC(127),&
              JCONG(54001),JD1(54001),JSTATE(54001),JDIST(54001),&
              JPARTY(54001),JCONGUNIQUE(54001),JD1UNIQUE(54001),&
              JSTATEUNIQUE(54001),JDISTUNIQUE(54001),&
              JPARTYUNIQUE(54001),PRIN,&
              ICONGRC(111991),INUMRC(111991,2),LMISS(3800),&
              K010RC(111991,10),NXDATA(127),KBIGLOG(54001,2)
      DIMENSION XDATA(54001,4),DYN(111991,4),ZMID(111991,4),WEIGHT(5),&
                YY(150000),CUMNML(150000),ZDF(150000,4),&
                XDATAUNIQUE(54001,4),XBIGLOG(54001),X(20)
!
      CHARACTER*1 KSTA(54001,7),LNAME(54001,11),&
                  JSTA(54001,7),JNAME(54001,11),&
                  JSTAUNIQUE(54001,7),JNAMEUNIQUE(54001,11)
!
      LOGICAL*1 RCVOTE1(54001,3800),RCVOTE9(54001,3800),&
                LWHERE(99999,127),&
                RCVOTE1TEMP(54001,3800),RCVOTE9TEMP(54001,3800),&
                RCVOTE0TEMP(54001,3800)
!
!
      OPEN(20,FILE='HANDSL01113C20.DAT',&
                   STATUS='OLD')
      OPEN(21,FILE='DWNOM21.DAT')
      OPEN(23,FILE='HS01114_BEGIN_2015.VOT',STATUS='OLD')
      OPEN(24,FILE='HANDSC01113C20.DAT',STATUS='OLD')
      OPEN(19,FILE='DWNOM19.DAT')
      OPEN(26,FILE='DWNOM26.DAT')
      OPEN(27,FILE='DWNOM27.DAT')
      OPEN(28,FILE='DWNOM28.DAT')
      OPEN(29,FILE='DWNOM29.DAT')
      OPEN(30,FILE='DWNOM30.DAT')
      OPEN(40,FILE='DWNOM40.DAT')
      OPEN(41,FILE='DWNOM41.DAT')
      OPEN(42,FILE='DWNOM42.DAT')
      OPEN(43,FILE='HANDSL01114A20_STAND_ALONE.DAT')
      OPEN(49,FILE='HANDSC01114A20_STAND_ALONE.DAT')
!
!
  175 FORMAT(I3,I5,4F7.3)
  176 FORMAT(I7,I3,I5,I2,4F7.3)
  200 FORMAT(I4,I6,I3,I2,1X,7A1,1X,I4,1X,11A1,3800I1)
  201 FORMAT(I4,I6,I3,I2,1X,7A1,1X,I4,1X,11A1,2F7.3,3I6)
  202 FORMAT(I9,4F15.8)
  203 FORMAT(2I6,I4,I6,I3,I2,1X,7A1,1X,I4,1X,11A1,10F7.3)
  250 FORMAT(I4,I6,I7,I6,I7)
  300 FORMAT(' NUMBER OF UNIQUE LEGISLATORS              ',3I7)
  301 FORMAT(' TOTAL LEGISLATORS ROLL CALL FILE          ',3I7)
  302 FORMAT(' TOTAL LEGISLATORS LAST RELEASE            ',3I7)
  303 FORMAT(' TOTAL AND UNSCALED RCS LAST RELEASE       ',3I7)
  304 FORMAT(' NUMBER OF ROLL CALLS LAST CONGRESS        ',3I7)
  305 FORMAT(' TOTAL SCALABLE LEGISLATORS AFTER CLEAN    ',3I7)
  306 FORMAT(' TOTAL ROLL CALLS INCLUDING NEW CONGRESS   ',3I7)
  307 FORMAT(' NUMBER OF ROLL CALLS LAST HOUSE AND SENATE',3I7)
  318 FORMAT(5I6)
  319 FORMAT(' FIRST ESTIMATE OF BETA',4F16.7)
  320 FORMAT(' LAST ESTIMATE OF BETA ',4F16.7)
  321 FORMAT(' FIRST ESTIMATE OF WEIGHT_2',4F16.7)
  322 FORMAT(' LAST ESTIMATE OF WEIGHT_2',4F16.7)
  701 FORMAT(' CORRECT CLASSIFICATION AND APRE',3I10,3F10.4)
 1000 FORMAT(' START TIME OF JOB    ',3F16.7)
 1001 FORMAT(' END TIME OF JOB      ',3F16.7)
 1002 FORMAT(' ELAPSED TIME OF JOB  ',3F16.7)
!
      call cpu_time(result)
      write(*,1000)result
      write(21,1000)result
!
!  HOUSE 113 HAS 640 ROLL CALLS -- 1202 OVERALL
!  SENATE 113 HAS 657 ROLL CALLS
!
      NDIMENSIONS=2
      NFIRST=1
      NLAST=114
      KVMIN=25
      XVMIN=0.025
!
!  INITIALIZE BETA AND WEIGHTS
!
      DO 48 K=1,NDIMENSIONS
!
!  D-NOMINATE MODEL
!
      WEIGHT(K)=.500
  48  CONTINUE
      WEIGHT(1)=1.0
      WEIGHT(2)=0.41130
!  HOUSE AND SENATE
      WEIGHT(NDIMENSIONS+1)=7.83344
      NPMAX=54001
      NQMAX=3800
      NSMAX=4
      NQTOTALMAX=111991
      NCONGMAX=127
      NSESSIONMAX=600
!
      DO 46 J=1,127
      NUMCONG(J)=0
      NUMCONGRC(J)=0
      DO 47 I=1,99999
      LWHERE(I,J)=.FALSE.
      KWHERE(I,J)=0
  47  CONTINUE
  46  CONTINUE
!
!
!  READ IN CUMULATIVE NORMAL DISTRIBUTION FUNCTION
!
!
      ndevit=50001
      XDEVIT=10000.0
!
      DO 501 I=1,NDEVIT
      YY(I)=FLOAT(I-1)/XDEVIT
      X2=(FLOAT(I-1)/XDEVIT)/SQRT(2.0)
      XX=DERF(X2)
      XX=XX/2.0 +.5
      CUMNML(I)=XX
 501  CONTINUE
!
      TWOPI=1.0/SQRT(2.0*3.1415926536)
      DO 901 I=1,NDEVIT
      ZDF(I,1)=YY(NDEVIT+1-I)*(-1.0)
      ZDF(I,2)=1.0-CUMNML(NDEVIT+1-I)
      ZDF(I,3)=LOG(ZDF(I,2))
 901  CONTINUE
      DO 902 I=2,NDEVIT
      ZDF(I-1+NDEVIT,1)=YY(I)
      ZDF(I-1+NDEVIT,2)=CUMNML(I)
      ZDF(I-1+NDEVIT,3)=LOG(ZDF(I-1+NDEVIT,2))
 902  CONTINUE
      DO 903 I=1,2*NDEVIT-1
      ZDF(I,4)=(TWOPI*EXP((-ZDF(I,1)**2)/2.0))/ZDF(I,2)
!      WRITE(33,202)I,(ZDF(I,J),J=1,4)
 903  CONTINUE    
      NOISE=9966
!      IF(NOISE.EQ.9966)STOP
!
!  READ ROLL CALL FILE LAST RELEASE
!
      KK=0
      I=0
  552 READ(24,175,END=452)ICONGRC(I+1),INUMRC(I+1,1),&
                (DYN(I+1,K),ZMID(I+1,K),K=1,NDIMENSIONS)
      I=I+1
      NUMCONGRC(ICONGRC(I))=NUMCONGRC(ICONGRC(I))+1
!
!  CHECK TO SEE IF THE ROLL CALL WAS SCALED
!
!             INUMRC(.,2) = 0  SCALED
!             INUMRC(.,2) = 1  NOT SCALED
!
      INUMRC(I,2)=0
      SUM=0.0
      DO 13 JJ=1,NDIMENSIONS
      SUM=SUM+DABS(DYN(I,JJ))+DABS(ZMID(I,JJ))
  13  CONTINUE
      IF(SUM.LT.0.001)THEN
         INUMRC(I,2)=1
         KK=KK+1
      ENDIF
!
!
      GO TO 552
  452 WRITE(*,303)I,KK
      WRITE(21,303)I,KK
      NROLLCALLSOLD=I
      NTOHERE=0
      DO 2 I=1,NLAST-1
      NTOHERE=NTOHERE+NUMCONGRC(I)
      WRITE(21,250)I,NUMCONGRC(I),NTOHERE
  2   CONTINUE
      NUMCONGRC(NLAST)=NQMAX
      NTOHERE=NTOHERE+NUMCONGRC(NLAST)
      WRITE(21,250)NLAST,NUMCONGRC(NLAST),NTOHERE
!
!  READ STACKED ROLL CALL MATRIX -- 
!
      KMARKER=0
!
      I=0
  550 READ(23,200,END=450)ICONG(I+1),ID1(I+1),ISTATE(I+1),IDIST(I+1),&
                  (KSTA(I+1,JJ),JJ=1,7),IPARTY(I+1),&
                  (LNAME(I+1,JJ),JJ=1,11),&
                  (LVOTE(JJ),JJ=1,NUMCONGRC(ICONG(I+1)))
!
      I=I+1
!      WRITE(*,3000)I
! 3000 FORMAT(I10)
!
!  RECORD POSITIONS OF ALL LEGISLATORS
!
      LWHERE(ID1(I),ICONG(I))=.TRUE.
      KWHERE(ID1(I),ICONG(I))=I
!
!  FIND FIRST LEGISLATOR LAST CONGRESS
!
      IF(KMARKER.GT.0)GO TO 5
      IF(ICONG(I).EQ.NLAST)THEN
         KMARKER=I
      ENDIF
  5   CONTINUE
!
!  GET COUNTS OF THE NUMBER OF LEGISLATORS IN EACH CONGRESS
!
      NUMCONG(ICONG(I))=NUMCONG(ICONG(I))+1
!
      KMISS=0
!
      DO 1 JJ=1,NUMCONGRC(ICONG(I))
      RCVOTE0TEMP(I,JJ)=.FALSE.
      RCVOTE1TEMP(I,JJ)=.FALSE.
      RCVOTE9TEMP(I,JJ)=.FALSE.
      IF(LVOTE(JJ).GE.1.AND.LVOTE(JJ).LE.3)THEN
         RCVOTE1TEMP(I,JJ)=.TRUE.
      ENDIF
      IF(LVOTE(JJ).EQ.0.OR.LVOTE(JJ).GT.6)THEN
         RCVOTE9TEMP(I,JJ)=.TRUE.
         KMISS=KMISS+1
      ENDIF
      IF(ICONG(I).EQ.NLAST)THEN
         IF(LVOTE(JJ).EQ.0)THEN
            RCVOTE0TEMP(I,JJ)=.TRUE.
         ENDIF
      ENDIF
  1   CONTINUE
!
!  SET UP TO SEE HOW MANY ROLL CALLS ARE IN THE LAST CONGRESS
!
!      WRITE(*,318)I,ICONG(I),NUMCONGRC(ICONG(I)),KMISS,KMARKER
!
      GO TO 550
  450 WRITE(*,301)I
      WRITE(21,301)I
      NLEGISLATORS=I
!
!  GET UNIQUE NUMBER OF LEGISLATORS
! 
      NPUNIQUE=0
      DO 9999 IIII=1,99999
      KK=0
      DO 409 J=NFIRST,NLAST
      IF(LWHERE(IIII,J).EQV..TRUE.)THEN
         KK=KK+1
      ENDIF
  409 CONTINUE
      IF(KK.GT.0)NPUNIQUE=NPUNIQUE+1
 9999 CONTINUE
      WRITE(*,300)NLEGISLATORS,NPUNIQUE
      WRITE(21,300)NLEGISLATORS,NPUNIQUE
!
!  FIGURE OUT HOW MANY ROLL CALLS ARE IN THE LAST CONGRESS
!
      KKMISS=0
      DO 7 JJ=1,NQMAX
      KMISS=0
      DO 6 I=KMARKER,NLEGISLATORS
      IF(RCVOTE9TEMP(I,JJ).EQV..TRUE.)THEN
         KMISS=KMISS+1
!         IF(JJ.EQ.739)WRITE(26,1400)JJ,I,KMISS,KMARKER,NLEGISLATORS
      ENDIF
  6   CONTINUE
!      WRITE(26,1400)JJ,KMISS,NLEGISLATORS-KMARKER+1
      IF(KMISS.EQ.(NLEGISLATORS-KMARKER+1))THEN
         KKMISS=KKMISS+1
!         WRITE(26,1400)JJ,KKMISS
 1400 FORMAT(5I7)
      ENDIF
  7   CONTINUE
      NUMCONGRC(NLAST)=NQMAX-KKMISS
      WRITE(*,304)KKMISS,NQMAX-KKMISS,NUMCONGRC(NLAST)
      WRITE(21,304)KKMISS,NQMAX-KKMISS,NUMCONGRC(NLAST)
!
!  GET NUMBER OF ROLL CALLS FOR NLAST HOUSE AND NLAST SENATE
!
      KMISSHBIG=10000
      KMISSSBIG=10000
      DO 606 I=KMARKER,NLEGISLATORS
      KMISSH=0
      KMISSS=0
      DO 607 JJ=1,NUMCONGRC(NLAST)
      IF(IDIST(I).NE.0)THEN
         IF(RCVOTE0TEMP(I,JJ).EQV..TRUE.)THEN
            KMISSH=KMISSH+1
         ENDIF
      ENDIF
      IF(IDIST(I).EQ.0.AND.ISTATE(I).NE.99)THEN
         IF(RCVOTE0TEMP(I,JJ).EQV..TRUE.)THEN
            KMISSS=KMISSS+1
         ENDIF
      ENDIF
  607 CONTINUE
      IF(KMISSH.GT.0)THEN
         KMISSHMIN=MIN0(KMISSH,KMISSHBIG)
         KMISSHBIG=KMISSHMIN
      ENDIF
      IF(KMISSS.GT.0)THEN
         KMISSSMIN=MIN0(KMISSS,KMISSSBIG)
         KMISSSBIG=KMISSSMIN
       ENDIF
      WRITE(21,1400)I,KMISSHMIN,KMISSSMIN
!      WRITE(21,1400)I,KMISSH,KMISSS
  606 CONTINUE
      NRCSENLAST=KMISSHMIN
      NRCHOUSELAST=KMISSSMIN
      WRITE(21,307)NRCHOUSELAST,NRCSENLAST,NUMCONGRC(NLAST)
      WRITE(*,307)NRCHOUSELAST,NRCSENLAST,NUMCONGRC(NLAST)
!
      NTOHERE=0
      NTOHERE2=0
      DO 3 I=1,NLAST
      NTOHERE=NTOHERE+NUMCONG(I)
      NTOHERE2=NTOHERE2+NUMCONGRC(I)
      WRITE(21,250)I,NUMCONG(I),NTOHERE,NUMCONGRC(I),NTOHERE2
  3   CONTINUE
      NQTOTAL=NTOHERE2
!  
!  PUT CONGRESS AND COUNTER NUMBERS AT THE BOTTOM OF ICONGRC(.) AND INUMRC(.)
!
      DO 1999 JJ=1,NRCHOUSELAST
      ICONGRC(NROLLCALLSOLD+JJ)=NLAST
      INUMRC(NROLLCALLSOLD+JJ,1)=JJ
 1999 CONTINUE
      DO 2999 JJ=1,NRCSENLAST
      ICONGRC(NROLLCALLSOLD+JJ+NRCHOUSELAST)=NLAST
      INUMRC(NROLLCALLSOLD+JJ+NRCHOUSELAST,1)=JJ
 2999 CONTINUE
!
!  READ LAST RELEASE OF DW-NOMINATE -- LEGISLATORS 
!
!  READ LEGISLATOR FILE LAST RELEASE
!
      I=0
  551 READ(20,201,END=451)JCONG(I+1),JD1(I+1),JSTATE(I+1),JDIST(I+1),&
              (JSTA(I+1,JJ),JJ=1,7),JPARTY(I+1),&
              (JNAME(I+1,JJ),JJ=1,11),(XDATA(I+1,JJ),JJ=1,NDIMENSIONS) 
!
      I=I+1
!
      GO TO 551
  451 WRITE(*,302)I
      WRITE(21,302)I
!
      NLEGSCALED=I
!
      CALL CLEAN(NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,NQMAX,&
                 NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC,&
                 KVMIN,XVMIN,ICONG,JCONG,ID1,JD1,ISTATE,JSTATE,&
                 IDIST,JDIST,IPARTY,JPARTY,KSTA,LNAME,JSTA,JNAME,&
                 RCVOTE1,RCVOTE9,LWHERE,&
                 RCVOTE1TEMP,RCVOTE9TEMP,K010RC,INUMRC,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA)
!
!  WRITE OUT NUMBER OF SCALABLE LEGISLATORS FROM CLEAN SUBROUTINE
!
      WRITE(*,305)NLEGSCALED
      WRITE(21,305)NLEGSCALED
!
!  RESET LWHERE(:) AND KWHERE(:)
!
      DO 1949 J=1,127
      DO 1950 I=1,99999
      LWHERE(I,J)=.FALSE.
      KWHERE(I,J)=0
 1950 CONTINUE
 1949 CONTINUE
!
!  WRITE OUT SCALABLE LEGISLATORS
!
      DO 1947 I=1,NLEGSCALED
!
!  IF NOT MISSING DATA
!
      KCHOICE=0
      KYES=0
      KNO=0
      DO 1948 J=1,NUMCONGRC(JCONG(I))
!
! IF NOT MISSING DATA
!
      IF(RCVOTE9(I,J).EQV..FALSE.)THEN
         KCHOICE=KCHOICE+1
!
!  YEA VOTE
!
         IF(RCVOTE1(I,J).EQV..TRUE.)THEN
            KYES=KYES+1
         ENDIF
!
!  NAY VOTE
!
         IF(RCVOTE1(I,J).EQV..FALSE.)THEN
            KNO=KNO+1
         ENDIF
      ENDIF
1948 CONTINUE
!
      WRITE(40,201)JCONG(I),JD1(I),JSTATE(I),JDIST(I),&
                  (JSTA(I,JJ),JJ=1,7),JPARTY(I),&
                  (JNAME(I,JJ),JJ=1,11),&
                  (XDATA(I,K),K=1,NDIMENSIONS),&
                   KCHOICE,KYES,KNO
!
!  RECORD POSITIONS OF ALL SCALABLE LEGISLATORS
!
      LWHERE(JD1(I),JCONG(I))=.TRUE.
      KWHERE(JD1(I),JCONG(I))=I
 1947 CONTINUE
!
!  FIND UNIQUE LEGISLATORS THAT ARE SCALABLE
!
      NPUNIQUE=0
      DO 9998 IIII=1,99999
      KK=0
      KLOCATE=0
      DO 49 J=NFIRST,NLAST
      IF(LWHERE(IIII,J).EQV..TRUE.)THEN
         KK=KK+1
         KLOCATE=KWHERE(IIII,J)
      ENDIF
  49  CONTINUE
      IF(KK.EQ.0)GO TO 9998
         WRITE(48,203)KLOCATE,KK,JCONG(KLOCATE),JD1(KLOCATE),&
                  JSTATE(KLOCATE),JDIST(KLOCATE),&
                  (JSTA(KLOCATE,JJ),JJ=1,7),JPARTY(KLOCATE),&
                  (JNAME(KLOCATE,JJ),JJ=1,11),&
                  (XDATA(KLOCATE,K),K=1,NDIMENSIONS)
      NPUNIQUE=NPUNIQUE+1
!
!  TRANSFER INTO UNIQUE VECTOR LOCATIONS
!
      JCONGUNIQUE(NPUNIQUE)=JCONG(KLOCATE)
      JD1UNIQUE(NPUNIQUE)=JD1(KLOCATE)
      JSTATEUNIQUE(NPUNIQUE)=JSTATE(KLOCATE)
      JDISTUNIQUE(NPUNIQUE)=JDIST(KLOCATE)
      DO 67 K=1,7
      JSTAUNIQUE(NPUNIQUE,K)=JSTA(KLOCATE,K)
  67  CONTINUE
      JPARTYUNIQUE(NPUNIQUE)=JPARTY(KLOCATE)
      DO 68 K=1,11
      JNAMEUNIQUE(NPUNIQUE,K)=JNAME(KLOCATE,K)
  68  CONTINUE
      DO 69 K=1,NDIMENSIONS
      XDATAUNIQUE(NPUNIQUE,K)=XDATA(KLOCATE,K)
  69  CONTINUE
 9998 CONTINUE
      DO 70 I=1,NPUNIQUE      
      WRITE(48,201)JCONGUNIQUE(I),JD1UNIQUE(I),&
                  JSTATEUNIQUE(I),JDISTUNIQUE(I),&
                  (JSTAUNIQUE(I,JJ),JJ=1,7),JPARTYUNIQUE(I),&
                  (JNAMEUNIQUE(I,JJ),JJ=1,11),&
                  (XDATAUNIQUE(I,K),K=1,NDIMENSIONS)
  70  CONTINUE
      WRITE(*,300)NLEGISLATORS,NPUNIQUE
      WRITE(21,300)NLEGISLATORS,NPUNIQUE
!
!  TOTAL ROLL CALLS AFTER CLEAN
!
      WRITE(*,306)NQTOTAL
      WRITE(21,306)NQTOTAL
      DO 1789 IZ=NFIRST,NLAST
      WRITE(21,1790)IZ,NXDATA(IZ),NUMCONGRC(IZ)
!      WRITE(*,1790)IZ,NXDATA(IZ),NUMCONGRC(IZ)
 1789 CONTINUE
 1790 FORMAT(I4,5I6)
!
      DO 1914 I=1,NQTOTAL
      WRITE(50,176)I,ICONGRC(I),INUMRC(I,1),INUMRC(I,2),&
                (DYN(I,K),ZMID(I,K),K=1,NDIMENSIONS)
 1914 CONTINUE
!
!
!  TEST MINIMUM FINDER -- BETA
!

      a =  1.0
      b = 50.0
      t = 0.001
      KICK=NDIMENSIONS+1
!      call test_local_min_rc_one ( a, b, t, g_05, 'g_05(x) test' )
      call test_local_min_rc_one (KICK,KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a, b, t, 'Find BETA' )
      WRITE(21,319)XPLOG,&
                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)

!
!  TEST MINIMUM FINDER -- WEIGHT
!

      a =  .1
      b = 2.0
      t = 0.001
      KICK=2
!      call test_local_min_rc_one ( a, b, t, g_05, 'g_05(x) test' )
      call test_local_min_rc_one (KICK,KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a, b, t, 'Find WEIGHT' )
      WRITE(21,321)XPLOG,&
                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)
!
      CALL CPU_TIME(RESULT1)
      WRITE(*,1001)RESULT1
      WRITE(21,1001)RESULT1
      TOTALTIME=RESULT1-RESULT
      WRITE(*,1002)TOTALTIME
      WRITE(21,1002)TOTALTIME
!      k7moa=1947
!      if(k7moa.eq.1947)stop
!
!
!  CALL SUBROUTINE THAT IN TURN DRIVES LBFGS
!
!  THE LBFGS FORTRAN CODE WAS WRITTEN BY JORGE NOCEDAL.  SEE THE FOLLOWING REFERENCES FOR THE THEORY
!    J. Nocedal. Updating Quasi-Newton Matrices with Limited Storage (1980), Mathematics of Computation 35, pp. 773-782.
!    D.C. Liu and J. Nocedal. On the Limited Memory Method for Large Scale Optimization (1989), Mathematical Programming B, 45, 3, pp. 503-528.
!
      CALL JORGE_NOCEDAL(KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,&
                 NLEGSCALED,NPUNIQUE,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 ICONG,ID1,ISTATE,IDIST,KSTA,IPARTY,LNAME, &
                 JCONG,JD1,JSTATE,JDIST,JSTA,JPARTY,JNAME,       &
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,&
                 JCONGUNIQUE,JD1UNIQUE,JSTATEUNIQUE,JPARTYUNIQUE,&
                 JDISTUNIQUE,JSTAUNIQUE,JNAMEUNIQUE,XDATAUNIQUE,&
                 ZDF,NDEVIT,XDEVIT)
!
!
!  CALCULATE BETA AFTER ZIG-ZAG
!
!      
!
!  TEST MINIMUM FINDER
!

      a =  WEIGHT(NDIMENSIONS+1) - 1.5
      b =  WEIGHT(NDIMENSIONS+1) + 1.5
      t = 0.005
!      a =  1.0
!      b = 50.0
!      t = 0.001
      KICK=NDIMENSIONS+1
!      call test_local_min_rc_one ( a, b, t, g_05, 'g_05(x) test' )
      call test_local_min_rc_one (KICK,KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a, b, t, 'Find BETA' )
      WRITE(21,320)XPLOG,&
                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)
!
!  TEST MINIMUM FINDER -- WEIGHT
!

      a =  WEIGHT(NDIMENSIONS) - .1
      b =  WEIGHT(NDIMENSIONS) + .1
      t = 0.005
!      a =  .1
!      b = 2.0
!      t = 0.001
      KICK=2
!      call test_local_min_rc_one ( a, b, t, g_05, 'g_05(x) test' )
      call test_local_min_rc_one (KICK,KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a, b, t, 'Find WEIGHT' )
      WRITE(21,322)XPLOG,&
                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)
!
!  Calculate APRE
!
      LPRE=0
      KK=0
      DO 996 III=NFIRST,NLAST
      NQNOW=NUMCONGRC(III)
      DO 995 I=1,NQNOW
      KK=KK+1
         KYES=K010RC(KK,2)+K010RC(KK,3)+K010RC(KK,4)
         KNO =K010RC(KK,5)+K010RC(KK,6)+K010RC(KK,7)
         KRCTOT=KYES+KNO
         KRCMIN=MIN0(KYES,KNO)
         IF(KRCTOT.GT.0)THEN
            XMARG=FLOAT(KRCMIN)/FLOAT(KRCTOT)
            IF(XMARG.GE..025)THEN
               LPRE=LPRE+KRCMIN
            ENDIF
         ENDIF
!
  995 CONTINUE
!
  996 CONTINUE
!
!
!  CALCULATE THE FINAL LOG-LIKELIHOOD TO GET THE APRE
!
      CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,NLEGISLATORS,&
                 NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
!
      XCLASS=FLOAT(KLASS)/FLOAT(KTOT)
      XPRE=FLOAT(LPRE-(KTOT-KLASS))/FLOAT(LPRE)
      WRITE(21,701)KTOT,KLASS,LPRE,XCLASS,XPRE,GMP
!
      CALL CPU_TIME(RESULT1)
      WRITE(*,1001)RESULT1
      WRITE(21,1001)RESULT1
      TOTALTIME=RESULT1-RESULT
      WRITE(*,1002)TOTALTIME
      WRITE(21,1002)TOTALTIME
!
      stop
      end

!
!  **********************************************************************
!    SUBROUTINE CLEAN---CALLED BY MAIN.  PROCESSES *.VOT MATRIX AND
!     SETS UP THE UNTRANSPOSED AND TRANSPOSED ROLL CALL MATRICES
!     
!  **********************************************************************
!
      SUBROUTINE CLEAN(NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 KVMIN,XVMIN,ICONG,JCONG,ID1,JD1,ISTATE,JSTATE,&
                 IDIST,JDIST,IPARTY,JPARTY,KSTA,LNAME,JSTA,JNAME,&
                 RCVOTE1,RCVOTE9,LWHERE,&
                 RCVOTE1TEMP,RCVOTE9TEMP,K010RC,INUMRC,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER ICONG(NPMAX),ID1(NPMAX),ISTATE(NPMAX),IDIST(NPMAX),&
              IPARTY(NPMAX),KWHERE(99999,NCONGMAX),&
              NUMCONG(NCONGMAX),NUMCONGRC(NCONGMAX),&
              JCONG(NPMAX),JD1(NPMAX),JSTATE(NPMAX),JDIST(NPMAX),&
              JPARTY(NPMAX),K010RC(NQTOTALMAX,10),&
              INUMRC(NQTOTALMAX,2),&
              ICONGRC(NQTOTALMAX),NXDATA(NCONGMAX)
!
      DIMENSION XDATA(NPMAX,NSMAX),DYN(NQTOTALMAX,NSMAX),&
                ZMID(NQTOTALMAX,NSMAX)
!
      CHARACTER*1 KSTA(NPMAX,7),LNAME(NPMAX,11),&
                  JSTA(NPMAX,7),JNAME(NPMAX,11)
!
      LOGICAL*1 RCVOTE1(NPMAX,NQMAX),RCVOTE9(NPMAX,NQMAX),&
                LWHERE(99999,NCONGMAX),&
                RCVOTE1TEMP(NPMAX,NQMAX),RCVOTE9TEMP(NPMAX,NQMAX)
!
!
      INTEGER, ALLOCATABLE :: LVOTE(:)
      DOUBLE PRECISION, ALLOCATABLE :: XDEM(:)
      DOUBLE PRECISION, ALLOCATABLE :: XREP(:)
      DOUBLE PRECISION, ALLOCATABLE :: SUMYES(:)
      DOUBLE PRECISION, ALLOCATABLE :: SUMNO(:)
!
!      DOUBLE PRECISION, ALLOCATABLE :: SA(:,:)
!
      ALLOCATE(LVOTE(NQMAX))
      ALLOCATE(XDEM(NSMAX))
      ALLOCATE(XREP(NSMAX))
      ALLOCATE(SUMYES(NSMAX))
      ALLOCATE(SUMNO(NSMAX))
!
!      ALLOCATE(SA(NISSUES,NISSUES))
!
  200 FORMAT(I4,I6,I3,I2,1X,7A1,1X,I4,1X,11A1,3800I1)
  201 FORMAT(I4,2I5,I7,I3,10I5)
  202 FORMAT(I4,I6,I3,I2,1X,7A1,1X,I4,1X,11A1,2I6,4F7.3)
  203 FORMAT(I4,I6,I3,I2,1X,7A1,1X,I4,1X,11A1,12X,4F7.3)
  204 FORMAT(I4,2I5,I7,I3,10I5,2I4,12F7.3)
  300 FORMAT(' NO MATCH!!',3I6)
  301 FORMAT(' NUMBER OF MATCHES',I6)
  305 FORMAT(' SCALED 1 - NLAST-1 AND NLAST',2I7)
  306 FORMAT(' NUMBER MATCHED IN NLAST CONGRESS',2I7)
  307 FORMAT(' DEMOCRAT MEANS  ',I5,2F7.3)
  308 FORMAT(' REPUBLICAN MEANS',I5,2F7.3)
!
!
      DO 1947 I=1,NCONGMAX
      NXDATA(I)=0
 1947 CONTINUE
!
!  GO THROUGH ROLL CALL FILE AND FLAG LEGISLATORS WITHOUT ENOUGH VOTES TO SCALE (KVMIN)
!
      DO 1 I=1,NLEGISLATORS
      KMISS=0
      DO 2 J=1,NUMCONGRC(ICONG(I))
!
!  IF NOT MISSING DATA
!
      IF(RCVOTE9TEMP(I,J).EQV..FALSE.)THEN
         KMISS=KMISS+1
      ENDIF
  2   CONTINUE
!      IF(KMISS.GE.10.AND.ID1(I).GE.99869)GO TO 1
      IF(ID1(I).LT.99869.AND.KMISS.LT.KVMIN)THEN
            WRITE(26,202)ICONG(I),ID1(I),ISTATE(I),IDIST(I),&
                  (KSTA(I,JJ),JJ=1,7),IPARTY(I),&
                  (LNAME(I,JJ),JJ=1,11),KMISS,&
                  NUMCONGRC(ICONG(I))
            GO TO 1
      ENDIF
      IF(ID1(I).GE.99869.AND.KMISS.LT.10)THEN
            WRITE(26,202)ICONG(I),ID1(I),ISTATE(I),IDIST(I),&
                  (KSTA(I,JJ),JJ=1,7),IPARTY(I),&
                  (LNAME(I,JJ),JJ=1,11),KMISS,&
                  NUMCONGRC(ICONG(I))
            GO TO 1
      ENDIF
!
!  WRITE OUT POTENTIALLY SCALABLE LEGISLATORS
!
      WRITE(29,202)ICONG(I),ID1(I),ISTATE(I),IDIST(I),&
                  (KSTA(I,JJ),JJ=1,7),IPARTY(I),&
                  (LNAME(I,JJ),JJ=1,11),KMISS,&
                  NUMCONGRC(ICONG(I))

  1   CONTINUE
!
!  MATCH WITH PREVIOUS RELEASE
!
!
      IMATCH=0
      DO 3 I=1,NLEGSCALED
      DO 4 II=1,NLEGISLATORS
      IF(JCONG(I).EQ.ICONG(II).AND.JD1(I).EQ.ID1(II))THEN
         IMATCH=IMATCH+1
         KEEPII=II
         GO TO 44
      ENDIF
  4   CONTINUE
!
!  NO MATCH
!
      WRITE(*,300)I,JCONG(I),JD1(I)
      STOP
  44  CONTINUE
      KMISS=0
      DO 22 J=1,NUMCONGRC(ICONG(II))
!
!  IF NOT MISSING DATA
!
      IF(RCVOTE9TEMP(KEEPII,J).EQV..FALSE.)THEN
         KMISS=KMISS+1
      ENDIF
!
!  TRANSFER VOTE DATA
!
      RCVOTE9(I,J)=RCVOTE9TEMP(KEEPII,J)
      RCVOTE1(I,J)=RCVOTE1TEMP(KEEPII,J)
  22  CONTINUE
!
      IF(KMISS.LT.KVMIN)THEN
         WRITE(27,202)JCONG(I),JD1(I),JSTATE(I),JDIST(I),&
                  (JSTA(I,JJ),JJ=1,7),IPARTY(I),&
                  (JNAME(I,JJ),JJ=1,11),KMISS

      ENDIF
!
      IF(KMISS.GE.KVMIN)THEN
         WRITE(28,202)JCONG(I),JD1(I),JSTATE(I),JDIST(I),&
                  (JSTA(I,JJ),JJ=1,7),JPARTY(I),&
                  (JNAME(I,JJ),JJ=1,11),KMISS

      ENDIF
!
!  SPECIAL KLUDGE FOR THE PRESIDENTS
!         IF(KHIT.LT.25.AND.KHIT.GE.10.AND.ID1(I).GE.99869)THEN
!
      IF(KMISS.LT.KVMIN.AND.KMISS.GE.10.AND.JD1(I).GE.99869)THEN
         WRITE(28,202)JCONG(I),JD1(I),JSTATE(I),JDIST(I),&
                  (JSTA(I,JJ),JJ=1,7),JPARTY(I),&
                  (JNAME(I,JJ),JJ=1,11),KMISS

      ENDIF
!
      
  3   CONTINUE
      WRITE(*,301)IMATCH
!
!  GO THROUGH ROLL CALL FILE AND ADD LAST CONGRESS ON THE BOTTOM
!
      KK=0
      DO 6 I=1,NLEGISLATORS
      KMISS=0
      IF(ICONG(I).EQ.NLAST)THEN
         DO 7 J=1,NUMCONGRC(ICONG(I))
!
!  IF NOT MISSING DATA
!
         IF(RCVOTE9TEMP(I,J).EQV..FALSE.)THEN
            KMISS=KMISS+1
         ENDIF
  7      CONTINUE
!
         IF(ID1(I).LT.99869.AND.KMISS.LT.KVMIN)GO TO 6
         IF(ID1(I).GE.99869.AND.KMISS.LT.10)GO TO 6
!
!  TRANSFER INTO BOTTOM OF CURRENT RELEASE
!
         KK=KK+1
         JCONG(NLEGSCALED+KK)=ICONG(I)
         JD1(NLEGSCALED+KK)=ID1(I)
         JSTATE(NLEGSCALED+KK)=ISTATE(I)
         JDIST(NLEGSCALED+KK)=IDIST(I)
         JPARTY(NLEGSCALED+KK)=IPARTY(I)
         DO 8 JJ=1,7
         JSTA(NLEGSCALED+KK,JJ)=KSTA(I,JJ)
  8      CONTINUE
         DO 9 JJ=1,11
         JNAME(NLEGSCALED+KK,JJ)=LNAME(I,JJ)
  9      CONTINUE
!
!  TRANSFER VOTE DATA
!
         DO 10 JJ=1,NUMCONGRC(NLAST)
         RCVOTE9(NLEGSCALED+KK,JJ)=RCVOTE9TEMP(I,JJ)
         RCVOTE1(NLEGSCALED+KK,JJ)=RCVOTE1TEMP(I,JJ)
  10     CONTINUE
      ENDIF
  6   CONTINUE
!
       WRITE(21,305)NLEGSCALED,KK
       WRITE(*,305)NLEGSCALED,KK
!
!  GO THROUGH AND CHECK FINAL LEGISLATOR FILE
!
      KK1=NLEGSCALED
      KK2=KK
      DO 11 I=1,NLEGSCALED+KK
      KMISS=0
      DO 12 J=1,NUMCONGRC(JCONG(I))
!
!  IF NOT MISSING DATA
!
      IF(RCVOTE9(I,J).EQV..FALSE.)THEN
         KMISS=KMISS+1
      ENDIF
  12  CONTINUE
!
!  FILL NLAST CONGRESS LEGISLATOR COORDINATES WITH ZEROES
!
      IF(JCONG(I).EQ.NLAST)THEN
         DO 13 K=1,NDIMENSIONS
         XDATA(I,K)=0.0
  13     CONTINUE
      ENDIF
      WRITE(30,202)JCONG(I),JD1(I),JSTATE(I),JDIST(I),&
                  (JSTA(I,JJ),JJ=1,7),JPARTY(I),&
                  (JNAME(I,JJ),JJ=1,11),KMISS,&
                  NUMCONGRC(JCONG(I)),(XDATA(I,K),K=1,NDIMENSIONS)
      NXDATA(JCONG(I))=NXDATA(JCONG(I))+1
  11  CONTINUE
!      DO 1948 I=1,NLAST
!      WRITE(99,1949)I,NXDATA(I)
! 1948 CONTINUE
! 1949 FORMAT(I4,I7)
!
      NLEGSCALED=NLEGSCALED+KK
!
!  GET COUNTS FOR EACH ROLL CALL VOTE -- 0,1,2,3,4,5,6,7,8,9
!
      REWIND(23)
      NTOTAL=0
      DO 21 J=1,10
      DO 20 I=1,NQTOTAL
      K010RC(I,J)=0
  20  CONTINUE
  21  CONTINUE
      DO 999 III=NFIRST,NLAST
      NQNOW=NUMCONGRC(III)
      NPNOW=NUMCONG(III)
      DO 998 I=1,NPNOW
      READ(23,200)ICONG(I),ID1(I),ISTATE(I),IDIST(I),&
                  (KSTA(I,JJ),JJ=1,7),IPARTY(I),&
                  (LNAME(I,JJ),JJ=1,11),&
                  (LVOTE(JJ),JJ=1,NQNOW)
      DO 997 J=1,NQNOW
      K010RC(J+NTOTAL,LVOTE(J)+1)=K010RC(J+NTOTAL,LVOTE(J)+1)+1
  997 CONTINUE
  998 CONTINUE
      NTOTAL=NTOTAL+NQNOW
!
  999 CONTINUE
!
      KK=0
      DO 996 III=NFIRST,NLAST
      NQNOW=NUMCONGRC(III)
      DO 995 I=1,NQNOW
      KK=KK+1
      IF(III.EQ.NLAST)THEN
!         INUMRC(KK,1)=I
         KYES=K010RC(KK,2)+K010RC(KK,3)+K010RC(KK,4)
         KNO =K010RC(KK,5)+K010RC(KK,6)+K010RC(KK,7)
         KRCTOT=KYES+KNO
         KRCMIN=MIN0(KYES,KNO)
         INUMRC(KK,2)=0
         IF(KRCTOT.GT.0)THEN
            XMARG=FLOAT(KRCMIN)/FLOAT(KRCTOT)
            IF(XMARG.LT..025)THEN
               INUMRC(KK,2)=1
            ENDIF
         ENDIF
      ENDIF
!      WRITE(41,201)III,I,INUMRC(KK,1),KK,INUMRC(KK,2),&
!                   (K010RC(KK,JJ),JJ=1,10)
  995 CONTINUE
!
  996 CONTINUE
!
!  GET COORDINATES FOR LEGISLATORS LAST CONGRESS
!
      KMATCH=0
      DO 14 II=KK1+1,NLEGSCALED
      DO 15 III=KK1-1800,KK1
      IF(JD1(II).EQ.JD1(III))THEN
         DO 16 K=1,NDIMENSIONS
         XDATA(II,K)=XDATA(III,K)
  16     CONTINUE
         KMATCH=KMATCH+1
         GO TO 14
      ENDIF
  15  CONTINUE
  14  CONTINUE
      WRITE(21,306)KMATCH
      WRITE(*,306)KMATCH
!
      KDEM=0
      KREP=0
      DO 25 K=1,NDIMENSIONS
      XREP(K)=0.0
      XDEM(K)=0.0
  25  CONTINUE
      DO 17 I=KK1+1,NLEGSCALED
      WRITE(30,203)JCONG(I),JD1(I),JSTATE(I),JDIST(I),&
                  (JSTA(I,JJ),JJ=1,7),JPARTY(I),&
                  (JNAME(I,JJ),JJ=1,11),&
                  (XDATA(I,K),K=1,NDIMENSIONS)
      SUM=0.0
      DO 24 K=1,NDIMENSIONS
      SUM=SUM+DABS(XDATA(I,K))
  24  CONTINUE
!
!  GET PARTY MEANS NON-MISSING MEMBERS
!
      IF(SUM.GT.0.001)THEN
         IF(JPARTY(I).EQ.100)THEN
            DO 18 K=1,NDIMENSIONS
            XDEM(K)=XDEM(K)+XDATA(I,K)
  18        CONTINUE
            KDEM=KDEM+1
         ENDIF
         IF(JPARTY(I).EQ.200)THEN
            DO 19 K=1,NDIMENSIONS
            XREP(K)=XREP(K)+XDATA(I,K)
  19        CONTINUE
            KREP=KREP+1
         ENDIF
      ENDIF
  17  CONTINUE
!
!  COMPUTE PARTY MEANS
!
      DO 23 K=1,NDIMENSIONS
      XDEM(K)=XDEM(K)/FLOAT(KDEM)
      XREP(K)=XREP(K)/FLOAT(KREP)
  23  CONTINUE
!
!  INSERT PARTY MEANS INTO MISSING FIELDS LAST CONGRESS
!
      DO 29 I=KK1+1,NLEGSCALED
      SUM=0.0
      DO 28 K=1,NDIMENSIONS
      SUM=SUM+DABS(XDATA(I,K))
  28  CONTINUE
      IF(SUM.LE.0.001)THEN
         IF(JPARTY(I).EQ.100)THEN
            DO 26 K=1,NDIMENSIONS
            XDATA(I,K)=XDEM(K)
  26        CONTINUE
         ENDIF
         IF(JPARTY(I).EQ.200)THEN
            DO 27 K=1,NDIMENSIONS
            XDATA(I,K)=XREP(K)
  27        CONTINUE
         ENDIF
      ENDIF
      WRITE(30,203)JCONG(I),JD1(I),JSTATE(I),JDIST(I),&
                  (JSTA(I,JJ),JJ=1,7),JPARTY(I),&
                  (JNAME(I,JJ),JJ=1,11),&
                  (XDATA(I,K),K=1,NDIMENSIONS)

  29  CONTINUE
!
      WRITE(21,307)KDEM,(XDEM(K),K=1,NDIMENSIONS)
      WRITE(21,308)KREP,(XREP(K),K=1,NDIMENSIONS)
      WRITE(*,307)KDEM,(XDEM(K),K=1,NDIMENSIONS)
      WRITE(*,308)KREP,(XREP(K),K=1,NDIMENSIONS)
!
!  GET COORDINATES FOR ROLL CALLS LAST CONGRESS
!
      KK=0
      NPUP=0
      NPDOWN=0
      DO 30 III=NFIRST,NLAST
      NQNOW=NUMCONGRC(III)
      NPDOWN=NPUP+1
      NPUP=NPUP+NXDATA(III)
!
      DO 31 J=1,NQNOW
      DO 36 K=1,NDIMENSIONS
      SUMYES(K)=0.0
      SUMNO(K)=0.0
  36  CONTINUE
      KK=KK+1
      IF(III.EQ.NLAST)THEN
         DO 33 K=1,NDIMENSIONS
         DYN(KK,K)=0.0
         ZMID(KK,K)=0.0
  33     CONTINUE
      ENDIF
      KYES=0
      KNO=0
      DO 32 I=NPDOWN,NPUP
!
!  IF NOT MISSING DATA
!
      IF(RCVOTE9(I,J).EQV..FALSE.)THEN
!
!  YEA VOTE
!
         IF(RCVOTE1(I,J).EQV..TRUE.)THEN
            KYES=KYES+1
            DO 34 K=1,NDIMENSIONS
            SUMYES(K)=SUMYES(K)+XDATA(I,K)
  34        CONTINUE  
         ENDIF
!
!  NAY VOTE
!
         IF(RCVOTE1(I,J).EQV..FALSE.)THEN
            KNO=KNO+1
            DO 35 K=1,NDIMENSIONS
            SUMNO(K)=SUMNO(K)+XDATA(I,K)
  35        CONTINUE  
         ENDIF
      ENDIF
  32  CONTINUE
!
      DO 37 K=1,NDIMENSIONS
      IF(KYES.GT.0)SUMYES(K)=SUMYES(K)/FLOAT(KYES)
      IF(KNO.GT.0) SUMNO(K) = SUMNO(K)/FLOAT(KNO)
  37  CONTINUE
!
!      WRITE(41,204)III,J,INUMRC(KK,1),KK,INUMRC(KK,2),&
!                   (K010RC(KK,JJ),JJ=1,10),&
!                   KYES,KNO,&
!                   (DYN(KK,K),ZMID(KK,K),K=1,NDIMENSIONS),&
!                   (SUMYES(K),K=1,NDIMENSIONS),&                   
!                   (SUMNO(K),K=1,NDIMENSIONS)
!
!  COMPUTE ZMID AND DL VALUES FROM THE YEA AND NAY MEANS
!
      IF(III.EQ.NLAST.AND.INUMRC(KK,2).EQ.0)THEN
         DO 38 K=1,NDIMENSIONS
         ZMID(KK,K)=(SUMYES(K)+SUMNO(K))/2.0
         DYN(KK,K)=(SUMYES(K)-SUMNO(K))/(-2.0)
  38     CONTINUE
!         WRITE(41,204)III,J,INUMRC(KK,1),KK,INUMRC(KK,2),&
!                   (K010RC(KK,JJ),JJ=1,10),&
!                   KYES,KNO,&
!                   (DYN(KK,K),ZMID(KK,K),K=1,NDIMENSIONS),&
!                   (SUMYES(K),K=1,NDIMENSIONS),&                   
!                   (SUMNO(K),K=1,NDIMENSIONS)
      ENDIF
  31  CONTINUE
  30  CONTINUE
!
      DEALLOCATE(LVOTE)
      DEALLOCATE(XDEM)
      DEALLOCATE(XREP)
      DEALLOCATE(SUMYES)
      DEALLOCATE(SUMNO)
!
!      DEALLOCATE(SA)
!
      RETURN
      END
!
!  ***************************************************************************
!    SUBROUTINE PLOG---CALCULATES LOG-LIKELIHOOD FOR THE NOMINATE PROBIT   
!                      MODEL      
!  ***************************************************************************
!
      SUBROUTINE PLOG(KTOT,KLASS,GMP,XPLOG,&
                      NDIMENSIONS,NLEGISLATORS,&
                      NLEGSCALED,&
                      NFIRST,NLAST,NPMAX,&
                      NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                      NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                      RCVOTE1,RCVOTE9,LWHERE,&
                      INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                      XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER KWHERE(99999,NCONGMAX),&
              NUMCONG(NCONGMAX),NUMCONGRC(NCONGMAX),&
              INUMRC(NQTOTALMAX,2),&
              ICONGRC(NQTOTALMAX),NXDATA(NCONGMAX),KBIGLOG(NPMAX,2)
!
      DIMENSION XDATA(NPMAX,NSMAX),DYN(NQTOTALMAX,NSMAX),&
                ZMID(NQTOTALMAX,NSMAX),ZDF(150000,4),&
                WEIGHT(NSMAX+1),&
                XBIGLOG(NPMAX)
!
      LOGICAL*1 RCVOTE1(NPMAX,NQMAX),RCVOTE9(NPMAX,NQMAX),&
                LWHERE(99999,NCONGMAX)
!
      DOUBLE PRECISION, ALLOCATABLE :: DYES(:)
      DOUBLE PRECISION, ALLOCATABLE :: DNO(:)
      DOUBLE PRECISION, ALLOCATABLE :: DCC(:)
      DOUBLE PRECISION, ALLOCATABLE :: DBB(:)
!
      ALLOCATE(DYES(NSMAX))
      ALLOCATE(DNO(NSMAX))
      ALLOCATE(DCC(NSMAX))
      ALLOCATE(DBB(NSMAX))
!
  201 FORMAT(' HELLO DOLLY',I4,5I6)
 1002 FORMAT(' LOG-L',2I10,2F8.5,F20.6,10F10.5)
!  
      XPLOG=0.0
      XXPLOG=0.0
!
      KTOT=0
      KTOTP=0
      KTOTQ=0
      KLASS=0
      KLASSYY=0
      KLASSNY=0
      KLASSYN=0
      KLASSNN=0
      KLASS2=0
      DO 1 II=NFIRST,NLAST
      NPC=NXDATA(II)
      NQC=NUMCONGRC(II)
!      WRITE(21,201)II,NPC,NQC,KTOTP,KTOTQ
!      WRITE(*,201)II,NPC,NQC,KTOTP,KTOTQ
      DO 2 I=1,NPC
      YPLOG=0.0
      KPLOG=0
      KPWRONG=0
      DO 33 J=1,NQC
      IF(INUMRC(J+KTOTQ,2).EQ.0)THEN
         DO 3 K=1,NDIMENSIONS
         DYES(K)=0.0
         DNO(K)=0.0
         DYES(K)=(XDATA(I+KTOTP,K)-ZMID(J+KTOTQ,K)+DYN(J+KTOTQ,K))**2
         DNO(K) =(XDATA(I+KTOTP,K)-ZMID(J+KTOTQ,K)-DYN(J+KTOTQ,K))**2
  3      CONTINUE
!
!  IF NOT MISSING DATA
!
         IF(RCVOTE9(I+KTOTP,J).EQV..FALSE.)THEN
            KTOT=KTOT+1
!
!
!  IF YES
!
            IF(RCVOTE1(I+KTOTP,J).EQV..TRUE.)THEN
               DC=0.0
               DB=0.0
               DO 4 K=1,NDIMENSIONS
               DC=DC+(-WEIGHT(K)*WEIGHT(K)*DYES(K))
               DB=DB+(-WEIGHT(K)*WEIGHT(K)*DNO(K))
               DCC(K)=DYES(K)
               DBB(K)=DNO(K)
  4            CONTINUE
               XCC=+1.0
            ENDIF
!
!  IF NO
!
            IF(RCVOTE1(I+KTOTP,J).EQV..FALSE.)THEN
               DC=0.0
               DB=0.0
               DO 5 K=1,NDIMENSIONS
               DC=DC+(-WEIGHT(K)*WEIGHT(K)*DNO(K))
               DB=DB+(-WEIGHT(K)*WEIGHT(K)*DYES(K))
               DCC(K)=DNO(K)
               DBB(K)=DYES(K)
  5            CONTINUE
               XCC=-1.0
            ENDIF
!
            ZS=WEIGHT(NDIMENSIONS+1)*(DEXP(DC)-DEXP(DB))
!            ZS=WEIGHT(NDIMENSIONS+1)*(DEXP(DC/2.0)-DEXP(DB/2.0))
!            WRITE(99,9191)II,I,I+KTOTP,J,KTOT,XCC,ZS,DC,DB
! 9191 FORMAT(4I7,I10,F8.3,3F15.6)
!
!
!  CROSS CLASSIFICATIONS
!
            IF(DABS(DC).LE.DABS(DB))KLASS=KLASS+1
            IF(DABS(DC).LE.DABS(DB).AND.XCC.EQ.+1.0)THEN
               KLASSYY=KLASSYY+1
            ENDIF
            IF(DABS(DC).GT.DABS(DB).AND.XCC.EQ.+1.0)THEN
               KLASSNY=KLASSNY+1
            ENDIF
            IF(DABS(DC).GT.DABS(DB).AND.XCC.EQ.-1.0)THEN
               KLASSYN=KLASSYN+1
            ENDIF
            IF(DABS(DC).LE.DABS(DB).AND.XCC.EQ.-1.0)THEN
               KLASSNN=KLASSNN+1
            ENDIF
            IF(ZS.GT.0.0)KLASS2=KLASS2+1
!
            WWIMJ=ZS*XDEVIT
            KWIMJ=IDINT(DABS(WWIMJ)+.5)
            IF(KWIMJ.GT.NDEVIT-2)KWIMJ=NDEVIT-2
            IF(ZS.GE.0.0)THEN
               ZDISTF=ZDF(NDEVIT-1+KWIMJ+1,2)
               CDFLOG=ZDF(NDEVIT-1+KWIMJ+1,3)
            ENDIF
            IF(ZS.LT.0.0)THEN
               ZDISTF=ZDF(NDEVIT-KWIMJ,2)
               CDFLOG=ZDF(NDEVIT-KWIMJ,3)
               KPWRONG=KPWRONG+1
            ENDIF
!
            XPLOG=XPLOG+CDFLOG
            YPLOG=YPLOG+CDFLOG
            KPLOG=KPLOG+1
         ENDIF
      ENDIF
!
  33  CONTINUE
      XBIGLOG(I+KTOTP)=YPLOG
      KBIGLOG(I+KTOTP,1)=KPLOG
      KBIGLOG(I+KTOTP,2)=KPWRONG
      XXPLOG=XXPLOG+YPLOG
  2   CONTINUE
      KTOTP=KTOTP+NPC
      KTOTQ=KTOTQ+NQC
!      WRITE(21,201)II,NPC,NQC,KTOTP,KTOTQ
  1   CONTINUE 
      GMP=EXP(XXPLOG/FLOAT(KTOT))
      XCLASS=FLOAT(KLASS)/FLOAT(KTOT)
!      WRITE(21,1002)KTOT,KLASS,XCLASS,GMP,XXPLOG,&
!                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)
!      WRITE(99,1002)KTOT,KLASS,XCLASS,GMP,XXPLOG,&
!                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)
      WRITE(*,1002)KTOT,KLASS,XCLASS,GMP,XXPLOG,&
                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)
!
      DEALLOCATE(DYES)
      DEALLOCATE(DNO)
      DEALLOCATE(DCC)
      DEALLOCATE(DBB)
!
      RETURN
      END
!
!
!************************************************************************
!
!  ADJUST BLOCK DATA STRUCTURE BELOW:
!
!      DATA MP,LP,GTOL,STPMIN,STPMAX/41,42,9.0D-01,1.0D-20,1.0D+20/
!
!
!
      SUBROUTINE JORGE_NOCEDAL(KTOT,KLASS,GMP,XPLOG,                &
                 NDIMENSIONS,NLEGISLATORS,                          &
                 NLEGSCALED,NPUNIQUE,NFIRST,NLAST,NPMAX,            &
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,                   &
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC,                  &
                 RCVOTE1,RCVOTE9,LWHERE,                            &
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,                     &
                 ICONG,ID1,ISTATE,IDIST,KSTA,IPARTY,LNAME,          &
                 JCONG,JD1,JSTATE,JDIST,JSTA,JPARTY,JNAME,          &
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,                     &
                 JCONGUNIQUE,JD1UNIQUE,JSTATEUNIQUE,JPARTYUNIQUE,   &
                 JDISTUNIQUE,JSTAUNIQUE,JNAMEUNIQUE,XDATAUNIQUE,    &
                 ZDF,NDEVIT,XDEVIT)
!
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER KWHERE(99999,NCONGMAX),                               &
              NUMCONG(NCONGMAX),NUMCONGRC(NCONGMAX),                 &
              INUMRC(NQTOTALMAX,2),                                  &
              ICONGRC(NQTOTALMAX),NXDATA(NCONGMAX),ICONG(NPMAX),     &
              ID1(NPMAX),ISTATE(NPMAX),IDIST(NPMAX),IPARTY(NPMAX),   &
              JCONG(NPMAX),JD1(NPMAX),JSTATE(NPMAX),JDIST(NPMAX),    &
              JPARTY(NPMAX),KBIGLOG(NPMAX,2),                        &
              JCONGUNIQUE(NPMAX),JD1UNIQUE(NPMAX),JSTATEUNIQUE(NPMAX),&
              JDISTUNIQUE(NPMAX),JPARTYUNIQUE(NPMAX)
!
      DIMENSION XDATA(NPMAX,NSMAX),DYN(NQTOTALMAX,NSMAX),            &
                ZMID(NQTOTALMAX,NSMAX),ZDF(150000,4),                &
                WEIGHT(NSMAX+1),                                     &
                XDATA2(NLEGISLATORS*NDIMENSIONS),                    &
                ZDATAY2(NQTOTALMAX*NDIMENSIONS),                     &
                ZDATAN2(NQTOTALMAX*NDIMENSIONS),IPRINT(2),           &
                XDATAUNIQUE(NPMAX,NSMAX),XBIGLOG(NPMAX)
!
      LOGICAL*1 RCVOTE1(NPMAX,NQMAX),RCVOTE9(NPMAX,NQMAX),&
                LWHERE(99999,NCONGMAX)
!
      CHARACTER*1 KSTA(NPMAX,7),LNAME(NPMAX,11),&
                  JSTAUNIQUE(NPMAX,7),JNAMEUNIQUE(NPMAX,11),&
                  JSTA(NPMAX,7),JNAME(NPMAX,11)
!
      INTEGER EXITSTATUS
      LOGICAL DIAGCO
!
!     The driver for LBFGS must always declare LB2 as EXTERNAL
!
      EXTERNAL LB2
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
!
!
      DOUBLE PRECISION, ALLOCATABLE :: ZMAT(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XTRUE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: VCOV(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XDATASAVE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XDATASAVESQ(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ZDATAY(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ZDATAN(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ZDATAYSAVE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ZDATANSAVE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: WVEC(:)
      DOUBLE PRECISION, ALLOCATABLE :: X(:)
      DOUBLE PRECISION, ALLOCATABLE :: G(:)
      DOUBLE PRECISION, ALLOCATABLE :: DIAG(:)
      DOUBLE PRECISION, ALLOCATABLE :: W(:)
      DOUBLE PRECISION, ALLOCATABLE :: THETA(:)
!
      ALLOCATE(VCOV(NDIMENSIONS,NDIMENSIONS))
      ALLOCATE(XDATASAVE(NLEGISLATORS,NDIMENSIONS))
      ALLOCATE(XDATASAVESQ(NLEGISLATORS,NDIMENSIONS))
      ALLOCATE(ZDATAY(NQTOTALMAX,NDIMENSIONS))
      ALLOCATE(ZDATAN(NQTOTALMAX,NDIMENSIONS))
      ALLOCATE(ZDATAYSAVE(NQTOTALMAX,NDIMENSIONS))
      ALLOCATE(ZDATANSAVE(NQTOTALMAX,NDIMENSIONS))
      ALLOCATE(WVEC(NLEGISLATORS))
      ALLOCATE(X((NLEGISLATORS+2*NQTOTALMAX)*NDIMENSIONS))
      ALLOCATE(G((NLEGISLATORS+2*NQTOTALMAX)*NDIMENSIONS))
      ALLOCATE(DIAG((NLEGISLATORS+2*NQTOTALMAX)*NDIMENSIONS))
      ALLOCATE(W(((NLEGISLATORS+2*NQTOTALMAX)*NDIMENSIONS*(2*5+1)+2*5)))
!
!
      call cpu_time(result)
!      write(*,1000)result
!
      KSEED=INT(1000000*RESULT)
      ARAND=RAND(KSEED)
!
  100 FORMAT(4I7)
  201 FORMAT(2I7,20F10.4)
  202 FORMAT(4I7,10F10.4)
  301 FORMAT(' NUMBER DIMENSIONS, SCALABLE LEGISLATORS, AND SCALABLE RO&
               LL CALLS BEFORE LBFGS',3I7)
  302 FORMAT(' PRELIMINARY LOG-LIKELIHOOD',I7,F20.6)
  303 FORMAT(' LOG-LIKELIHOOD CHECK',4I7,F20.6)
  304 FORMAT(' NS, NQ, NP, N',4I10)
  305 FORMAT(' MAXIMUM SQUARED',F7.3)
  311 FORMAT(' ITERATION = ',I6,F20.6)
  320 FORMAT(' ESTIMATE OF BETA BETWEEN ZIGZAG    ',4F16.7)
  322 FORMAT(' ESTIMATE OF WEIGHT_2 BETWEEN ZIGZAG',4F16.7)
!
!
!  CALCULATE THE INITIAL LOG-LIKELIHOOD
!
      CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,NLEGISLATORS,&
                 NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
!
!  HOUSE DATA:
!    NQTOTAL=52968      -- TOTAL ROLL CALLS (THE UNANIMOUS ARE NOT FILTERED OUT)
!    NQ=46378           -- TOTAL SCALABLE ROLL CALLS
!    NLEGSCALED=37512   -- TOTAL SCALABLE LEGISLATORS -- ALSO TOTAL IN RCVOTE1(:) AND RCVOTE9(:)
!    NPUNIQUE=10724     -- TOTAL UNIQUE LEGISLATORS IN RCVOTE1(:), RCVOTE9(:), LWHERE(:), KWHERE(:)
!    NLEGISLATORS=38331 -- TOTAL LEGISLATORS IN THE STACKED ROLL CALL FILE
!    
     
!  GET NQ VALUE = NUMBER OF SCALABLE ROLL CALLS -- SET NP, NS, AND NQ HERE
! 
      NS=NDIMENSIONS
      NP=NLEGSCALED
      NQ=0
      DO 1 I=1,NQTOTAL 
      IF(INUMRC(I,2).EQ.0)NQ=NQ+1
   1  CONTINUE
      WRITE(*,301)NS,NP,NQ
!
      M=5
      ALPHA=0.95
      BETA=WEIGHT(NDIMENSIONS+1)
!     
! 80309.033641 from NOMINATE_LBFGSC.FOR
!
      IPRINT(1)= 1
      IPRINT(2)= 0
      JPRINT=1
!

      IF(JPRINT.EQ.1)WRITE(21,100)NLEGISLATORS,NQTOTALMAX,NP,NQ
!
!
      IF(JPRINT.EQ.1)WRITE(21,304)NS,NQ,NP
      IF(JPRINT.EQ.1)WRITE(*,304)NS,NQ,NP
!
!  CALCULATE YEA AND NAY OUTCOME COORDINATES
!
      KK=0
      DO 85 J=1,NQTOTAL
      DO 87 K=1,NS
      ZDATAY(J,K)=0.0
      ZDATAN(J,K)=0.0
  87  CONTINUE
      IF(INUMRC(J,2).EQ.0)THEN
         KK=KK+1
         KMARKY=0
         KMARKN=0
         DO 89 K=1,NS
         ZDATAY(J,K)=ZMID(J,K)-DYN(J,K)
         ZDATAN(J,K)=ZMID(J,K)+DYN(J,K)
  89     CONTINUE
         WRITE(19,201)J,KK,(ZDATAY(J,K),K=1,NS),&
                           (ZDATAN(J,K),K=1,NS),&
                      (DYN(J,K),ZMID(J,K),K=1,NS)
      ENDIF
  85  CONTINUE
      WRITE(*,309)KK
!      WRITE(41,309)KK
  309 FORMAT(' NUMBER SCALABLE ROLL CALLS',I7)
!
!  DO ZIZ-ZAG BETWEEN THE ROLLCALL AND LEGISLATOR COORDINATES USING LBFGS
!
      DO 999 IIII=1,5
      NIIII=IIII
!
!  CALCULATE YEA AND NAY OUTCOME COORDINATES
!
      KK=0
      DO 885 J=1,NQTOTAL
      IF(INUMRC(J,2).EQ.1)THEN
         DO 887 K=1,NS
         ZDATAY(J,K)=0.0
         ZDATAN(J,K)=0.0
 887     CONTINUE
      ENDIF
      IF(INUMRC(J,2).EQ.0)THEN
!
         KK=KK+1
         KMARKY=0
         KMARKN=0
         DO 889 K=1,NS
         ZDATAY(J,K)=ZMID(J,K)-DYN(J,K)
         ZDATAN(J,K)=ZMID(J,K)+DYN(J,K)
 889     CONTINUE
         WRITE(19,201)J,KK,(ZDATAY(J,K),K=1,NS),&
                           (ZDATAN(J,K),K=1,NS),&
                      (DYN(J,K),ZMID(J,K),K=1,NS)
      ENDIF
 885  CONTINUE
      WRITE(*,309)KK
      CALL ROLLCALL_ZIGZAG(NIIII,NFIRST,NLAST,NPUNIQUE,&
                           NLEGISLATORS,NQTOTAL,NSMAX,&
                           NQTOTALMAX,NDIMENSIONS,&
                           NP,NQ,NS,M,JPRINT,IPRINT,&
                           XDATA,ZDATAY,ZDATAN,ALPHA,NCONGMAX,&
                           NPMAX,NQMAX,NXDATA,NUMCONG,NUMCONGRC,&
                           INUMRC,DYN,ZMID,ICONGRC,&
                           RCVOTE1,RCVOTE9,LWHERE,KWHERE,&
                           WEIGHT,&
                 JCONGUNIQUE,JD1UNIQUE,JSTATEUNIQUE,JPARTYUNIQUE,   &
                 JDISTUNIQUE,JSTAUNIQUE,JNAMEUNIQUE,XDATAUNIQUE,    &
                 ZDF,NDEVIT,XDEVIT,EXITSTATUS)
!

      a =  WEIGHT(NDIMENSIONS+1) - 1.5
      b =  WEIGHT(NDIMENSIONS+1) + 1.5
      t = 0.005
      KICK=NDIMENSIONS+1
!      call test_local_min_rc_one ( a, b, t, g_05, 'g_05(x) test' )
      call test_local_min_rc_one (KICK,KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a, b, t, 'Find BETA' )
      WRITE(21,320)XPLOG,&
                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)
!
!  TEST MINIMUM FINDER -- WEIGHT
!

      a =  WEIGHT(NDIMENSIONS) - .1
      b =  WEIGHT(NDIMENSIONS) + .1
      t = 0.005
      KICK=2
!      call test_local_min_rc_one ( a, b, t, g_05, 'g_05(x) test' )
      call test_local_min_rc_one (KICK,KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a, b, t, 'Find WEIGHT' )
      WRITE(21,322)XPLOG,&
                    (WEIGHT(JJ),JJ=1,NDIMENSIONS+1)
!
!
      CALL X_ZIGZAG(NIIII,KTOT,KLASS,GMP,NFIRST,NLAST,NPUNIQUE,&
                           NLEGSCALED,NLEGISLATORS,NQTOTAL,NSMAX,&
                           NQTOTALMAX,NDIMENSIONS,&
                           NP,NQ,NS,M,JPRINT,IPRINT,&
                           XDATA,ZDATAY,ZDATAN,ALPHA,NCONGMAX,&
                           NPMAX,NQMAX,NXDATA,NUMCONG,NUMCONGRC,&
                           INUMRC,DYN,ZMID,ICONGRC,&
                           RCVOTE1,RCVOTE9,LWHERE,KWHERE,&
                           ICONG,ID1,ISTATE,IDIST,KSTA,IPARTY,&
                           JCONG,JD1,JSTATE,JDIST,JSTA,JPARTY,&
                           JNAME,LNAME,WEIGHT,KBIGLOG,XBIGLOG,&
                 JCONGUNIQUE,JD1UNIQUE,JSTATEUNIQUE,JPARTYUNIQUE,   &
                 JDISTUNIQUE,JSTAUNIQUE,JNAMEUNIQUE,XDATAUNIQUE,    &
                          ZDF,NDEVIT,XDEVIT,EXITSTATUS)
!
 999  CONTINUE
!
!
      EXITSTATUS=0
!
      DEALLOCATE(VCOV) 
      DEALLOCATE(XDATASAVE) 
      DEALLOCATE(XDATASAVESQ) 
      DEALLOCATE(ZDATAY) 
      DEALLOCATE(ZDATAN) 
      DEALLOCATE(WVEC)
      DEALLOCATE(X)
      DEALLOCATE(G)
      DEALLOCATE(DIAG)
      DEALLOCATE(W)
!
      RETURN
      END
!
! 
!
!     ----------------------------------------------------------------------
!     This file contains the LBFGS algorithm and supporting routines
!
!     ****************
!     LBFGS SUBROUTINE
!     ****************
!
      SUBROUTINE LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
!
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER N,M,IPRINT(2),IFLAG
      DOUBLE PRECISION X(N),G(N),DIAG(N),W(N*(2*M+1)+2*M)
      DOUBLE PRECISION F,EPS,XTOL
      LOGICAL DIAGCO
!
!        LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
!                          JORGE NOCEDAL
!                        *** July 1990 ***
!
! 
!     This subroutine solves the unconstrained minimization problem
! 
!                      min F(x),    x= (x1,x2,...,xN),
!
!      using the limited memory BFGS method. The routine is especially
!      effective on problems involving a large number of variables. In
!      a typical iteration of this method an approximation Hk to the
!      inverse of the Hessian is obtained by applying M BFGS updates to
!      a diagonal matrix Hk0, using information from the previous M steps.
!      The user specifies the number M, which determines the amount of
!      storage required by the routine. The user may also provide the
!      diagonal matrices Hk0 if not satisfied with the default choice.
!      The algorithm is described in "On the limited memory BFGS method
!      for large scale optimization", by D. Liu and J. Nocedal,
!      Mathematical Programming B 45 (1989) 503-528.
! 
!      The user is required to calculate the function value F and its
!      gradient G. In order to allow the user complete control over
!      these computations, reverse  communication is used. The routine
!      must be called repeatedly under the control of the parameter
!      IFLAG. 
!
!      The steplength is determined at each iteration by means of the
!      line search routine MCVSRCH, which is a slight modification of
!      the routine CSRCH written by More' and Thuente.
! 
!      The calling statement is 
! 
!          CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
! 
!      where
! 
!     N       is an INTEGER variable that must be set by the user to the
!             number of variables. It is not altered by the routine.
!             Restriction: N>0.
!
!     M       is an INTEGER variable that must be set by the user to
!             the number of corrections used in the BFGS update. It
!             is not altered by the routine. Values of M less than 3 are
!             not recommended; large values of M will result in excessive
!             computing time. 3<= M <=7 is recommended. Restriction: M>0.
! 
!     X       is a DOUBLE PRECISION array of length N. On initial entry
!             it must be set by the user to the values of the initial
!             estimate of the solution vector. On exit with IFLAG=0, it
!             contains the values of the variables at the best point
!             found (usually a solution).
! 
!     F       is a DOUBLE PRECISION variable. Before initial entry and on
!             a re-entry with IFLAG=1, it must be set by the user to
!             contain the value of the function F at the point X.
! 
!     G       is a DOUBLE PRECISION array of length N. Before initial
!             entry and on a re-entry with IFLAG=1, it must be set by
!             the user to contain the components of the gradient G at
!             the point X.
! 
!     DIAGCO  is a LOGICAL variable that must be set to .TRUE. if the
!             user  wishes to provide the diagonal matrix Hk0 at each
!             iteration. Otherwise it should be set to .FALSE., in which
!             case  LBFGS will use a default value described below. If
!             DIAGCO is set to .TRUE. the routine will return at each
!             iteration of the algorithm with IFLAG=2, and the diagonal
!              matrix Hk0  must be provided in the array DIAG.
! 
! 
!     DIAG    is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE.,
!             then on initial entry or on re-entry with IFLAG=2, DIAG
!             it must be set by the user to contain the values of the 
!             diagonal matrix Hk0.  Restriction: all elements of DIAG
!             must be positive.
! 
!     IPRINT  is an INTEGER array of length two which must be set by the
!             user.
! 
!             IPRINT(1) specifies the frequency of the output:
!                IPRINT(1) < 0 : no output is generated,
!                IPRINT(1) = 0 : output only at first and last iteration,
!                IPRINT(1) > 0 : output every IPRINT(1) iterations.
! 
!             IPRINT(2) specifies the type of output generated:
!                IPRINT(2) = 0 : iteration count, number of function 
!                                evaluations, function value, norm of the
!                                gradient, and steplength,
!                IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of
!                                variables and  gradient vector at the
!                                initial point,
!                IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of
!                                variables,
!                IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector.
! 
! 
!     EPS     is a positive DOUBLE PRECISION variable that must be set by
!             the user, and determines the accuracy with which the solution
!             is to be found. The subroutine terminates when
!
!                         ||G|| < EPS max(1,||X||),
!
!             where ||.|| denotes the Euclidean norm.
! 
!     XTOL    is a  positive DOUBLE PRECISION variable that must be set by
!             the user to an estimate of the machine precision (e.g.
!             10**(-16) on a SUN station 3/60). The line search routine will
!             terminate if the relative width of the interval of uncertainty
!             is less than XTOL.
! 
!     W       is a DOUBLE PRECISION array of length N(2M+1)+2M used as
!             workspace for LBFGS. This array must not be altered by the
!             user.
! 
!     IFLAG   is an INTEGER variable that must be set to 0 on initial entry
!             to the subroutine. A return with IFLAG<0 indicates an error,
!             and IFLAG=0 indicates that the routine has terminated without
!             detecting errors. On a return with IFLAG=1, the user must
!             evaluate the function F and gradient G. On a return with
!             IFLAG=2, the user must provide the diagonal matrix Hk0.
! 
!             The following negative values of IFLAG, detecting an error,
!             are possible:
! 
!              IFLAG=-1  The line search routine MCSRCH failed. The
!                        parameter INFO provides more detailed information
!                        (see also the documentation of MCSRCH):
!
!                       INFO = 0  IMPROPER INPUT PARAMETERS.
!
!                       INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF
!                                 UNCERTAINTY IS AT MOST XTOL.
!
!                       INFO = 3  MORE THAN 20 FUNCTION EVALUATIONS WERE
!                                 REQUIRED AT THE PRESENT ITERATION.
!
!                       INFO = 4  THE STEP IS TOO SMALL.
!
!                       INFO = 5  THE STEP IS TOO LARGE.
!
!                       INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS. 
!                                 THERE MAY NOT BE A STEP WHICH SATISFIES
!                                 THE SUFFICIENT DECREASE AND CURVATURE
!                                 CONDITIONS. TOLERANCES MAY BE TOO SMALL.
!
! 
!              IFLAG=-2  The i-th diagonal element of the diagonal inverse
!                        Hessian approximation, given in DIAG, is not
!                        positive.
!           
!              IFLAG=-3  Improper input parameters for LBFGS (N or M are
!                        not positive).
! 
!
!
!    ON THE DRIVER:
!
!    The program that calls LBFGS must contain the declaration:
!
!                       EXTERNAL LB2
!
!    LB2 is a BLOCK DATA that defines the default values of several
!    parameters described in the COMMON section. 
!
! 
! 
!    COMMON:
! 
!     The subroutine contains one common area, which the user may wish to
!    reference:
! 
         COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
! 
!    MP  is an INTEGER variable with default value 6. It is used as the
!        unit number for the printing of the monitoring information
!        controlled by IPRINT.
! 
!    LP  is an INTEGER variable with default value 6. It is used as the
!        unit number for the printing of error messages. This printing
!        may be suppressed by setting LP to a non-positive value.
! 
!    GTOL is a DOUBLE PRECISION variable with default value 0.9, which
!        controls the accuracy of the line search routine MCSRCH. If the
!        function and gradient evaluations are inexpensive with respect
!        to the cost of the iteration (which is sometimes the case when
!        solving very large problems) it may be advantageous to set GTOL
!        to a small value. A typical small value is 0.1.  Restriction:
!        GTOL should be greater than 1.D-04.
! 
!    STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which
!        specify lower and uper bounds for the step in the line search.
!        Their default values are 1.D-20 and 1.D+20, respectively. These
!        values need not be modified unless the exponents are too large
!        for the machine being used, or unless the problem is extremely
!        badly scaled (in which case the exponents should be increased).
! 
!
!  MACHINE DEPENDENCIES
!
!        The only variables that are machine-dependent are XTOL,
!        STPMIN and STPMAX.
! 
!
!  GENERAL INFORMATION
! 
!    Other routines called directly:  DAXPY, DDOT, LB1, MCSRCH
! 
!    Input/Output  :  No input; diagnostic messages on unit MP and
!                     error messages on unit LP.
! 
! 
!     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
      DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,DDOT,STP1,FTOL,STPMIN,&
                       STPMAX,STP,YS,YY,SQ,YR,BETA,XNORM
      INTEGER MP,LP,ITER,NFUN,POINT,ISPT,IYPT,MAXFEV,INFO,&
              BOUND,NPT,CP,I,NFEV,INMC,IYCN,ISCN
      LOGICAL FINISH
!
      SAVE
      DATA ONE,ZERO/1.0D+0,0.0D+0/
!
!     INITIALIZE
!     ----------
!
      IF(IFLAG.EQ.0) GO TO 10
      GO TO (172,100) IFLAG
  10  ITER= 0
      IF(N.LE.0.OR.M.LE.0) GO TO 196
      IF(GTOL.LE.1.D-04) THEN
        IF(LP.GT.0) WRITE(LP,245)
        GTOL=9.D-01
      ENDIF
      NFUN= 1
      POINT= 0
      FINISH= .FALSE.
      IF(DIAGCO) THEN
         DO 30 I=1,N
 30      IF (DIAG(I).LE.ZERO) GO TO 195
      ELSE
         DO 40 I=1,N
 40      DIAG(I)= 1.0D0
      ENDIF
!
!     THE WORK VECTOR W IS DIVIDED AS FOLLOWS:
!     ---------------------------------------
!     THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND
!         OTHER TEMPORARY INFORMATION.
!     LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO.
!     LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED
!         IN THE FORMULA THAT COMPUTES H*G.
!     LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH
!         STEPS.
!     LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M
!         GRADIENT DIFFERENCES.
!
!     THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A
!     CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT.
!
      ISPT= N+2*M
      IYPT= ISPT+N*M     
      DO 50 I=1,N
 50   W(ISPT+I)= -G(I)*DIAG(I)
      GNORM= DSQRT(DDOT(N,G,1,G,1))
      STP1= ONE/GNORM
!
!     PARAMETERS FOR LINE SEARCH ROUTINE
!     
!      FTOL= 1.1
      FTOL= 1.0D-4
      MAXFEV= 20
      XNORM=1.0
!
      IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,&
                       XNORM,GNORM,N,M,X,F,G,STP,FINISH)
!
!    --------------------
!     MAIN ITERATION LOOP
!    --------------------
!
 80   ITER= ITER+1
      INFO=0
      BOUND=ITER-1
      IF(ITER.EQ.1) GO TO 165
      IF (ITER .GT. M)BOUND=M
!
         YS= DDOT(N,W(IYPT+NPT+1),1,W(ISPT+NPT+1),1)
      IF(.NOT.DIAGCO) THEN
         YY= DDOT(N,W(IYPT+NPT+1),1,W(IYPT+NPT+1),1)
         DO 90 I=1,N
   90    DIAG(I)= YS/YY
      ELSE
         IFLAG=2
         RETURN
      ENDIF
 100  CONTINUE
      IF(DIAGCO) THEN
        DO 110 I=1,N
 110    IF (DIAG(I).LE.ZERO) GO TO 195
      ENDIF
!
!     COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980,
!     "Updating quasi-Newton matrices with limited storage",
!     Mathematics of Computation, Vol.24, No.151, pp. 773-782.
!     ---------------------------------------------------------
!
      CP= POINT
      IF (POINT.EQ.0) CP=M
      W(N+CP)= ONE/YS
      DO 112 I=1,N
 112  W(I)= -G(I)
      CP= POINT
      DO 125 I= 1,BOUND
         CP=CP-1
         IF (CP.EQ. -1)CP=M-1
         SQ= DDOT(N,W(ISPT+CP*N+1),1,W,1)
         INMC=N+M+CP+1
         IYCN=IYPT+CP*N
         W(INMC)= W(N+CP+1)*SQ
         CALL DAXPY(N,-W(INMC),W(IYCN+1),1,W,1)
 125  CONTINUE
!
      DO 130 I=1,N
 130  W(I)=DIAG(I)*W(I)
!
      DO 145 I=1,BOUND
         YR= DDOT(N,W(IYPT+CP*N+1),1,W,1)
         BETA= W(N+CP+1)*YR
         INMC=N+M+CP+1
         BETA= W(INMC)-BETA
         ISCN=ISPT+CP*N
         CALL DAXPY(N,BETA,W(ISCN+1),1,W,1)
         CP=CP+1
         IF (CP.EQ.M)CP=0
 145  CONTINUE
!
!     STORE THE NEW SEARCH DIRECTION
!     ------------------------------
!
       DO 160 I=1,N
 160   W(ISPT+POINT*N+I)= W(I)
!
!     OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION 
!     BY USING THE LINE SEARCH ROUTINE MCSRCH
!     ----------------------------------------------------
 165  NFEV=0
      STP=ONE
      IF (ITER.EQ.1) STP=STP1
      DO 170 I=1,N
 170  W(I)=G(I)
 172  CONTINUE
      CALL MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,&
                  XTOL,MAXFEV,INFO,NFEV,DIAG)
      IF (INFO .EQ. -1) THEN
        IFLAG=1
        RETURN
      ENDIF
!
!  ERROR RETURN
!
      IF (INFO .NE. 1) GO TO 190
      NFUN= NFUN + NFEV
!
!     COMPUTE THE NEW STEP AND GRADIENT CHANGE 
!     -----------------------------------------
!
      NPT=POINT*N
      DO 175 I=1,N
      W(ISPT+NPT+I)= STP*W(ISPT+NPT+I)
 175  W(IYPT+NPT+I)= G(I)-W(I)
      POINT=POINT+1
      IF (POINT.EQ.M)POINT=0
!
!     TERMINATION TEST
!     ----------------
!
      GNORM= DSQRT(DDOT(N,G,1,G,1))
      XNORM= DSQRT(DDOT(N,X,1,X,1))
      XNORM= DMAX1(1.0D0,XNORM)
!
!     
      WRITE(*,1493)ITER,NFUN,XNORM,GNORM,GNORM/XNORM,EPS,F
 1493 FORMAT(2I6,4F12.6,F15.6)
      IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE.

      IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,&
                 XNORM,GNORM,N,M,X,F,G,STP,FINISH)
      IF (FINISH) THEN
         IFLAG=0
         RETURN
      ENDIF
      GO TO 80
!
!     ------------------------------------------------------------
!     END OF MAIN ITERATION LOOP. ERROR EXITS.
!     ------------------------------------------------------------
!
 190  IFLAG=-1
      IF(LP.GT.0) WRITE(LP,200) INFO
      
      WRITE(*,1981)ITER,NFUN,&
                 XNORM,GNORM,GNORM/XNORM
 1981 FORMAT(2I5,3F15.6)
      RETURN
 195  IFLAG=-2
      IF(LP.GT.0) WRITE(LP,235) I
      RETURN
 196  IFLAG= -3
      IF(LP.GT.0) WRITE(LP,240)
!
!     FORMATS
!     -------
!
 200  FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED. SEE'&
               ' DOCUMENTATION OF ROUTINE MCSRCH',/' ERROR RETURN'&
               ' OF LINE SEARCH: INFO= ',I2,/&
               ' POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT',/,&
               ' OR INCORRECT TOLERANCES')
 235  FORMAT(/' IFLAG= -2',/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',/,&
             ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE')
 240  FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N OR M',&
             ' ARE NOT POSITIVE)')
 245  FORMAT(/'  GTOL IS LESS THAN OR EQUAL TO 1.D-04',&
             / ' IT HAS BEEN RESET TO 9.D-01')
      RETURN
      END
!
!     LAST LINE OF SUBROUTINE LBFGS
! *********************************************************************
!
      SUBROUTINE LB1(IPRINT,ITER,NFUN,XNORM,&
                           GNORM,N,M,X,F,G,STP,FINISH)
!
!     -------------------------------------------------------------
!     THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND
!     AMOUNT OF OUTPUT ARE CONTROLLED BY IPRINT.
!     -------------------------------------------------------------
!
      INTEGER IPRINT(2),ITER,NFUN,LP,MP,N,M
      DOUBLE PRECISION X(N),G(N),F,GNORM,STP,GTOL,STPMIN,STPMAX
      DOUBLE PRECISION XNORM
      LOGICAL FINISH
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
!
      IF (ITER.EQ.0)THEN
           WRITE(MP,10)
           WRITE(MP,20) N,M
           WRITE(MP,30)F,GNORM
                 IF (IPRINT(2).GE.1)THEN
                     WRITE(MP,40)
                     WRITE(MP,50) (X(I),I=1,N)
                     WRITE(MP,60)
                     WRITE(MP,50) (G(I),I=1,N)
                  ENDIF
           WRITE(MP,10)
           WRITE(MP,70)
      ELSE
          IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN
              IF (IPRINT(1).NE.0)THEN
                   IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN
                         IF(IPRINT(2).GT.1.AND.ITER.GT.1) WRITE(MP,70)
                   WRITE(MP,80)ITER,NFUN,F,XNORM,GNORM,GNORM/XNORM,STP
                   ELSE
                         RETURN
                   ENDIF
              ELSE
                   IF( IPRINT(2).GT.1.AND.FINISH) WRITE(MP,70)
                   WRITE(MP,80)ITER,NFUN,F,XNORM,GNORM,GNORM/XNORM,STP
              ENDIF
              IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN
                    IF (FINISH)THEN
                        WRITE(MP,90)
                    ELSE
                        WRITE(MP,40)
                    ENDIF
                      WRITE(MP,50)(X(I),I=1,N)
                  IF (IPRINT(2).EQ.3)THEN
                      WRITE(MP,60)
                      WRITE(MP,50)(G(I),I=1,N)
                  ENDIF
              ENDIF
            IF (FINISH) WRITE(MP,100)
      ENDIF
!
 10   FORMAT(' ************************************************')
 20   FORMAT('  N=',I15,'   NUMBER OF CORRECTIONS=',I2,&
             /,  '       INITIAL VALUES')
 30   FORMAT(' F= ',1PD10.3,'   GNORM= ',1PD10.3)
 40   FORMAT(' VECTOR X= ')
 50   FORMAT(6(2X,1PD10.3))
 60   FORMAT(' GRADIENT VECTOR G= ')
 70   FORMAT(/'   I   NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/)
! 80   FORMAT(2(I4,1X),3X,3(1PD10.3,2X))
 80   FORMAT(2I6,2X,F16.4,5F12.4)
 90   FORMAT(' FINAL POINT X= ')
 100  FORMAT(/' THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.',&
             /' IFLAG = 0')
!
      RETURN
      END
!     ***********************************************************
!
!
!   ----------------------------------------------------------
!     DATA 
!
!    MP  is an INTEGER variable with default value 41. It is used as the
!        unit number for the printing of the monitoring information
!        controlled by IPRINT.
! 
!    LP  is an INTEGER variable with default value 42. It is used as the
!        unit number for the printing of error messages. This printing
!        may be suppressed by setting LP to a non-positive value.
! 
!    GTOL is a DOUBLE PRECISION variable with default value 0.9, which
!        controls the accuracy of the line search routine MCSRCH. If the
!        function and gradient evaluations are inexpensive with respect
!        to the cost of the iteration (which is sometimes the case when
!        solving very large problems) it may be advantageous to set GTOL
!        to a small value. A typical small value is 0.1.  Restriction:
!        GTOL should be greater than 1.D-04.
! 
!    STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which
!        specify lower and uper bounds for the step in the line search.
!        Their default values are 1.D-20 and 1.D+20, respectively. These
!        values need not be modified unless the exponents are too large
!        for the machine being used, or unless the problem is extremely
!        badly scaled (in which case the exponents should be increased).
! 
!   ----------------------------------------------------------
!
      BLOCK DATA LB2
      INTEGER LP,MP
      DOUBLE PRECISION GTOL,STPMIN,STPMAX
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
!      DATA MP,LP,GTOL,STPMIN,STPMAX/6,6,1.0D-01,1.0D-20,1.0D+20/
!      DATA MP,LP,GTOL,STPMIN,STPMAX/6,6,0.5,1.0D-20,1.0D+20/
      DATA MP,LP,GTOL,STPMIN,STPMAX/41,42,9.0D-01,1.0D-20,1.0D+20/
      END
!
!
!   ----------------------------------------------------------
!
      subroutine daxpy(n,da,dx,incx,dy,incy)
!
!     constant times a vector plus a vector.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!
      double precision dx(1),dy(1),da
      integer i,incx,incy,ix,iy,m,mp1,n
!
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end
!
!
!   ----------------------------------------------------------
!
      double precision function ddot(n,dx,incx,dy,incy)
!
!     forms the dot product of two vectors.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!
      double precision dx(1),dy(1),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
!
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +&
         dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end
!    ------------------------------------------------------------------
!
!     **************************
!     LINE SEARCH ROUTINE MCSRCH
!     **************************
!
      SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL,MAXFEV,INFO,NFEV,WA)
      INTEGER N,MAXFEV,INFO,NFEV
      DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX
      DOUBLE PRECISION X(N),G(N),S(N),WA(N)
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
      SAVE
!
!                     SUBROUTINE MCSRCH
!                
!     A slight modification of the subroutine CSRCH of More' and Thuente.
!     The changes are to allow reverse communication, and do not affect
!     the performance of the routine. 
!
!     THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES
!     A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.
!
!     AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF
!     UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF
!     UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A
!     MINIMIZER OF THE MODIFIED FUNCTION
!
!          F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
!
!     IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION
!     HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE,
!     THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT
!     CONTAINS A MINIMIZER OF F(X+STP*S).
!
!     THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES
!     THE SUFFICIENT DECREASE CONDITION
!
!           F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
!
!     AND THE CURVATURE CONDITION
!
!           ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
!
!     IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION
!     IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES
!     BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH
!     CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING
!     ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY
!     SATISFIES THE SUFFICIENT DECREASE CONDITION.
!
!     THE SUBROUTINE STATEMENT IS
!
!        SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA)
!     WHERE
!
!       N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
!         OF VARIABLES.
!
!       X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
!         BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS
!         X + STP*S.
!
!       F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F
!         AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.
!
!       G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
!         GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT
!         OF F AT X + STP*S.
!
!       S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE
!         SEARCH DIRECTION.
!
!       STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN
!         INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT
!         STP CONTAINS THE FINAL ESTIMATE.
!
!       FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverse
!         communication implementation GTOL is defined in a COMMON
!         statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE
!         CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
!         SATISFIED.
!
!       XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
!         WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
!         IS AT MOST XTOL.
!
!       STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
!         SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverse
!         communication implementatin they are defined in a COMMON
!         statement).
!
!       MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
!         OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
!         MAXFEV BY THE END OF AN ITERATION.
!
!       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
!
!         INFO = 0  IMPROPER INPUT PARAMETERS.
!
!         INFO =-1  A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.
!
!         INFO = 1  THE SUFFICIENT DECREASE CONDITION AND THE
!                   DIRECTIONAL DERIVATIVE CONDITION HOLD.
!
!         INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
!                   IS AT MOST XTOL.
!
!         INFO = 3  NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
!
!         INFO = 4  THE STEP IS AT THE LOWER BOUND STPMIN.
!
!         INFO = 5  THE STEP IS AT THE UPPER BOUND STPMAX.
!
!         INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS.
!                   THERE MAY NOT BE A STEP WHICH SATISFIES THE
!                   SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
!                   TOLERANCES MAY BE TOO SMALL.
!
!       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
!         CALLS TO FCN.
!
!       WA IS A WORK ARRAY OF LENGTH N.
!
!     SUBPROGRAMS CALLED
!
!       MCSTEP
!
!       FORTRAN-SUPPLIED...ABS,MAX,MIN
!
!     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
!     JORGE J. MORE', DAVID J. THUENTE
!
!     **********
      INTEGER INFOC,J
      LOGICAL BRACKT,STAGE1
      DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,&
             FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,&
             STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO
      DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/
      IF(INFO.EQ.-1) GO TO 45
      INFOC = 1
!
!     CHECK THE INPUT PARAMETERS FOR ERRORS.
!
      IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR.&
          GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO&
          .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) RETURN
!
!     COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
!     AND CHECK THAT S IS A DESCENT DIRECTION.
!
      DGINIT = ZERO
      DO 10 J = 1, N
         DGINIT = DGINIT + G(J)*S(J)
   10    CONTINUE
      IF (DGINIT .GE. ZERO) then
         write(LP,15)
   15    FORMAT(/'  THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION')
         RETURN
         ENDIF
!
!     INITIALIZE LOCAL VARIABLES.
!
      BRACKT = .FALSE.
      STAGE1 = .TRUE.
      NFEV = 0
      FINIT = F
      DGTEST = FTOL*DGINIT
      WIDTH = STPMAX - STPMIN
      WIDTH1 = WIDTH/P5
      DO 20 J = 1, N
         WA(J) = X(J)
   20    CONTINUE
!
!     THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
!     FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
!     THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
!     FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
!     THE INTERVAL OF UNCERTAINTY.
!     THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
!     FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
!
      STX = ZERO
      FX = FINIT
      DGX = DGINIT
      STY = ZERO
      FY = FINIT
      DGY = DGINIT
!
!     START OF ITERATION.
!
   30 CONTINUE
!
!        SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
!        TO THE PRESENT INTERVAL OF UNCERTAINTY.
!
         IF (BRACKT) THEN
            STMIN = MIN(STX,STY)
            STMAX = MAX(STX,STY)
         ELSE
            STMIN = STX
            STMAX = STP + XTRAPF*(STP - STX)
            END IF
!
!        FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
!
         STP = MAX(STP,STPMIN)
         STP = MIN(STP,STPMAX)
!
!        IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
!        STP BE THE LOWEST POINT OBTAINED SO FAR.
!
         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))&
            .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0&
            .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX
!
!        EVALUATE THE FUNCTION AND GRADIENT AT STP
!        AND COMPUTE THE DIRECTIONAL DERIVATIVE.
!        We return to main program to obtain F and G.
!
         DO 40 J = 1, N
            X(J) = WA(J) + STP*S(J)
   40       CONTINUE
         INFO=-1
         RETURN
!
   45    INFO=0
         NFEV = NFEV + 1
         DG = ZERO
         DO 50 J = 1, N
            DG = DG + G(J)*S(J)
   50       CONTINUE
         FTEST1 = FINIT + STP*DGTEST
!
!        TEST FOR CONVERGENCE.
!
         IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))&
            .OR. INFOC .EQ. 0) INFO = 6
         IF (STP .EQ. STPMAX .AND.&
             F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5
         IF (STP .EQ. STPMIN .AND.&
             (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4
         IF (NFEV .GE. MAXFEV) INFO = 3
         IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2
         IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1
!
!        CHECK FOR TERMINATION.
!
         IF (INFO .NE. 0) RETURN
!
!        IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
!        FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
!
         IF (STAGE1 .AND. F .LE. FTEST1 .AND.&
             DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
!
!        A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
!        WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
!        FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
!        DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
!        OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
!
         IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN
!
!           DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
!
            FM = F - STP*DGTEST
            FXM = FX - STX*DGTEST
            FYM = FY - STY*DGTEST
            DGM = DG - DGTEST
            DGXM = DGX - DGTEST
            DGYM = DGY - DGTEST
!
!           CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
!           AND TO COMPUTE THE NEW STEP.
!
            CALL MCSTEP(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,&
                       BRACKT,STMIN,STMAX,INFOC)
!
!           RESET THE FUNCTION AND GRADIENT VALUES FOR F.
!
            FX = FXM + STX*DGTEST
            FY = FYM + STY*DGTEST
            DGX = DGXM + DGTEST
            DGY = DGYM + DGTEST
         ELSE
!
!           CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
!           AND TO COMPUTE THE NEW STEP.
!
            CALL MCSTEP(STX,FX,DGX,STY,FY,DGY,STP,F,DG,&
                       BRACKT,STMIN,STMAX,INFOC)
            END IF
!
!        FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
!        INTERVAL OF UNCERTAINTY.
!
         IF (BRACKT) THEN
            IF (ABS(STY-STX) .GE. P66*WIDTH1)&
               STP = STX + P5*(STY - STX)
            WIDTH1 = WIDTH
            WIDTH = ABS(STY-STX)
            END IF
!
!        END OF ITERATION.
!
         GO TO 30
!
!     LAST LINE OF SUBROUTINE MCSRCH.
!
      END
!  **********************************************************************
!
      SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,&
                       STPMIN,STPMAX,INFO)
      INTEGER INFO
      DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX
      LOGICAL BRACKT,BOUND
!
!     SUBROUTINE MCSTEP
!
!     THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR
!     A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR
!     A MINIMIZER OF THE FUNCTION.
!
!     THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION
!     VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS
!     ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE
!     DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A
!     MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY
!     WITH ENDPOINTS STX AND STY.
!
!     THE SUBROUTINE STATEMENT IS
!
!       SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
!                        STPMIN,STPMAX,INFO)
!
!     WHERE
!
!       STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP,
!         THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED
!         SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION
!         OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE
!         SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY.
!
!       STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP,
!         THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF
!         THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE
!         UPDATED APPROPRIATELY.
!
!       STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP,
!         THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP.
!         IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE
!         BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP.
!
!       BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER
!         HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED
!         THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER
!         IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE.
!
!       STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER
!         AND UPPER BOUNDS FOR THE STEP.
!
!       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
!         IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED
!         ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE
!         INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS.
!
!     SUBPROGRAMS CALLED
!
!       FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT
!
!     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
!     JORGE J. MORE', DAVID J. THUENTE
!
      DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA
      INFO = 0
!
!     CHECK THE INPUT PARAMETERS FOR ERRORS.
!
      IF ((BRACKT .AND. (STP .LE. MIN(STX,STY) .OR.&
          STP .GE. MAX(STX,STY))) .OR.&
          DX*(STP-STX) .GE. 0.0 .OR. STPMAX .LT. STPMIN) RETURN
!
!     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
!
      SGND = DP*(DX/ABS(DX))
!
!     FIRST CASE. A HIGHER FUNCTION VALUE.
!     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER
!     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN,
!     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN.
!
      IF (FP .GT. FX) THEN
         INFO = 1
         BOUND = .TRUE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
         IF (STP .LT. STX) GAMMA = -GAMMA
         P = (GAMMA - DX) + THETA
         Q = ((GAMMA - DX) + GAMMA) + DP
         R = P/Q
         STPC = STX + R*(STP - STX)
         STPQ = STX + ((DX/((FX-FP)/(STP-STX)+DX))/2)*(STP - STX)
         IF (ABS(STPC-STX) .LT. ABS(STPQ-STX)) THEN
            STPF = STPC
         ELSE
           STPF = STPC + (STPQ - STPC)/2
           END IF
         BRACKT = .TRUE.
!
!     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF
!     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC
!     STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP,
!     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN.
!
      ELSE IF (SGND .LT. 0.0) THEN
         INFO = 2
         BOUND = .FALSE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
         GAMMA = S*SQRT((THETA/S)**2 - (DX/S)*(DP/S))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = ((GAMMA - DP) + GAMMA) + DX
         R = P/Q
         STPC = STP + R*(STX - STP)
         STPQ = STP + (DP/(DP-DX))*(STX - STP)
         IF (ABS(STPC-STP) .GT. ABS(STPQ-STP)) THEN
            STPF = STPC
         ELSE
            STPF = STPQ
            END IF
         BRACKT = .TRUE.
!
!     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
!     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES.
!     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY
!     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC
!     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE
!     EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO
!     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP
!     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
!
      ELSE IF (ABS(DP) .LT. ABS(DX)) THEN
         INFO = 3
         BOUND = .TRUE.
         THETA = 3*(FX - FP)/(STP - STX) + DX + DP
         S = MAX(ABS(THETA),ABS(DX),ABS(DP))
!
!        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND
!        TO INFINITY IN THE DIRECTION OF THE STEP.
!
         GAMMA = S*SQRT(MAX(0.0D0,(THETA/S)**2 - (DX/S)*(DP/S)))
         IF (STP .GT. STX) GAMMA = -GAMMA
         P = (GAMMA - DP) + THETA
         Q = (GAMMA + (DX - DP)) + GAMMA
         R = P/Q
         IF (R .LT. 0.0 .AND. GAMMA .NE. 0.0) THEN
            STPC = STP + R*(STX - STP)
         ELSE IF (STP .GT. STX) THEN
            STPC = STPMAX
         ELSE
            STPC = STPMIN
            END IF
         STPQ = STP + (DP/(DP-DX))*(STX - STP)
         IF (BRACKT) THEN
            IF (ABS(STP-STPC) .LT. ABS(STP-STPQ)) THEN
               STPF = STPC
            ELSE
               STPF = STPQ
               END IF
         ELSE
            IF (ABS(STP-STPC) .GT. ABS(STP-STPQ)) THEN
               STPF = STPC
            ELSE
               STPF = STPQ
               END IF
            END IF
!
!     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE
!     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES
!     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP
!     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN.
!
      ELSE
         INFO = 4
         BOUND = .FALSE.
         IF (BRACKT) THEN
            THETA = 3*(FP - FY)/(STY - STP) + DY + DP
            S = MAX(ABS(THETA),ABS(DY),ABS(DP))
            GAMMA = S*SQRT((THETA/S)**2 - (DY/S)*(DP/S))
            IF (STP .GT. STY) GAMMA = -GAMMA
            P = (GAMMA - DP) + THETA
            Q = ((GAMMA - DP) + GAMMA) + DY
            R = P/Q
            STPC = STP + R*(STY - STP)
            STPF = STPC
         ELSE IF (STP .GT. STX) THEN
            STPF = STPMAX
         ELSE
            STPF = STPMIN
            END IF
         END IF
!
!     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT
!     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE.
!
      IF (FP .GT. FX) THEN
         STY = STP
         FY = FP
         DY = DP
      ELSE
         IF (SGND .LT. 0.0) THEN
            STY = STX
            FY = FX
            DY = DX
            END IF
         STX = STP
         FX = FP
         DX = DP
         END IF
!
!     COMPUTE THE NEW STEP AND SAFEGUARD IT.
!
      STPF = MIN(STPMAX,STPF)
      STPF = MAX(STPMIN,STPF)
      STP = STPF
      IF (BRACKT .AND. BOUND) THEN
         IF (STY .GT. STX) THEN
            STP = MIN(STX+0.66*(STY-STX),STP)
         ELSE
            STP = MAX(STX+0.66*(STY-STX),STP)
            END IF
         END IF
      RETURN
!
!     LAST LINE OF SUBROUTINE MCSTEP.
!
      END
! 
!  *************************************************************************
!
      SUBROUTINE ROLLCALL_ZIGZAG(NIIII,NFIRST,NLAST,NPUNIQUE,       &
                           NLEGISLATORS,NQTOTAL,NSMAX,              &
                           NQTOTALMAX,NDIMENSIONS,                  &
                           NP,NQ,NS,M,JPRINT,IPRINT,                &
                           XDATA,ZDATAY,ZDATAN,ALPHA,NCONGMAX,      &
                           NPMAX,NQMAX,NXDATA,NUMCONG,NUMCONGRC,    &
                           INUMRC,DYN,ZMID,ICONGRC,                 &
                           RCVOTE1,RCVOTE9,LWHERE,KWHERE,           &
                           WEIGHT,                                  &
                 JCONGUNIQUE,JD1UNIQUE,JSTATEUNIQUE,JPARTYUNIQUE,   &
                 JDISTUNIQUE,JSTAUNIQUE,JNAMEUNIQUE,XDATAUNIQUE,    &
                           ZDF,NDEVIT,XDEVIT,EXITSTATUS)
!
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION XDATA(NPMAX,NDIMENSIONS),WEIGHT(NDIMENSIONS+1),     &
                ZDATAY(NQTOTALMAX,NDIMENSIONS),NUMCONG(NCONGMAX),   &
                ZDATAN(NQTOTALMAX,NDIMENSIONS),                     &
                X((NPMAX+2*NQTOTALMAX)*NDIMENSIONS),                &
                G((NPMAX+2*NQTOTALMAX)*NDIMENSIONS),                &
                DIAG((NPMAX+2*NQTOTALMAX)*NDIMENSIONS),             &
           W(((NPMAX+2*NQTOTALMAX)*NDIMENSIONS*(2*5+1)+2*5)),       &
                NXDATA(NCONGMAX),NUMCONGRC(NCONGMAX),               &
                INUMRC(NQTOTALMAX,2),ZDF(150000,4),                 &
                ZMID(NQTOTALMAX,NSMAX),DYN(NQTOTALMAX,NSMAX),       &
                XDATAUNIQUE(NPMAX,NSMAX),XBIGLOG(NPMAX)
!
      INTEGER EXITSTATUS,IPRINT(2),KWHERE(99999,NCONGMAX),&
              ICONGRC(NQTOTALMAX),&
              JCONGUNIQUE(NPMAX),JD1UNIQUE(NPMAX),JSTATEUNIQUE(NPMAX),&
              JDISTUNIQUE(NPMAX),JPARTYUNIQUE(NPMAX),KBIGLOG(NPMAX,2)
      LOGICAL DIAGCO
      LOGICAL*1 RCVOTE1(NPMAX,NQMAX),RCVOTE9(NPMAX,NQMAX),&
                LWHERE(99999,NCONGMAX)

      CHARACTER*1 JSTAUNIQUE(NPMAX,7),JNAMEUNIQUE(NPMAX,11)
!
!     The driver for LBFGS must always declare LB2 as EXTERNAL
!
!  HOUSE DATA:
!    NQTOTAL=52968      -- TOTAL ROLL CALLS (THE UNANIMOUS ARE NOT FILTERED OUT)
!    NQ=46378           -- TOTAL SCALABLE ROLL CALLS
!    NLEGSCALED=37512   -- TOTAL SCALABLE LEGISLATORS -- ALSO TOTAL IN RCVOTE1(:) AND RCVOTE9(:)
!    NPUNIQUE=10724     -- TOTAL UNIQUE LEGISLATORS IN RCVOTE1(:), RCVOTE9(:), LWHERE(:), KWHERE(:)
!    NLEGISLATORS=38331 -- TOTAL LEGISLATORS IN THE STACKED ROLL CALL FILE
!    
!
      EXTERNAL LB2
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
!
      INTEGER, ALLOCATABLE :: KK1(:)
      INTEGER, ALLOCATABLE :: KK2(:)
      DOUBLE PRECISION, ALLOCATABLE :: XSAVE(:)
      DOUBLE PRECISION, ALLOCATABLE :: XLOGJ(:)
      DOUBLE PRECISION, ALLOCATABLE :: OLDZMID(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: OLDDL(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XNEWZMID(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XNEWDL(:,:)
!
      ALLOCATE(KK1(NDIMENSIONS))
      ALLOCATE(KK2(NDIMENSIONS))
      ALLOCATE(XSAVE((NPMAX+2*NQTOTALMAX)*NDIMENSIONS))
      ALLOCATE(XLOGJ(NQTOTAL+1))
      ALLOCATE(OLDZMID(NQTOTALMAX,NDIMENSIONS))
      ALLOCATE(OLDDL(NQTOTALMAX,NDIMENSIONS))
      ALLOCATE(XNEWZMID(NQTOTALMAX,NDIMENSIONS))
      ALLOCATE(XNEWDL(NQTOTALMAX,NDIMENSIONS))
!
  189 FORMAT(I6,2F12.6)
  201 FORMAT(I3,I5,F12.6,20F7.3)
!  175 FORMAT(I3,I5,4F7.3)
  302 FORMAT(' PRELIMINARY LOG-LIKELIHOOD O PHASE',2I12,F20.6)
  303 FORMAT(' LOG-LIKELIHOOD Z-PHASE',2I5,6I12,F20.6)
  304 FORMAT(' NS, NQ, NP, N',4I7)
  305 FORMAT(' NS, NQ, NP, N, BETA',4I7,F10.4)
  306 FORMAT(' R-SQUARE CURRENT VS. PREVIOUS ROLL CALL PHASE',I12,F9.4)
  307 FORMAT(' R-SQUARE CURRENT VS. PREVIOUS MIDPOINT',I4,I12,F9.4)
  308 FORMAT(' R-SQUARE CURRENT VS. PREVIOUS SPREAD  ',I4,I12,F9.4)
 1198 FORMAT(4I7,7X,12F8.3)
 1199 FORMAT(5I7,4F8.3)
!
      N=2*NS*NQ
      TWOPI=1.0/SQRT(2.0*3.1415926536)
!
      IF(JPRINT.EQ.1)WRITE(*,304)NS,NQ,NP,N
!
!  CALCULATE PRELIMINARY LOG-LIKELIHOOD
!
      BETA=WEIGHT(NS+1)
      IF(JPRINT.EQ.1)WRITE(*,305)NS,NQ,NP,N,BETA
!
!
!  CALCULATE THE INITIAL LOG-LIKELIHOOD
!
      CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,NLEGISLATORS,&
                 NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
!
      IF(JPRINT.EQ.1)WRITE(*,302)NIIII,KTOT,XPLOG
      IF(JPRINT.EQ.1)WRITE(21,302)NIIII,KTOT,XPLOG
!
!
      SUMSQUARED=0.0
      KCHOICE=0
      KTOTP=0
      KTOTQ=0
      DO 111 II=NFIRST,NLAST
      NPC=NXDATA(II)
      NQC=NUMCONGRC(II)
      DO 1 I=1,NPC
      KK=0
      DO 2 J=1,NQC
      IF(INUMRC(J+KTOTQ,2).EQ.0)THEN
         DISTYEA=0.0
         DISTNAY=0.0
         KK=KK+1
         DO 3 K=1,NS
         DISTYEA=DISTYEA+(WEIGHT(K)*(XDATA(I+KTOTP,K)-&
                          ZDATAY(J+KTOTQ,K)))**2
         DISTNAY=DISTNAY+(WEIGHT(K)*(XDATA(I+KTOTP,K)-&
                          ZDATAN(J+KTOTQ,K)))**2
  3      CONTINUE
         PSIYEA=DEXP((-1.0)*(DISTYEA/1.0)) - DEXP((-1.0)*(DISTNAY/1.0))
         PSINAY=DEXP((-1.0)*(DISTNAY/1.0)) - DEXP((-1.0)*(DISTYEA/1.0))
         SMALLX = (BETA*PSIYEA)/SQRT(2.0)
         VALUE = DERF(SMALLX)       
         PHIYEA = (VALUE/2.0) + 0.5
         IF(PHIYEA.LE..000001)PHIYEA=.00001
         SMALLX = (BETA*PSINAY)/SQRT(2.0)
         VALUE = DERF(SMALLX)       
         PHINAY = (VALUE/2.0) + 0.5
         IF(PHINAY.LE..000001)PHINAY=.00001
!
!  IF NOT MISSING DATA
!
         IF(RCVOTE9(I+KTOTP,J).EQV..FALSE.)THEN
            KCHOICE=KCHOICE+1
!
!  IF YES
!
            IF(RCVOTE1(I+KTOTP,J).EQV..TRUE.)THEN
               SUMSQUARED=SUMSQUARED+LOG(PHIYEA)
            ENDIF
!
!  IF NO
!
            IF(RCVOTE1(I+KTOTP,J).EQV..FALSE.)THEN
               SUMSQUARED=SUMSQUARED+LOG(PHINAY)
            ENDIF
         ENDIF
      ENDIF
  2   CONTINUE
  1   CONTINUE
      KTOTP=KTOTP+NPC
      KTOTQ=KTOTQ+NQC
  111 CONTINUE
!
      IF(JPRINT.EQ.1)WRITE(*,302)NIIII,KCHOICE,SUMSQUARED
      IF(JPRINT.EQ.1)WRITE(21,302)NIIII,KCHOICE,SUMSQUARED
!
!
!  TRANSFER STARTING COORDINATES INTO X(.)
!
      KKK=0
      KTOTQ=0
      DO 113 II=NFIRST,NLAST
      NQC=NUMCONGRC(II)
      DO 91 I=1,NQC
      IF(INUMRC(I+KTOTQ,2).EQ.0)THEN
         DO 92 K=1,NS
         KKK=KKK+1
         X(KKK)=ZDATAY(I+KTOTQ,K)
         XSAVE(KKK)=X(KKK)
!         WRITE(99,1199)II,I,K,KKK,KTOTQ,X(KKK),ZDATAY(I+KTOTQ,K),&
!                    DYN(I+KTOTQ,K),ZMID(I+KTOTQ,K)
  92     CONTINUE
         DO 93 K=1,NS
         KKK=KKK+1
         X(KKK)=ZDATAN(I+KTOTQ,K)
         XSAVE(KKK)=X(KKK)
!         WRITE(99,1199)II,I,K,KKK,KTOTQ,X(KKK),ZDATAN(I+KTOTQ,K),&
!                    DYN(I+KTOTQ,K),ZMID(I+KTOTQ,K)
  93     CONTINUE
      ENDIF
  91  CONTINUE
      KTOTQ=KTOTQ+NQC
  113 CONTINUE
!
!
!     We do not wish to provide the diagonal matrices Hk0, and 
!     therefore set DIAGCO to FALSE.
!
      DIAGCO= .FALSE.
      EPS= 1.0
!
!  XTOL IS USED IN THE LINE SEARCH ROUTINE
!
!      XTOL= 0.0001
      XTOL= 0.01
      TWOPI=1.0/SQRT(2.0*3.1415926536)
      ICALL=0
      IFLAG=0
!
      N=2*NS*NQ
!
!  BEGIN LBFGS LOOP
!
!  X(.) IS THE VECTOR OF PARAMETERS
!  G(.) IS THE VECTOR OF DERIVATIVES
!  F IS THE VALUE OF THE OBJECTIVE FUNCTION GIVEN X(.)
!
 20   CONTINUE
      F= 0.0
      DO 59 I=1,N
      G(I)=0.0
  59  CONTINUE
!
!  TRANSFER X(.) PARAMETER VECTOR BACK INTO COORDINATE MATRICES
!
      
      KKK=0
      KTOTQ=0
      DO 114 II=NFIRST,NLAST
      NQC=NUMCONGRC(II)
      DO 11 I=1,NQC
      IF(INUMRC(I+KTOTQ,2).EQ.0)THEN
         DO 12 K=1,NS
         KKK=KKK+1
         ZDATAY(I+KTOTQ,K)=X(KKK)
  12     CONTINUE
         DO 13 K=1,NS
         KKK=KKK+1
         ZDATAN(I+KTOTQ,K)=X(KKK)
  13     CONTINUE
      ENDIF
!      WRITE(99,1198)II,I,KKK,KTOTQ,X(KKK),&
!              (ZDATAY(I+KTOTQ,K),K=1,NS),(ZDATAN(I+KTOTQ,K),K=1,NS),&
!              (DYN(I+KTOTQ,K),(ZMID(I+KTOTQ,K)),K=1,NS)
  11  CONTINUE
      KTOTQ=KTOTQ+NQC
  114 CONTINUE
!
!  CALCULATE THE OBJECTIVE FUNCTION, F
!
      KCHOICE=0
      KLASSIFY=0
      KTOTP=0
      KTOTQ=0
      KKK=0
      JJJJ=0
      DO 112 II=NFIRST,NLAST
      NPC=NXDATA(II)
      NQC=NUMCONGRC(II)
      DO 17 J=1,NQC
      JJJJ=JJJJ+1
      XLOGJ(JJJJ)=0.0
      IF(INUMRC(J+KTOTQ,2).EQ.0)THEN
!
         DO 1591 K=1,NS
         KKK=KKK+1
         KK1(K)=KKK
 1591    CONTINUE
         DO 1592 K=1,NS
         KKK=KKK+1
         KK2(K)=KKK
 1592    CONTINUE
!        
         DO 16 I=1,NPC
         DISTYEA=0.0
         DISTNAY=0.0
         KK=KK+1
         DO 18 K=1,NS
         DISTYEA=DISTYEA+(WEIGHT(K)*(XDATA(I+KTOTP,K)-&
                          ZDATAY(J+KTOTQ,K)))**2
         DISTNAY=DISTNAY+(WEIGHT(K)*(XDATA(I+KTOTP,K)-&
                          ZDATAN(J+KTOTQ,K)))**2
  18     CONTINUE
         KK=KK+1
         PSIYEA=DEXP((-1.0)*(DISTYEA/1.0)) - DEXP((-1.0)*(DISTNAY/1.0))
         PSINAY=DEXP((-1.0)*(DISTNAY/1.0)) - DEXP((-1.0)*(DISTYEA/1.0))
         SMALLX = (BETA*PSIYEA)/SQRT(2.0)
         VALUE = DERF(SMALLX)       
         PHIYEA = (VALUE/2.0) + 0.5
         IF(PHIYEA.LE..000001)PHIYEA=.00001
         SMALLX = (BETA*PSINAY)/SQRT(2.0)
         VALUE = DERF(SMALLX)       
         PHINAY = (VALUE/2.0) + 0.5
         IF(PHINAY.LE..000001)PHINAY=.00001
         PSIYEA2 = TWOPI*DEXP((-.5)*BETA*BETA*PSIYEA*PSIYEA)
         PSINAY2 = TWOPI*DEXP((-.5)*BETA*BETA*PSINAY*PSINAY)
         XINVMILLSY = PSIYEA2/PHIYEA
         XINVMILLSN = PSINAY2/PHINAY
!
!  IF NOT MISSING DATA
!
         IF(RCVOTE9(I+KTOTP,J).EQV..FALSE.)THEN
            KCHOICE=KCHOICE+1
!
!  IF YES
!
            IF(RCVOTE1(I+KTOTP,J).EQV..TRUE.)THEN
               F=F+LOG(PHIYEA)
               XLOGJ(JJJJ)=XLOGJ(JJJJ)+LOG(PHIYEA)
               IF(DISTYEA.LT.DISTNAY)KLASSIFY=KLASSIFY+1
            ENDIF
!
!  IF NO
!
            IF(RCVOTE1(I+KTOTP,J).EQV..FALSE.)THEN
               F=F+LOG(PHINAY)
               XLOGJ(JJJJ)=XLOGJ(JJJJ)+LOG(PHINAY)
               IF(DISTYEA.GT.DISTNAY)KLASSIFY=KLASSIFY+1
            ENDIF
!
!  CALCULATE THE VECTOR OF PARTIAL DERIVATIVES, G(.)
!
! DERIVATIVES FOR ROLL CALLS
! 
!         DO 23 K=1,NS
!
!  IF YES
!
            IF(RCVOTE1(I+KTOTP,J).EQV..TRUE.)THEN
               DO 24 K=1,NS
               G(KK1(K))=G(KK1(K))-XINVMILLSY*&
                 ((DEXP((-.5)*WEIGHT(K)*WEIGHT(K)*DISTYEA))*&
                 (XDATA(I+KTOTP,K)-ZDATAY(J+KTOTQ,K)))
  24           CONTINUE
               DO 25 K=1,NS
               G(KK2(K))=G(KK2(K))&
                 -XINVMILLSY*&
                 ((-1.0)*(DEXP((-.5)*WEIGHT(K)*WEIGHT(K)*DISTNAY))*&
                 (XDATA(I+KTOTP,K)-ZDATAN(J+KTOTQ,K)))
  25           CONTINUE               
            ENDIF
!
!  IF NO
!
            IF(RCVOTE1(I+KTOTP,J).EQV..FALSE.)THEN
               DO 26 K=1,NS
               G(KK1(K))=G(KK1(K))-XINVMILLSN*&
                 ((-1.0)*(DEXP((-.5)*WEIGHT(K)*WEIGHT(K)*DISTYEA))*&
                 (XDATA(I+KTOTP,K)-ZDATAY(J+KTOTQ,K)))
  26           CONTINUE
               DO 27 K=1,NS
               G(KK2(K))=G(KK2(K))-&
                  XINVMILLSN*&
                 ((DEXP((-.5)*WEIGHT(K)*WEIGHT(K)*DISTNAY))*&
                 (XDATA(I+KTOTP,K)-ZDATAN(J+KTOTQ,K)))
  27           CONTINUE
            ENDIF
!  23     CONTINUE
         ENDIF
  16     CONTINUE
      ENDIF
  17  CONTINUE
      KTOTP=KTOTP+NPC
      KTOTQ=KTOTQ+NQC
  112 CONTINUE
!
      F=-F
!
!      DO 188 I=1,N
!      WRITE(99,189)I,X(I),G(I)
!  188 CONTINUE
      IF(JPRINT.EQ.1)WRITE(21,303)NIIII,ICALL,KKK,N,KLASSIFY,KCHOICE,&
                                  KTOTP,KTOTQ,F
      WRITE(*,303)NIIII,ICALL,KKK,N,KLASSIFY,KCHOICE,KTOTP,KTOTQ,F
!
!
!  N = # OF PARAMETERS
!  M = 5 FOR LINE SEARCH
!  X(N) IS THE PARAMETER VECTOR
!  F = FUNCTION VALUE
!  G(N) IS THE VECTOR OF DERIVATIVES
!  DIAGCO = .FALSE.
!  DIAG(N) = DIAGONAL OF HESSION IF DIAGCO = .TRUE.
!  IPRINT(2) = WRITES OUTPUT
!  EPS = 0.01
!  XTOL = 0.0001
!  W(N*(2*M + 1) + 2*M)
!  IFLAG = 0, NO ERRORS
!
      CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
      IF(IFLAG.LE.0) GO TO 50
      ICALL=ICALL + 1
!     We allow at most 2000 evaluations of F and G
      IF(ICALL.GT.100) GO TO 50
      GO TO 20
  50  CONTINUE
!
!  TRANSFER X(.) PARAMETER VECTOR BACK INTO COORDINATE MATRICES
!
!    
      KKK=0
      KK=0
      DO 61 I=1,NQ
      KK=KK+1
      DO 62 K=1,NS
      KKK=KKK+1
      ZDATAY(I,K)=X(KKK)
  62  CONTINUE
      KK=KK+1
      DO 63 K=1,NS
      KKK=KKK+1
      ZDATAN(I,K)=X(KKK)
  63  CONTINUE
  61  CONTINUE
!
!  TRANSFER OUTCOME COORDINATES BACK INTO MIDPOINTS AND SPREADS
!
      KK=0
      DO 85 J=1,NQTOTAL
      IF(INUMRC(J,2).EQ.1)THEN
         DO 887 K=1,NS
         ZMID(J,K)=0.0
         DYN(J,K)=0.0
  887    CONTINUE
      ENDIF
      IF(INUMRC(J,2).EQ.0)THEN
         KK=KK+1
         SUM=0.0
         DO 89 K=1,NS
         OLDZMID(KK,K)=ZMID(J,K)
         OLDDL(KK,K)=DYN(J,K)
         XNEWZMID(KK,K)=(ZDATAY(KK,K)+ZDATAN(KK,K))/2.0
         XNEWDL(KK,K)=(ZDATAY(KK,K)-ZDATAN(KK,K))/2.0
         ZMID(J,K)=(ZDATAY(KK,K)+ZDATAN(KK,K))/2.0
         SUM=SUM+ZMID(J,K)**2
         DYN(J,K)=(ZDATAN(KK,K)-ZDATAY(KK,K))/2.0
!         ZDATAY(KK,K)=ZMID(KK,K)-DYN(KK,K)
!         ZDATAN(KK,K)=ZMID(KK,K)+DYN(KK,K)
  89     CONTINUE
!
!  DO CHECK ON ROLL CALL MIDPOINT -- IF IT IS OUTSIDE THE
!   HYPERSPHERE PULL IT BACK INSIDE
!
         IF(SUM.GT.1.0)THEN
            DO 55 K=1,NS
            ZMID(J,K)=(ZMID(J,K)/SQRT(SUM))
            XNEWZMID(KK,K)=ZMID(J,K)
  55        CONTINUE
         ENDIF
!
      ENDIF
  85  CONTINUE
!
!  CORRELATE CURRENT ESTIMATE WITH LAST ESTIMATE
!
      CALL CORR2(KKK,NPMAX,NQTOTALMAX,NDIMENSIONS,X,XSAVE,R)
      WRITE(42,306)KKK,R
      WRITE(*,306)KKK,R
!
!  DO MIDPOINTS
!
      DO 96 K=1,NS
      DO 97 I=1,KK
      XSAVE(I)=OLDZMID(I,K)
      X(I)=XNEWZMID(I,K)
   97 CONTINUE
      CALL CORR2(KK,NPMAX,NQTOTALMAX,NDIMENSIONS,X,XSAVE,R)
      WRITE(42,307)K,KK,R
      WRITE(*,307)K,KK,R
   96 CONTINUE
!
!  DO SPREADS
!
      DO 94 K=1,NS
      DO 95 I=1,KK
      XSAVE(I)=OLDDL(I,K)
      X(I)=XNEWDL(I,K)
   95 CONTINUE
      CALL CORR2(KK,NPMAX,NQTOTALMAX,NDIMENSIONS,X,XSAVE,R)
      WRITE(42,308)K,KK,R
      WRITE(*,308)K,KK,R
   94 CONTINUE
!
!
!  TRANSFER OUTCOME COORDINATES BACK INTO OVERALL MATRIX WITH THE MISSING VOTES
!
      KK=0
      DO 885 I=1,NQTOTAL
      DO 899 K=1,NS
      ZDATAY(I,K)=ZMID(I,K)-DYN(I,K)
      ZDATAN(I,K)=ZMID(I,K)+DYN(I,K)
  899 CONTINUE
      IF(NIIII.EQ.5)THEN
      WRITE(49,201)ICONGRC(I),INUMRC(I,1),XLOGJ(I),&
                  ((ZDATAY(I,K)*WEIGHT(K)),K=1,NS),&
                  ((ZDATAN(I,K)*WEIGHT(K)),K=1,NS),&
                  ((DYN(I,K)*WEIGHT(K)),(ZMID(I,K)*WEIGHT(K)),K=1,NS)
      ENDIF
!      WRITE(*,201)ICONGRC(I),INUMRC(I,1),XLOGJ(I),&
!                           (ZDATAY(I,K),K=1,NS),&
!                           (ZDATAN(I,K),K=1,NS),&
!                      (DYN(I,K),ZMID(I,K),K=1,NS)
  885 CONTINUE
!
      DEALLOCATE(KK1)
      DEALLOCATE(KK2)
      DEALLOCATE(XSAVE)
      DEALLOCATE(XLOGJ)
      DEALLOCATE(OLDZMID)
      DEALLOCATE(XNEWZMID)
      DEALLOCATE(OLDDL)
      DEALLOCATE(XNEWDL)
!
      RETURN
      END
!
!  **********************************************************************
!    SUBROUTINE CORR2---CALLED BY ROLLCALL_ZIGZAG.  COMPUTES CORRELATION 
!     FOR CURRENT ESTIMATED ROLL CALL PARAMETERS WITH INPUT ROLL CALL 
!     PARAMETERS
!  **********************************************************************
!
      SUBROUTINE CORR2(KKK,NPMAX,NQTOTALMAX,NDIMENSIONS,X,XSAVE,R)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      
      DIMENSION X((NPMAX+2*NQTOTALMAX)*NDIMENSIONS),                &
                XSAVE((NPMAX+2*NQTOTALMAX)*NDIMENSIONS)                
!
      ASUM=0.0
      BSUM=0.0
      CSUM=0.0
      DSUM=0.0
      ESUM=0.0
!
      DO 1 I=1,KKK
      ASUM=ASUM+X(I)
      BSUM=BSUM+XSAVE(I)
      CSUM=CSUM+X(I)*X(I)
      DSUM=DSUM+XSAVE(I)*XSAVE(I)
      ESUM=ESUM+X(I)*XSAVE(I)
  1   CONTINUE
!
      AA=FLOAT(KKK)*CSUM-ASUM*ASUM
      BB=FLOAT(KKK)*DSUM-BSUM*BSUM
      CC=FLOAT(KKK)*ESUM-ASUM*BSUM
      R=(CC*CC)/(AA*BB)
!
      RETURN
      END
! 
!  *************************************************************************
!
      SUBROUTINE X_ZIGZAG(NIIII,KTOT,KLASS,GMP,NFIRST,NLAST,NPUNIQUE,&
                           NLEGSCALED,NLEGISLATORS,NQTOTAL,NSMAX,   &
                           NQTOTALMAX,NDIMENSIONS,                  &
                           NP,NQ,NS,M,JPRINT,IPRINT,                &
                           XDATA,ZDATAY,ZDATAN,ALPHA,NCONGMAX,      &
                           NPMAX,NQMAX,NXDATA,NUMCONG,NUMCONGRC,    &
                           INUMRC,DYN,ZMID,ICONGRC,                 &
                           RCVOTE1,RCVOTE9,LWHERE,KWHERE,           &
                           ICONG,ID1,ISTATE,IDIST,KSTA,IPARTY,      &
                           JCONG,JD1,JSTATE,JDIST,JSTA,JPARTY,      &
                           JNAME,LNAME,WEIGHT,KBIGLOG,XBIGLOG,      &
                 JCONGUNIQUE,JD1UNIQUE,JSTATEUNIQUE,JPARTYUNIQUE,   &
                 JDISTUNIQUE,JSTAUNIQUE,JNAMEUNIQUE,XDATAUNIQUE,    &
                           ZDF,NDEVIT,XDEVIT,EXITSTATUS)
!
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION XDATA(NPMAX,NDIMENSIONS),WEIGHT(NDIMENSIONS+1),     &
                ZDATAY(NQTOTALMAX,NDIMENSIONS),NUMCONG(NCONGMAX),   &
                ZDATAN(NQTOTALMAX,NDIMENSIONS),                     &
                X((NPMAX+2*NQTOTALMAX)*NDIMENSIONS),                &
                G((NPMAX+2*NQTOTALMAX)*NDIMENSIONS),                &
                DIAG((NPMAX+2*NQTOTALMAX)*NDIMENSIONS),             &
           W(((NPMAX+2*NQTOTALMAX)*NDIMENSIONS*(2*5+1)+2*5)),       &
                NXDATA(NCONGMAX),NUMCONGRC(NCONGMAX),               &
                INUMRC(NQTOTALMAX,2),ZDF(150000,4),                 &
                ZMID(NQTOTALMAX,NSMAX),DYN(NQTOTALMAX,NSMAX),       &
                XDATAUNIQUE(NPMAX,NSMAX),XBIGLOG(NPMAX)
!
      INTEGER EXITSTATUS,IPRINT(2),KWHERE(99999,NCONGMAX),&
              ICONGRC(NQTOTALMAX),ICONG(NPMAX),ID1(NPMAX),&
              ISTATE(NPMAX),IDIST(NPMAX),IPARTY(NPMAX),&
              JCONG(NPMAX),JD1(NPMAX),&
              JSTATE(NPMAX),JDIST(NPMAX),JPARTY(NPMAX),&
              JCONGUNIQUE(NPMAX),JD1UNIQUE(NPMAX),JSTATEUNIQUE(NPMAX),&
              JDISTUNIQUE(NPMAX),JPARTYUNIQUE(NPMAX),KBIGLOG(NPMAX,2)
      CHARACTER*1 KSTA(NPMAX,7),LNAME(NPMAX,11),&
                  JSTA(NPMAX,7),JNAME(NPMAX,11),&
                  JSTAUNIQUE(NPMAX,7),JNAMEUNIQUE(NPMAX,11)
!
      LOGICAL DIAGCO
      LOGICAL*1 RCVOTE1(NPMAX,NQMAX),RCVOTE9(NPMAX,NQMAX),&
                LWHERE(99999,NCONGMAX)
!
!     The driver for LBFGS must always declare LB2 as EXTERNAL
!
!  HOUSE DATA:
!    NQTOTAL=52968      -- TOTAL ROLL CALLS (THE UNANIMOUS ARE NOT FILTERED OUT)
!    NQ=46378           -- TOTAL SCALABLE ROLL CALLS
!    NLEGSCALED=37512   -- TOTAL SCALABLE LEGISLATORS -- ALSO TOTAL IN RCVOTE1(:) AND RCVOTE9(:)
!    NPUNIQUE=10724     -- TOTAL UNIQUE LEGISLATORS IN RCVOTE1(:), RCVOTE9(:), LWHERE(:), KWHERE(:)
!    NLEGISLATORS=38331 -- TOTAL LEGISLATORS IN THE STACKED ROLL CALL FILE
!    
!
      EXTERNAL LB2
      COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
!
      INTEGER, ALLOCATABLE :: KLOCATE(:)
      INTEGER, ALLOCATABLE :: LLOCATE(:)
      DOUBLE PRECISION, ALLOCATABLE :: XSAVE(:)
      DOUBLE PRECISION, ALLOCATABLE :: XSAVE2(:)
!
      ALLOCATE(XSAVE((NPMAX+2*NQTOTALMAX)*NDIMENSIONS))
      ALLOCATE(XSAVE2((NPMAX+2*NQTOTALMAX)*NDIMENSIONS))
!
      ALLOCATE(KLOCATE(NCONGMAX))
      ALLOCATE(LLOCATE(NCONGMAX))
!
  200 FORMAT(2I6,I4,I6,I3,I2,1X,7A1,1X,I4,1X,11A1,3800I1)
  201 FORMAT(I4,I6,I3,I2,1X,7A1,1X,I4,1X,11A1,2F7.3,1X,F12.5,2I5,F7.3)
  202 FORMAT(1X,F12.5,2I5,F7.3)
  302 FORMAT(' PRELIMINARY LOG-LIKELIHOOD X PHASE',2I10,F20.6,I10)
  303 FORMAT(' LOG-LIKELIHOOD X',4I10,F20.6)
  304 FORMAT(' X_ZIGZAG -- NS, NQ, NP, N',4I10)
  305 FORMAT(' NS, NQ, NP, N, BETA',4I7,F10.4)
  306 FORMAT(' LOG-LIKELIHOOD X PHASE',2I5,3I10,F20.6,I10)
  307 FORMAT(' R-SQUARE CURRENT VS. PREVIOUS X PHASE (UNIQUE)',&
                I12,F9.4)
  308 FORMAT(' R-SQUARE CURRENT VS. PREVIOUS X PHASE (ALL)   ',&
                I12,F9.4)
!
      N=NS*NPUNIQUE
      TWOPI=1.0/SQRT(2.0*3.1415926536)
!
      IF(JPRINT.EQ.1)WRITE(21,304)NS,NQ,NP,N
      IF(JPRINT.EQ.1)WRITE(*,304)NS,NQ,NP,N
!
!  CALCULATE PRELIMINARY LOG-LIKELIHOOD
!
      BETA=WEIGHT(NS+1)
      IF(JPRINT.EQ.1)WRITE(*,305)NS,NQ,NP,N,BETA
!
      KK=0
      DO 649 I=1,NLEGSCALED
      DO 648 K=1,NS
      KK=KK+1
      XSAVE2(KK)=XDATA(I,K)
  648 CONTINUE
  649 CONTINUE
!
!  CALCULATE THE INITIAL LOG-LIKELIHOOD
!
      CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,NLEGISLATORS,&
                 NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
!
      IF(JPRINT.EQ.1)WRITE(*,302)NIIII,KTOT,XPLOG
      IF(JPRINT.EQ.1)WRITE(21,302)NIIII,KTOT,XPLOG
!
!
      SUMSQUARED=0.0
      KCHOICE=0
      KTOTP=0
      KTOTQ=0
      DO 111 II=NFIRST,NLAST
      NPC=NXDATA(II)
      NQC=NUMCONGRC(II)
      DO 1 I=1,NPC
      DO 2 J=1,NQC
      IF(INUMRC(J+KTOTQ,2).EQ.0)THEN
         DISTYEA=0.0
         DISTNAY=0.0
         DO 3 K=1,NS
         DISTYEA=DISTYEA+(WEIGHT(K)*(XDATA(I+KTOTP,K)-&
                          ZDATAY(J+KTOTQ,K)))**2
         DISTNAY=DISTNAY+(WEIGHT(K)*(XDATA(I+KTOTP,K)-&
                          ZDATAN(J+KTOTQ,K)))**2
  3      CONTINUE
         PSIYEA=DEXP((-1.0)*(DISTYEA/1.0)) - DEXP((-1.0)*(DISTNAY/1.0))
         PSINAY=DEXP((-1.0)*(DISTNAY/1.0)) - DEXP((-1.0)*(DISTYEA/1.0))
         SMALLX = (BETA*PSIYEA)/SQRT(2.0)
         VALUE = DERF(SMALLX)       
         PHIYEA = (VALUE/2.0) + 0.5
         IF(PHIYEA.LE..000001)PHIYEA=.00001
         SMALLX = (BETA*PSINAY)/SQRT(2.0)
         VALUE = DERF(SMALLX)       
         PHINAY = (VALUE/2.0) + 0.5
         IF(PHINAY.LE..000001)PHINAY=.00001
!
!  IF NOT MISSING DATA
!
         IF(RCVOTE9(I+KTOTP,J).EQV..FALSE.)THEN
            KCHOICE=KCHOICE+1
!
!  IF YES
!
            IF(RCVOTE1(I+KTOTP,J).EQV..TRUE.)THEN
               SUMSQUARED=SUMSQUARED+LOG(PHIYEA)
            ENDIF
!
!  IF NO
!
            IF(RCVOTE1(I+KTOTP,J).EQV..FALSE.)THEN
               SUMSQUARED=SUMSQUARED+LOG(PHINAY)
            ENDIF
         ENDIF
      ENDIF
  2   CONTINUE
  1   CONTINUE
      KTOTP=KTOTP+NPC
      KTOTQ=KTOTQ+NQC
  111 CONTINUE
!
      IF(JPRINT.EQ.1)WRITE(*,302)NIIII,KCHOICE,SUMSQUARED
      IF(JPRINT.EQ.1)WRITE(21,302)NIIII,KCHOICE,SUMSQUARED
!
!  TRANSFER STARTING COORDINATES INTO X(.)
!
      KKK=0
      DO 613 I=1,NPUNIQUE
      DO 614 K=1,NS
      KKK=KKK+1
      X(KKK)=XDATAUNIQUE(I,K)
      XSAVE(KKK)=X(KKK)
  614 CONTINUE
  613 CONTINUE
!
!     We do not wish to provide the diagonal matrices Hk0, and 
!     therefore set DIAGCO to FALSE.
!
      DIAGCO= .FALSE.
      EPS= 2.0
!      XTOL= 0.0001
      XTOL= 0.01
      TWOPI=1.0/SQRT(2.0*3.1415926536)
      ICALL=0
      IFLAG=0
!
      N=NS*NPUNIQUE
!
!  BEGIN LBFGS LOOP
!
!  X(.) IS THE VECTOR OF PARAMETERS
!  G(.) IS THE VECTOR OF DERIVATIVES
!  F IS THE VALUE OF THE OBJECTIVE FUNCTION GIVEN X(.)
!
 20   CONTINUE
      F= 0.0
      DO 59 I=1,N
      G(I)=0.0
  59  CONTINUE
!
!  TRANSFER X(.) PARAMETER VECTOR BACK INTO COORDINATE MATRICES
!
      KKK=0
      DO 615 I=1,NPUNIQUE
      DO 616 K=1,NS
      KKK=KKK+1
      XDATAUNIQUE(I,K)=X(KKK)
  616 CONTINUE
!      WRITE(48,201)JCONGUNIQUE(I),JD1UNIQUE(I),&
!                  JSTATEUNIQUE(I),JDISTUNIQUE(I),&
!                  (JSTAUNIQUE(I,JJ),JJ=1,7),JPARTYUNIQUE(I),&
!                  (JNAMEUNIQUE(I,JJ),JJ=1,11),&
!                  (XDATAUNIQUE(I,K),K=1,NDIMENSIONS)
  615 CONTINUE
      
!
!  CALCULATE CONSTANT MODEL COORDINATES
!
!
      KCHOICE=0
      KORRECT=0
      F=0.0
      KKKK=0
      DO 9999 IIII=1,99999
      DO 4 I=1,NCONGMAX
      KLOCATE(I)=0
      LLOCATE(I)=0
  4   CONTINUE
      KK=0
!
      DO 449 JJ=NFIRST,NLAST
      IF(LWHERE(IIII,JJ).EQV..TRUE.)THEN
         KK=KK+1
         KLOCATE(KK)=KWHERE(IIII,JJ)
         LLOCATE(KK)=JJ
      ENDIF
  449 CONTINUE
      IF(KK.EQ.0)GO TO 9999
      KKKK=KKKK+1
!
!  LEGISLATOR IN KK CONGRESSES
!
      DO 8888 JJJJ=1,KK
      KTOTQ=0
      IF(LLOCATE(KK).GT.1)THEN
         DO 34 II=1,LLOCATE(JJJJ)-1
         KTOTQ=KTOTQ+NUMCONGRC(II)
  34     CONTINUE
      ENDIF
      NQC=NUMCONGRC(LLOCATE(JJJJ))
!      WRITE(43,1844)IIII,JJJJ,KTOTQ,NQC
 1844 FORMAT(2I6,2I8)
      DO 35 J=1,NQC
         IF(INUMRC(J+KTOTQ,2).EQ.0)THEN
            DISTYEA=0.0
            DISTNAY=0.0
!
!  CONSTANT MODEL SO USE THE FIRST APPEARANCE OF THE LEGISLATOR
!
            DO 33 K=1,NS
            DISTYEA=DISTYEA+(WEIGHT(K)*(XDATAUNIQUE(KKKK,K)-&
                          ZDATAY(J+KTOTQ,K)))**2
            DISTNAY=DISTNAY+(WEIGHT(K)*(XDATAUNIQUE(KKKK,K)-&
                          ZDATAN(J+KTOTQ,K)))**2
  33        CONTINUE
            PSIYEA=DEXP((-1.0)*(DISTYEA/1.0)) - &
                   DEXP((-1.0)*(DISTNAY/1.0))
            PSINAY=DEXP((-1.0)*(DISTNAY/1.0)) - &
                   DEXP((-1.0)*(DISTYEA/1.0))
            SMALLX = (BETA*PSIYEA)/SQRT(2.0)
            VALUE = DERF(SMALLX)       
            PHIYEA = (VALUE/2.0) + 0.5
            IF(PHIYEA.LE..000001)PHIYEA=.00001
            SMALLX = (BETA*PSINAY)/SQRT(2.0)
            VALUE = DERF(SMALLX)       
            PHINAY = (VALUE/2.0) + 0.5
            IF(PHINAY.LE..000001)PHINAY=.00001
            PSIYEA2 = TWOPI*DEXP((-.5)*BETA*BETA*PSIYEA*PSIYEA)
            PSINAY2 = TWOPI*DEXP((-.5)*BETA*BETA*PSINAY*PSINAY)
            XINVMILLSY = PSIYEA2/PHIYEA
            XINVMILLSN = PSINAY2/PHINAY
!
!  IF NOT MISSING DATA
!
            IF(RCVOTE9(KLOCATE(JJJJ),J).EQV..FALSE.)THEN
               KCHOICE=KCHOICE+1
!
!  IF YES
!
               IF(RCVOTE1(KLOCATE(JJJJ),J).EQV..TRUE.)THEN
                  F=F+LOG(PHIYEA)
                  IF(PHIYEA.GT.PHINAY)KORRECT=KORRECT+1
!
! DERIVATIVES FOR LEGISLATORS
!
                  DO 67 K=1,NS
       G(NS*(KKKK-1)+K)=G(NS*(KKKK-1)+K)-XINVMILLSY*&
        (((-1.0)*(DEXP((-1.0)*DISTYEA))*&
        (XDATAUNIQUE(KKKK,K)-&
        ZDATAY(J+KTOTQ,K)))+((DEXP((-1.0)*DISTNAY))*&
        (XDATAUNIQUE(KKKK,K)-ZDATAN(J+KTOTQ,K))))
  67              CONTINUE 
               ENDIF
!
!  IF NO
!
               IF(RCVOTE1(KLOCATE(JJJJ),J).EQV..FALSE.)THEN
                  F=F+LOG(PHINAY)
                  IF(PHIYEA.LT.PHINAY)KORRECT=KORRECT+1
!
! DERIVATIVES FOR LEGISLATORS
!
                  DO 68 K=1,NS
       G(NS*(KKKK-1)+K)=G(NS*(KKKK-1)+K)-XINVMILLSN*&
        (((-1.0)*(DEXP((-1.0)*WEIGHT(K)*WEIGHT(K)*DISTNAY))*&
        (XDATAUNIQUE(KKKK,K)-&
        ZDATAN(J+KTOTQ,K)))+((DEXP((-1.0)*WEIGHT(K)*WEIGHT(K)*DISTYEA))*&
        (XDATAUNIQUE(KKKK,K)-ZDATAY(J+KTOTQ,K))))
  68              CONTINUE
               ENDIF
            ENDIF
         ENDIF
  35  CONTINUE
!
!
 8888 CONTINUE
 9999 CONTINUE
!
      F=-F
!
      IF(JPRINT.EQ.1)WRITE(*,306)NIIII,ICALL,KCHOICE,KORRECT,KKKK,F
      IF(JPRINT.EQ.1)WRITE(21,306)NIIII,ICALL,KCHOICE,KORRECT,KKKK,F
!
!  N = # OF PARAMETERS
!  M = 5 FOR LINE SEARCH
!  X(N) IS THE PARAMETER VECTOR
!  F = FUNCTION VALUE
!  G(N) IS THE VECTOR OF DERIVATIVES
!  DIAGCO = .FALSE.
!  DIAG(N) = DIAGNOAL OF HESSION IF DIAGCO = .TRUE.
!  IPRINT(2) = WRITES OUTPUT
!  EPS = 0.01
!  XTOL = 0.0001
!  W(N*(2*M + 1) + 2*M)
!  IFLAG = 0, NO ERRORS
!
!      DO 69 I=1,N
!      WRITE(99,699)I,X(I),G(I)
!  69  CONTINUE
! 699  FORMAT(I7,2F15.6)
      CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
      WRITE(*,303)NIIII,ICALL,KKKK,N,F
      IF(IFLAG.LE.0) GO TO 50
      ICALL=ICALL + 1
!     We allow at most 2000 evaluations of F and G
      IF(ICALL.GT.100) GO TO 50
      GO TO 20
  50  CONTINUE
!
      CALL CORR2(N,NPMAX,NQTOTALMAX,NDIMENSIONS,X,XSAVE,R)
      WRITE(42,307)N,R
      WRITE(*,307)N,R
!
!      WRITE(*,202)
!  202 FORMAT(' HELLO DOLLY',I4,5I6)
! 2022 FORMAT(' HELLO DOLLY',I10,5I6)
! 2023 FORMAT(' HELLO POLLY',I10,5I6)
!
!  TRANSFER X(.) PARAMETER VECTOR BACK INTO COORDINATE MATRICES
!
      KKK=0
      DO 617 I=1,NPUNIQUE
      SUM=0.0
      DO 618 K=1,NS
      KKK=KKK+1
      XDATAUNIQUE(I,K)=X(KKK)
      SUM=SUM+XDATAUNIQUE(I,K)**2
  618 CONTINUE
!
!  DO CHECK ON LEGISLATOR POINT -- IF IT IS OUTSIDE THE
!   HYPERSPHERE PULL IT BACK INSIDE
!
      IF(SUM.GT.1.0)THEN
         DO 55 K=1,NS
         XDATAUNIQUE(I,K)=(XDATAUNIQUE(I,K)/SQRT(SUM))
  55     CONTINUE
      ENDIF
!
!      WRITE(*,2022)I
!      WRITE(48,201)JCONGUNIQUE(I),JD1UNIQUE(I),&
!                  JSTATEUNIQUE(I),JDISTUNIQUE(I),&
!                  (JSTAUNIQUE(I,JJ),JJ=1,7),JPARTYUNIQUE(I),&
!                  (JNAMEUNIQUE(I,JJ),JJ=1,11),&
!                  (XDATAUNIQUE(I,K),K=1,NDIMENSIONS)
  617 CONTINUE
!
!  TRANSFER X(.) BACK INTO XDATA(:)
!
      KKKK=0
      DO 9998 IIII=1,99999
      DO 44 I=1,NCONGMAX
      KLOCATE(I)=0
      LLOCATE(I)=0
  44  CONTINUE
      KK=0
      DO 49 J=NFIRST,NLAST
      IF(LWHERE(IIII,J).EQV..TRUE.)THEN
         KK=KK+1
         IF(KK.EQ.1)KKKK=KKKK+1
         KLOCATE(KK)=KWHERE(IIII,J)
         LLOCATE(KK)=J
         DO 69 K=1,NDIMENSIONS
         XDATA(KLOCATE(KK),K)=XDATAUNIQUE(KKKK,K)
  69     CONTINUE
      ENDIF
  49  CONTINUE
      IF(KK.EQ.0)GO TO 9998
!      WRITE(*,2023)kkkk
 9998 CONTINUE
!
!  WRITE OUT SCALABLE LEGISLATORS
!
!  CALCULATE THE OUTPUT LOG-LIKELIHOOD
!
      CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,NLEGISLATORS,&
                 NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
!
!
      KK=0
      DO 1947 I=1,NLEGSCALED
!
      GMPX=EXP(XBIGLOG(I)/FLOAT(KBIGLOG(I,1)))
      IF(NIIII.EQ.5)THEN
      WRITE(43,201)JCONG(I),JD1(I),JSTATE(I),JDIST(I),&
                  (JSTA(I,JJ),JJ=1,7),JPARTY(I),&
                  (JNAME(I,JJ),JJ=1,11),&
                  ((XDATA(I,K)*WEIGHT(K)),K=1,NDIMENSIONS),&
                   XBIGLOG(I),(KBIGLOG(I,JJ),JJ=1,2),&
                  GMPX
!      WRITE(43,202)XBIGLOG(I),(KBIGLOG(I,JJ),JJ=1,2),&
!                  GMPX
      ENDIF
      DO 647 K=1,NS
      KK=KK+1
      XSAVE(KK)=XDATA(I,K)
  647 CONTINUE
!
 1947 CONTINUE
!
      CALL CORR2(KK,NPMAX,NQTOTALMAX,NDIMENSIONS,XSAVE2,XSAVE,R)
      WRITE(42,308)KK,R
      WRITE(*,308)KK,R
!
!
      DEALLOCATE(XSAVE)
      DEALLOCATE(XSAVE2)
      DEALLOCATE(KLOCATE)
      DEALLOCATE(LLOCATE)
!
      RETURN
      END
!*********************************************************************

!      subroutine test_local_min_rc_one ( a, b, t, f, title )
      subroutine test_local_min_rc_one (KICK,KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a, b, t, TITLE )

!*********************************************************************72
!
!c TEST_LOCAL_MIN_RC_ONE tests LOCAL_MIN_RC on one test function.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    15 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, double precision A, B, the endpoints of the interval.
!
!    Input, double precision T, a positive absolute error tolerance.
!
!!    function, of the form "FUNCTION F ( X )", which evaluates the
!    function whose local minimum is being sought.
!
!    Input, character*(*) TITLE, a title for the problem.
!
!      implicit none
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER KWHERE(99999,NCONGMAX),&
              NUMCONG(NCONGMAX),NUMCONGRC(NCONGMAX),&
              INUMRC(NQTOTALMAX,2),&
              ICONGRC(NQTOTALMAX),NXDATA(NCONGMAX),KBIGLOG(NPMAX,2)
!
      DIMENSION XDATA(NPMAX,NSMAX),DYN(NQTOTALMAX,NSMAX),&
                ZMID(NQTOTALMAX,NSMAX),ZDF(150000,4),&
                WEIGHT(NSMAX+1),&
                XBIGLOG(NPMAX)
!
      LOGICAL*1 RCVOTE1(NPMAX,NQMAX),RCVOTE9(NPMAX,NQMAX),&
                LWHERE(99999,NCONGMAX)
!
      double precision a
      double precision a2
      double precision arg
      double precision b
      double precision b2
!      double precision f
!      external f
      double precision fa
      double precision fb
      double precision fx
      integer status
      integer step
      double precision t
      character*(*) title
      double precision value

      write (21, '(a)' ) ' '
      write (21, '(a)' ) title
      write (21, '(a)' ) ' '
      write (21, '(a)' ) '  Step      X                          F(X)'
      write (21, '(a)' ) ' '
      step = 0

      arg = a
!      value = f ( arg )
       WEIGHT(KICK)=ARG      
        CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,&
                 NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
      value=-XPLOG
      write (21, '(2x,i4,2x,g24.16,2x,2g24.16)' ) step, arg, value
!
      arg = b
!      value = f ( arg )
       WEIGHT(KICK)=ARG      
       CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,&
                 NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
      value=-XPLOG
      write (21, '(2x,i4,2x,g24.16,2x,2g24.16)' ) step, arg, value
!
      a2 = a
      b2 = b
      status = 0
      WEIGHT(KICK)=(A+B)/2.0
10    continue

!        call local_min_rc ( a2, b2, arg, status, value )
         ARG=WEIGHT(KICK)
        CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,&
                 NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
        value=-XPLOG  
         call local_min_rc (KTOT,KLASS,GMP,XPLOG,&
                 NDIMENSIONS,NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a2, b2, arg, status, value )

        if ( status .lt. 0 ) then
          write (21, '(a)' ) ' '
          write (21, '(a)' ) 'TEST_LOCAL_MIN_RC_ONE - Fatal error!'
          write (21, '(a)' )&
          '  LOCAL_MIN_RC returned negative status.'
          go to 20
        end if

!        value = f ( arg )
        WEIGHT(KICK)=ARG      
        CALL PLOG(KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,&
                 NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT)
        value=-XPLOG  
        step = step + 1
        write (21, '(2x,i4,2x,g24.16,2x,2g24.16)' ) step, arg, value

        if ( status .eq. 0 ) then
          go to 20
        end if

      go to 10

20    continue

      return
      end
!*********************************************************************

      subroutine local_min_rc (KTOT,KLASS,GMP,XPLOG,NDIMENSIONS,&
                 NLEGISLATORS,NLEGSCALED,&
                 NFIRST,NLAST,NPMAX,&
                 NQMAX,NSMAX,NQTOTALMAX,NCONGMAX,&
                 NQTOTAL,KWHERE,NUMCONG,NUMCONGRC, &
                 RCVOTE1,RCVOTE9,LWHERE,&
                 INUMRC,WEIGHT,KBIGLOG,XBIGLOG,&
                 XDATA,DYN,ZMID,ICONGRC,NXDATA,ZDF,NDEVIT,XDEVIT,&
                 a, b, arg, status, value )
!      subroutine local_min_rc ( a, b, arg, status, value )

!*********************************************************************72
!
!c LOCAL_MIN_RC seeks a minimizer of a scalar function of a scalar variable.
!
!  Discussion:
!
!    This routine seeks an approximation to the point where a function
!    F attains a minimum on the interval (A,B).
!
!    The method used is a combination of golden section search and
!    successive parabolic interpolation.  Convergence is never much
!    slower than that for a Fibonacci search.  If F has a continuous
!    second derivative which is positive at the minimum (which is not
!    at A or B), then convergence is superlinear, and usually of the
!    order of about 1.324...
!
!    The routine is a revised version of the Brent local minimization 
!    algorithm, using reverse communication.
!
!    It is worth stating explicitly that this routine will NOT be
!    able to detect a minimizer that occurs at either initial endpoint
!    A or B.  If this is a concern to the user, then the user must
!    either ensure that the initial interval is larger, or to check
!    the function value at the returned minimizer against the values
!    at either endpoint.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    16 April 2008
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization Without Derivatives,
!    Dover, 2002,
!    ISBN: 0-486-41998-3,
!    LC: QA402.5.B74.
!
!    David Kahaner, Cleve Moler, Steven Nash,
!    Numerical Methods and Software,
!    Prentice Hall, 1989,
!    ISBN: 0-13-627258-4,
!    LC: TA345.K34.
!
!  Parameters
!
!    Input/output, double precision A, B.  On input, the left and right
!    endpoints of the initial interval.  On output, the lower and upper
!    bounds for an interval containing the minimizer.  It is required
!    that A < B.
!
!    Output, double precision ARG, the currently considered point.  The user
!    does not need to initialize this value.  On return with STATUS positive,
!    the user is requested to evaluate the function at ARG, and return
!    the value in VALUE.  On return with STATUS zero, ARG is the routine's
!    estimate for the function minimizer.
!
!    Input/output, integer STATUS, used to communicate between the user
!    and the routine.  The user only sets STATUS to zero on the first call,
!    to indicate that this is a startup call.  The routine returns STATUS
!    positive to request that the function be evaluated at ARG, or returns
!    STATUS as 0, to indicate that the iteration is complete and that
!    ARG is the estimated minimizer.
!
!    Input, double precision VALUE, the function value at ARG, as requested
!    by the routine on the previous call.
!
!  Local parameters:
!
!    C is the squared inverse of the golden ratio.
!
!    EPS is the square root of the relative machine precision.
!
!      implicit none
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER KWHERE(99999,NCONGMAX),&
              NUMCONG(NCONGMAX),NUMCONGRC(NCONGMAX),&
              INUMRC(NQTOTALMAX,2),&
              ICONGRC(NQTOTALMAX),NXDATA(NCONGMAX),KBIGLOG(NPMAX,2)
!
      DIMENSION XDATA(NPMAX,NSMAX),DYN(NQTOTALMAX,NSMAX),&
                ZMID(NQTOTALMAX,NSMAX),ZDF(150000,4),&
                WEIGHT(NSMAX+1),&
                XBIGLOG(NPMAX)
!
      LOGICAL*1 RCVOTE1(NPMAX,NQMAX),RCVOTE9(NPMAX,NQMAX),&
                LWHERE(99999,NCONGMAX)
!
      double precision a
      double precision arg
      double precision b
      double precision c
      double precision d
      double precision e
      double precision eps
      double precision fu
      double precision fv
      double precision fw
      double precision fx
      double precision midpoint
      double precision p
      double precision q
      double precision r
!      double precision r8_epsilon
      integer status
      double precision tol
      double precision tol1
      double precision tol2
      double precision u
      double precision v
      double precision value
      double precision w
      double precision x

      save c
      save d
      save e
      save eps
      save fu
      save fv
      save fw
      save fx
      save midpoint
      save p
      save q
      save r
      save tol
      save tol1
      save tol2
      save u
      save v
      save w
      save x
!
!  STATUS (INPUT) = 0, startup.
!
      if ( status .eq. 0 ) then

        if ( b .le. a ) then
          write (21, '(a)' ) ' '
          write (21, '(a)' ) 'LOCAL_MIN_RC - Fatal error!'
          write (21, '(a)' ) '  A < B is required, but'
          write (21, '(a,g14.6)' ) '  A = ', a
          write (21, '(a,g14.6)' ) '  B = ', b
          status = -1
          stop
        end if

        c = 0.5D+00 * ( 3.0D+00 - sqrt ( 5.0D+00 ) )
!
!
        XPRECISION=22.5
        T=EPSILON(XPRECISION)
        eps = sqrt (T)
        tol = T
!        eps = sqrt ( r8_epsilon ( ) )
!        tol = r8_epsilon ( )
!        eps = sqrt (0.0001)
!        tol = 0.0001
!
        v = a + c * ( b - a )
        w = v
        x = v
        e = 0.0D+00

        status = 1
        arg = x

        return
!
!  STATUS (INPUT) = 1, return with initial function value of FX.
!
      else if ( status .eq. 1 ) then

        fx = value
        fv = fx
        fw = fx
!
!  STATUS (INPUT) = 2 or more, update the data.
!
      else if ( 2 .le. status ) then

        fu = value

        if ( fu .le. fx ) then

          if ( x .le. u ) then
            a = x
          else
            b = x
          end if

          v = w
          fv = fw
          w = x
          fw = fx
          x = u
          fx = fu

        else

          if ( u .lt. x ) then
            a = u
          else
            b = u
          end if

          if ( fu .le. fw .or. w .eq. x ) then
            v = w
            fv = fw
            w = u
            fw = fu
          else if ( fu .le. fv .or. v .eq. x .or. v .eq. w ) then
            v = u
            fv = fu
          end if

        end if

      end if
!
!  Take the next step.
!
      midpoint = 0.5D+00 * ( a + b )
      tol1 = eps * abs ( x ) + tol / 3.0D+00
      tol2 = 2.0D+00 * tol1
!
!  If the stopping criterion is satisfied, we can exit.
!
      if ( abs ( x - midpoint ) .le. &
         ( tol2 - 0.5D+00 * ( b - a ) ) ) then
        status = 0
        return
      end if
!
!  Is golden-section necessary?
!
      if ( abs ( e ) .le. tol1 ) then
        if ( midpoint .le. x ) then
          e = a - x
        else
          e = b - x
        end if

        d = c * e
!
!  Consider fitting a parabola.
!
      else

        r = ( x - w ) * ( fx - fv )
        q = ( x - v ) * ( fx - fw )
        p = ( x - v ) * q - ( x - w ) * r
        q = 2.0D+00 * ( q - r )
        if ( 0.0D+00 .le. q ) then
          p = - p
        end if
        q = abs ( q )
        r = e
        e = d
!
!  Choose a golden-section step if the parabola is not advised.
!
        if ( &
          ( abs ( 0.5D+00 * q * r ) .le. abs ( p ) ) .or. &
          ( p .le. q * ( a - x ) ) .or. &
          ( q * ( b - x ) .le. p ) ) then

          if ( midpoint .le. x ) then
            e = a - x
          else
            e = b - x
          end if

          d = c * e
!
!  Choose a parabolic interpolation step.
!
        else

          d = p / q
          u = x + d

          if ( ( u - a ) .lt. tol2 ) then
            d = sign ( tol1, midpoint - x )
          end if

          if ( ( b - u ) .lt. tol2 ) then
            d = sign ( tol1, midpoint - x )
          end if

        end if

      end if
!
!  F must not be evaluated too close to X.
!
      if ( tol1 .le. abs ( d ) ) then
        u = x + d
      end if

      if ( abs ( d ) .lt. tol1 ) then
        u = x + sign ( tol1, d )
      end if
!
!  Request value of F(U).
!
      arg = u
      status = status + 1

      return
      end
