!
!
      INTEGER NRESPONDENTS,NISSUES,NDIMENSIONS,NMISSING,POLARITY,    &
                EXITSTATUS,KMISS(6*3),MINSCALE,KISSUE(1614*6),       &
                MID(1614),LRESPONDENTS(6+1614),LMARK(1614)
      DOUBLE PRECISION FITS(5*7),PSIMATRIX(1614*((5*6)/2)+2*1614*5), &
                       WMATRIX(6*((5*6)/2)+2*6*5),FITS2(6),          &
                       STIMCOORDS(5),EIGENVALUES(5)
      CHARACTER*21 CAND(6)
!
!
!  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
!  
!    IBLACKWHAT=1  ALDRICH-MCKELVEY SCALING
!    IBLACKWHAT=2  BLACKBOX
!    IBLACKWHAT=3  BLACKBOX_TRANSPOSE
!
      OPEN(22,FILE='C:\INETPUB\FTPROOT\WF1\NES1980.DAT',STATUS='OLD')
      OPEN(23,FILE='BLACKT23F95.DAT')
      OPEN(24,FILE='BLACKT24F95.DAT')
      OPEN(28,FILE='BLACKT28F95.DAT')
!
!  FORMAT FOR 1980 DATA
!
  100 FORMAT(8X,I4,1810X,4I1,11X,2I1)
  113 FORMAT(1X,I5,8F15.4)
  201 FORMAT(2I7,35F7.3)
  225 FORMAT(1X,A21,I7,35F7.3)
  300 FORMAT(1X,'NUMBER OF ROWS               =',F8.1/               &
             1X,'NUMBER OF COLUMNS            =',F8.1/               &
             1X,'TOTAL NUMBER OF DATA ENTRIES =',F8.1/               &
             1X,'NUMBER MISSING ENTRIES       =',F8.1/               &
             1X,'PERCENT MISSING DATA         =',F15.5/              &
             1X,'SUM OF SQUARES GRAND MEAN    =',F15.5)              
  301 FORMAT(I5,2F10.4,1X,A21)
  302 FORMAT(2I6,4F12.4)
  303 FORMAT(' FITS',5F10.4)
!
! 1980 LIBERAL-CONSERVATIVE SCALE EXAMPLE
!
      NRESPONDENTS=1614
      NDIMENSIONS=3
      NISSUES=6
      NMISSING=3
!
      MINSCALE=NDIMENSIONS+2
      EXITSTATUS=0
!
      IBLACKWHAT=3
!
      DO 1 I=1,NISSUES
      DO 2 J=1,NMISSING
      KMISS(J+(I-1)*NMISSING)=0
  2   CONTINUE
      KMISS(2+(I-1)*NMISSING)=8
      KMISS(3+(I-1)*NMISSING)=9
  1   CONTINUE
!
!  1980 7-POINT SCALE
!
      CAND(1)= 'JIMMY CARTER         '
      CAND(2)= 'RONALD REAGAN        '
      CAND(3)= 'TED KENNEDY          '
      CAND(4)= 'JOHN ANDERSON        '
      CAND(5)= 'REPUBLICAN PARTY     '
      CAND(6)= 'DEMOCRATIC PARTY     '
!
      DO 111 I=1,NRESPONDENTS
      READ(22,100)MID(I),(KISSUE(J+(I-1)*NISSUES),J=1,NISSUES)
  111 CONTINUE
!
      IF(IBLACKWHAT.EQ.1)THEN
         CALL MCKALNEW(NRESPONDENTS,NISSUES,NSELFPOS,NMISSING,      &
           KMISS,POLARITY,MID,KISSUE,CAND,FITS,PSIMATRIX,STIMCOORDS,&
           EIGENVALUES,EXITSTATUS)
      ENDIF
!
      IF(IBLACKWHAT.EQ.2)THEN
         CALL BLACKBOX(NRESPONDENTS,NISSUES,NDIMENSIONS,NMISSING,    &
           KMISS,MINSCALE,MID,KISSUE,CAND,FITS,PSIMATRIX,WMATRIX,    &
           LRESPONDENTS,LMARK,FITS2,EXITSTATUS)
      ENDIF
!
      IF(IBLACKWHAT.EQ.3)THEN
         CALL BLACKBOXT(NRESPONDENTS,NISSUES,NDIMENSIONS,NMISSING,   &
           KMISS,MINSCALE,MID,KISSUE,CAND,FITS,PSIMATRIX,WMATRIX,    &
           LRESPONDENTS,LMARK,FITS2,EXITSTATUS)
      ENDIF
!
!
!  WRITE OUT ESTIMATES AS A CHECK
!
!  CHECK OUTPUT
!
!
! &&&&&&&&&&&&&&&&&&&  ALDRICH-MCKELVEY OUTPUT  &&&&&&&&&&&&&&&&&&&&&&&&&&&
!
      IF(IBLACKWHAT.EQ.1)THEN
!
         WRITE(23,303)(FITS(J),J=1,5)
         DO 52 I=1,NISSUES-1
         WRITE(23,301)I,EIGENVALUES(I),STIMCOORDS(I),CAND(I)
         WRITE(*,301)I,EIGENVALUES(I),STIMCOORDS(I),CAND(I)
  52     CONTINUE
         DO 53 I=1,NRESPONDENTS
         KK=(I-1)*4
         WRITE(24,302)I,MID(I),(PSIMATRIX(KK+J),J=1,4)
  53     CONTINUE
      ENDIF
!
! &&&&&&&&&&&&&&&&&&&&&&  BLACKBOX OUTPUT  &&&&&&&&&&&&&&&&&&&&&&&&&&&
!
      IF(IBLACKWHAT.EQ.2)THEN
!
!  WRITE OUT W MATRIX
!
         KKKK=0
         DO 621 II=1,NDIMENSIONS
         DO 622 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
  622    CONTINUE
         KKKK=KKKK+II+2
  621    CONTINUE
!
!
!  WRITE OUT PSI MATRIX
!
         KKKK=0
         DO 721 II=1,NDIMENSIONS
         DO 722 I=1,NRESPONDENTS
!
!  WRITE OUT ONLY SCALED RESPONDENTS
!
         IF(LMARK(I).EQ.0)GO TO 722
!
!  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
  722    CONTINUE
         KKKK=KKKK+II
  721    CONTINUE
      ENDIF
!
!
! &&&&&&&&&&&&&&&&&&&&&&  BLACKBOX_TRANSPOSE OUTPUT  &&&&&&&&&&&&&&&&&&&&&&&&&&&
!
      IF(IBLACKWHAT.EQ.3)THEN
!
!  WRITE OUT PSI MATRIX
!
         KKKK=0
         DO 221 II=1,NDIMENSIONS
         DO 222 I=1,NISSUES
!
!  TRANSFER WMATRIX(PSI MATRIX) PARAMETERS
!
         IF(KKKK.EQ.0)THEN
            WRITE(23,225)CAND(I),LRESPONDENTS(I+NRESPONDENTS),&
           (WMATRIX((I-1)*(II+1)+JJ),JJ=1,II+1)
         ENDIF
         IF(KKKK.GT.0)THEN
            WRITE(23,225)CAND(I),LRESPONDENTS(I+NRESPONDENTS),&
           (WMATRIX((KKKK)*NISSUES+(I-1)*(II+1)+JJ),JJ=1,II+1)
         ENDIF
  222    CONTINUE
         KKKK=KKKK+II+1
  221    CONTINUE
!
!  WRITE OUT W MATRIX PARAMETERS
!
         KKKK=0
         DO 421 II=1,NDIMENSIONS
         KXKK=0
         DO 422 I=1,NRESPONDENTS
!
!  WRITE OUT ONLY SCALED RESPONDENTS
!
         IF(LMARK(I).EQ.0)GO TO 422
         KXKK=KXKK+1
!
!  TRANSFER PSIMATRIX PARAMETERS
!
         IF(KKKK.EQ.0)THEN
            WRITE(28,201)MID(I),LRESPONDENTS(KXKK),&
            (PSIMATRIX((I-1)*(II+2)+JJ),JJ=1,II+2)
         ENDIF
         IF(KKKK.GT.0)THEN
          WRITE(28,201)MID(I),LRESPONDENTS(KXKK),&
           (PSIMATRIX(KKKK*NRESPONDENTS+(I-1)*(II+2)+JJ),JJ=1,II+2)
         ENDIF
  422    CONTINUE
         KKKK=KKKK+II+2
  421    CONTINUE
!
      ENDIF
!
      IF(IBLACKWHAT.NE.1)THEN
         DO 69 J=1,NDIMENSIONS
         WRITE(23,113)J,(FITS(JJ+(J-1)*7),JJ=1,7)
  69     CONTINUE
         WRITE(23,300)(FITS2(J),J=1,6)
      ENDIF
!
      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,      &
                          FITS2,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(7*NDIMENSIONS),FITS2(6), &
        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-1)*NISSUES+J).EQ.KMISS((J-1)*NMISSING+K))THEN 
         KISSUE((I-1)*NISSUES+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-1)*NISSUES+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,FITS2,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)*7)=RSAVE(J,1)
      FITS(2+(J-1)*7)=CCC
      FITS(3+(J-1)*7)=SUM*100.0
      FITS(4+(J-1)*7)=AA*100.0
      FITS(5+(J-1)*7)=RSAVE(J,2)
      FITS(6+(J-1)*7)=RSAVE(J,3)
      FITS(7+(J-1)*7)=YHAT(J)
!
      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(WORK)
      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,FITS2,IPRNT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XBIGONE(NRESPONDENTS,NISSUES),                       &
                XDATA(NRESPONDENTS,NDIMENSIONS),                     &
                W(NISSUES,NDIMENSIONS+2),FITS2(6) 
!                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
         FITS2(1)=NP
         FITS2(2)=NY
         FITS2(3)=KTOT
         FITS2(4)=LTOT
         FITS2(5)=PMISS*100.0
         FITS2(6)=SVSUM
         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
!
!************************************************************************
!************************************************************************
!*                                                                     **
!*  BLACKTBOX_TRANSPOSE.FOR                                            **
!*                                                                     **
!*           FORTRAN 95 CONVERSION APRIL 2009                          **
!*                                                                     **
!*  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.            **
!*                                                                     **
!*                                                                     **
!*                                                                     **
!*  TECHNIQUE GENERALIZED TO HANDLE TRANSPOSED MATRICES, JULY 1996.    **
!*                                                                     **
!*                                                                     **
!*                                                                     **
!*    THIS VERSION CURRENT AS OF:   28 MAY 2008 (UPDATED 1998 VERSION) **
!*                                                                     **
!*                                                                     **
!************************************************************************
!************************************************************************
!
      SUBROUTINE BLACKBOXT(NRESPONDENTS,NISSUES,NDIMENSIONS,NMISSING,&
                          KMISS,MINSCALE,MID,KISSUE,CAND,FITS,       &
                          PSIMATRIX,WMATRIX,LRESPONDENTS,LMARK,      &
                          FITS2,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+NRESPONDENTS),LMARK(NRESPONDENTS)
      DIMENSION FITS(7*NDIMENSIONS),FITS2(6),                        &
        PSIMATRIX(NRESPONDENTS*((NDIMENSIONS*(NDIMENSIONS+1))/2)+2*&
      NRESPONDENTS*NDIMENSIONS),                                     &
        WMATRIX(NISSUES*((NDIMENSIONS*(NDIMENSIONS+1))/2)+2*NISSUES*&
      NDIMENSIONS)  
!
      CHARACTER*21 CAND(NISSUES)
!
      INTEGER, ALLOCATABLE :: KID(:)
      INTEGER, ALLOCATABLE :: LID(:)
      DOUBLE PRECISION, ALLOCATABLE :: YHAT(:)
      DOUBLE PRECISION, ALLOCATABLE :: WORK(:)
      DOUBLE PRECISION, ALLOCATABLE :: WORK4(:)
      DOUBLE PRECISION, ALLOCATABLE :: WORK5(:)
      DOUBLE PRECISION, ALLOCATABLE :: XBIGONE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: W(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: W2(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: WSAVE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XDATA(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XDATA2(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: PSISAVE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XT(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: RSAVE(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: UUU(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: U4(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: VVV(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: V4(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: WORK2(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: WORK3(:,:)
!
      ALLOCATE(KID(2*NRESPONDENTS))
      ALLOCATE(LID(2*NRESPONDENTS))
      ALLOCATE(YHAT(NRESPONDENTS))
      ALLOCATE(WORK(NRESPONDENTS+NISSUES))
      ALLOCATE(WORK4(NDIMENSIONS))
      ALLOCATE(WORK5(NRESPONDENTS))
      ALLOCATE(XBIGONE(NISSUES,NRESPONDENTS))
      ALLOCATE(W(NRESPONDENTS,NDIMENSIONS+2))
      ALLOCATE(W2(NRESPONDENTS,NDIMENSIONS+2))
      ALLOCATE(WSAVE(NRESPONDENTS,NDIMENSIONS+2))
      ALLOCATE(XDATA(NISSUES,NDIMENSIONS))
      ALLOCATE(XDATA2(NISSUES,NDIMENSIONS))
      ALLOCATE(PSISAVE(NISSUES,NDIMENSIONS))
      ALLOCATE(XT(NRESPONDENTS,NISSUES))
      ALLOCATE(RSAVE(NDIMENSIONS,3))
      ALLOCATE(UUU(NRESPONDENTS,NISSUES))
      ALLOCATE(U4(NDIMENSIONS,NDIMENSIONS))
      ALLOCATE(VVV(NRESPONDENTS,NISSUES))
      ALLOCATE(V4(NDIMENSIONS,NDIMENSIONS))
      ALLOCATE(WORK2(NDIMENSIONS,NDIMENSIONS))
      ALLOCATE(WORK3(NDIMENSIONS,NDIMENSIONS))
!
  70  FORMAT(' ITERATION RECORD'/3X,'DIM',5X,'ERROR',9X,'EXPLAINED', &
               7X,'PERCENT',3X,'CUM PERCENT',3X,'R-SQUARE',2X,       &
                                                     'STD ERR EST')
  113 FORMAT(1X,I3,2F16.4,7F12.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(1X,a21,I5,1X,35F7.3)
  202 FORMAT(5X,'W-TRANSPOSE*W')
  225 FORMAT(I7,I5,35F7.3)
  301 FORMAT(' NUMBER OF CASES',3I6)
  315 FORMAT(' R-SQUARE CHECK',I8,F7.3)
  330 FORMAT(I5,7F10.4)
  331 FORMAT(' ROTATION MATRIX FIRST COLUMN',I4,10F12.6)
  332 FORMAT(' WORK4',I5,7F12.6)
  333 FORMAT(' U4   ',I5,7F12.6)
  334 FORMAT(' V4   ',I5,7F12.6)
  335 FORMAT(' WORK3',I5,7F12.6)
  337 FORMAT(' RANK CHECK OF PSI*W',I5)
  338 FORMAT(' CHECK OF VVV * XDATA')
 1002 FORMAT(8I5)
 1003 FORMAT(A10)
 1082 format(A120)
 2222 FORMAT(1X,78('*'))
!
!
!  PRINTER SWITCH -- IPRNT=1 NO WRITE TO DISK
!                    IPRNT=0 WRITE TO DISK
      IPRNT=1
!
      LWORK=NRESPONDENTS+NISSUES
      NS=NDIMENSIONS
      NOBS=NISSUES
      NMISS=NMISSING
!
! MINIMUM RESPONSE IN A COLUMN
!
      KVMIN=NS+2
!
! MINIMUM RESPONSES IN A ROW
!
      KVMIN2=8
!
!      READ(4,1002)NS,NOBS,NMISS
      II=0
      DO 8888 I=1,NRESPONDENTS
      LMARK(I)=0
!  500 READ(22,FMT,END=400)MID,(KISSUE(J),J=1,NOBS)
      IMARK=0
      DO 1 J=1,NOBS
      DO 112 K=1,NMISS
      IF(KISSUE((I-1)*NISSUES+J).EQ.KMISS((J-1)*NMISSING+K))THEN 
         KISSUE((I-1)*NISSUES+J)=-999
         IMARK=IMARK+1
         GO TO 1
      ENDIF
  112 CONTINUE
  1   CONTINUE
      IF((NOBS-IMARK).LT.KVMIN)GO TO 8888
      II=II+1
      LMARK(I)=1
      KID(II)=MID(I)
      DO 11 J=1,NOBS    
      XBIGONE(J,II)=FLOAT(KISSUE((I-1)*NISSUES+J))
  11  CONTINUE
 8888 CONTINUE
      IF(IPRNT.EQ.0)WRITE(23,301)I,II
      NY=II
!      WRITE(23,301)NY
!      WRITE(*,301)NY
!
!  SUBROUTINE BLACKB--GENERATES MULTI-DIMENSIONAL STARTING VALUES 
!    FOR LEGISLATORS BY PERFORMING A PSEUDO-SINGULAR VALUE 
!    DECOMPOSITION OF NY ROLL CALLS
!
!
      KKKK=0
      LLLL=0
!      WRITE(23,2222)
!      WRITE(23,2222)
      DO 999 KKK=1,NS
!      WRITE(23,2222)
!      WRITE(23,196)KKK
!      WRITE(*,196)KKK
      call BLACKBT(NRESPONDENTS,NISSUES,NDIMENSIONS,NOBS,NY,KKK,     &
                   KVMIN2,XBIGONE,XDATA,W,SVSUM,FITS2)
!
!  WRITE OUT ROW PARAMETERS
!
!      DO 2 I=1,NOBS
!      WRITE(24,201)CAND(I),(XDATA(I,J),J=1,KKK)
!  2   CONTINUE
!  
!  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
!      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
!      WRITE(23,197)
      DO 229 I=1,KKK
!      WRITE(23,200)I,(VVV(I,J),J=1,KKK)
  229 CONTINUE
!      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
!      WRITE(23,200)I,(VVV(I,J),J=1,KKK)
  232 CONTINUE
!
!  CALCULATE R-SQUARES FOR COLUMNS OF  PSI*W + Jp'C 
!
      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(J,I)=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
!      WRITE(23,315)KK,RRR
      RSAVE(KKK,1)=SUME
      RSAVE(KKK,2)=RRR
      RSAVE(KKK,3)=SQRT(SUME/FLOAT(KK - KKK*(NOBS+NY) + NY))
!
!  WRITE OUT COLUMN PARAMETERS
!
      DO 222 I=1,NY
      LRESPONDENTS(I)=LID(I)
!      WRITE(28,225)KID(I),LID(I),(W(I,JJ),JJ=1,KKK+2)
  222 CONTINUE
!
!  *****patched in code*****
!
!      XTOL=.001
!      CALL LSVRR(NY,NOBS,XT,NRESPONDENTS,21,XTOL,IRANK,YHAT,UUU,
!     C           NRESPONDENTS,VVV,NRESPONDENTS)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
      CALL DGESVD('S','S',NY,NOBS,XT,NRESPONDENTS,YHAT,UUU,          &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
!      WRITE(23,337)INFO
!
!  COMPARE VVV(NS,NOBS) WITH XDATA(NOBS,NS) AS A SAFETY CHECK
!
      DO 671 JP=1,KKK
      DO 672 KP=1,KKK
      SUM=0.0
      DO 670 I=1,NOBS
      SUM=SUM+VVV(JP,I)*XDATA(I,KP)
  670 CONTINUE
      WORK2(JP,KP)=SUM
  672 CONTINUE
  671 CONTINUE
!      WRITE(23,338)
!      WRITE(*,338)
      DO 673 KP=1,KKK
!      WRITE(23,330)KP,(WORK2(KP,JP),JP=1,KKK)
!      WRITE(*,330)KP,(WORK2(KP,JP),JP=1,KKK)
  673 CONTINUE
!
!
!  CALCULATE ROTATION MATRIX HERE
!
      DO 765 JJ=1,KKK
      WORK4(JJ)=0.0
      DO 765 JP=1,KKK
      WORK2(JJ,JP)=0.0
      WORK3(JJ,JP)=0.0
  765 CONTINUE
      DO 464 I=1,NOBS
      IF(KKK.EQ.2)THEN
         DO 763 JP=1,KKK-1
         DO 764 JJ=1,KKK
         WORK2(JJ,JP)=WORK2(JJ,JP)+(VVV(JJ,I)*PSISAVE(I,JP))
  764    CONTINUE
  763    CONTINUE
      ENDIF
!      WRITE(23,330)I,YHAT(I)
!      WRITE(*,330)I,YHAT(I)
      WORK5(I)=YHAT(I)
  464 CONTINUE
      IF(KKK.EQ.2)THEN
         DO 743 JP=1,KKK-1
         SUM=0.0
         DO 842 JJ=1,KKK
         SUM=SUM+WORK2(JJ,JP)**2
  842    CONTINUE
         DO 741 JJ=1,KKK
         WORK2(JJ,JP)=WORK2(JJ,JP)/SQRT(SUM)
  741    CONTINUE
  743    CONTINUE
!
         DO 744 JJ=1,KKK
!         WRITE(23,331)JJ,(WORK2(JJ,JP),JP=1,KKK-1)
!         WRITE(*,331)JJ,(WORK2(JJ,JP),JP=1,KKK-1)
  744    CONTINUE
!
         WORK2(1,2)=-WORK2(2,1)
         WORK2(2,2)=WORK2(1,1)
!
         DO 761 KP=1,KKK
         DO 762 JP=1,KKK
         SUM=0.0
         DO 766 JJ=1,KKK
         SUM=SUM+WORK2(JJ,JP)*WORK2(JJ,KP)
  766    CONTINUE
         WORK3(JP,KP)=SUM
  762    CONTINUE
  761    CONTINUE
         DO 748 JJ=1,KKK
!         WRITE(23,335)JJ,(WORK3(JJ,JP),JP=1,KKK)
!         WRITE(*,335)JJ,(WORK3(JJ,JP),JP=1,KKK)
  748    CONTINUE
         DO 760 JJ=1,KKK
!         WRITE(23,331)JJ,(WORK2(JJ,JP),JP=1,KKK)
!         WRITE(*,331)JJ,(WORK2(JJ,JP),JP=1,KKK)
  760    CONTINUE
!
!  ROTATE W
!
         DO 758 J=1,NY
         DO 757 IJJ=1,KKK
         SUM=0.0
         DO 756 JJ=1,KKK
         SUM=SUM+WORK2(JJ,IJJ)*W(J,JJ+1)
  756    CONTINUE
         W2(J,IJJ+1)=SUM
         W2(J,1)=W(J,1)
         W2(J,IJJ+2)=W(J,IJJ+2)
  757    CONTINUE
  758    CONTINUE
!
!  ROTATE N BY S PSI-MATRIX 
!
         DO 801 I=1,NOBS 
         DO 701 IJJ=1,KKK
         SUM=0.0
         SUM2=0.0
         DO 723 JJ=1,KKK
         SUM=SUM+VVV(JJ,I)*WORK2(JJ,IJJ)
         SUM2=SUM2+XDATA(I,JJ)*WORK2(JJ,IJJ)
  723    CONTINUE
         XDATA2(I,IJJ)=SUM2
         UUU(I,IJJ)=SUM
  701    CONTINUE
  801    CONTINUE
      ENDIF
!
!
!  CALCULATE R-SQUARES FOR ROWS OF PSI*W + Jp'C 
!
      ASUM=0.0
      BSUM=0.0
      CSUM=0.0
      DSUM=0.0
      ESUM=0.0
      SUME=0.0
      KK=0
      DO 623 I=1,NOBS
      AASUM=0.0
      BBSUM=0.0
      CCSUM=0.0
      DDSUM=0.0
      EESUM=0.0
      KJJ=0
      DO 624 J=1,NY
      SUM=0.0
      DO 626 K=1,KKK
      SUM=SUM+XDATA(I,K)*W(J,K+1)
  626 CONTINUE
!      XT(J,I)=SUM
      AA=SUM+W(J,1)
      IF(XBIGONE(I,J).EQ.-999.0)GO TO 624
      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
  624 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)
      YHAT(I)=RRR
      LID(I)=KJJ
      ASUM=ASUM+AASUM
      BSUM=BSUM+BBSUM
      CSUM=CSUM+CCSUM
      DSUM=DSUM+DDSUM
      ESUM=ESUM+EESUM
      KK=KK+KJJ
  623 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
!      WRITE(23,315)KK,RRR
!
!
!  CALCULATE R-SQUARES FOR ROTATED ROWS OF PSI*W + Jp'C 
!
      ASUM=0.0
      BSUM=0.0
      CSUM=0.0
      DSUM=0.0
      ESUM=0.0
      SUME=0.0
      ASUM2=0.0
      BSUM2=0.0
      CSUM2=0.0
      DSUM2=0.0
      ESUM2=0.0
      SUME2=0.0
      KK=0
      DO 823 I=1,NOBS
      DO 923 JJ=1,KKK
      IF(KKK.EQ.1)PSISAVE(I,JJ)=VVV(JJ,I)
      IF(KKK.EQ.2)THEN
         VVV(JJ,I)=UUU(I,JJ)
      ENDIF
  923 CONTINUE
      AASUM=0.0
      BBSUM=0.0
      CCSUM=0.0
      DDSUM=0.0
      EESUM=0.0
      AASUM2=0.0
      BBSUM2=0.0
      CCSUM2=0.0
      DDSUM2=0.0
      EESUM2=0.0
      KJJ=0
      DO 824 J=1,NY
      SUM=0.0
      SUM2=0.0
      DO 826 K=1,KKK
      SUM=SUM+XDATA(I,K)*W(J,K+1)
      IF(KKK.EQ.2)THEN
         SUM2=SUM2+XDATA2(I,K)*W2(J,K+1)
      ENDIF
  826 CONTINUE
!      XT(J,I)=SUM
      AA=SUM+W(J,1)
      IF(XBIGONE(I,J).EQ.-999.0)GO TO 824
      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
      IF(KKK.EQ.2)THEN
         AA2=SUM2+W2(J,1)
         SUME2=SUME2+(AA2-BB)**2
         AASUM2=AASUM2+AA2
         BBSUM2=BBSUM2+BB
         CCSUM2=CCSUM2+AA2*AA2
         DDSUM2=DDSUM2+BB*BB
         EESUM2=EESUM2+AA2*BB
      ENDIF
      KJJ=KJJ+1
  824 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)
      YHAT(I)=RRR
      LID(I)=KJJ
      ASUM=ASUM+AASUM
      BSUM=BSUM+BBSUM
      CSUM=CSUM+CCSUM
      DSUM=DSUM+DDSUM
      ESUM=ESUM+EESUM
      IF(KKK.EQ.2)THEN
         AAA2=FLOAT(KJJ)*EESUM2-AASUM2*BBSUM2
         BBB2=FLOAT(KJJ)*CCSUM2-AASUM2*AASUM2
         CCC2=FLOAT(KJJ)*DDSUM2-BBSUM2*BBSUM2
         RRR2=0.0
         IF(ABS(BBB2*CCC2).GT.0.0)RRR2=(AAA2*AAA2)/(BBB2*CCC2)
         YHAT(I)=RRR
         LID(I)=KJJ
         ASUM2=ASUM2+AASUM2
         BSUM2=BSUM2+BBSUM2
         CSUM2=CSUM2+CCSUM2
         DSUM2=DSUM2+DDSUM2
         ESUM2=ESUM2+EESUM2
      ENDIF
      KK=KK+KJJ
  823 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
!      WRITE(23,315)KK,RRR
      IF(KKK.EQ.2)THEN
         AA=FLOAT(KK)*ESUM2-ASUM2*BSUM2
         BB=FLOAT(KK)*CSUM2-ASUM2*ASUM2
         CC=FLOAT(KK)*DSUM2-BSUM2*BSUM2
         RRR=(AA*AA)/(BB*CC)
!         WRITE(*,315)KK,RRR
!         WRITE(23,315)KK,RRR
      ENDIF
!
!  TRANSFER ROTATED MATRICES
!
      IF(KKK.EQ.2)THEN
         DO 806 K=1,KKK
         DO 807 I=1,NOBS
         XDATA(I,K)=XDATA2(I,K)
  807    CONTINUE
         DO 808 J=1,NY
         W(J,K+1)=W2(J,K+1)
  808    CONTINUE
  806    CONTINUE
      ENDIF
!
!  WRITE OUT COLUMN PARAMETERS
!
      KI=0
      DO 622 I=1,NRESPONDENTS
!
!  TRANSFER W PARAMETERS TO PSIMATRIX(.,.)
!
!  STORE ONLY SCALED RESPONDENTS
!
      IF(LMARK(I).EQ.0)GO TO 633
      KI=KI+1
!
      DO 632 JJ=1,KKK+2
      IF(LLLL.EQ.0)THEN
         PSIMATRIX((I-1)*(KKK+2)+JJ)=W(KI,JJ)
      ENDIF
      IF(LLLL.GT.0)THEN
         PSIMATRIX((LLLL)*NRESPONDENTS+(I-1)*(KKK+2)+JJ)=W(KI,JJ)
      ENDIF
  632 CONTINUE
  633 CONTINUE
  622 CONTINUE
!
!  WRITE OUT ROW PARAMETERS
!
!
      DO 742 I=1,NOBS
!
!  TRANSFER PSIMATRIX PARAMETERS
!
      DO 532 JJ=1,KKK
      IF(KKKK.EQ.0)THEN
         WMATRIX((I-1)*(KKK+1)+JJ)=VVV(JJ,I)
         IF(JJ.EQ.KKK)THEN
            WMATRIX((I-1)*(KKK+1)+KKK+1)=YHAT(I)
         ENDIF
      ENDIF
      IF(KKKK.GT.0)THEN
         WMATRIX((KKKK)*NOBS+(I-1)*(KKK+1)+JJ)=VVV(JJ,I)
         IF(JJ.EQ.KKK)THEN
            WMATRIX((KKKK)*NOBS+(I-1)*(KKK+1)+KKK+1)=YHAT(I)
         ENDIF
      ENDIF
  532 CONTINUE
      LRESPONDENTS(I+NRESPONDENTS)=LID(I)
!      WRITE(24,201)CAND(I),LID(I),(VVV(J,I),J=1,KKK),YHAT(I)
!      WRITE(*,201)CAND(I), LID(I),(VVV(J,I),J=1,KKK),YHAT(I)
  742 CONTINUE
      LLLL=LLLL+KKK+2
      KKKK=KKKK+KKK+1
  999 CONTINUE
!      WRITE(23,2222)
!      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)*7)=RSAVE(J,1)
      FITS(2+(J-1)*7)=CCC
      FITS(3+(J-1)*7)=SUM*100.0
      FITS(4+(J-1)*7)=AA*100.0
      FITS(5+(J-1)*7)=RSAVE(J,2)
      FITS(6+(J-1)*7)=RSAVE(J,3)
      FITS(7+(J-1)*7)=WORK5(J)
!
!      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)
  69  CONTINUE
!      WRITE(23,2222)
!
      DEALLOCATE(KID)
      DEALLOCATE(LID)
      DEALLOCATE(YHAT)
      DEALLOCATE(WORK)
      DEALLOCATE(WORK4)
      DEALLOCATE(WORK5)
      DEALLOCATE(XBIGONE)
      DEALLOCATE(W)
      DEALLOCATE(W2)
      DEALLOCATE(WSAVE)
      DEALLOCATE(XDATA)
      DEALLOCATE(XDATA2)
      DEALLOCATE(PSISAVE)
      DEALLOCATE(XT)
      DEALLOCATE(RSAVE)
      DEALLOCATE(UUU)
      DEALLOCATE(U4)
      DEALLOCATE(VVV)
      DEALLOCATE(V4)
      DEALLOCATE(WORK2)
      DEALLOCATE(WORK3)
!
      RETURN
      END
!
!
!  ***********************************************************************
!    SUBROUTINE BLACKBT
!      
!      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 BLACKBT(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NY,NF,  &
                         NFX,XBIGONE,XDATA,W,SVSUM,FITS2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION W(NRESPONDENTS,NDIMENSIONS+2),                       &
                XBIGONE(NISSUES,NRESPONDENTS),                       &
                XDATA(NISSUES,NDIMENSIONS),FITS2(6)
!
      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 :: YHAT(:)
      DOUBLE PRECISION, ALLOCATABLE :: WK(:)
      DOUBLE PRECISION, ALLOCATABLE :: WORK(:)
!
      DOUBLE PRECISION, ALLOCATABLE :: X(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: PSIX(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XS(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ROOTC(:,:)
      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(NRESPONDENTS))
      ALLOCATE(MPOS(NRESPONDENTS))
!
      ALLOCATE(XX(2*NRESPONDENTS))
      ALLOCATE(D(NRESPONDENTS))
      ALLOCATE(DD(NISSUES))
      ALLOCATE(DC(NRESPONDENTS))
      ALLOCATE(CC(NRESPONDENTS))
      ALLOCATE(TSUM(2*NRESPONDENTS))
      ALLOCATE(DDD(NISSUES))
      ALLOCATE(YHAT(NRESPONDENTS))
      ALLOCATE(WK(2*NRESPONDENTS))
      ALLOCATE(WORK(NRESPONDENTS+NISSUES))
!
      ALLOCATE(X(NISSUES,NRESPONDENTS))
      ALLOCATE(PSIX(NISSUES,NDIMENSIONS+1))
      ALLOCATE(XS(NISSUES,NRESPONDENTS))
      ALLOCATE(ROOTC(NRESPONDENTS,NISSUES))
      ALLOCATE(XT(NISSUES,NRESPONDENTS))
      ALLOCATE(XSS(NISSUES,NRESPONDENTS))
      ALLOCATE(UUU(NRESPONDENTS,NISSUES))
      ALLOCATE(VVV(NRESPONDENTS,NISSUES))
      ALLOCATE(WSAVE(NRESPONDENTS,NDIMENSIONS+1))
      ALLOCATE(XSAVE(NISSUES,NDIMENSIONS+2))
!
  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:  PSI*W')
  113 FORMAT(1X,I5,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               =',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=NISSUES+NRESPONDENTS
      IF(NY.LT.NP)NYMIN=NY
      IF(NP.LE.NY)NYMIN=NP
!      NYMIN=AMIN0(NY,NP)
      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
!      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
!      IF(X(I,J).EQ.0.0)X(I,J)=9.0
!      IF(X(I,J).NE.9.0)SUM=SUM+X(I,J)*X(I,J)
!      IF(X(I,J).NE.9.0)SUMA=SUMA+X(I,J)
!      IF(X(I,J).NE.9.0)KK=KK+1
      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
!
      SVSUM=SVSUM-(SWSUM**2)/FLOAT(KTOT)
      LTOT=NP*NY
      PMISS=FLOAT(LTOT-KTOT)/FLOAT(LTOT)
      LTOT=LTOT-KTOT
      IF(NF.EQ.1)THEN
         FITS2(1)=NP
         FITS2(2)=NY
         FITS2(3)=KTOT
         FITS2(4)=LTOT
         FITS2(5)=PMISS*100.0
         FITS2(6)=SVSUM
!         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   
      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
!
!  COMPUTE CORRELATION NY BY NY CORRELATION MATRIX AND VECTOR OF SIGN
!    CHANGES FOR COLUMNS OF X/XS/XSS
!
      call CORR3(NRESPONDENTS,NISSUES,NP,NY,X,LL,MPOS,KS,KPOS)
!      DO 667 J=1,NY
!      WRITE(23,301)J,LL(J)
!  667 CONTINUE
!
!  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
!
      XT(I,1)=-999.0
      IF(KK.LT.NFX)GO TO 3
      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
      IF(XT(I,1).EQ.-999.0)GO TO 40
      XT(I,1)=(XT(I,1)-TXB)
      XT(I,2)=1.0
!      WRITE(23,338)I,XT(I,1),XT(I,2)
  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 REGT(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,1,NY,TSUM,W,XS,X,&
                                                     XT,1,1,KKK,AREG)
      call REG2T(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,1,NY,W,XS,X,XT, &
                                                  PXB,PXS,KKK,1,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)MM,XXK,QXS,PXS,XCOR,PXB,TXB
      DO 444 I=1,NP
      IF(XT(I,1).EQ.-999.0)PSIX(I,JJJ)=-999.0
      IF(PSIX(I,JJJ).EQ.-999.0)GO TO 444
!
!  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 CORR3(NRESPONDENTS,NISSUES,NP,NY,X,LL,MPOS,KS,KPOS)
 9999 CONTINUE
!
!  WRITE OUT CORRELATION AND CROSS PRODUCT MATRIX OF P
!
!      WRITE(23,23)
!      call CORR2(NRESPONDENTS,NISSUES,NP,NF,PSIX,R,LL,MPOS,KS,KPOS,1)
!      WRITE(23,521)
!      call PSIPRM(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF1,PSIX)
!
!  PERFORM A.L.S. WITH FULL P MATRIX UNTIL CONVERGENCE
!
      DO 222 NN=1,5
      call REGT(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF,NY,TSUM,W,XS,&
                                                   X,PSIX,1,1,NF,AREG)
      call REG2T(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF,NY,W,XS,X,&
                                              PSIX,PXB,PXS,NF,NF,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
      IF(PSIX(I,1).EQ.-999.0)GO TO 60
      DO 61 K=1,NY
      WK(K)=0.0
      WK(K+NY)=0.0
      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
!
      IF(NY.GT.NP)THEN
!         DO 710 I=1,NP
!         DO 710 J=1,NP
!         SUM=0.0
!         DO 730 K=1,NY
!  730    SUM=SUM+XT(I,K)*XT(J,K)
!  710    ROOTC(I,J)=SUM
!         CALL RS(NRESPONDENTS,NP,ROOTC,D,1,CROOT,FAL1,FAL2,IERR)
         DO 729 I=1,NP
         DO 729 J=1,NY
         ROOTC(J,I)=XT(I,J)
  729    CONTINUE
         XTOL=.001
!         CALL LSVRR(NY,NP,ROOTC,NRESPONDENTS,21,XTOL,IRANK,YHAT,UUU,
!     C           NRESPONDENTS,VVV,NRESPONDENTS)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
      CALL DGESVD('S','S',NY,NP,ROOTC,NRESPONDENTS,YHAT,UUU,         &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
!      WRITE(23,337)INFO
!         WRITE(23,337)IRANK
         DO 464 I=1,NP
!         WRITE(23,338)I,YHAT(I)
         D(I)=YHAT(I)
  464    CONTINUE
      ENDIF
!      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'
!
      IF(NY.GT.NP)THEN
         DO 750 I=1,NP
         DO 750 J=1,NP
         SUM=0.0
         DO 760 K=1,NY
         SUM=SUM+XSS(I,K)*XSS(J,K)
  760    CONTINUE
  750    ROOTC(I,J)=SUM
!         CALL RS(NRESPONDENTS,NP,ROOTC,DDD,1,CROOT,FAL1,FAL2,IERR)
         DO 728 I=1,NP
         DO 728 J=1,NY
         ROOTC(J,I)=XSS(I,J)
  728    CONTINUE
         XTOL=.001
!         CALL LSVRR(NY,NP,ROOTC,NRESPONDENTS,21,XTOL,IRANK,YHAT,UUU,
!     C           NRESPONDENTS,VVV,NRESPONDENTS)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
!      LWORK=3539
      CALL DGESVD('S','S',NY,NP,ROOTC,NRESPONDENTS,YHAT,UUU,         &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
!      WRITE(23,337)INFO
!         WRITE(23,337)IRANK
         DO 463 I=1,NP
!         WRITE(23,338)I,YHAT(I)
         DDD(I)=YHAT(I)
  463    CONTINUE
      ENDIF
!      WRITE(23,7200)I
!      WRITE(*,7200)IERR
!
!  PERFORM SINGULAR VALUE DECOMPOSITION OF PW'
!
      IF(NY.GT.NP)THEN
         DO 770 I=1,NP
         DO 770 J=1,NP
         SUM=0.0
         DO 780 K=1,NY
         SUM=SUM+XT(I,K)*XT(J,K)
  780    CONTINUE
  770    ROOTC(I,J)=SUM
!         CALL RS(NRESPONDENTS,NP,ROOTC,DD,1,CROOT,FAL1,FAL2,IERR)
         DO 727 I=1,NP
         DO 727 J=1,NY
         ROOTC(J,I)=XT(I,J)
  727    CONTINUE
         XTOL=.001
!         CALL LSVRR(NY,NP,ROOTC,NRESPONDENTS,21,XTOL,IRANK,YHAT,UUU,
!     C           NRESPONDENTS,VVV,NRESPONDENTS)
!
!  *****NOTE THAT DGESVD RETURNS U, LAMBDA, AND V_transpose !!!  *****
!
!      LWORK=3539
      CALL DGESVD('S','S',NY,NP,ROOTC,NRESPONDENTS,YHAT,UUU,         &
                  NRESPONDENTS,VVV,NRESPONDENTS,WORK,LWORK,INFO)
!      WRITE(23,337)INFO
!         WRITE(23,337)IRANK
         DO 462 I=1,NP
         DD(I)=YHAT(I)
         DO 462 JJ=1,NF
!         PSIX(I,JJ)=VVV(I,JJ)*SQRT(YHAT(JJ))
         PSIX(I,JJ)=VVV(JJ,I)*SQRT(YHAT(JJ))
  462    CONTINUE
         DO 460 I=1,NY
         DO 460 JJ=1,NF
         W(I,JJ)=UUU(I,JJ)*SQRT(YHAT(JJ))
  460    CONTINUE
!  
!  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,NFPS)
         DO 227 J=1,NFPS
         DO 227 K=1,NFPS
         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,NFPS
!         WRITE(23,200)I,(VVV(I,J),J=1,NFPS)
  229    CONTINUE
         DO 230 J=1,NFPS
         DO 230 K=1,NFPS
         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,NFPS
!         WRITE(23,200)I,(VVV(I,J),J=1,NFPS)
  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        
!
      ENDIF
!      WRITE(23,86)
      DO 87 J=1,NYMIN
      IF(J.GT.NF+3)GO TO 87
!      WRITE(23,1113)J,D(J),DDD(J),DD(J)
!      WRITE(*,1113)J,D(J),DDD(J),DD(J)
  87  CONTINUE
!
!  COMPUTE STARTING ESTIMATE OF LEGISLATOR COORDINATES AS
!      XDATA(NOBS,NF) = P*(NOBS,NF) , WHERE
!      P*P*' = PW'WP'
!
!    ****NOTE THAT P* IS NOBSxNY, BUT ONLY FIRST NF COLUMNS USED
!        AS THE STARTS****
!
!      IF(NY.GT.NP)THEN
!         DO 900 K=1,NP
!         SUM1=0.0
!         SUM2=0.0
!         DO 890 I=1,NP
!         PSIX(I,K)=CROOT(I,NP+1-K)*SQRT(DD(K))
!         XDATA(I,K)=PSIX(I,K)
!         SUM1=SUM1+PSIX(I,K)
!         SUM2=SUM2+PSIX(I,K)**2
!  890    CONTINUE
!  900    CONTINUE
!      ENDIF
!
!      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(YHAT)
      DEALLOCATE(WK)
      DEALLOCATE(WORK)
!
      DEALLOCATE(X)
      DEALLOCATE(PSIX)
      DEALLOCATE(XS)
      DEALLOCATE(ROOTC)
      DEALLOCATE(XT)
      DEALLOCATE(XSS)
      DEALLOCATE(UUU)
      DEALLOCATE(VVV)
      DEALLOCATE(WSAVE)
      DEALLOCATE(XSAVE)
!
      RETURN
      END
!
!  **********************************************************************
!    SUBROUTINE CORR3---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 CORR3(NRESPONDENTS,NISSUES,NP,NY,X,LL,MPOS,KS,KPOS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(NISSUES,NRESPONDENTS),LL(NRESPONDENTS),            &
                MPOS(NRESPONDENTS)
!                
      DOUBLE PRECISION, ALLOCATABLE :: R(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: SA(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: SB(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: SC(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: SD(:,:)
!
      ALLOCATE(R(NRESPONDENTS,NRESPONDENTS))
      ALLOCATE(SA(NRESPONDENTS,NRESPONDENTS))
      ALLOCATE(SB(NRESPONDENTS,NRESPONDENTS))
      ALLOCATE(SC(NRESPONDENTS,NRESPONDENTS))
      ALLOCATE(SD(NRESPONDENTS,NRESPONDENTS))
!
  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)
  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
      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
!
      DEALLOCATE(R)
      DEALLOCATE(SA)
      DEALLOCATE(SB)
      DEALLOCATE(SC)
      DEALLOCATE(SD)
!
      RETURN
      END
!
!
!  ************************************************************************
!    SUBROUTINE REGT---CALLED BY BLACKBT.  PERFORMS REGRESSION TO ESTIMATE
!      W AND c.
!  ************************************************************************
!
      SUBROUTINE REGT(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF,NY,TSUM,&
                                    W,XS,X,PSI,IPRNT,ILAST,KKK,AREG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TSUM(2*NRESPONDENTS),W(NRESPONDENTS,NDIMENSIONS+2),  &
                XS(NISSUES,NRESPONDENTS),X(NISSUES,NRESPONDENTS),    &
                PSI(NISSUES,NRESPONDENTS)
!
      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(NRESPONDENTS))
      ALLOCATE(LL(NRESPONDENTS))
      ALLOCATE(MPOS(NRESPONDENTS))
!
      ALLOCATE(C(NRESPONDENTS))
!      ALLOCATE(RSUM(NRESPONDENTS))
      ALLOCATE(WVEC(4*NDIMENSIONS+8))
      ALLOCATE(WK(4*NDIMENSIONS+8))
      ALLOCATE(A(NDIMENSIONS+2,NDIMENSIONS+2))
      ALLOCATE(B(NDIMENSIONS+2,NDIMENSIONS+2))
      ALLOCATE(XT(NISSUES,NRESPONDENTS))
      ALLOCATE(R(NISSUES,NISSUES))
      ALLOCATE(ZMAT(NDIMENSIONS+2,NDIMENSIONS+2))
!
  100 FORMAT(' PERFORMANCE INDEX MATRIX INVERSE=',2I5)
  101 FORMAT(' PERFORMANCE INDEX EIGENVALUE/VECTOR ROUTINE=',3I5)
  102 FORMAT(2I5,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=',F16.4)
  300 FORMAT(' SSE AND FINAL WEIGHT MATRIX')
      NF1=NF+1
      LWORK=4*(NDIMENSIONS+2)
      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(J.EQ.1.AND.JJ.EQ.1)WRITE(23,665)I,PSI(I,1),PSI(I,2)
 665  FORMAT(I5,2F10.4)
      IF(XS(I,K).EQ.-999.0)GO TO 1
      IF(PSI(I,1).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,NDIMENSIONS+2,WVEC,WK,LWORK,INFO)
!      call rs(127,nf1,a,wvec,1,ZMAT,fv1,fv2,ier)
!      WRITE(*,667)NF1,INFO,(WVEC(JC),JC=1,NF1),ZMAT(1,1),ZMAT(1,2),&
!                   ZMAT(2,1),ZMAT(2,2),A(1,1),A(1,2),A(2,1),A(2,2)
!  667 FORMAT(' REG1 ',2I5,20F10.4)
!      KEITH=999
!      IF(KEITH.EQ.999)STOP
      DO 60 I=1,NF1
      DO 60 KX=1,NF1
      SUM=0.0
      DO 61 J=1,NF1
      IF(ABS(WVEC(J)).GT..001)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
      IF(PSI(I,1).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
      XT(I,1)=-999.0
      XT(I,2)=-999.0
      IF(XS(I,K).EQ.-999.0)GO TO 8
      IF(PSI(I,1).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 887 K=1,NY
  887 XT(I,K)=-999.0
      IF(PSI(I,1).EQ.-999.0)GO TO 888
      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
!
!  TRANSPOSE CALL
!
      call RSQUR(NISSUES,NRESPONDENTS,NP,NY,RR,XT,XS,IPRNT)
!      call RSQUR(NRESPONDENTS,NISSUES,NP,NY,RR,XT,XS,IPRNT)
      TSUM(KK+NY)=RR
  90  CONTINUE
!      IF(ILAST.EQ.1)WRITE(23,250)KKK,TSUM(KK)
!      IF(ILAST.EQ.0)WRITE(23,250)KK,TSUM(KK)
      AREG=TSUM(KK)
  11  CONTINUE
!
!  END OF DIMENSION ITERATION LOOP--ONLY IF ILAST = 0
!
!
      DEALLOCATE(C)
!      DEALLOCATE(RSUM)
      DEALLOCATE(WVEC)
      DEALLOCATE(WK)
      DEALLOCATE(A)
      DEALLOCATE(B)
      DEALLOCATE(XT)
      DEALLOCATE(R)
      DEALLOCATE(ZMAT)
!
      RETURN
      END
!
!  **************************************************************************
!    SUBROUTINE REG2T---CALLED BY BLACKBT.  PERFORMS REGRESSION TO ESTIMATE P.
!  **************************************************************************
!
      SUBROUTINE REG2T(NRESPONDENTS,NISSUES,NDIMENSIONS,NP,NF,NY,W,XS,&
                                     X,PSI,PXB,PXS,KKK,NWHO,BREG)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION W(NRESPONDENTS,NDIMENSIONS+2),XS(NISSUES,NRESPONDENTS),&
                X(NISSUES,NRESPONDENTS),PSI(NISSUES,NRESPONDENTS)
!
      DOUBLE PRECISION, ALLOCATABLE :: V(:)
      DOUBLE PRECISION, ALLOCATABLE :: Y(:)
      DOUBLE PRECISION, ALLOCATABLE :: SQU(:)
      DOUBLE PRECISION, ALLOCATABLE :: A(:,:)
!
      ALLOCATE(V(NRESPONDENTS))
      ALLOCATE(Y(NRESPONDENTS))
      ALLOCATE(SQU(2*NDIMENSIONS))
      ALLOCATE(A(NRESPONDENTS,NDIMENSIONS+2))
!
      NF1=NWHO+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
      IF(PSI(I,1).EQ.-999.0)GO TO 1
      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 REGAT(NRESPONDENTS,NDIMENSIONS,NS,NF,A,Y,V)
!
!  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
!      WRITE(23,350)KKK,ESUM
  350 FORMAT(' DIMENSION=',I3,1X,'TOTAL SSE REG2=',F16.4)
      BREG=ESUM
!
      DEALLOCATE(V)
      DEALLOCATE(Y)
      DEALLOCATE(SQU)
      DEALLOCATE(A)
!
      RETURN
      END
!
!
!  **************************************************************************
!    SUBROUTINE REGAT---CALLED BY REG2T.  PERFORMS THE REGRESSION:
!
!         [W'W]-1*W'[X(I,NY) - c']
!
!      TO GET THE ROW ENTRIES OF P.
!  **************************************************************************
!
      SUBROUTINE REGAT(NRESPONDENTS,NDIMENSIONS,NS,NF,A,Y,V)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NRESPONDENTS,NDIMENSIONS+2),Y(NRESPONDENTS),       &
                V(NRESPONDENTS)
!
      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(4*NDIMENSIONS+8))
      ALLOCATE(WK(4*NDIMENSIONS+8))
      ALLOCATE(B(NDIMENSIONS+2,NDIMENSIONS+2))
      ALLOCATE(C(NDIMENSIONS+2,NDIMENSIONS+2))
      ALLOCATE(BB(NDIMENSIONS+2,NRESPONDENTS))
      ALLOCATE(ZMAT(NDIMENSIONS+2,NDIMENSIONS+2))
!
  101 FORMAT(' PERFORMANCE INDEX EIGENVALUE/VECTOR ROUTINE=',3I5)
  102 FORMAT(2I5,30F10.4)
      LWORK=4*(NDIMENSIONS+2)
      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,NDIMENSIONS+2,WVEC,WK,LWORK,INFO)
!      call rs(127,nf,B,wvec,1,ZMAT,fv1,fv2,ier)
!      WRITE(23,101)NS,NF,IER
!      WRITE(23,102)NS,NF,(WVEC(J),J=1,NF)
      DO 60 I=1,NF
      DO 60 K=1,NF
      SUM=0.0
      DO 61 J=1,NF
      IF(ABS(WVEC(J)).GT..00001)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
!
!
!  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!
!  **************************************************************************
!
!   MCKALNEW.F95 -- 27 MARCH 2009                            
!
!             LAST UPDATE OF FORTRAN IV VERSION:  10 JULY 1996
!
!   ALDRICH-MCKELVEY SCALING PROGRAM
!
!
!   WRITTEN BY KEITH T. POOLE AT THE UNIVERSITY OF OREGON
!   OCTOBER 1978.  FOR THE THEORY OF THE PROGRAM SEE
!   ALDRICH AND MCKELVEY, "A METHOD OF SCALING WITH APPLICATIONS
!   TO THE 1968 AND 1972 PRESIDENTIAL ELECTIONS," AMERICAN
!   POLITICAL SCIENCE REVIEW, 1977, PP. 111-130
!
!   MODIFIED FOR USE ON IBM PC'S BY KEITH T. POOLE AT THE
!   GRADUATE SCHOOL OF INDUSTRIAL ADMINISTRATION, NOVEMBER, 1986.
!                
!                REVISED, JULY, 1996.
!
!   NQ=# OF STIMULI
!   NRESP= THE POSITION THE RESPONDENT IS IN THE DATA
!          IF NRESP=0, ONLY THE STIMULI CONFIGURATION WILL BE
!          COMPUTED
!   NMISS=# OF MISSING DATA VALUES
!   NLEFT=# OF THE STIMULUS WHICH SHOULD BE RECOVERED FURTHEST TO THE
!         LEFT OF THE DIMENSION.
!   XMISS IS AN ARRAY CONTAINING THE MISSING VALUES
!   KFMT IS AN ARRAY CONTAINING THE FORMAT OF THE INPUT DATA
!   KNAM IS AN ARRAY CONTAINING TITLES OF THE STIMULI.  UP TO 13 LETTERS
!        MAY BE USED.  ALWAYS PLACE "INDIV" AT END.
!   
!
!   UP TO 5 MISSING DATA VALUES MAY BE USED
!
!   
!   SAMPLE CONTROL CARDS
!
!   1011 4 7
!   (2X,I4,70X,F2.0,9(1X,F2.0),8X,F2.0)
!    00.0 97.0 98.0 99.0
!   SFIO
!   ALLRP
!   CENDM
!   PSU
!   REPIN
!   MRP
!   PCF
!   RAD
!   UNR
!   FGDS
!   INDIV
!
!
!   UNIT 21 IS USED TO READ IN THE CONTROL CARDS
!   UNIT 1 IS USED TO READ IN THE DATA SET
!   UNIT 20 IS USED FOR SCRATCH (NO LONGER USED)
!   UNIT 23 IS USED TO WRITE OUT THE INDIVIDUAL COORDINATES
!   UNIT 24 IS USED TO WRITE OUT THE OLD-FASHIONED PRINTER OUTPUT OF THE RUN
!
!
      SUBROUTINE MCKALNEW(NRESPONDENTS,NISSUES,NSELFPOS,NMISSING,   &
           KMISS,POLARITY,MID,KISSUE,CAND,FITS,PSIMATRIX,STIMCOORDS,&
           EIGENVALUES,EXITSTATUS)
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER NRESPONDENTS,NISSUES,NSELFPOS,NMISSING,               &
              KMISS(NISSUES*NMISSING),POLARITY,MID(NRESPONDENTS),   &
              KISSUE(NRESPONDENTS*NISSUES),EXITSTATUS               
!                
      DOUBLE PRECISION FITS(5),PSIMATRIX(NRESPONDENTS*4),           &
                       STIMCOORDS(NISSUES-1),EIGENVALUES(NISSUES-1)
!
      CHARACTER*21 CAND(NISSUES)
      CHARACTER*21 KNAM(NISSUES)
      INTEGER, ALLOCATABLE :: KDATA(:,:)
      INTEGER, ALLOCATABLE :: KDEGO(:)
      DOUBLE PRECISION, ALLOCATABLE :: A(:)
      DOUBLE PRECISION, ALLOCATABLE :: ZZ(:)
      DOUBLE PRECISION, ALLOCATABLE :: XMISS(:)
      DOUBLE PRECISION, ALLOCATABLE :: B(:)
      DOUBLE PRECISION, ALLOCATABLE :: D(:)
      DOUBLE PRECISION, ALLOCATABLE :: WK(:)
      DOUBLE PRECISION, ALLOCATABLE :: ZSUM(:)
      DOUBLE PRECISION, ALLOCATABLE :: AWV(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: XDATA(:,:)
!
      ALLOCATE(KDATA(NRESPONDENTS,2))
      ALLOCATE(KDEGO(3*NISSUES))
      ALLOCATE(A(NISSUES-1*NISSUES-1))
      ALLOCATE(ZZ(3*NISSUES))
      ALLOCATE(XMISS(NMISSING))
      ALLOCATE(B(3*NISSUES))
      ALLOCATE(D(4*NISSUES))
      ALLOCATE(WK(4*NISSUES))
      ALLOCATE(ZSUM(3*NISSUES))
      ALLOCATE(AWV(NISSUES-1,NISSUES-1))
      ALLOCATE(XDATA(NRESPONDENTS,NISSUES))

  301 FORMAT(8F9.4/8F9.4/8F9.4)
  303 FORMAT(4X,8(A13,4X)/4x,8(A13,4x)/4x,8(A13,4x))
!
      LWORK=4*NISSUES
!
      NQ=NISSUES-1
      NRESP=NSELFPOS
      NMISS=NMISSING
      NLEFT=POLARITY
      IF(NLEFT.EQ.0)NLEFT=1
      IF(NMISS.GT.0)THEN
         DO 89 JJ=1,NMISS
         XMISS(JJ)=KMISS(JJ)
  89     CONTINUE
      ENDIF
      NNQ=NQ+1
      IF(NRESP.EQ.0)NNQ=NQ
      DO 1 I=1,NNQ
      KDEGO(I)=0
  1   KNAM(I)=CAND(I)
!
!
!      open(21,file='VIET68.CTL',status='old')
!      call gettim(itim1,itim2,itim3,itim4)
!      write(24,3)
  3   FORMAT(1X,78('*'))
!
      NP=0
      XNQ=NQ
      XKNTL=.01
      KK=1
      KDEGN=0
      KDEGM=0
      DO 27 K=1,NQ
      ZSUM(K)=0.0
      DO 27 L=1,NQ
!      WV(K,L)=0.0
      awv(k,l)=0.0
!      A(KK)=0.0
!      Q(KK)=0.0
  27  KK=KK+1
      I=0
      ZZSUM=0.0
  9   CONTINUE
      I=I+1
      IF(I.GT.NRESPONDENTS)GO TO 10
      DO 88 JJ=1,NNQ
      A(JJ)=KISSUE(JJ+(I-1)*NNQ)
      I1=MID(I)
  88  CONTINUE
!      READ(1,100)IDRESP(I),                                          
!          (RESPDATA(J+(I-1)*(NSTIMULI+1)),J=1,NSTIMULI+1)
!      READ(1,KFMT,END=10)I1,(A(J),J=1,NNQ)
      SUM=0.0
      KXX=1
      DO 22 JJ=1,NNQ
      B(JJ)=A(JJ)
      IF(NRESP.EQ.1)KXX=2
      IF(NRESP.EQ.JJ)GO TO 22
      SUM=SUM+ABS(A(JJ)-A(KXX))
  22  CONTINUE
  678 FORMAT(2I5,F9.4,10F4.1)
      IF(SUM.GT.XKNTL)GO TO 99
!      WRITE(24,678)I,KXX,SUM,(A(JX),JX=1,NNQ)
      IF(NMISS.EQ.0)KDEGN=KDEGN+1
      IF(NMISS.EQ.0)GO TO 9
      DO 98 II=1,NMISS
      IF(A(KXX).EQ.XMISS(II))KDEGM=KDEGM+1
      IF(A(KXX).EQ.XMISS(II))GO TO 9
  98  CONTINUE
      KDEGN=KDEGN+1
      GO TO 9
  99  IF(NMISS.EQ.0)GO TO 997
      KERR=0
      DO 100 II=1,NMISS
      KK=0
      DO 100 K=1,NNQ
      IF(K.EQ.NRESP)GO TO 100
      KK=KK+1
      AA=ABS(A(K)-XMISS(II))
      IF(AA.LE.XKNTL)KDEGO(KK)=KDEGO(KK)+1
      IF(AA.LE.XKNTL)KERR=1
  100 CONTINUE
  997 CONTINUE
      IF(KERR.EQ.1)GO TO 9
!      WRITE(20)I,I1,(A(JJ),JJ=1,NNQ)
      NP=NP+1
      KDATA(NP,1)=I
      KDATA(NP,2)=I1
      DO 617 JJ=1,NNQ
      XDATA(NP,JJ)=A(JJ)
  617 CONTINUE
      KK=0
      DO 60 II=1,NNQ
      IF(II.EQ.NRESP)GO TO 60
      KK=KK+1
      A(KK)=B(II)
  60  CONTINUE
      ASUM=0.0
      BSUM=0.0
      DO 4 II=1,NQ
      ASUM=ASUM+A(II)
      BSUM=BSUM+A(II)*A(II)
      ZSUM(II)=ZSUM(II)+A(II)
   4  CONTINUE
      ZZSUM=ZZSUM+ASUM
      DET=1.0/(XNQ*BSUM-ASUM*ASUM)
      DO 7 II=1,NQ
      DO 8 J=1,II
      aWV(II,J)=aWV(II,J)+DET*(BSUM-A(II)*ASUM+XNQ*A(II)*            &
                       a(j)-a(j)*asum)
   8  CONTINUE
  7   CONTINUE
      GO TO 9
  10  CONTINUE
      I=I-1
      XNP=NP
!      write(24,12)
!      write(*,12)
  12  FORMAT(' NUMBER OF CASES')
!      write(*,13) NP,I
!      write(24,13) NP,I
  13  FORMAT(2I10)
!      write(24,97)
!      write(*,97)
  97  FORMAT(' MISSING CASES')
!      write(24,96)
!      write(*,96)
  96  FORMAT(' DEGENERATE--MISS  NOMISS')
!      write(24,13)KDEGM,KDEGN
!      write(*,13)KDEGM,KDEGN
!      write(*,333)(KNAM(I),I=1,NQ)
!      write(24,333)(KNAM(I),I=1,NQ)
!      write(*,94)(KDEGO(J),J=1,NQ)
!      write(24,94)(KDEGO(J),J=1,NQ)
  94  FORMAT(1X,20I7)
  333 FORMAT(3X,20(A13,2X))
!      write(24,3)
      NN=I
      KK=1
      BTSUM=0.0
      ZZSUM=ZZSUM/(XNQ*XNP)
      DO 14 I=1,NQ
      aWV(I,I)=aWV(I,I)-XNP
      ZSUM(I)=(ZSUM(I)/XNP)-ZZSUM
      BTSUM=BTSUM+ZSUM(I)**2
      DO 15 J=1,I
      A(KK)=aWV(I,J)
      KK=KK+1
  15  CONTINUE
  14  CONTINUE
      BTSUM=SQRT(BTSUM)
!      write(24,30)
  30  FORMAT(' X(XPRIMEX)-1XPRIME-NI MATRIX')
      DO 24 I=1,NQ
      ZSUM(I)=ZSUM(I)/BTSUM
      K7=I+(I-1)*(I-2)/2
      K8=I+I*(I-1)/2
!      write(24,25)(A(II),II=K7,K8)
  25  FORMAT(60F12.4)
  24  CONTINUE
!      write(24,3)
      CALL DSYEV('V','L', NQ, AWV, NQ, D, WK, LWORK, INFO )
!      call rs(20,nq,awv,d,1,wv,wk,a,ier)
!      write(*,833)INFO
!      write(24,833)INFO
  833 FORMAT(' PERFORMANCE INDEX=',I5)
!      write(*,16)
!      write(24,16)
  16  FORMAT(' EIGENVALUES')
      DO 17 I=1,NQ
!
!  TRANSFER EIGENVALUES HERE
!
      EIGENVALUES(I)=D(I)
!      write(24,18) D(I)
!      write(*,18) D(I)
  18  FORMAT(F12.4)
  17  CONTINUE
      XROOT=D(NQ-1)
!      write(24,3)
      A1=XNP*XROOT*(-1.0)
      A2=XNQ*((XNP+XROOT)**2)
      XSIGM=A1/A2
!
!  TRANSFER CORRECTED GOODNESS OF FIT HERE
!
      FITS(1)=XSIGM
      XWV=AWV(NLEFT,NQ-1)
      DO 42 J=1,NQ
  42  IF(XWV.GT.0.0)AWV(J,NQ-1)=-AWV(J,NQ-1)
!      write(24,3)
!      write(*,19)
!      write(24,19)
  19  FORMAT(' STIMULUS COORDINATES')
!      write(*,303)(KNAM(I),I=1,NQ)
!      write(24,303)(KNAM(I),I=1,NQ)
!      write(24,301)(AWV(J,NQ-1),J=1,NQ)
!      write(*,301)(AWV(J,NQ-1),J=1,NQ)
!      write(*,199)
!      write(24,199)
  199 FORMAT(' STIMULUS COORDINATES RAW DATA')
!      write(24,301)(ZSUM(J),J=1,NQ)
!      write(*,301)(ZSUM(J),J=1,NQ)
!      write(24,3)
      DO 40 I=1,NQ
!
!  TRANSFER STIMULUS COORDINATES HERE
!
      STIMCOORDS(I)=AWV(I,NQ-1)
      ZZ(I)=AWV(I,NQ-1)*(XNP/(XNP+XROOT))
  40  CONTINUE
!      REWIND 20
      CALL AMREG(NRESPONDENTS,NP,NQ,NRESP,NMISS,ZZ,XMISS,KNAM,BTSUM, &
               ZZSUM,ZSUM,XSIGM,KDATA,XDATA,FITS,PSIMATRIX)
  41  CONTINUE
        DEALLOCATE(KDATA)
        DEALLOCATE(KDEGO)
        DEALLOCATE(A)
        DEALLOCATE(ZZ)
        DEALLOCATE(XMISS)
        DEALLOCATE(B)
        DEALLOCATE(D)
        DEALLOCATE(WK)
        DEALLOCATE(ZSUM)
        DEALLOCATE(AWV)
        DEALLOCATE(XDATA)
!
      EXITSTATUS=1
!
      RETURN
      END
!
!  *****************************************************************
!     SUBROUTINE AMREG -- DOES SIMPLE OLS TO RECOVER THE INDIVIDUAL
!                       TRANSFORMATION
!  *****************************************************************
!
!
      SUBROUTINE AMREG(NRESPONDENTS,NN,NQ,NRESP,NMISS,ZZ,XMISS,KNAM, &
                     BTSUM,ZZSUM,ZTSUM,XSIGM,KDATA,XDATA,FITS,       &
                     PSIMATRIX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ZZ(3*NQ),PSIMATRIX(NRESPONDENTS*4),                 &
                XMISS(NMISS),                                        &
                ZTSUM(3*NQ),FITS(5),                                 &
                KDATA(NRESPONDENTS,2),XDATA(NRESPONDENTS,NQ+1)
      CHARACTER*13 KNAM(NQ+1)
      DOUBLE PRECISION, ALLOCATABLE :: A(:)
      DOUBLE PRECISION, ALLOCATABLE :: XSUM(:)
      DOUBLE PRECISION, ALLOCATABLE :: YSUM(:)
      DOUBLE PRECISION, ALLOCATABLE :: ZSUM(:)
      DOUBLE PRECISION, ALLOCATABLE :: Q(:)
      DOUBLE PRECISION, ALLOCATABLE :: QQ(:)
      DOUBLE PRECISION, ALLOCATABLE :: QQQ(:)
      DOUBLE PRECISION, ALLOCATABLE :: XV(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: YV(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: ZV(:,:)
!
      ALLOCATE(A(NQ+1))
      ALLOCATE(XSUM(2*NQ))
      ALLOCATE(YSUM(NQ+1))
      ALLOCATE(ZSUM(NQ+1))
      ALLOCATE(Q(NN+NQ))
      ALLOCATE(QQ(NN+NQ))
      ALLOCATE(QQQ(NN+NQ))
      ALLOCATE(XV(NQ,2))
      ALLOCATE(YV(NQ,2))
      ALLOCATE(ZV(2,NQ))
!
      NP=0
      XNQ=NQ
      NNQ=NQ+1
      IF(NRESP.EQ.0)NNQ=NQ
      XKNTL=.001
      XMAX=-99.0
      TMAX=-99.0
      ASUM=0.0
      BSUM=0.0
      CSUM=0.0
      DSUM=0.0
      ESUM=0.0
      TSUM=0.0
      TTSUM=0.0
      DO 111 I=1,NQ
      TMAX=AMAX1(TMAX,ABS(ZZ(I)))
      Q(I)=ZZ(I)
      QQ(I)=ZZ(I)
      QQQ(I)=ZZ(I)
      XSUM(I)=0.0
      ZSUM(I)=0.0
  111 XSUM(I+NQ)=0.0
  3   FORMAT(1X,78('*'))
      IK=0
      IPOS=0
      INEG=0
      SSE=0.0
  2   CONTINUE
      IK=IK+1
      IF(IK.GT.NN)GO TO 10
!      READ(20)IC,I1,(A(J),J=1,NNQ)
      IC=KDATA(IK,1)
      I1=KDATA(IK,2)
      DO 617 J=1,NNQ
      A(J)=XDATA(IK,J)
  617 CONTINUE
      IF(NRESP.EQ.0)GO TO 101
      IF(NMISS.EQ.0)GO TO 1011
      DO 100 I=1,NMISS
      AA=ABS(A(NRESP)-XMISS(I))
  100 IF(AA.LE.XKNTL)GO TO 2
 1011 XRES1=A(NRESP)
  101 CONTINUE
      KK=0
      NP=NP+1
      DO 8 I=1,NNQ
      IF(I.EQ.NRESP)GO TO 8
      KK=KK+1
      XV(KK,1)=1.0
      XV(KK,2)=A(I)
      SSE=SSE+(((A(I)-ZZSUM)/BTSUM)-ZTSUM(KK))**2
  8   CONTINUE
      DO 9 K=1,2
      DO 11 J=1,2
      SUM=0.0
      DO 12 I=1,NQ
      SUM=SUM+XV(I,K)*XV(I,J)
  12  CONTINUE
      YV(K,J)=SUM
  11  CONTINUE
  9   CONTINUE
      DET=1.0/(YV(1,1)*YV(2,2)-YV(1,2)*YV(2,1))
      AA=YV(1,1)
      YV(1,1)=YV(2,2)*DET
      YV(2,2)=AA*DET
      YV(1,2)=YV(1,2)*(-1.0)*DET
      YV(2,1)=YV(2,1)*(-1.0)*DET
      DO 13 K=1,NQ
      DO 14 I=1,2
      SUM=0.0
      DO 15 J=1,2
      SUM=SUM+YV(I,J)*XV(K,J)
  15  CONTINUE
      ZV(I,K)=SUM
  14  CONTINUE
  13  CONTINUE
      DO 16 I=1,2
      SUM=0.0
      DO 17 J=1,NQ
      SUM=SUM+ZV(I,J)*ZZ(J)
  17  CONTINUE
      A(I)=SUM
  16  CONTINUE
      SUMA=0.0
      SUMB=0.0
      SUMC=0.0
      SUMD=0.0
      SUME=0.0
      DO 22 J=1,NQ
      AA=A(1)+A(2)*XV(J,2)
      SUMA=SUMA+AA*ZZ(J)
      SUMB=SUMB+AA
      SUMC=SUMC+ZZ(J)
      SUMD=SUMD+AA*AA
  22  SUME=SUME+ZZ(J)*ZZ(J)
      AA=XNQ*SUMA-SUMB*SUMC
      BB=XNQ*SUMD-SUMB*SUMB
      CC=XNQ*SUME-SUMC*SUMC
      rsqrt=aa/sqrt(abs(bb*cc))
      if(a(2).lt.0.0)rsqrt=-rsqrt
      RR=(AA*AA)/(BB*CC)
      IF(NP.GT.1)GO TO 63
!      if(nresp.ne.0)write(23,64)
  64  FORMAT('      LINE #   CASE #  R POS   ALPHA    BETA',4X,      &
               'SCALED POS',4X,'RSQ')
!      if(nresp.eq.0)write(23,644)
  644 FORMAT('      LINE #   CASE #    ALPHA      BETA',6X,'RSQ',5X, &
               'RAW CORR')
  63  CONTINUE
      AA=A(1)+A(2)*XRES1
      TSUM=TSUM+AA
      TTSUM=TTSUM+AA*AA
  31  FORMAT(2I10,F6.1,5f10.4)
  311 FORMAT(2I10,4F10.4)
!      IF(NRESP.NE.0)WRITE(23,31)IC,I1,XRES1,A(1),A(2),AA,RR,rsqrt
!      IF(NRESP.EQ.0)WRITE(23,311)IC,I1,A(1),A(2),RR,RSQRT
      IF(NRESP.NE.0)Q(NQ+NP)=AA
!
!  TRANSFER RESPONDENT COORDINATES
!
         PSIMATRIX(((KDATA(IK,1)-1)*4)+1)=A(1)
         PSIMATRIX(((KDATA(IK,1)-1)*4)+2)=A(2)
         PSIMATRIX(((KDATA(IK,1)-1)*4)+3)=AA
         PSIMATRIX(((KDATA(IK,1)-1)*4)+4)=RR
!
      IF(A(2).GE.0.0)IPOS=IPOS+1
      IF(A(2).GE.0.0)QQ(IPOS+NQ)=AA
      IF(A(2).LT.0.0)INEG=INEG+1
      IF(A(2).LT.0.0)QQQ(INEG+NQ)=AA
      DO 43 II=1,NQ
      AA=A(1)+XV(II,2)*A(2)
      XSUM(II)=XSUM(II)+AA
      XSUM(II+NQ)=XSUM(II+NQ)+AA*AA
      XCAND=ZZ(II)
      ASUM=ASUM+AA
      BSUM=BSUM+XCAND
      CSUM=CSUM+AA*AA
      DSUM=DSUM+XCAND*XCAND
      ESUM=ESUM+AA*XCAND
  43  CONTINUE
      GO TO 2
  10  CONTINUE
      XNP=NP
      SSE=SQRT(SSE/(XNP*XNQ))
      YSUM(NQ+1)=SQRT((TTSUM-(TSUM*TSUM)/XNP)/(XNP))
      DO 200 I=1,NQ
      ZSUM(I)=XSUM(I)/XNP
  200 YSUM(I)=SQRT((XSUM(I+NQ)-(XSUM(I)*XSUM(I))/XNP)/(XNP))
      ZSUM(NQ+1)=TSUM/XNP
!      write(*,50)
!      write(24,50)
  50  FORMAT(' CORRECTED GOODNESS OF FIT AND RAW FIT')
!      write(*,18)XSIGM,SSE
!      write(24,18)XSIGM,SSE
  18  FORMAT(2F12.4)
!      write(24,3)
!      write(*,19)
!      write(24,19)
  19  FORMAT(' TOTAL CASES--NUMBER POSITIVE  NUMBER NEGATIVE')
!      write(*,20) NP,IPOS,INEG
!      write(24,20) NP,IPOS,INEG
  20  FORMAT(1X,I8,7X,I8,7X,I10)
!
!  TRANSFER CASE NUMBERS HERE
!
      FITS(3)=NP
      FITS(4)=IPOS
      FITS(5)=INEG
!
!      write(24,3)
      RA=XNP*ESUM-ASUM*BSUM
      RB=XNP*CSUM-ASUM*ASUM
      RC=XNP*DSUM-BSUM*BSUM
      R1=RA/(SQRT(RB*RC))
      R1S=R1*R1
!
!  TRANSFER R-SQUARE HERE
!
      FITS(2)=R1S
!
!      write(*,40)
!      write(24,40)
  40  FORMAT(' R AND R-SQUARE')
!      write(*,41) R1,R1S
!      write(24,41) R1,R1S
  41  FORMAT(2F10.4)
!      write(*,302)
!      write(24,302)
  302 FORMAT(' STIMULUS COORDINATES FROM COMPUTED BETA WEIGHTS')
!      write(*,303)(KNAM(I),I=1,NNQ)
!      write(24,303)(KNAM(I),I=1,NNQ)
  303 FORMAT(4X,8(A13,4X)/4x,8(a13,4x)/4x,8(a13,4x))
!      write(*,301)(ZSUM(J),J=1,NNQ)
!      write(24,301)(ZSUM(J),J=1,NNQ)
!      write(*,300)
!      write(24,300)
  300 FORMAT(' STANDARD DEVIATIONS')
!      write(*,301)(YSUM(J),J=1,NNQ)
!      write(24,301)(YSUM(J),J=1,NNQ)
  301 FORMAT(8F9.4/8F9.4/8F9.4)
!      write(24,3)
!
!
        DEALLOCATE(A)
        DEALLOCATE(XSUM)
        DEALLOCATE(YSUM)
        DEALLOCATE(ZSUM)
        DEALLOCATE(Q)
        DEALLOCATE(QQ)
        DEALLOCATE(QQQ)
        DEALLOCATE(XV)
        DEALLOCATE(YV)
        DEALLOCATE(ZV)
      RETURN
      END
