!
!  
!
      INTEGER NRESPONDENTS,NISSUES,NDIMENSIONS,NMISSING,             &
                EXITSTATUS,KMISS(14,4),MINSCALE,KISSUE(1614,14),     &
                MID(1614),LRESPONDENTS(14),LMARK(1614)
      DOUBLE PRECISION FITS(3*6),PSIMATRIX(1614*((3*4)/2)),          &
                       WMATRIX(14*((3*4)/2)+2*14*3)
      CHARACTER*21 CAND(14)
!
!
!  NRESPONDENTS = NUMBER OF RESPONDENTS IN DATA FILE
!  NISSUES = NUMBER OF ISSUE SCALES IN ANALYSIS
!  NDIMENSIONS = NUMBER OF LATENT DIMENSIONS TO ESTIMATE
!  NMISSING = MAXIMUM NUMBER OF MISSING VALUES FOR ANY ISSUE SCALE
!  
!
      OPEN(22,FILE='NES1980.DAT',STATUS='OLD')
      OPEN(23,FILE='BLACK23F95.DAT')
      OPEN(24,FILE='BLACK24F95.DAT')
      OPEN(28,FILE='BLACK28F95.DAT')
!
!  FORMAT FOR 1980 DATA
!
  100 FORMAT(8X,I4,527X,I1,13X,I1,11X,I1,11X,I1,11X,I1,13X,I1,1217X,&
                   I1,36X,I1,17X,I1,17X,I1,17X,I1,18X,I1,5X,I1,2X,I1)
  113 FORMAT(1X,I5,8F12.4)
  201 FORMAT(2I7,35F7.3)
  225 FORMAT(1X,A21,I7,35F7.3)
!
      NRESPONDENTS=1614
      NDIMENSIONS=3
      NISSUES=14
      NMISSING=4
      MINSCALE=8
      EXITSTATUS=0
!
      DO 1 I=1,NISSUES
      DO 2 J=1,NMISSING
      KMISS(I,J)=0
  2   CONTINUE
      KMISS(I,2)=8
      KMISS(I,3)=9
      IF(I.EQ.5.OR.I.EQ.14)THEN
         KMISS(I,4)=7
      ENDIF
  1   CONTINUE
!
      CAND(1)= 'LIBERAL/CONSERVATIVE '
      CAND(2)= 'DEFENSE              '
      CAND(3)= 'GOVT SERVICES        '
      CAND(4)= 'INFLATION            '
      CAND(5)= 'ABORTION             '
      CAND(6)= 'TAX CUTS             '
      CAND(7)= 'LIBERAL/CONSERVATIVE '
      CAND(8)= 'GOVT HELP MINORITIES '
      CAND(9)= 'RUSSIA               '
      CAND(10)='WOMENS EQUAL ROLE    '
      CAND(11)='GOVT JOBS            '
      CAND(12)='EQUAL RIGHTS AMEND   '
      CAND(13)='BUSING               '
      CAND(14)='ABORTION             '
!
      DO 111 I=1,NRESPONDENTS
      READ(22,100)MID(I),(KISSUE(I,J),J=1,NISSUES)
  111 CONTINUE
!
      CALL BLACKBOX(NRESPONDENTS,NISSUES,NDIMENSIONS,NMISSING,       &
           KMISS,MINSCALE,MID,KISSUE,CAND,FITS,PSIMATRIX,WMATRIX,    &
           LRESPONDENTS,LMARK,EXITSTATUS)
!
!  WRITE OUT ESTIMATES AS A CHECK
!
!
!  WRITE OUT W MATRIX
!
      KKKK=0
      DO 221 II=1,NDIMENSIONS
      DO 222 I=1,NISSUES
!
!  TRANSFER WMATRIX PARAMETERS
!
      IF(KKKK.EQ.0)THEN
         WRITE(23,225)CAND(I),LRESPONDENTS(I),&
           (WMATRIX((I-1)*(II+2)+JJ),JJ=1,II+2)
      ENDIF
      IF(KKKK.GT.0)THEN
         WRITE(23,225)CAND(I),LRESPONDENTS(I),&
           (WMATRIX((KKKK)*NISSUES+(I-1)*(II+2)+JJ),JJ=1,II+2)
      ENDIF
  222 CONTINUE
      KKKK=KKKK+II+2
  221 CONTINUE
!
!
!  WRITE OUT PSI MATRIX
!
      KKKK=0
      DO 421 II=1,NDIMENSIONS
      DO 422 I=1,NRESPONDENTS
!
!  WRITE OUT ONLY SCALED RESPONDENTS
!
      IF(LMARK(I).EQ.0)GO TO 422
!
!  TRANSFER PSIMATRIX PARAMETERS
!
      IF(KKKK.EQ.0)THEN
         WRITE(24,201)I,MID(I),&
           (PSIMATRIX((I-1)*II+JJ),JJ=1,II)
      ENDIF
      IF(KKKK.GT.0)THEN
         WRITE(24,201)I,MID(I),&
           (PSIMATRIX(KKKK*NRESPONDENTS+(I-1)*II+JJ),JJ=1,II)
      ENDIF
  422 CONTINUE
      KKKK=KKKK+II
  421 CONTINUE
!
      DO 69 J=1,NDIMENSIONS
      WRITE(23,113)J,(FITS(JJ+(J-1)*6),JJ=1,6)
  69  CONTINUE
!
      STOP
      END
!
!************************************************************************
!************************************************************************
!*                                                                     **
!* MODIFIED FOR FORTRAN 95 APRIL 2009                                  **
!*                                                                     **
!*                                                                     **
!*  BLACKBOX.FOR          23 MAY 1998                                  **
!*                                                                     **
!*  THE BLACK BOX -- PERFORMS SINGULAR VALUE DECOMPOSITION ON A MATRIX **
!*    WITH MISSING ENTRIES.  DEVELOPED BY KEITH T. POOLE, JULY-OCTOBER **
!*    1982.  FURTHER DEVELOPMENT IN NOVEMBER-DECEMBER 1992.            **
!*                                                                     **
!*                                                                     **
!************************************************************************
!************************************************************************
!
      SUBROUTINE BLACKBOX(NRESPONDENTS,NISSUES,NDIMENSIONS,NMISSING, &
                          KMISS,MINSCALE,MID,KISSUE,CAND,FITS,       &
                          PSIMATRIX,WMATRIX,LRESPONDENTS,LMARK,      &
                          EXITSTATUS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER NRESPONDENTS,NISSUES,NDIMENSIONS,NMISSING,             &
                EXITSTATUS,KMISS(NISSUES,NMISSING),MINSCALE,         &
                KISSUE(NRESPONDENTS,NISSUES),MID(NRESPONDENTS),      &
                LRESPONDENTS(NISSUES),LMARK(NRESPONDENTS)
      DIMENSION FITS(6*NDIMENSIONS), &
        PSIMATRIX(NRESPONDENTS*((NDIMENSIONS*(NDIMENSIONS+1))/2)),   &
        WMATRIX(NISSUES*((NDIMENSIONS*(NDIMENSIONS+1))/2)+2*NISSUES*&
      NDIMENSIONS)  
!               XBIGONE(3539,60),KID(3539),       &
!               W(152,127),XDATA(3539,125),XT(3539,25),   &
!               RSAVE(25,25),YHAT(3539),UUU(3539,25),VVV(3539,25),    &
!               LID(3539),WORK(3539)
! 
      CHARACTER*21 CAND(NISSUES)
!
      INTEGER, ALLOCATABLE :: KID(:)
      INTEGER, ALLOCATABLE :: LID(:)
      DOUBLE PRECISION, ALLOCATABLE :: YHAT(:)
      DOUBLE PRECISION, ALLOCATABLE :: WORK(:)
      DOUBLE PRECISION, ALLOCATABLE :: XBIGONE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: W(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XDATA(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XT(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: RSAVE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: UUU(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: VVV(:,:)
!
      ALLOCATE(KID(NRESPONDENTS))
      ALLOCATE(LID(NRESPONDENTS))
      ALLOCATE(YHAT(NRESPONDENTS))
      ALLOCATE(WORK(NRESPONDENTS+NISSUES))
      ALLOCATE(XBIGONE(NRESPONDENTS,NISSUES))
      ALLOCATE(W(NISSUES,NDIMENSIONS+2))
      ALLOCATE(XDATA(NRESPONDENTS,NDIMENSIONS))
      ALLOCATE(XT(NRESPONDENTS,NISSUES))
      ALLOCATE(RSAVE(NDIMENSIONS,3))
      ALLOCATE(UUU(NRESPONDENTS,NISSUES))
      ALLOCATE(VVV(NRESPONDENTS,NISSUES))
!
  70  FORMAT(' ITERATION RECORD'                                     &
      /3X,'DIM',5X,'ERROR',4X,'EXPLAINED',                           &
      6X,'PERCENT',2X,'CUM PERCENT',3X,'R-SQUARE',2X,'STD ERR EST')
  113 FORMAT(1X,I5,8F12.4)
  195 FORMAT(' ITERATION SUMMARY')
  196 FORMAT(' NUMBER OF DIMENSIONS='I4)
  197 FORMAT(5X,'PSI-TRANSPOSE*PSI')
  198 FORMAT(I6,20F7.3)
  199 FORMAT(1X,'CONSTRAINT CHECKS ON PSI AND W'/                    &
             5X,'SUM OF COLUMNS OF PSI'/                             &
             6X,35F10.4)
  200 FORMAT(I6,35F10.4)
  201 FORMAT(2I7,35F7.3)
  202 FORMAT(5X,'W-TRANSPOSE*W')
  225 FORMAT(1X,A21,I7,35F7.3)
  301 FORMAT(' NUMBER OF CASES',3I6)
  315 FORMAT(' R-SQUARE ',I8,F7.3)
  330 FORMAT(I5,7F10.4)
  337 FORMAT(' RANK CHECK OF PSI*W',I5)
 2222 FORMAT(1X,78('*'))
!
!
!
!  PRINTER SWITCH -- IPRNT=1 NO WRITE TO DISK
!                    IPRNT=0 WRITE TO DISK
      IPRNT=1
!
      LWORK=NRESPONDENTS+NISSUES
      NS=NDIMENSIONS
      NY=NISSUES
      NMISS=NMISSING
      KVMIN=MINSCALE
!
!      READ(4,1002)NS,NY,NMISS,KVMIN
!
      II=0
      DO 8888 I=1,NRESPONDENTS
      LMARK(I)=0
!  500 READ(22,FMT,END=400)MID,(KISSUE(J),J=1,NY)
      IMARK=0
      DO 1 J=1,NY
      DO 112 K=1,NMISS
      IF(KISSUE(I,J).EQ.KMISS(J,K))THEN 
         KISSUE(I,J)=-999
         IMARK=IMARK+1
         GO TO 1
      ENDIF
  112 CONTINUE
  1   CONTINUE
      IF((NY-IMARK).LT.KVMIN)GO TO 8888
      II=II+1
      LMARK(I)=1
      KID(II)=MID(I)
      DO 11 J=1,NY
      XBIGONE(II,J)=FLOAT(KISSUE(I,J))
  11  CONTINUE
 8888 CONTINUE
      IF(IPRNT.EQ.0)WRITE(23,301)I,II
      NOBS=II
      IF(IPRNT.EQ.0)WRITE(23,301)NOBS
!
!  SUBROUTINE BLACKB--GENERATES MULTI-DIMENSIONAL STARTING VALUES 
!    FOR LEGISLATORS BY PERFORMING A PSEUDO-SINGULAR VALUE 
!    DECOMPOSITION OF NY ROLL CALLS
!
!
      KKKK=0
      KKKKK=0
      IF(IPRNT.EQ.0)WRITE(23,2222)
      IF(IPRNT.EQ.0)WRITE(23,2222)
      DO 999 KKK=1,NS
      IF(IPRNT.EQ.0)WRITE(23,2222)
      IF(IPRNT.EQ.0)WRITE(23,196)KKK
!      WRITE(*,196)KKK
      call BLACKB(NOBS,NRESPONDENTS,NISSUES,NDIMENSIONS,KKK,         &
                  XBIGONE,XDATA,W,SVSUM,IPRNT)
!  
!  CONSTRAINT CHECKS
!
      DO 233 K=1,KKK
      SUM=0.0
      DO 234 I=1,NOBS
      SUM=SUM+XDATA(I,K)
  234 CONTINUE
      VVV(1,K)=SUM/FLOAT(NOBS)
  233 CONTINUE
      IF(IPRNT.EQ.0)WRITE(23,199)(VVV(1,K),K=1,KKK)
      DO 227 J=1,KKK
      DO 227 K=1,KKK
      SUM=0.0
      DO 228 I=1,NOBS
      SUM=SUM+XDATA(I,J)*XDATA(I,K)
  228 CONTINUE
  227 VVV(J,K)=SUM
      IF(IPRNT.EQ.0)WRITE(23,197)
      DO 229 I=1,KKK
      IF(IPRNT.EQ.0)WRITE(23,200)I,(VVV(I,J),J=1,KKK)
  229 CONTINUE
      IF(IPRNT.EQ.0)WRITE(23,202)
      DO 230 J=1,KKK
      DO 230 K=1,KKK
      SUM=0.0
      DO 231 I=1,NY
      SUM=SUM+W(I,J+1)*W(I,K+1)
  231 CONTINUE
  230 VVV(J,K)=SUM
      DO 232 I=1,KKK
      IF(IPRNT.EQ.0)WRITE(23,200)I,(VVV(I,J),J=1,KKK)
  232 CONTINUE
      ASUM=0.0
      BSUM=0.0
      CSUM=0.0
      DSUM=0.0
      ESUM=0.0
      SUME=0.0
      KK=0
      DO 223 J=1,NY
      AASUM=0.0
      BBSUM=0.0
      CCSUM=0.0
      DDSUM=0.0
      EESUM=0.0
      KJJ=0
      DO 224 I=1,NOBS
      SUM=0.0
      DO 726 K=1,KKK
      SUM=SUM+XDATA(I,K)*W(J,K+1)
  726 CONTINUE
      XT(I,J)=SUM
      AA=SUM+W(J,1)
      IF(XBIGONE(I,J).EQ.-999.0)GO TO 224
      BB=XBIGONE(I,J)
      SUME=SUME+(AA-BB)**2
      AASUM=AASUM+AA
      BBSUM=BBSUM+BB
      CCSUM=CCSUM+AA*AA
      DDSUM=DDSUM+BB*BB
      EESUM=EESUM+AA*BB
      KJJ=KJJ+1
  224 CONTINUE
      AAA=FLOAT(KJJ)*EESUM-AASUM*BBSUM
      BBB=FLOAT(KJJ)*CCSUM-AASUM*AASUM
      CCC=FLOAT(KJJ)*DDSUM-BBSUM*BBSUM
      RRR=0.0
      IF(ABS(BBB*CCC).GT.0.0)RRR=(AAA*AAA)/(BBB*CCC)
      W(J,KKK+2)=RRR
      LID(J)=KJJ
      ASUM=ASUM+AASUM
      BSUM=BSUM+BBSUM
      CSUM=CSUM+CCSUM
      DSUM=DSUM+DDSUM
      ESUM=ESUM+EESUM
      KK=KK+KJJ
  223 CONTINUE
      AA=FLOAT(KK)*ESUM-ASUM*BSUM
      BB=FLOAT(KK)*CSUM-ASUM*ASUM
      CC=FLOAT(KK)*DSUM-BSUM*BSUM
      RRR=(AA*AA)/(BB*CC)
!      WRITE(*,315)KK,RRR
      IF(IPRNT.EQ.0)WRITE(23,315)KK,RRR
!
      RSAVE(KKK,1)=SUME
      RSAVE(KKK,2)=RRR
      RSAVE(KKK,3)=SQRT(SUME/FLOAT(KK - KKK*(NOBS+NY) + NY))
!
!  *****patched in code*****
!
      XTOL=.001
!      CALL LSVRR(NOBS,NY,XT,3539,21,XTOL,IRANK,YHAT,UUU,             &
!                 3539,VVV,3539)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
      CALL DGESVD('S','S',NOBS,NY,XT,NRESPONDENTS,YHAT,UUU,           &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
      IF(IPRNT.EQ.0)WRITE(23,337)INFO
!      WRITE(*,337)INFO
      DO 464 I=1,NS+3
      IF(IPRNT.EQ.0)WRITE(23,330)I,YHAT(I)
!      WRITE(*,330)I,YHAT(I)
  464 CONTINUE
!
!
!  WRITE OUT COLUMN PARAMETERS
!
      DO 222 I=1,NY
!
!  TRANSFER WMATRIX PARAMETERS
!
      DO 532 JJ=1,KKK+2
      IF(KKKK.EQ.0)THEN
         WMATRIX((I-1)*(KKK+2)+JJ)=W(I,JJ)
      ENDIF
      IF(KKKK.GT.0)THEN
         WMATRIX((KKKK)*NY+(I-1)*(KKK+2)+JJ)=W(I,JJ)
      ENDIF
  532 CONTINUE
      LRESPONDENTS(I)=LID(I)
!      WRITE(28,225)CAND(I),LID(I),(W(I,JJ),JJ=1,KKK+2)
  222 CONTINUE
!
!  WRITE OUT ROW PARAMETERS
!
      II=0
      DO 742 I=1,NRESPONDENTS
      IF(LMARK(I).EQ.0)GO TO 742
      II=II+1
      DO 741 JJ=1,KKK
      IF(KKKKK.EQ.0)THEN
         PSIMATRIX((I-1)*KKK+JJ)=XDATA(II,JJ)
      ENDIF
      IF(KKKKK.GT.0)THEN
         PSIMATRIX(KKKKK*NRESPONDENTS+(I-1)*KKK+JJ)=XDATA(II,JJ)
      ENDIF
!
  741 CONTINUE
!      WRITE(24,201)II,KID(II),(XDATA(II,J),J=1,KKK)
!      WRITE(*,201)II,KID(I), (UUU(II,J),J=1,KKK)
  742 CONTINUE
      KKKK=KKKK+KKK+2
      KKKKK=KKKKK+KKK
  999 CONTINUE
      IF(IPRNT.EQ.0)WRITE(23,2222)
      IF(IPRNT.EQ.0)WRITE(23,70)
!      WRITE(*,70)
      DO 69 J=1,NS
      AA=(SVSUM-RSAVE(J,1))/SVSUM
      CCC=SVSUM-RSAVE(J,1)
      IF(J.EQ.1)BB=SVSUM
      IF(J.GT.1)BB=RSAVE(J-1,1)
      SUM=(BB-RSAVE(J,1))/SVSUM
!
!  TRANSFER FITS
!
      FITS(1+(J-1)*6)=RSAVE(J,1)
      FITS(2+(J-1)*6)=CCC
      FITS(3+(J-1)*6)=SUM*100.0
      FITS(4+(J-1)*6)=AA*100.0
      FITS(5+(J-1)*6)=RSAVE(J,2)
      FITS(6+(J-1)*6)=RSAVE(J,3)
!
      IF(IPRNT.EQ.0)THEN
      WRITE(23,113)J,RSAVE(J,1),CCC,SUM*100.0,AA*100.0,RSAVE(J,2),   &
                    RSAVE(J,3)
!      WRITE(* ,113)J,RSAVE(J,1),CCC,SUM*100.0,AA*100.0,RSAVE(J,2),   &
!                    RSAVE(J,3)
      ENDIF
  69  CONTINUE
      IF(IPRNT.EQ.0)WRITE(23,2222)
!
!      WRITE(23,2222)
!
      DEALLOCATE(KID)
      DEALLOCATE(LID)
      DEALLOCATE(YHAT)
      DEALLOCATE(XBIGONE)
      DEALLOCATE(W)
      DEALLOCATE(XDATA)
      DEALLOCATE(XT)
      DEALLOCATE(RSAVE)
      DEALLOCATE(UUU)
      DEALLOCATE(VVV)
!
      RETURN
      END
!
!  ***********************************************************************
!    SUBROUTINE BLACKB
!      
!      PERFORMS THE FOLLOWING DECOMPOSITION:
!
!           X = [PW' + Jc'] + E
!         
!          WHERE X IS NPxNY, P IS NPxNF, W IS NYxNF, J IS A NP LENGTH
!            VECTORS OF 1'S, c IS A NY VECTOR OF CONSTANTS,
!            AND E IS NPxNY MATRIX OF ERROR WHICH IS MINIMIZED
!
!  ***********************************************************************
!  
!
      SUBROUTINE BLACKB(NP,NRESPONDENTS,NISSUES,NDIMENSIONS,         &
                        NF,XBIGONE,XDATA,W,SVSUM,IPRNT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XBIGONE(NRESPONDENTS,NISSUES),                       &
                XDATA(NRESPONDENTS,NDIMENSIONS),                     &
                W(NISSUES,NDIMENSIONS+2) 
!                X(3539,152),XX(408),FAL1(152),FAL2(152),             &
!                PSIX(3539,152),XS(3539,152),LL(152),                 &
!                R(152,152),D(152),ROOTC(3539,152),                   &
!                CROOT(3539,152),DD(152),DC(152),CC(152),XT(3539,152),&
!                XSS(3539,152),MPOS(152),TSUM(408),DDD(152),          &
!                DX(152),            &
!                UUU(3539,152),VVV(3539,152),WSAVE(127,127),          &
!                XSAVE(3539,60),YHAT(3539),WORK(3539)
!
      INTEGER, ALLOCATABLE :: LL(:)
      INTEGER, ALLOCATABLE :: MPOS(:)
      DOUBLE PRECISION, ALLOCATABLE :: XX(:)
      DOUBLE PRECISION, ALLOCATABLE :: D(:)
      DOUBLE PRECISION, ALLOCATABLE :: DD(:)
      DOUBLE PRECISION, ALLOCATABLE :: DC(:)
      DOUBLE PRECISION, ALLOCATABLE :: CC(:)
      DOUBLE PRECISION, ALLOCATABLE :: TSUM(:)
      DOUBLE PRECISION, ALLOCATABLE :: DDD(:)
      DOUBLE PRECISION, ALLOCATABLE :: DX(:)
      DOUBLE PRECISION, ALLOCATABLE :: YHAT(:)
      DOUBLE PRECISION, ALLOCATABLE :: WORK(:)
      DOUBLE PRECISION, ALLOCATABLE :: X(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: PSIX(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XS(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: R(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ROOTC(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: CROOT(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XT(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XSS(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: UUU(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: VVV(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: WSAVE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XSAVE(:,:)
!
      ALLOCATE(LL(NISSUES))
      ALLOCATE(MPOS(NISSUES))
      ALLOCATE(XX(2*NISSUES))
      ALLOCATE(D(NISSUES))
      ALLOCATE(DD(NISSUES))
      ALLOCATE(DC(NISSUES))
      ALLOCATE(CC(NISSUES))
      ALLOCATE(TSUM(2*NISSUES))
      ALLOCATE(DDD(NISSUES))
      ALLOCATE(DX(NISSUES))
      ALLOCATE(YHAT(NISSUES))
      ALLOCATE(WORK(NRESPONDENTS+NISSUES))
!
      ALLOCATE(X(NRESPONDENTS,NISSUES))
      ALLOCATE(PSIX(NRESPONDENTS,NISSUES))
      ALLOCATE(XS(NRESPONDENTS,NISSUES))
      ALLOCATE(R(NISSUES,NISSUES))
      ALLOCATE(ROOTC(NRESPONDENTS,NISSUES))
      ALLOCATE(CROOT(NRESPONDENTS,NISSUES))
      ALLOCATE(XT(NRESPONDENTS,NISSUES))
      ALLOCATE(XSS(NRESPONDENTS,NISSUES))
      ALLOCATE(UUU(NRESPONDENTS,NISSUES))
      ALLOCATE(VVV(NRESPONDENTS,NISSUES))
      ALLOCATE(WSAVE(NISSUES,2*NDIMENSIONS))
      ALLOCATE(XSAVE(NRESPONDENTS,NISSUES))
!
  23  FORMAT(' CORRELATION MATRIX FOR PSIX')
  70  FORMAT(' ITERATION RECORD'/3X,'DIM',5X,'ERROR',5X,'EXPLAINED', &
      6X,'PERCENT',2X,'CUM PERCENT',3X,'R-SQUARE')
  86  FORMAT(1X,'SINGULAR VALUES OF ESTIMATED MATRICES'/             &
             1X,'FIRST COLUMN:  ORIGINAL MATRIX WITH FILLED IN MISSING E&
      NTRIES'/                                                       &
             1X,'SECOND COLUMN: REPRODUCED MATRIX -- PSI*W + Jc'/    &
             1X,'THIRD COLUMN:  ORIGINAL MATRIX WITH FILLED IN MISSING E&
      NTRIES MINUS THE ORIGINAL COLUMN MEANS'/                       &
             1X,'FOURTH COLUMN:  PSI*W')
  113 FORMAT(1X,2I5,8F12.4)
  199 FORMAT(' CONSTRAINT CHECKS ON PSI AND W'/6X,35F7.3)
  200 FORMAT(I6,35F10.4)
  300 FORMAT(1X,'NUMBER OF ROWS               =',I6/                 &
             1X,'NUMBER OF COLUMNS            =',I6/                 &
             1X,'TOTAL NUMBER OF DATA ENTRIES =',I6/                 &
             1X,'NUMBER MISSING ENTRIES       =',I6/                 &
             1X,'PERCENT MISSING DATA         =',F15.5/              &
             1X,'SUM OF SQUARES GRAND MEAN    =',F15.5)              
  301 FORMAT(2I5,F7.3)
  335 FORMAT(' CORRECTION ROW=',I3,1X,' # POSITIVE=',I5)
  336 FORMAT(1X,25I3/1X,25I3)
  337 FORMAT(' RANK CHECK',I5)
  338 FORMAT(I5,7F10.4)
  521 FORMAT(' CROSS PRODUCT MATRIX FOR PSIX')
 1113 FORMAT(I3,12F10.3)
 2222 FORMAT(1X,78('*'))
 7200 FORMAT(' ERROR PARAMETER EIGENVALUE ROUTINE=',I4)
 9011 FORMAT(I5,10X,3F10.4)
 9012 FORMAT(3I5,3F10.4)
!
      LWORK=NRESPONDENTS+NISSUES
      NY=NISSUES
      NF1=NF+1
      DO 349 J=1,NY
      DC(J)=0.0
      CC(J)=0.0
  349 LL(J)=0.0
      KTOT=0
      SVSUM=0.0
      SWSUM=0.0
      IF(IPRNT.EQ.0)WRITE(23,2222)
!
!
      DO 600 I=1,NP
      SUM=0.0
      SUMA=0.0
      KK=0
      DO 30 J=1,NY
      X(I,J)=XBIGONE(I,J)
      IF(X(I,J).NE.-999.0)THEN
         SUM=SUM+X(I,J)*X(I,J)
         SUMA=SUMA+X(I,J)
         KK=KK+1
      ENDIF
      XS(I,J)=X(I,J)
  30  CONTINUE
      KTOT=KTOT+KK
      SVSUM=SVSUM+SUM
      SWSUM=SWSUM+SUMA
      DO 31 J=1,NY
      XSS(I,J)=X(I,J)
      IF(X(I,J).EQ.-999.0)GO TO 31
      LL(J)=LL(J)+1
      DC(J)=DC(J)+X(I,J)
      CC(J)=CC(J)+X(I,J)*X(I,J)
  31  CONTINUE
  600 CONTINUE
!
!   COMPUTE NY COLUMN MEANS PLACE IN DC(NY), COMPUTE SUM OF SQUARES
!     AND PERCENT MISSING DATA
!
      DO 32 J=1,NY
      D(J)=CC(J)-(DC(J)*DC(J))/FLOAT(LL(J))
      DC(J)=DC(J)/FLOAT(LL(J))
!      WRITE(*,301)J,LL(J),DC(J)
!      WRITE(23,301)J,LL(J),DC(J)
  32  CONTINUE
      SVSUM=SVSUM-(SWSUM**2)/FLOAT(KTOT)
      LTOT=NP*NY
      PMISS=FLOAT(LTOT-KTOT)/FLOAT(LTOT)
      LTOT=LTOT-KTOT
      IF(NF.EQ.1)THEN
         IF(IPRNT.EQ.0)THEN
            WRITE(23,300)NP,NY,KTOT,LTOT,PMISS*100.0,SVSUM
!         WRITE(* ,300)NP,NY,KTOT,LTOT,PMISS*100.0,SVSUM
            WRITE(23,2222)
!         WRITE(*,2222)
         ENDIF
      ENDIF   
!
!  COMPUTE CORRELATION NY BY NY CORRELATION MATRIX AND VECTOR OF SIGN
!    CHANGES FOR COLUMNS OF X/XS/XSS
!
      call CORR2(NRESPONDENTS,NISSUES,NP,NY,X,R,LL,MPOS,KS,KPOS,1)
!
!  MASTER LOOP FOR NF DIMENSIONS--ONE DIMENSION IS EXTRACTED AT A TIME.
!     THE RESIDUAL MATRIX IS COMPUTED AFTER EACH DIMENSION IS EXTRACTED
!     AND IT IS USED AS X IN THE FORMULA BELOW.  AT THE END OF THIS
!     PROCESS:
!
!              P'P = I, WITH P BEING NOBSxNF, W NYxNF, c NYx1
!
!           X = [PW' + Jc'] + E
!
      DO 9999 JJJ=1,NF
      KKK=JJJ
      DO 10 J=1,NY
      XX(J)=1.0
      XX(J+NY)=-DC(J)
      IF(JJJ.GT.1)XX(J+NY)=0.0
  10  CONTINUE
      XXK=0.0
      TXB=0.0
      KKT=0
!
!  AT FIRST ITERATION JJJ = 1 THE COLUMN MEANS ARE SUBTRACTED
!     FOR SUBSEQUENT ITERATIONS, COLUMN MEANS = 0.  NOTE USE OF
!     COLUMN SIGN CHANGES (LL(NY)) FROM CORR2
!
      DO 3 I=1,NP
      KK=0
      SUM=0.0
      DO 1 J=1,NY
      IF(X(I,J).EQ.-999.0)GO TO 1
      SUM=SUM+(X(I,J)*XX(J)+XX(J+NY))*FLOAT(LL(J))
      KK=KK+1
  1   CONTINUE
!
!  XT(NOBS,1) = (X - COL. MEAN)xSIGN CHANGE, SUMMED OVER ITH ROW--THIS
!     IS THE STARTING ESTIMATE OF THE FIRST COLUMN OF P (PSIX BELOW).
!     THE MEAN OF XT(NOBS,1) IS THEN SUBTRACTED OUT BELOW SO THAT THE
!     STARTING ESTIMATE SUMS TO ZERO.
!     XT(NOBS,2) = 1 FOR THE LINEAR REGRESSIONS BELOW---TO PICK UP THE
!       INTERCEPT VECTOR, c
!
      WXB=SUM/(FLOAT(KK))
!
      KKT=KKT+1
      TXB=TXB+WXB
      XT(I,1)=WXB
      XXK=XXK+XT(I,1)*XT(I,1)
  3   CONTINUE
      TXB=TXB/FLOAT(KKT)
      XXK=XXK-FLOAT(KKT)*TXB*TXB
      DO 40 I=1,NP
      XT(I,1)=(XT(I,1)-TXB)
      XT(I,2)=1.0
  40  CONTINUE
!
!  ALTERNATING LEAST SQUARES IS NOW PERFORMED TO FIRST ESTIMATE W AND c--
!    SUBROUTINE REG--ESTIMATES OF W AND c ARE RETURNED IN W(NY,NF+1)
!  W AND c ARE THEN USED TO ESTIMATE P IN SUBROUTINE REG2--
!  RESIDUALS ARE RETURNED IN X(NOBS,NY)
!
      DO 444 MM=1,4
      call REG(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,1,NY,TSUM,        &
                                   W,XS,X,XT,IPRNT,1,KKK,AREG,BREG)
      call REG2(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,1,NY,W,XS,       &
                                  X,XT,PXB,PXS,KKK,IPRNT,AREG,BREG)
!
!  PXB IS THE MEAN OF THE ESTIMATED COLUMN OF P
!  PXS IS THE VARIANCE OF THE ESTIMATED COLUMN OF P
!
      QXS=PXS
      XCOR=SQRT(XXK/PXS)
!      WRITE(23,113)KKK,MM,XXK,QXS,PXS,XCOR,PXB,TXB
!      WRITE(*,113)KKK,MM,XXK,QXS,PXS,XCOR,PXB,TXB
      DO 444 I=1,NP
!
!  THE ESTIMATED MEAN, PXB, IS SUBTRACTED OFF AND THE COLUMN IS
!    ADJUSTED SO IT HAS THE ORIGINAL SUM OF SQUARES WITH MEAN ZERO
!
      PSIX(I,JJJ)=(XT(I,1)-PXB)*XCOR
      XT(I,1)=(XT(I,1)-PXB)*XCOR
  444 CONTINUE
!
!  END OF A.L.S. LOOP
!
!  IF THIS IS THE LAST ITERATION, I.E. JJJ = NF, THEN PUT A VECTOR
!    OF ONES INTO THE NF+1 COLUMN OF THE ESTIMATED P MATRIX
!    AND RESET X(NOBS,NY)--WHICH STARTED AS THE ORIGINAL DATA AND
!    THEN BECAME THE RESIDUAL MATRIX AT EACH STAGE--BACK TO THE
!    ORIGINAL DATA WHICH WAS SAVED IN XSS(NOBS,NY)
!
!  IF THIS IS NOT THE LAST ITERATION, PUT THE RESIDUALS X(NOBS,NY) INTO
!    XS(NOBS,NY), AND CALL CORR2 TO GET THE SIGN CHANGES FOR THE COLUMNS
!    WHICH ARE USED ABOVE
!
      DO 461 I=1,NP
      IF(KKK.EQ.NF)PSIX(I,NF+1)=1.0
      DO 461 J=1,NY
      IF(KKK.EQ.NF)X(I,J)=XSS(I,J)
  461 XS(I,J)=X(I,J)
      IF(KKK.EQ.NF)GO TO 9999
      call CORR2(NRESPONDENTS,NISSUES,NP,NY,X,R,LL,MPOS,KS,KPOS,1)
 9999 CONTINUE
!
!  WRITE OUT CORRELATION AND CROSS PRODUCT MATRIX OF P
!
!      WRITE(23,23)
!      call CORR2(NP,NF,PSIX,R,LL,MPOS,KS,KPOS,1)
!      WRITE(23,521)
!      call PSIPRM(NP,NF1,PSIX,IPRNT)
!
!  PERFORM A.L.S. WITH FULL P MATRIX UNTIL CONVERGENCE
!
      DO 222 NN=1,5
      call REG(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF,NY,TSUM,       &
                                      W,XS,X,PSIX,1,1,NF,AREG,BREG)
      call REG2(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF,NY,W,XS,      &
                                 X,PSIX,PXB,PXS,NF,IPRNT,AREG,BREG)
! 
!  CENTER THE PSI MATRIX AT ZERO
!
      DO 42 K=1,NF
      SUM=0.0
      DO 41 I=1,NP
      SUM=SUM+PSIX(I,K)
  41  CONTINUE
      SUM=SUM/FLOAT(NP)
      DO 43 I=1,NP
      PSIX(I,K)=PSIX(I,K)-SUM
  43  CONTINUE
  42  CONTINUE
      AA=ABS(AREG-BREG)
      IF(AA.LT..01)GO TO 223
  222 CONTINUE
  223 NFPS=NF+1
!
!  COMPUTE PW' + Jc' AND PUT INTO X(NOBS,NY)
!  FOR MISSING DATA, INSERT ESTIMATED VALUES IN X INTO THE 
!    ORIGINAL DATA MATRIX USING XT(NOBS,NY) SO THAT SINGULAR VALUES
!    CAN BE ESTRACTED FROM ORIGINAL DATA
!
      DO 60 I=1,NP
      DO 601 K=1,NY
  601 XT(I,K)=0.0
      DO 61 K=1,NY
      SUM=0.0
      DO 62 J=1,NFPS
  62  SUM=SUM+PSIX(I,J)*W(K,J)
  61  X(I,K)=SUM
      DO 59 JJ=1,NY
      IF(XS(I,JJ).NE.-999.0)XT(I,JJ)=XS(I,JJ)
      IF(XS(I,JJ).EQ.-999.0)XT(I,JJ)=X(I,JJ)
  59  CONTINUE
  60  CONTINUE
!
!  EXTRACT SINGULAR VALUES FROM ORIGINAL DATA MATRIX WITH ESTIMATES
!   INSERTED FOR MISSING DATA
!
      DO 71 I=1,NY
      DO 71 J=1,NY
      SUM=0.0
      DO 73 K=1,NP
  73  SUM=SUM+XT(K,I)*XT(K,J)
  71  ROOTC(I,J)=SUM
!      CALL RS(3539,NY,ROOTC,D,1,CROOT,FAL1,FAL2,IERR)
      DO 729 I=1,NP
      DO 729 J=1,NY
      ROOTC(I,J)=XT(I,J)
      CROOT(I,J)=XT(I,J)-DC(J)
  729 CONTINUE
      XTOL=.001
!      CALL LSVRR(NP,NY,ROOTC,3539,21,XTOL,IRANK,YHAT,UUU,            &
!                 3539,VVV,3539)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
!      LWORK=3539
      CALL DGESVD('S','S',NP,NY,ROOTC,NRESPONDENTS,YHAT,UUU,         &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
!      WRITE(23,337)INFO
!      WRITE(23,337)IRANK
      DO 464 I=1,NY
!      WRITE(23,338)I,YHAT(I)
      D(I)=YHAT(I)
  464 CONTINUE
      XTOL=.001
!      CALL LSVRR(NP,NY,CROOT,3539,21,XTOL,IRANK,YHAT,UUU,            &
!                 3539,VVV,3539)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
!      LWORK=3539
      CALL DGESVD('S','S',NP,NY,CROOT,NRESPONDENTS,YHAT,UUU,         &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
!      WRITE(23,337)INFO
!      WRITE(23,337)IRANK
      DO 465 I=1,NY
!      WRITE(23,338)I,YHAT(I)
      DX(I)=YHAT(I)
  465 CONTINUE
!      WRITE(23,7200)IERR
!      WRITE(*,7200)IERR
      ESUM=0.0
      SUMM=0.0
      SUMM1=0.0
      SUMM2=0.0
      DO 83 I=1,NP
      DO 83 J=1,NY
      SUMM=SUMM+X(I,J)**2
      IF(XS(I,J).NE.-999.0)ESUM=ESUM+(X(I,J)-XS(I,J))**2
      IF(XS(I,J).NE.-999.0)SUMM1=SUMM1+X(I,J)**2
      IF(XS(I,J).EQ.-999.0)SUMM2=SUMM2+X(I,J)**2
!
!  STORE PW' + Jc' IN XSS(NOBS,NY) AND COMPUTE PW' AND STORE IN XT(NOBS,NY)
!
      XSS(I,J)=X(I,J)
  83  XT(I,J)=X(I,J)-W(J,NF+1)
!
!  PERFORM SINGULAR VALUE DECOMPOSITION OF PW' + Jc'
!
      DO 75 I=1,NY
      DO 75 J=1,NY
      SUM=0.0
      DO 76 K=1,NP
      SUM=SUM+XSS(K,I)*XSS(K,J)
  76  CONTINUE
  75  ROOTC(I,J)=SUM
!      CALL RS(3539,NY,ROOTC,DDD,1,CROOT,FAL1,FAL2,IERR)
      DO 728 I=1,NP
      DO 728 J=1,NY
      ROOTC(I,J)=XSS(I,J)
  728 CONTINUE
      XTOL=.001
!      CALL LSVRR(NP,NY,ROOTC,3539,21,XTOL,IRANK,YHAT,UUU,            &
!                 3539,VVV,3539)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
!      LWORK=3539
      CALL DGESVD('S','S',NP,NY,ROOTC,NRESPONDENTS,YHAT,UUU,         &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
!      WRITE(23,337)INFO
!      WRITE(23,337)IRANK
      DO 463 I=1,NY
!      WRITE(23,338)I,YHAT(I)
      DDD(I)=YHAT(I)
  463 CONTINUE
!      WRITE(23,7200)IERR
!      WRITE(*,7200)IERR
!
!  PERFORM SINGULAR VALUE DECOMPOSITION OF PW'
!
      DO 77 I=1,NY
      DO 77 J=1,NY
      SUM=0.0
      DO 78 K=1,NP
      SUM=SUM+XT(K,I)*XT(K,J)
  78  CONTINUE
  77  ROOTC(I,J)=SUM
!      CALL RS(3539,NY,ROOTC,DD,1,CROOT,FAL1,FAL2,IERR)
      DO 727 I=1,NP
      DO 727 J=1,NY
      ROOTC(I,J)=XT(I,J)
  727 CONTINUE
      XTOL=.001
!      CALL LSVRR(NP,NY,ROOTC,3539,21,XTOL,IRANK,YHAT,UUU,            &
!                 3539,VVV,3539)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
!      LWORK=3539
      CALL DGESVD('S','S',NP,NY,ROOTC,NRESPONDENTS,YHAT,UUU,         &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
!      WRITE(23,337)INFO
!      WRITE(23,337)IRANK
      DO 462 I=1,NP
      IF(I.LE.NY)THEN
         DD(I)=YHAT(I)
!         WRITE(23,338)I,YHAT(I)
      ENDIF
      DO 462 JJ=1,NF
      PSIX(I,JJ)=UUU(I,JJ)*SQRT(YHAT(JJ))
  462 CONTINUE
      DO 460 I=1,NY
      DO 460 JJ=1,NF
!      W(I,JJ)=VVV(I,JJ)*SQRT(YHAT(JJ))
      W(I,JJ)=VVV(JJ,I)*SQRT(YHAT(JJ))
  460 CONTINUE
!      WRITE(23,7200)IERR
!      WRITE(*,7200)IERR
!  
!  CONSTRAINT CHECKS
!
      DO 233 K=1,NFPS
      SUM=0.0
      DO 234 I=1,NP
      SUM=SUM+PSIX(I,K)
  234 CONTINUE
      VVV(1,K)=SUM/FLOAT(NP)
  233 CONTINUE
!      WRITE(23,199)(VVV(1,K),K=1,NF)
      DO 227 J=1,NF
      DO 227 K=1,NF
      SUM=0.0
      DO 228 I=1,NP
      SUM=SUM+PSIX(I,J)*PSIX(I,K)
  228 CONTINUE
  227 VVV(J,K)=SUM
      DO 229 I=1,NF
!      WRITE(23,200)I,(VVV(I,J),J=1,NF)
  229 CONTINUE
      DO 230 J=1,NF
      DO 230 K=1,NF
      SUM=0.0
      DO 231 I=1,NY
      SUM=SUM+W(I,J)*W(I,K)
  231 CONTINUE
  230 VVV(J,K)=SUM
      DO 232 I=1,NF
!      WRITE(23,200)I,(VVV(I,J),J=1,NF)
  232 CONTINUE
!
!  SAVE PSI AND W HERE
!
      DO 266 J=1,NF
      DO 266 I=1,NP
      XSAVE(I,J)=PSIX(I,J)
  266 CONTINUE
      DO 265 J=1,NF
      DO 265 I=1,NY
      WSAVE(I,J+1)=W(I,J)
      WSAVE(I,1)=W(I,NF+1)
  265 CONTINUE        
!
!
!  WRITE OUT SINGULAR VALUES FOR ORIGINAL DATA (USING ESTIMATES FOR
!   MISSING ENTRIES), PW' + Jc', AND PW'
!
      IF(IPRNT.EQ.0)WRITE(23,86)
      DO 87 J=1,NY
      IF(J.GT.NF+3)GO TO 87
      IF(IPRNT.EQ.0)WRITE(23,1113)J,D(J),DDD(J),DX(J),DD(J)
!      WRITE(*,1113)J,D(J),DDD(J),DX(J),DD(J)
  87  CONTINUE
!
!
      IF(IPRNT.EQ.0)WRITE(23,2222)
!
!  RESTORE PSI AND W HERE
!
      DO 264 J=1,NF
      DO 264 I=1,NP
      XDATA(I,J)=XSAVE(I,J)
  264 CONTINUE
      DO 263 J=1,NF+1
      DO 263 I=1,NY
      W(I,J)=WSAVE(I,J)
  263 CONTINUE        
!
      DEALLOCATE(LL)
      DEALLOCATE(MPOS)
      DEALLOCATE(XX)
      DEALLOCATE(D)
      DEALLOCATE(DD)
      DEALLOCATE(DC)
      DEALLOCATE(CC)
      DEALLOCATE(TSUM)
      DEALLOCATE(DDD)
      DEALLOCATE(DX)
      DEALLOCATE(YHAT)
      DEALLOCATE(WORK)
!
      DEALLOCATE(X)
      DEALLOCATE(PSIX)
      DEALLOCATE(XS)
      DEALLOCATE(R)
      DEALLOCATE(ROOTC)
      DEALLOCATE(CROOT)
      DEALLOCATE(XT)
      DEALLOCATE(XSS)
      DEALLOCATE(UUU)
      DEALLOCATE(VVV)
      DEALLOCATE(WSAVE)
      DEALLOCATE(XSAVE)
!
      RETURN
      END
!
!  **********************************************************************
!    SUBROUTINE CORR2---CALLED BY BLACKB.  COMPUTES CORRELATION MATRIX
!     FOR INPUT MATRIX X AND COMPUTES THE VECTOR OF COLUMN SIGN CHANGES
!     WHICH IS STORED IN VECTOR LL.
!  **********************************************************************
!
      SUBROUTINE CORR2(NRESPONDENTS,NISSUES,NP,NY,X,R,LL,MPOS,KS, &
                                                        KPOS,IPRNT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(NRESPONDENTS,NISSUES),R(NISSUES,NISSUES),       &
                LL(NISSUES),MPOS(NISSUES)
!
      DOUBLE PRECISION, ALLOCATABLE :: SA(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: SB(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: SC(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: SD(:,:)
!
      ALLOCATE(SA(NISSUES,NISSUES))
      ALLOCATE(SB(NISSUES,NISSUES))
      ALLOCATE(SC(NISSUES,NISSUES))
      ALLOCATE(SD(NISSUES,NISSUES))

  200 FORMAT(1X,50F7.3)
      DO 1 J=1,NY
      DO 1 JJ=1,NY
      SA(J,JJ)=0.0
      SB(J,JJ)=0.0
      SC(J,JJ)=0.0
  1   SD(J,JJ)=0.0
      DO 40 I=1,NP
      DO 31 J=1,NY
      DO 31 JJ=1,J
      IF(X(I,J).EQ.-999.0)GO TO 31
      IF(X(I,JJ).EQ.-999.0)GO TO 31
      SA(J,JJ)=SA(J,JJ)+X(I,J)
      IF(J.NE.JJ)SA(JJ,J)=SA(JJ,J)+X(I,JJ)
      SB(J,JJ)=SB(J,JJ)+X(I,J)*X(I,J)
      IF(J.NE.JJ)SB(JJ,J)=SB(JJ,J)+X(I,JJ)*X(I,JJ)
      SC(J,JJ)=SC(J,JJ)+X(I,J)*X(I,JJ)
      SC(JJ,J)=SC(J,JJ)
      SD(J,JJ)=SD(J,JJ)+1.0
  31  CONTINUE
  40  CONTINUE
      DO 32 J=1,NY
      DO 32 JJ=1,J
      AA=SD(J,JJ)*SC(J,JJ)-SA(J,JJ)*SA(JJ,J)
      BB=SD(J,JJ)*SB(J,JJ)-SA(J,JJ)*SA(J,JJ)
      CC=SD(J,JJ)*SB(JJ,J)-SA(JJ,J)*SA(JJ,J)
      BBCC=BB*CC
      IF(BBCC.LE..0)R(JJ,J)=0.0
      IF(BBCC.LE..0)GO TO 343
      R(JJ,J)=AA/SQRT(BB*CC)
  343 CONTINUE
  32  R(J,JJ)=R(JJ,J)
      IF(IPRNT.EQ.1)GO TO 62
      DO 34 J=1,NY
  34  WRITE(23,200)(R(J,JJ),JJ=1,NY)
  62  BB=-99.0
!
!  FIND ROW WITH LARGEST TOTAL SUM OF ABSOLUTE VALUED CORRELATIONS
!
      KS=0
      DO 50 J=1,NY
      SUM=0.0
      DO 51 JJ=1,NY
  51  SUM=SUM+ABS(R(J,JJ))
      IF(SUM.GT.BB)THEN
         BB=SUM
         KS=J
      ENDIF
!      BB=AMAX1(SUM,BB)
!      IF(BB.EQ.SUM)KS=J
  50  CONTINUE
      KSOLD=KS
      DO 52 J=1,NY
      IF(R(KS,J).LE.0.0)LL(J)=-1
  52  IF(R(KS,J).GT.0.0)LL(J)=1
!
!  ITERATIVELY CHANGE SIGNS OF COLUMNS TO MAXIMIZE NUMBER OF ENTRIES
!   IN THE CORRELATION MATRIX WHICH ARE POSITIVE.  THIS VECTOR IS THEN
!   RETURNED AS THE SIGN CHANGES FOR THE DATA MATRIX X.
!
      NYD2=(NY-1)/2
      KPOS=0
      DO 60 JK=1,NY
      DO 60 J=1,NY
      KK=0
      KSUM=0
      DO 61 JJ=1,NY
      AA=R(J,JJ)*FLOAT(LL(JJ))*FLOAT(LL(J))
      IF(JK.EQ.NY.AND.AA.GE.0.0)KPOS=KPOS+1
      IF(JK.EQ.NY.AND.AA.GE.0.0)KSUM=KSUM+1
  61  IF(AA.LT.0.0)KK=KK+1
      IF(KK.GT.NYD2)LL(J)=LL(J)*(-1)
      IF(JK.EQ.NY)MPOS(J)=KSUM
  60  IF(KK.GT.NYD2)KS=999
      IF(IPRNT.EQ.0)WRITE(23,1009)NY,NYD2,KPOS,KSUM,KS,KSOLD,        &
                                         (LL(JCJ),JCJ=1,NY)
 1009 FORMAT(' SIGN VECTOR',40I4)
!
      DEALLOCATE(SA)
      DEALLOCATE(SB)
      DEALLOCATE(SC)
      DEALLOCATE(SD)
!
      RETURN
      END
!
!  ************************************************************************
!    SUBROUTINE REG---CALLED BY BLACKB.  PERFORMS REGRESSION TO ESTIMATE
!      W AND c.
!  ************************************************************************
!
      SUBROUTINE REG(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF,NY,      &
                           TSUM,W,XS,X,PSI,IPRNT,ILAST,KKK,AREG,BREG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TSUM(2*NISSUES),W(NISSUES,NDIMENSIONS+2),            &
                XS(NRESPONDENTS,NISSUES),X(NRESPONDENTS,NISSUES),    &
                PSI(NRESPONDENTS,NISSUES)
!            A(127,127),B(127,127),    &
!                C(152),LLL(152),    &
!                XT(3539,152),R(152,152),LL(152),RSUM(152),  &
!                MPOS(152),WVEC(127),ZMAT(127,127),WK(527)
!
      INTEGER, ALLOCATABLE :: LLL(:)
      INTEGER, ALLOCATABLE :: LL(:)
      INTEGER, ALLOCATABLE :: MPOS(:)
      DOUBLE PRECISION, ALLOCATABLE :: C(:)
      DOUBLE PRECISION, ALLOCATABLE :: RSUM(:)
      DOUBLE PRECISION, ALLOCATABLE :: WVEC(:)
      DOUBLE PRECISION, ALLOCATABLE :: WK(:)
      DOUBLE PRECISION, ALLOCATABLE :: A(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: B(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XT(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: R(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ZMAT(:,:)
! 
      ALLOCATE(LLL(NISSUES))
      ALLOCATE(LL(NISSUES))
      ALLOCATE(MPOS(NISSUES))
      ALLOCATE(C(NISSUES))
      ALLOCATE(RSUM(NISSUES))
      ALLOCATE(WVEC(3*NDIMENSIONS))
      ALLOCATE(WK(4*NDIMENSIONS))
      ALLOCATE(A(3*NDIMENSIONS,3*NDIMENSIONS))
      ALLOCATE(B(3*NDIMENSIONS,3*NDIMENSIONS))
      ALLOCATE(XT(NRESPONDENTS,NISSUES))
      ALLOCATE(R(NISSUES,NISSUES))
      ALLOCATE(ZMAT(3*NDIMENSIONS,3*NDIMENSIONS))
!
  100 FORMAT(' PERFORMANCE INDEX MATRIX INVERSE=',2I5)
  101 FORMAT(' PERFORMANCE INDEX EIGENVALUE/VECTOR ROUTINE=',4I5)
  102 FORMAT(3I5,50F13.5)
  111 FORMAT(1X,8A1,13F10.4)
  200 FORMAT(3I5,15F10.4/15x,15f10.4/15x,15f10.4)
  250 FORMAT(' DIMENSION=',I3,1X,'TOTAL SSE REG1=',F15.4)
  300 FORMAT(' SSE AND FINAL WEIGHT MATRIX')
      NF1=NF+1
      LWORK=4*NDIMENSIONS
      DO 15 K=1,NY
      LLL(K)=0
      TSUM(K)=0.0
      C(K)=0.0
      DO 15 J=1,NF1
  15  W(K,J)=0.0
!
!  IF ILAST = 1 THEN ONLY PERFORM THE ESTIMATE OF W AND c FOR THE 
!   CURRENT COLUMN OF P.
!  IF ILAST = 0 THEN ESTIMATE W AND c ASSUMING NF=1, THEN ASSUMING
!   NF=2, UP TO NF=NF.  THIS ALLOWS THE COMPUTATION OF THE SS ACCOUNTED
!   FOR BY A 1 DIMENSION PW' + Jc' MODEL, A 2 DIMENSION PW' + Jc' MODEL,
!   AND SO ON.
!
      IF(ILAST.EQ.1)KKA=NF
      IF(ILAST.EQ.0)KKA=1
      DO 11 KK=KKA,NF
      DO 16 K=1,NY
      DO 16 J=1,NF+1
  16  W(K,J)=0.0
!
!  RUN REGRESSION TO ESTIMATE W AND c ONE COLUMN AT A TIME
!   BECAUSE OF MISSING DATA
!
      NF1=KK+1
      NF11=KK+2
      DO 7 K=1,NY
      LLL(K)=0
!
!  COMPUTE [P'P]-1
!
      DO 2 J=1,NF1
      DO 2 JJ=1,NF1
      SUM=0.0
      DO 1 I=1,NP
      IF(XS(I,K).EQ.-999.0)GO TO 1
      SUM=SUM+PSI(I,J)*PSI(I,JJ)
  1   CONTINUE
      A(J,JJ)=SUM
  2   ZMAT(J,JJ)=SUM
      CALL DSYEV('V','L',NF1,ZMAT,3*NDIMENSIONS,WVEC,WK,LWORK,INFO)
!      call rs(127,nf1,a,wvec,1,ZMAT,fv1,fv2,ier)
      DO 60 I=1,NF1
      DO 60 KX=1,NF1
      SUM=0.0
      DO 61 J=1,NF1
      IF(ABS(WVEC(J)).GT..0001)THEN
          SUM=SUM+ZMAT(KX,J)*(1.0/WVEC(J))*ZMAT(I,J)
      ENDIF
  61  CONTINUE
  60  B(I,KX)=SUM
!
!  COMPUTE [P'P]-1*P'X = W'/c'
!
      DO 3 I=1,NP
      IF(XS(I,K).EQ.-999.0)GO TO 3
      DO 5 J=1,NF1
      SUM=0.0
      DO 4 JJ=1,NF1
  4   SUM=SUM+B(J,JJ)*PSI(I,JJ)
  5   C(J)=SUM
      DO 6 J=1,NF1
  6   W(K,J)=W(K,J)+C(J)*XS(I,K)
  3   CONTINUE
!
!  CALCULATE R-SQUARE AND SSE FOR REGRESSION
!
      ESUM=0.0
      DO 8 I=1,NP
      IF(XS(I,K).EQ.-999.0)GO TO 8
      LLL(K)=LLL(K)+1
      SUM=0.0
      DO 9 J=1,NF1
  9   SUM=SUM+PSI(I,J)*W(K,J)
      ESUM=ESUM+(SUM-XS(I,K))**2
      X(I,K)=SUM-XS(I,K)
      XT(I,1)=SUM
      XT(I,2)=XS(I,K)
  8   CONTINUE
      W(K,NF11)=ESUM
      TSUM(KK)=TSUM(KK)+ESUM
      call CORR2(NRESPONDENTS,NISSUES,NP,2,XT,R,LL,MPOS,KS,KPOS,1)
      RSUM(K)=R(1,2)*R(1,2)
  7   CONTINUE
!
!  END OF REGRESSION LOOP FOR CURRENT DIMENSION
!
      IF(ILAST.EQ.1)GO TO 90
!
!  IF ILAST = 0 THEN CALCULATE THE R-SQUARE FOR THE CURRENT
!   NUMBER OF DIMENSIONS
!
      DO 888 I=1,NP
      DO 88 K=1,NY
      SUM=0.0
      DO 89 J=1,NF1
  89  SUM=SUM+PSI(I,J)*W(K,J)
  88  XT(I,K)=SUM
  888 CONTINUE
      call RSQUR(NRESPONDENTS,NISSUES,NP,NY,RR,XT,XS,1)
      TSUM(KK+NY)=RR
  90  CONTINUE
      IF(IPRNT.EQ.1)GO TO 333
!      WRITE(23,300)
!      DO 10 J=1,NY
!  10  WRITE(23,200)KK,J,LLL(J),W(J,NF11),RSUM(J),                    &
!                          (W(J,JJ),JJ=1,NF1)
  333 CONTINUE
      IF(ILAST.EQ.1)THEN
         IF(IPRNT.EQ.0)WRITE(23,250)KKK,TSUM(KK)
!         WRITE(*,250)KKK,TSUM(KK)
      ENDIF
      IF(ILAST.EQ.0)THEN
         IF(IPRNT.EQ.0)WRITE(23,250)KK,TSUM(KK)
!         WRITE(*,250)KK,TSUM(KK)
      ENDIF
      AREG=TSUM(KK)
  11  CONTINUE
!
!  END OF DIMENSION ITERATION LOOP--ONLY IF ILAST = 0
!
      DEALLOCATE(LLL)
      DEALLOCATE(LL)
      DEALLOCATE(MPOS)
      DEALLOCATE(C)
      DEALLOCATE(RSUM)
      DEALLOCATE(WVEC)
      DEALLOCATE(WK)
      DEALLOCATE(A)
      DEALLOCATE(B)
      DEALLOCATE(XT)
      DEALLOCATE(R)
      DEALLOCATE(ZMAT)
!
      RETURN
      END
!
!  **************************************************************************
!    SUBROUTINE REG2---CALLED BY BLACKB.  PERFORMS REGRESSION TO ESTIMATE P.
!  **************************************************************************
!
      SUBROUTINE REG2(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF,NY,W,XS,&
                      X,PSI,PXB,PXS,KKK,IPRNT,AREG,BREG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION W(NISSUES,NDIMENSIONS+2),                            &
                XS(NRESPONDENTS,NISSUES),X(NRESPONDENTS,NISSUES),    &
                PSI(NRESPONDENTS,NISSUES)
!                A(152,127),V(152),Y(152)
      DOUBLE PRECISION, ALLOCATABLE :: V(:)
      DOUBLE PRECISION, ALLOCATABLE :: Y(:)
      DOUBLE PRECISION, ALLOCATABLE :: A(:,:)
!
      ALLOCATE(V(NISSUES))
      ALLOCATE(Y(NISSUES))
      ALLOCATE(A(NISSUES,NDIMENSIONS+2))
!
      NF1=NF+1
      ESUM=0.0
      PXB=0.0
      PXS=0.0
      XNS=0.0
!
!  REGRESS ONE ROW AT A TIME BECAUSE OF MISSING DATA
!
      DO 1 I=1,NP
      KK=0
!
!  SET UP DEPENDENT VARIABLE-- Y(NY) = X[I,NY] - NY COL. MEANS
!    AND PUT W(NY,NF) INTO MATRIX A(NY,NF)
!
      DO 2 J=1,NY
      V(J)=0.0
      IF(XS(I,J).EQ.-999.0)GO TO 2
      KK=KK+1
      Y(KK)=XS(I,J)-W(J,NF1)
      DO 3 JJ=1,NF
  3   A(KK,JJ)=W(J,JJ)
  2   CONTINUE
!
!  CALL SUBROUTINE REGA TO PERFORM LEAST SQUARES
!
      NS=KK
      call REGA(NISSUES,NDIMENSIONS,NS,NF,A,Y,V,1)
!
!  STORE ESTIMATED ROW ENTRIES OF P AND COMPUTE SSE
!
      DO 5 K=1,NY
      SUM=0.0
      IF(XS(I,K).EQ.-999.0)GO TO 5
      DO 4 J=1,NF
      PSI(I,J)=V(J)
      SUM=SUM+PSI(I,J)*W(K,J)
  4   CONTINUE
      SUM=SUM+W(K,NF1)
      X(I,K)=SUM-XS(I,K)
      ESUM=ESUM+(SUM-XS(I,K))**2
  5   CONTINUE
      PXB=PXB+PSI(I,1)
      PXS=PXS+PSI(I,1)**2
      XNS=XNS+1.0
  1   CONTINUE
!
!  END OF REGRESSION LOOP
!
!  PXB IS THE SUM OF THE ESTIMATED COLUMN OF P AND
!  PXS IS THE VARIANCE OF THE ESTIMATED COLUMN OF P.
!  THESE ARE RETURNED TO BLACKB FOR ADJUSTMENT PURPOSES
!
      PXB=PXB/XNS
      PXS=PXS-XNS*PXB*PXB
!      SQU(KKK)=ESUM
      IF(IPRNT.EQ.0)WRITE(23,350)KKK,ESUM
!      IF(IPRNT.EQ.0)WRITE(*,350)KKK,ESUM
  350 FORMAT(' DIMENSION=',I3,1X,'TOTAL SSE REG2=',F15.4)
      BREG=ESUM
!
      DEALLOCATE(V)
      DEALLOCATE(Y)
      DEALLOCATE(A)
!
      RETURN
      END
!
!  **************************************************************************
!    SUBROUTINE REGA---CALLED BY REG2.  PERFORMS THE REGRESSION:
!
!         [W'W]-1*W'[X(I,NY) - c']
!
!      TO GET THE ROW ENTRIES OF P.
!  **************************************************************************
!
      SUBROUTINE REGA(NISSUES,NDIMENSIONS,NS,NF,A,Y,V,IPRNT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NISSUES,NDIMENSIONS+2),Y(NISSUES),V(NISSUES)
!
!             B(127,127),C(127,127), &
!               BB(127,152),ZMAT(127,127),WVEC(127),&
!               WK(527)
!
      DOUBLE PRECISION, ALLOCATABLE :: WVEC(:)
      DOUBLE PRECISION, ALLOCATABLE :: WK(:)
      DOUBLE PRECISION, ALLOCATABLE :: B(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: C(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: BB(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ZMAT(:,:)
!
      ALLOCATE(WVEC(3*NDIMENSIONS))
      ALLOCATE(WK(3*NDIMENSIONS))
      ALLOCATE(B(3*NDIMENSIONS,3*NDIMENSIONS))
      ALLOCATE(C(3*NDIMENSIONS,3*NDIMENSIONS))
      ALLOCATE(BB(3*NDIMENSIONS,NISSUES))
      ALLOCATE(ZMAT(3*NDIMENSIONS,3*NDIMENSIONS))
!
  101 FORMAT(' PERFORMANCE INDEX EIGENVALUE/VECTOR ROUTINE=',3I5)
  102 FORMAT(2I5,30F10.4)
      LWORK=3*NDIMENSIONS
      DO 1 J=1,NF
      DO 1 JJ=1,NF
      SUM=0.0
      DO 2 I=1,NS
  2   SUM=SUM+A(I,J)*A(I,JJ)
      B(J,JJ)=SUM
  1   ZMAT(J,JJ)=SUM
!
      CALL DSYEV('V','L', NF, ZMAT,3*NDIMENSIONS,WVEC,WK,LWORK,INFO)
!      call rs(127,nf,B,wvec,1,ZMAT,fv1,fv2,ier)
      DO 60 I=1,NF
      DO 60 K=1,NF
      SUM=0.0
      DO 61 J=1,NF
      IF(ABS(WVEC(J)).GT..0001)THEN
          SUM=SUM+ZMAT(K,J)*(1.0/WVEC(J))*ZMAT(I,J)
      ENDIF
  61  CONTINUE
  60  C(I,K)=SUM
      DO 3 I=1,NS
      DO 3 J=1,NF
      SUM=0.0
      DO 4 JJ=1,NF
  4   SUM=SUM+C(J,JJ)*A(I,JJ)
  3   BB(J,I)=SUM
      DO 5 JJ=1,NF
      SUM=0.0
      DO 6 J=1,NS
  6   SUM=SUM+BB(JJ,J)*Y(J)
  5   V(JJ)=SUM
!
      DEALLOCATE(WVEC)
      DEALLOCATE(WK)
      DEALLOCATE(B)
      DEALLOCATE(C)
      DEALLOCATE(BB)
      DEALLOCATE(ZMAT)
!
      RETURN
      END
!
!  ***************************************************************************
!    SUBROUTINE RSQUR---CALLED BY REG.  COMPUTES THE SIMPLE PEARSON R-SQUARE
!     BETWEEN TWO INPUT MATRICES---A(NOBS,NY) AND B(NOBS,NY)
!  ***************************************************************************
!
      SUBROUTINE RSQUR(NRESPONDENTS,NISSUES,NP,NY,R,A,B,IPRNT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NRESPONDENTS,NISSUES),B(NRESPONDENTS,NISSUES)
      ASUM=0.0
      BSUM=0.0
      CSUM=0.0
      DSUM=0.0
      ESUM=0.0
      XNT=0.0
      DO 1 I=1,NP
      DO 1 J=1,NY
      IF(A(I,J).EQ.-999.0)GO TO 1
      IF(B(I,J).EQ.-999.0)GO TO 1
      ASUM=ASUM+A(I,J)
      BSUM=BSUM+B(I,J)
      CSUM=CSUM+A(I,J)*A(I,J)
      DSUM=DSUM+B(I,J)*B(I,J)
      ESUM=ESUM+A(I,J)*B(I,J)
      XNT=XNT+1.0
  1   CONTINUE
      AA=XNT*ESUM-ASUM*BSUM
      BB=XNT*CSUM-ASUM*ASUM
      CC=XNT*DSUM-BSUM*BSUM
      R=(AA*AA)/(BB*CC)
      RETURN
      END
!
!  ***************************************************************************
!    SUBROUTINE PSIPRM---CALLED BY BLACKB.  CALCULATES P'P AND WRITES IT OUT.
!  ***************************************************************************
!
      SUBROUTINE PSIPRM(NP,NF,PSI,IPRNT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION PSI(3539,152),R(127,127)
  200 FORMAT(1X,50F10.4)
      DO 1 J=1,NF
      DO 1 K=1,NF
      SUM=0.0
      DO 2 I=1,NP
      SUM=SUM+PSI(I,J)*PSI(I,K)
  2   CONTINUE
  1   R(J,K)=SUM
      IF(IPRNT.EQ.0)THEN
         DO 3 I=1,NF
  3      WRITE(23,200)(R(I,J),J=1,NF)
      ENDIF
      RETURN
      END

!
!
!  ************************************************************************
!    SORT SUBROUTINE--SORTS A VECTOR 'A' OF REAL ELEMENTS INTO ASCENDING
!    ORDER.  'LA' IS THE NUMBER OF ELEMENTS TO BE SORTED AND 'IR' IS A
!    VECTOR OF INTEGERS THAT RECORDS THE PERMUTATIONS--USUALLY SET TO
!    1,2,3,4,...
!  ************************************************************************
!
!
      SUBROUTINE RSORT(A,LA,IR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LA),IU(21),IL(21),IR(LA)
      IF (LA.LE.0) RETURN
      M = 1
      I = 1
      J = LA
      R = .375
    5 IF (I.EQ.J) GO TO 45
      IF (R.GT..5898437) GO TO 10
      R = R+3.90625E-2
      GO TO 15 
   10 R = R-.21875 
   15 K = I 
!
! SELECT A CENTRAL ELEMENT OF THE  
! ARRAY AND SAVE IT IN LOCATION T  
!
      IJ = I+(J-I)*R 
      T = A(IJ)   
      IT = IR(IJ) 
!
! FIRST ELEMENT OF ARRAY IS GREATER
! THAN T, INTERCHANGE WITH T       
!
      IF (A(I).LE.T) GO TO 20  
      A(IJ) = A(I) 
      A(I) = T     
      T = A(IJ) 
      IR(IJ) = IR(I)  
      IR(I) = IT
      IT = IR(IJ) 
   20 L = J
!
! IF LAST ELEMENT OF ARRAY IS LESS THAN
! T, INTERCHANGE WITH T
!
      IF (A(J).GE.T) GO TO 30
      A(IJ) = A(J)
      A(J) = T
      T = A(IJ)
      IR(IJ) = IR(J)
      IR(J) = IT
      IT = IR(IJ)
!
! IF FIRST ELEMENT OF ARRAY IS GREATER
! THAN T, INTERCHANGE WITH T
!
      IF (A(I).LE.T) GO TO 30
      A(IJ) = A(I)
      A(I) = T
      T = A(IJ)
      IR(IJ) = IR(I)
      IR(I) = IT
      IT = IR(IJ)
      GO TO 30
   25 IF (A(L).EQ.A(K)) GO TO 30
      TT = A(L)
      A(L) = A(K)
      A(K) = TT
      ITT = IR(L)
      IR(L) = IR(K)
      IR(K) = ITT
!
! FIND AN ELEMENT IN THE SECOND HALF OF
! THE ARRAY WHICH IS SMALLER THAN T
!
   30 L = L-1
      IF (A(L).GT.T) GO TO 30
!
! FIND AN ELEMENT IN THE FIRST HALF OF
! THE ARRAY WHICH IS GREATER THAN T
!
   35 K = K+1
      IF (A(K).LT.T) GO TO 35
!
! INTERCHANGE THESE ELEMENTS
!
      IF (K.LE.L) GO TO 25
!
! SAVE UPPER AND LOWER SUBSCRIPTS OF
! THE ARRAY YET TO BE SORTED
!
      IF (L-I.LE.J-K) GO TO 40
      IL(M) = I
      IU(M) = L
      I = K
      M = M+1
      GO TO 50
   40 IL(M) = K
      IU(M) = J
      J = L
      M = M+1
      GO TO 50
!
! BEGIN AGAIN ON ANOTHER PORTION OF
! THE UNSORTED ARRAY
!
   45 M = M-1
      IF (M.EQ.0) RETURN
      I = IL(M)
      J = IU(M)
   50 IF (J-I.GE.11) GO TO 15
      IF (I.EQ.1) GO TO 5
      I = I-1
   55 I = I+1
      IF (I.EQ.J) GO TO 45
      T = A(I+1)
      IT = IR(I+1)
      IF (A(I).LE.T) GO TO 55
      K = I
   60 A(K+1) = A(K)
      IR(K+1) = IR(K)
      K = K-1
      IF (T.LT.A(K)) GO TO 60
      A(K+1) = T
      IR(K+1) = IT
      GO TO 55
      END
