!
! ***********************************************************************
!
!  MLSMU6_2010.F95----METRIC LEAST SQUARES MULTIDIMENSIONAL UNFOLDING
!                      
!                     PROGRAM WRITTEN AT THE UNIVERSITY OF OREGON, 1979
!
!                      FORTRAN 95: 8 NOVEMBER 2010
!                 
! ***********************************************************************
!
!
!  ***METRIC LEAST SQUARES MULTIDIMENSIONAL UNFOLDING***
!
!
!  *******WARNING: INPUT MUST BE TRANSFORMABLE TO SQUARED DISTANCES*******
!
!
!
!  ***CARD 1, UNIT 21, FORMAT(5I5)***
!  NQ= # OF STIMULI
!  NR= # OF DIMENSIONS(LOWER LIMIT)
!  NRMAX= # OF DIMENSIONS(UPPER LIMIT)
!  ******NOTE: PROGRAM COMPUTES SOLUTION FOR ALL DIMENSIONS
!              NR TO NRMAX INCLUSIVE
!  KJEP2=# OF ITERATIONS IN THE METRIC GRADIENT ROUTINE
!  NSYM=1 IF A SYMMETRIC DISTANCE MATRIX IS TO BE ANALYZED
!  NSYM=0 IF A NONSYMMETRIC DISTANCE MATRIX IS TO BE ANALYZED
!  NN=NP IF KREAD=3 AND THE STIMULI COORDINATES THAT ARE READ IN ARE TO
!        BE HELD FIXED
!  NN=NQ IF KREAD=2 AND THE INDIVIDUAL COORDINATES THAT ARE READ IN ARE
!        TO BE HELD FIXED
!  NN=0  OTHERWISE
!
!  ***CARD 2, UNIT 21, FORMAT(5I5)***
!  NPRR=1 IF NO OUTPUT OF INDIVIDUAL COORDINATES TO THE
!         TERMINAL IS DESIRED
!  NPRR=0 IF OUTPUT TO TERMINAL OF INDIVIDUAL COORDINATES IS DESIRED
!  KREAD=1  NO STARTING CONFIGURATION WILL BE READ IN
!  KREAD=0  STARTING CONFIGURATION WILL BE READ IN
!  KREAD=2  INDIVIDUAL COORDINATES ONLY ARE READ IN.  IF THESE COORD-
!           INATES ARE TO BE HELD FIXED, SET NN=NQ (SEE ABOVE)
!  KREAD=3  STIMULI COORDINATES ONLY ARE READ IN.  IF THESE COORDINATES
!           ARE TO BE HELD FIXED, SET NN=NP (SEE ABOVE)
!  NOTIT=1  NO TITLES FOR THE INDIVIDUALS WILL BE READ IN
!  NOTIT=0  TITLES FOR THE INDIVIDUALS WILL BE READ IN
!  KTIT= # OF LETTERS IN INDIVIDUAL TITLES, MAXIMUM OF 13 LETTERS,
!        WILL BE READ IN A1 FORMAT
!  NMISS= # OF MISSING DATA VALUES (MAXIMUM OF 6 IN INPUT DISTANCE MATRIX)
!  NMISS=0  NO MISSING DATA
! 
!  ***CARD 3, UNIT 21, FORMAT(7F8.3)
!  TOL=  TOLERANCE TO STOP THE ITERATION IN focus SUBROUTINE,
!        NORMALLY IS SET TO 0.001 OR 0.010
!  AMX,BMX,CMX ARE PARAMETERS TO TRANSFORM THE INPUT DATA INTO SQUARED
!              DISTANCES.  THE FOLLOWING EQUATION IS USED:  
!                     D-SQUARED=(AMX*T+BMX)**CMX
!              WHERE T=INPUT DATA.  THIS FORMULA TAKES A LINEAR TRANS-
!              FORMATION OF THE INPUT DATA TO THE POWER CMX.  FOR EXAMPLE,
!              IF THE DATA, T, ARE DISTANCES, SET AMX=1.0, BMX=0.0, CMX=2.0;
!              OR IF THE DATA ARE THERMOMETER RATINGS, SET AMX=-.02, BMX=2.0,
!              CMX=2.0, WHICH IS EQUIVALENT TO SUBTRACTING THE THERMOMETER
!              SCORE FROM 100 AND DIVIDING BY 5O AND THEN SQUARING.  THIS
!              CONVERTS T FROM A 0-100 SCALE TO A 4-0 SCALE.  IF THE DATA
!              ARE CORRELATIONS, SET AMX=-1.0, BMX=1.0, AND CMX=2.0 OR
!              1.0 IF THE CORRELATIONS ARE INITIALLY REGARDED AS UNSQUARED
!              OR SQUARED DISTANCES RESPECTIVELY.
!  XMAX= MAXIMUM ABSOLUTE EXPECTED CORRDINATE VALUE ON ANY DIMENSION.
!        IT IS USED FOR PLOTTING PURPOSES.  IF THE SQUARED DISTANCES
!        ARE CONFINED TO A 4-0 SCALE, XMAX=1.5 IS USUALLY SUFFICIENT.
!  DMIN= MINIMUM EXPECTED VALUE OF INPUT DATA
!  DMAX= MAXIMUM EXPECTED VALUE OF INPUT DATA; DMIN AND DMAX ARE USED
!        TO CATCH CODING ERRORS IN THE INPUT DATA.  ANYTHING OUT OF
!        RANGE IS TREATED AS MISSING DATA.
!
!  ***CARD 4, UNIT 21, FORMAT(20A5)***
!  KFMT= THE FORMAT OF THE INPUT SQUARED DISTANCE MATRIX
!
!  ***OPTIONAL READ OF STARTING COORDINATES***
!  ***CARD 5, UNIT 21, FORMAT(20A5)***       
!  KFMT1= THE FORMAT OF THE STARTING CONFIGURATION
!
!  ***CARD 6 OR 5, UNIT 21, FORMAT(NMISSF3.0)***
!  XMISS= THE ARRAY CONTAINING THE MISSING DATA VALUES
!         (FORMAT MUST BE F3.0)
!
!  ***CARDS 7 TO NQ+7 OR 6 TO NQ+6, UNIT 21, FORMAT(13A1)***
!  LABELS FOR THE STIMULI MUST BE READ IN--# LABELS=NQ,
!     MAXIMUM LENGTH OF 13 LETTERS
!
!  LABELS FOR THE INDIVIDUALS ARE OPTIONAL AND ARE LIMITED
!     TO A MAXIMUM OF 13 LETTERS--IF LESS THAN 13 LETTERS
!     KFMT MUST BE FORMATED APPROPRIATELY AND THIS MUST BE
!     PRIOR TO THE FORMAT FOR THE DISTANCES.
!     ALL THIS INFORMATION IS CONTAINED IN ONE FORMAT--KFMT
!
!    UNIT NUMBERS:  1 IS USED AS SCRATCH
!                  20 IS USED TO READ IN STARTING COORDINATES
!                  21 IS THE CONTROL FILE
!                  22 IS USED TO WRITE OUT THE FINAL COORDINATES
!                  23 IS USED TO READ IN THE SQUARED DISTANCE MATRIX
!                  24 IS USED AS SCRATCH
!
!
!
!
!  *********SAMPLE CONTROL INPUT*********
!     20    1    2   10    0    0
!      1    1    1    0    4    3
!      .001  -0.02    2.0     2.0     1.5     0.0   100.0
!  (8X,4A1,318X,16F3.0,1089X,4F3.0)
!  997998999
!  CARTER
!  REAGAN
!  KENNEDY
!  CONNALLY
!  FORD
!  BROWN
!  BAKER
!  MONDALE
!  BUSH
!  ANDERSON
!  LUCEY
!  MCGOVERN
!  WALLACE
!  NIXON
!  DEM PTY
!  REP PTY
!  CARTER
!  REAGAN
!  KENNEDY
!  ANDERSON
!
!
!
!   **********ARRAY SIZES**********
!   DIST=NP*NQ
!   XX AND XXX=NP X NRMAX
!   ZZ AND ZZZ=NQ X NRMAX
!   A=NQ*(NQ+1)/2
!   XROW=NP
!   XCOL=NQ
!   Q=NP+NQ+2
!   RSQ=NP+NQ+2
!   SSQ=NP+NQ+2
!   QQ=NQ*NQ
!   D=NQ
!   WV=NQ X NQ
!   DAT=NQ X 7
!
! ******************************** 
!            ARRAY SIZES BASED ON NQ = 143, NP = 85,000
!                          NP*NQ = 1,900,999
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DIST(1900999),A(10300),Q(85143),RSQ(85143),&
                ZZ(143,4),XX(85000,4),DAT(100,7),D(143),WV(143,143),&
                XMISS(6),XCOL(143),XROW(85000),QQ(85143),&
                SSQ(85143),TT(4,4),YCENTER(143,143),FV1(4*143),&
                FV2(4*143)
      REAL RESULT,RESULT1
      CHARACTER*1 LTIT,TIT(13)
      CHARACTER*5 KFMT(60),KFMT1(60)
      CHARACTER*64 FNAME
!
      OPEN(21,FILE='UNFOLD_1968.CTL',STATUS='OLD')
      OPEN(3,FILE='UNFOLD3.DAT')
      OPEN(5,FILE='UNFOLD5.DAT')
!      OPEN(1)
!      OPEN(24)
 1000 FORMAT(' START TIME OF JOB    ',3F12.7)
 1001 FORMAT(' END TIME OF JOB      ',3F12.7)
 1002 FORMAT(' ELAPSED TIME OF JOB  ',3F12.7)
!
      call cpu_time(result)
      write(*,1000)result
!
!
      LTIT=' '
      read(21,1082)FNAME 
      write(*,1082)FNAME
      OPEN(unit=23,FILE=FNAME,STATUS='OLD')
 1082 format(a64)
  801 FORMAT(10I5)
      READ(21,803) NQ,NR,NRMAX,KJEP2,NSYM,NN
  803 FORMAT(10I5)
      READ(21,803)NPRR,KREAD,NOTIT,KTIT,NMISS
  804 FORMAT(10I5)
      READ(21,806)TOL,AMX,BMX,CMX,XMAX,DMIN,DMAX
      READ(21,510) KFMT
      IF(KREAD.EQ.1)GO TO 290
      READ(21,510)KFMT1
  290 IF(NMISS.EQ.0)GO TO 260
      READ(21,805)(XMISS(I),I=1,NMISS)
  260 CONTINUE
  510 FORMAT(60A5)
      WRITE(3,3)
      WRITE(3,808)
      WRITE(5,809)
  808 FORMAT(46X,'INPUT PARAMETERS AS READ FROM UNIT 21')
  809 FORMAT(' INPUT PARAMETERS AS READ FROM UNIT 21')
      WRITE(3,810)
      WRITE(5,811)
  810 FORMAT(53X,'NQ',3X,'NR',3X,'NRMAX',1X,'KJEP2',1X,'NSYM',2X,&
                                         'NN')
  811 FORMAT(3X,'NQ',3X,'NR',3X,'NRMAX',1X,'KJEP2',1X,'NSYM',2X,'NN')
      WRITE(3,8034)NQ,NR,NRMAX,KJEP2,NSYM,NN
      WRITE(5,8035)NQ,NR,NRMAX,KJEP2,NSYM,NN
      WRITE(3,812)
      WRITE(5,813)
 8033 FORMAT(50X,8I5)
  812 FORMAT(52X,'NPRR',1X,'KREAD',1X,'NOTIT',1X,'KTIT',1X,'NMISS')
  813 FORMAT(2X,'NPRR',1X,'KREAD',1X,'NOTIT',1X,'KTIT',1X,'NMISS')
      WRITE(3,8034)NPRR,KREAD,NOTIT,KTIT,NMISS
      WRITE(5,8035)NPRR,KREAD,NOTIT,KTIT,NMISS
 8034 FORMAT(50X,I5,2I6,I5,2I6)
 8035 FORMAT(I5,2I6,I5,2I6)
      WRITE(3,814)
      WRITE(5,815)
  814 FORMAT(40X,'TOL',5X,'AMX',5X,'BMX',5X,'CMX',4X,'XMAX',4X,'DMIN',4X&
      ,'DMAX')
  815 FORMAT(3X,'TOL',5X,'AMX',5X,'BMX',5X,'CMX',4X,'XMAX',4X,'DMIN',4X,&
      'DMAX')
      WRITE(3,8036)TOL,AMX,BMX,CMX,XMAX,DMIN,DMAX
 8036 FORMAT(37X,7F8.3)
      WRITE(5,806)TOL,AMX,BMX,CMX,XMAX,DMIN,DMAX
  3   FORMAT(1X,130('*'))
  333 FORMAT('1')
  66  FORMAT(6F10.4)
  6   FORMAT(6F15.8)
  37  FORMAT(41X,4F15.8)
  377 FORMAT(4F10.4)
  96  FORMAT(I8,6F15.8)
  966 FORMAT(I8,6F10.4)
  59  FORMAT(1F15.8)
  155 FORMAT(2F15.8)
  250 FORMAT(13A1)
  251 FORMAT(1X,13A1,9F15.8)
 2551 FORMAT(42X,13A1,9F15.8)
 2511 FORMAT(11X,13A1,9F15.8)
  252 FORMAT(1X,13A1,9F10.4)
 2522 FORMAT(1X,13A1,9F9.4)
  255 FORMAT(1X,I13,9F15.8)
  256 FORMAT(1X,I13,9F10.4)
  805 FORMAT(6F3.0)
  806 FORMAT(10F8.3)
      INR=NR
      XNQ=NQ
      XNP=NP
      NR1=NR+1
      NR2=2*NR
      XNR=NR
      NNQ=NQ+1
      KK=0
      KKK=0
      I=1
  13  IF(NOTIT.EQ.1)READ(23,KFMT,END=1333)(DIST((I-1)*NQ+J),J=1,NQ)
      IF(NOTIT.EQ.0)READ(23,KFMT,END=1333)(TIT(J),J=1,KTIT),&
                      (DIST((I-1)*NQ+J),J=1,NQ)
      KK=KK+1
      WRITE(*,4444)KK
 4444 FORMAT(' HELLO',I6)
      KMISS=0
      DO 133 J=1,NQ
      IF(NMISS.EQ.0)GO TO 1234
      IF(DIST((I-1)*NQ+J).LT.DMIN.OR.DIST((I-1)*NQ+J).GT.DMAX)DIST((I-1)&
      *NQ+J)=XMISS(1)
      DO 134 JJ=1,NMISS
      IF(DIST((I-1)*NQ+J).EQ.XMISS(JJ))KMISS=KMISS+1
  134 IF(DIST((I-1)*NQ+J).EQ.XMISS(JJ))GO TO 133
 1234 DIST((I-1)*NQ+J)=(AMX*DIST((I-1)*NQ+J)+BMX)**CMX
  133 CONTINUE
  886 FORMAT(1X,5A1,3F5.1)
      KKK=KKK+1
      IF(KMISS.GE.(NQ-NRMAX))GO TO 13
      IF(NOTIT.EQ.0)WRITE(1)(TIT(J),J=1,KTIT)
      I=I+1
      GO TO 13
 1333 NP=I-1
      XNP=NP
      WRITE(3,3)
      WRITE(5,1334)
      WRITE(3,1335)
 1335 FORMAT(31X,'INITIAL READ OF CASES--TOTAL CASES SELECTED CASES&
       NONMISSING CASES')
 1334 FORMAT(' INITIAL READ OF CASES--TOTAL CASES  SELECTED CASES  NONMI&
      SSING CASES')
      WRITE(3,8134)KK,KKK,NP
      WRITE(*,8133)KK,KKK,NP
      WRITE(5,8133)KK,KKK,NP
 8133 FORMAT(27X,I5,9X,I5,12X,I5)
 8134 FORMAT(58X,I5,9X,I5,12X,I5)
      IF(NMISS.EQ.0)GO TO 280
      SUM=0.0
      DO 261 I=1,NP
      XROW(I)=0.0
      IF(I.GT.NQ)GO TO 270
      XCOL(I)=0.0
  270 DO 261 J=1,NQ
      DO 262 K=1,NMISS
  262 IF(DIST((I-1)*NQ+J).EQ.XMISS(K))DIST((I-1)*NQ+J)=999.0
  261 IF(DIST((I-1)*NQ+J).NE.999.0)SUM=SUM+DIST((I-1)*NQ+J)
      SUM=SUM/(XNP*XNQ)
      XNPQT=0.0
      NPQT=0
      DO 263 I=1,NP
      DO 263 J=1,NQ
      IF(DIST((I-1)*NQ+J).NE.999.0)XCOL(J)=XCOL(J)+1.0
      IF(DIST((I-1)*NQ+J).NE.999.0)XROW(I)=XROW(I)+1.0
      IF(DIST((I-1)*NQ+J).NE.999.0)NPQT=NPQT+1
  263 IF(DIST((I-1)*NQ+J).NE.999.0)XNPQT=XNPQT+1.0
      IF(NMISS.NE.0)GO TO 283
  280 CONTINUE
      DO 281 I=1,NP
  281 XROW(I)=XNQ
      DO 282 J=1,NQ
  282 XCOL(J)=XNP
      XNPQT=XNP*XNQ
      NPQT=NP*NQ
  283 CONTINUE
      IF(NQ.EQ.1)GO TO 350
      SUMB=0.0
      DO 205 I=1,NP
      SUMA=0.0
      DO 206 J=1,NQ
      IF(DIST((I-1)*NQ+J).EQ.999.0)GO TO 650
      SUMA=SUMA+DIST((I-1)*NQ+J)
      GO TO 206
  650 SUMA=SUMA+SUM
  206 CONTINUE
      Q(I)=SUMA/XNQ
  205 SUMB=SUMB+Q(I)
      SUMB=SUMB/XNP
      DO 207 J=1,NQ
      SUMA=0.0
      DO 208 I=1,NP
      IF(DIST((I-1)*NQ+J).EQ.999.0)GO TO 651
      SUMA=SUMA+DIST((I-1)*NQ+J)
      GO TO 208
  651 SUMA=SUMA+SUM
  208 CONTINUE
  207 Q(NP+J)=SUMA/XNP
      KK=1
      DO 211 I=1,NQ
      DO 211 J=1,I
      SUM=0.0
      DO 213 K=1,NP
      IF(DIST((K-1)*NQ+J).EQ.999.0)GO TO 652
      AA=(DIST((K-1)*NQ+J)-Q(K)-Q(J+NP)+SUMB)/(-2.0)
      GO TO 653
  652 AA=(2.0*SUMB-Q(K)-Q(J+NP))/(-2.0)
  653 IF(DIST((K-1)*NQ+I).EQ.999.0)GO TO 654
      BB=(DIST((K-1)*NQ+I)-Q(K)-Q(I+NP)+SUMB)/(-2.0)
      GO TO 213
  654 BB=(2.0*SUMB-Q(K)-Q(I+NP))/(-2.0)
  213 SUM=SUM+AA*BB
      A(KK)=SUM
      YCENTER(I,J)=SUM
      YCENTER(J,I)=SUM
  211 KK=KK+1
      WRITE(*,1947)KK
 1947 FORMAT(' HELLO BEAVIS',I8)
!
!  EIGENVECTOR-EIGENVALUE DECOMPOSITION DOUBLE-CENTERED AGREEMENT
!     SCORE MATRIX
!
      LWORK=4*143
      CALL DSYEV('V','L',NQ,YCENTER,143,FV1,FV2,LWORK,INFO)
!
      WRITE(*,1948)INFO
!      call rs(143,NQ,YCENTER,D,1,WV,FV1,FV2,IER)
      DO 26 I=1,NQ
      WRITE(*,1999)I,FV1(I)
 1999 FORMAT(I5,F12.4)
      A(I)=SQRT(ABS(FV1(NQ-I+1)))
  26  QQ(I)=SQRT(A(I))
      DO 500 J=1,NRMAX
      DO 500 I=1,NQ
      ZZ(I,J)=YCENTER(I,NQ-J+1)*QQ(J)
  500 IF(NSYM.EQ.1)XX(I,J)=ZZ(I,J)
 1948 FORMAT(' PERFORMANCE INDEX=',I5)
      IF(NSYM.EQ.1)GO TO 350
      DO 214 I=1,NP
      DO 214 K=1,NRMAX
      SUM=0.0
      DO 215 J=1,NQ
      IF(DIST((I-1)*NQ+J).EQ.999.0)GO TO 655
      AA=(DIST((I-1)*NQ+J)-Q(I)-Q(J+NP)+SUMB)/(-2.0)
      GO TO 215
  655 AA=SUMB
  215 SUM=SUM+AA*YCENTER(J,NQ-K+1)*(1.0/QQ(K))
      XX(I,K)=SUM
  214 CONTINUE
  350 IF(KREAD.EQ.1)GO TO 608
      IF(KREAD.EQ.2)GO TO 612
      OPEN(20,FILE='NES2004_10_THERM_Z.DAT',STATUS='OLD')
      DO 609 I=1,NQ
  614 READ(20,KFMT1,END=6122)(ZZ(I,J),J=1,NR)
      IF(ZZ(I,1).EQ.9.0)GO TO 614
  609 CONTINUE
 6122 IF(KREAD.EQ.0)GO TO 612
      DO 400 J=1,NR
      DO 400 K=1,NR
      SUM=0.0
      DO 401 I=1,NQ
  401 SUM=SUM+(1.0/QQ(J))*YCENTER(I,NQ-J+1)*ZZ(I,K)
  400 TT(J,K)=SUM
      DO 402 I=1,NP
      DO 403 J=1,NR
      SUM=0.0
      DO 404 K=1,NR
  404 SUM=SUM+XX(I,K)*TT(K,J)
  403 SSQ(J)=SUM
      DO 402 K=1,NR
  402 XX(I,K)=SSQ(K)
      IF(KREAD.EQ.3)GO TO 608
  612 CONTINUE
      DO 610 I=1,NP
  613 READ(20,KFMT1,END=6088)(XX(I,J),J=1,NR)
      IF(XX(I,1).EQ.9.0)GO TO 613
  610 CONTINUE
 6088 IF(KREAD.EQ.0)GO TO 608
      DO 405 KK=1,NR
      DO 405 K=1,NR
      ASUM=0.0
      DO 407 I=1,NP
      SUM=0.0
      DO 406 J=1,NQ
      IF(DIST((I-1)*NQ+J).EQ.999.0)GO TO 410
      AA=(DIST((I-1)*NQ+J)-Q(I)-Q(J+NP)+SUMB)/(-2.0)
      GO TO 406
  410 AA=SUMB
  406 SUM=SUM+AA*WV(J,NQ-K+1)*(1.0/A(K))
  407 ASUM=ASUM+(1.0/QQ(K))*SUM*XX(I,KK)
  405 TT(K,KK)=ASUM
      DO 408 I=1,NQ
      DO 409 J=1,NR
      SUM=0.0
      DO 411 K=1,NR
  411 SUM=SUM+ZZ(I,K)*TT(K,J)
  409 SSQ(J)=SUM
      DO 408 K=1,NR
  408 ZZ(I,K)=SSQ(K)
  608 CONTINUE
  611 CONTINUE
  999 REWIND 1
      REWIND 24
      WRITE(5,336)
      WRITE(3,3336)
 3336 FORMAT(31X,'NUMBER OF INDIVIDUALS  NUMBER OF STIMULI NUMBER OF&
       DIMENSIONS')
  336 FORMAT(' NUMBER OF INDIVIDUALS  NUMBER OF STIMULI  NUMBER OF DIMEN&
      SIONS')
      WRITE(5,337) NP,NQ,NR
      WRITE(3,3337)NP,NQ,NR
  337 FORMAT(I13,I19,I21)
 3337 FORMAT(39X,I5,14X,I5,17X,I5)
      WRITE(3,3)
      WRITE(3,331)
      WRITE(5,31)
  331 FORMAT(41X,'STARTING ESTIMATE OF STIMULUS COORDINATES')
  31  FORMAT(' STARTING ESTIMATE OF STIMULUS COORDINATES')
      DO 32 I=1,NQ
      IF(NR.GT.INR)READ(24)(TIT(J),J=1,13)
      IF(NR.GT.INR)GO TO 998
      READ(21,250)(TIT(J),J=1,13)
      WRITE(24)(TIT(J),J=1,13)
  998 CONTINUE
      WRITE(3,2551)(TIT(J),J=1,13),(ZZ(I,J),J=1,NR)
      WRITE(5,252)(TIT(J),J=1,13),(ZZ(I,J),J=1,NR)
  32  CONTINUE
      WRITE(3,3)
      IF(NQ.EQ.1)GO TO 220
      IF(NR.GT.INR)GO TO 220
      SUM=0.0
      DO 33 I=1,NQ
  33  SUM=SUM+A(I)
      SUM1=0.0
      DO 34 I=1,NQ
      Q(I)=(A(I)/SUM)*100.0
      SUM1=SUM1+Q(I)
  34  Q(I+NQ)=SUM1
      WRITE(3,1100)INFO
 1100 FORMAT(41X,'PERFORMANCE INDEX=',I5)
      WRITE(5,1100)INFO
      WRITE(5,35)
      WRITE(3,335)
  35  FORMAT(' S.V.D. ROOTS OF DOUBLE CENTERED DISTANCE MATRIX')
 335  FORMAT(41X,'S.V.D. ROOTS OF DOUBLE CENTERED DISTANCE MATRIX')
      DO 36 I=1,NQ
      WRITE(3,37) A(I),Q(I),Q(I+NQ)
  36  WRITE(5,377)A(I),Q(I),Q(I+NQ)
      WRITE(3,3)
  220 CALL focus(NP,NQ,NR,NN,KJEP2,NSYM,TOL,XNPQT,XCOL,XROW,DIST,XX,ZZ,&
                        AMX,BMX,CMX)
      CALL STATKP(NP,NQ,NR,99,XNPQT,XCOL,XROW,RSQ,SSQ,DIST,XX,ZZ,DAT,&
                            AMX,BMX,CMX)
      WRITE(3,54)
      WRITE(5,544)
      REWIND 24
  54  FORMAT(10X,'STIMULI',13X,'SSED',11X,'ALPHA',12X,'BETA',12X,&
      'RSQ',11X,'MAED',10X,'STRESS',8X,'STRESS/ID')
  544 FORMAT('STIMULI',10X,'SSED',5X,'ALPHA',4X,'BETA',6X,'RSQ',&
      5X,'MAED',4X,'STRESS',2X,'STRESS/ID')
      DO 55 I=1,NQ
      READ(24)(TIT(J),J=1,13)
      WRITE(3,2511)(TIT(J),J=1,13),(DAT(I,J),J=1,7)
      WRITE(5,2522)(TIT(J),J=1,13),(DAT(I,J),J=1,7)
  55  CONTINUE
      WRITE(3,3)
      WRITE(3,200)
      WRITE(5,200)
  200 FORMAT(' FINAL ESTIMATE OF STIMULUS COORDINATES, SSE, R-SQUARE, AN&
      D # OF CASES')
      REWIND 24
      DO 201 I=1,NQ
      Q(I)=ZZ(I,1)
!
!  SIGN KLUDGE
!
!      ZZ(I,1)=-ZZ(I,1)
!      ZZ(I,2)=-ZZ(I,2)
!
      READ(24)(TIT(J),J=1,13)
      WRITE(3,251)(TIT(J),J=1,13),(ZZ(I,J),J=1,NR),SSQ(I),RSQ(I),XCOL(I)
      WRITE(5,252)(TIT(J),J=1,13),(ZZ(I,J),J=1,NR),SSQ(I),RSQ(I),XCOL(I)
      WRITE(22,252)(TIT(J),J=1,13),(ZZ(I,J),J=1,NR),SSQ(I),RSQ(I),&
                                                          XCOL(I)
  201 CONTINUE
      WRITE(3,3)
      WRITE(3,202)
      REWIND 1
      IF(NPRR.EQ.1)GO TO 203
      WRITE(5,202)
  203 CONTINUE
  202 FORMAT(' FINAL ESTIMATE OF INDIVIDUAL COORDINATES, SSE, R-SQUARE, &
      AND # OF CASES')
      DO 204 I=1,NP
      Q(I+NQ)=XX(I,1)
      IF(NOTIT.EQ.1)GO TO 253
      READ(1)(TIT(J),J=1,KTIT)
      IF(KTIT.GE.13.OR.I.GT.1)GO TO 2044
      KMJ=KTIT+1
      DO 2045 J=KMJ,13
 2045 TIT(J)=LTIT
 2044 CONTINUE
!
!  SIGN KLUDGE
!
!      XX(I,1)=-XX(I,1)
!      XX(I,2)=-XX(I,2)
!
      WRITE(3,251)(TIT(J),J=1,13),(XX(I,J),J=1,NR),SSQ(I+NQ),RSQ(I+NQ),&
                                                      XROW(I)
      WRITE(22,252)(TIT(J),J=1,13),(XX(I,J),J=1,NR),SSQ(I+NQ),RSQ(I+NQ),&
                                                      XROW(I)
      IF(NPRR.EQ.1)GO TO 204
      WRITE(5,252)(TIT(J),J=1,13),(XX(I,J),J=1,NR),SSQ(I+NQ),RSQ(I+NQ),&
                                                      XROW(I)
      GO TO 204
  253 WRITE(3,255)I,(XX(I,J),J=1,NR),SSQ(I+NQ),RSQ(I+NQ),XROW(I)
      WRITE(22,256)I,(XX(I,J),J=1,NR),SSQ(I+NQ),RSQ(I+NQ),XROW(I)
      IF(NPRR.EQ.1)GO TO 204
      WRITE(5,256)I,(XX(I,J),J=1,NR),SSQ(I+NQ),RSQ(I+NQ),XROW(I)
  204 CONTINUE
      WRITE(3,3)
!      IF(NR.EQ.1)CALL ONEPLT(NP,NQ,Q,XMAX)
!      IF(NR.GT.1)CALL PLOTIT(NP,NQ,NR,XX,ZZ,Q,RSQ,XMAX)
!      CALL PLOT3(NP,NQ,NR,XX,ZZ,DIST)
      NR=NR+1
      IF(NR.LE.NRMAX)WRITE(3,333)
      IF(NR.LE.NRMAX)GO TO 999
!
      call cpu_time(result1)
      write(*,1001)result1
      TOTALTIME=RESULT1-RESULT
      WRITE(*,1002)TOTALTIME
!
      STOP
      END
!
!
! *******************************************************************
!   SUBROUTINE STAT COMPUTES VARIOUS MEASURES OF FIT BETWEEN THE
!   INPUT DISTANCES (D*'s) AND THE REPRODUCED DISTANCES
! *******************************************************************
!
!
      SUBROUTINE STATKP(NP,NQ,NR,KP,XNPQT,XCOL,XROW,RSQ,SSQ,DIST,XX,ZZ,&
                      DAT,AMX,BMX,CMX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DIST(1900999),XX(85000,4),ZZ(143,4),DAT(100,7),&
                XCOL(143),RSQ(85143),XROW(85000),SSQ(85143)
      XNP=NP
      XNQ=NQ
      XNR=NR
      SUMA=0.0
      SUMB=0.0
      SUMC=0.0
      SUMCC=0.0
      SUMD=0.0
      SUME=0.0
      SUMF=0.0
      SUMG=0.0
      SUMH=0.0
      SUMI=0.0
      SUMJ=0.0
      SUMK=0.0
      SUML=0.0
      DO 13 I=1,NP
      IF(KP.NE.99)GO TO 21
      SUMA=0.0
      SUME=0.0
      SUMG=0.0
      SUMH=0.0
      SUMI=0.0
      SUMJ=0.0
  21  DO 14 J=1,NQ
      IF(DIST((I-1)*NQ+J).EQ.999.0)GO TO 14
      SUM=0.0
      DO 15 K=1,NR
  15  SUM=SUM+(XX(I,K)-ZZ(J,K))**2
      TDH=((SUM**(1.0/CMX))-BMX)/AMX
      TDS=((DIST((I-1)*NQ+J)**(1.0/CMX))-BMX)/AMX
      SUM=SQRT(SUM)
      SUMX=SQRT(DIST((I-1)*NQ+J))
      SUMA=SUMA+(SUMX-SUM)**2
      SUMC=SUMC+ABS(SUMX-SUM)
      SUME=SUME+(SUMX)*(SUM)
      SUMG=SUMG+SUMX
      SUMH=SUMH+SUM
      SUMI=SUMI+SUMX*SUMX
      SUMJ=SUMJ+SUM*SUM
      SUMB=SUMB+(TDH-TDS)**2
      SUMCC=SUMCC+TDH*TDH
  14  CONTINUE
      IF(KP.NE.99)GO TO 13
      AA=XROW(I)*SUME-SUMG*SUMH
      BB=XROW(I)*SUMI-SUMG*SUMG
      CC=XROW(I)*SUMJ-SUMH*SUMH
      IF(BB.EQ.0.0.OR.CC.EQ.0.0)RSQ(I+NQ)=99.0
      IF(BB.EQ.0.0.OR.CC.EQ.0.0)GO TO 13
      RSQ(I+NQ)=(AA/(SQRT(BB*CC)))**2
      SSQ(I+NQ)=SUMA
  13  CONTINUE
      IF(KP.EQ.99)GO TO 16
      DAT(KP,1)=SUMA
      AA=XNPQT*SUME-SUMG*SUMH
      BB=XNPQT*SUMI-SUMG*SUMG
      CC=XNPQT*SUMJ-SUMH*SUMH
      DAT(KP,3)=AA/CC
      DAT(KP,2)=SUMG/XNPQT-DAT(KP,3)*SUMH/XNPQT
      DAT(KP,4)=(AA*AA)/(BB*CC)
      DAT(KP,5)=(SUMC/SUMG)*100.0
      DAT(KP,6)=SQRT(SUMA/SUMJ)
      DAT(KP,7)=SQRT(SUMB/SUMCC)
      RETURN
  16  CONTINUE
      DO 18 J=1,NQ
      SUMA=0.0
      SUMB=0.0
      SUMC=0.0
      SUMD=0.0
      SUME=0.0
      SUMF=0.0
      SUMG=0.0
      SUMH=0.0
      SUMI=0.0
      DO 19 I=1,NP
      IF(DIST((I-1)*NQ+J).EQ.999.0)GO TO 19
      SUM=0.0
      DO 20 K=1,NR
  20  SUM=SUM+(XX(I,K)-ZZ(J,K))**2
      TDH=((SUM**(1.0/CMX))-BMX)/AMX
      TDS=((DIST((I-1)*NQ+J)**(1.0/CMX))-BMX)/AMX
      SUM=SQRT(SUM)
      SUMX=SQRT(DIST((I-1)*NQ+J))
      SUMA=SUMA+(SUMX-SUM)**2
      SUMB=SUMB+ABS(SUMX-SUM)
      SUMC=SUMC+SUMX*SUM
      SUMD=SUMD+SUMX
      SUME=SUME+SUM
      SUMF=SUMF+SUMX*SUMX
      SUMG=SUMG+SUM*SUM
      SUMH=SUMH+(TDH-TDS)**2
      SUMI=SUMI+TDH*TDH
  19  CONTINUE
      DAT(J,1)=SUMA
      AA=XCOL(J)*SUMC-SUMD*SUME
      BB=XCOL(J)*SUMF-SUMD*SUMD
      CC=XCOL(J)*SUMG-SUME*SUME
      DAT(J,3)=AA/CC
      DAT(J,2)=SUMD/XCOL(J)-DAT(J,3)*SUME/XCOL(J)
      DAT(J,4)=(AA*AA)/(BB*CC)
      RSQ(J)=DAT(J,4)
      SSQ(J)=SUMA
      DAT(J,5)=(SUMB/SUMD)*100.0
      DAT(J,6)=SQRT(SUMA/SUMG)
      DAT(J,7)=SQRT(SUMH/SUMI)
  18  CONTINUE
      RETURN
      END
!
!
!  *********************************************************************
!    SUBROUTINE FOCUS IS A QUASI-GRADIENT ALGORITHM.  IT COMPUTES THE
!    COORDINATES
!  *********************************************************************
!
!
      SUBROUTINE focus(NP,NQ,NR,NN,KJEP2,NSYM,TOL,XNPQT,XCOL,XROW,DIST,&
                                      XX,ZZ,AMX,BMX,CMX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DIST(1900999),XX(85000,4),ZZ(143,4),ZZZ(143,4),&
                  XXX(85000,4),DAT(100,7),XCOL(143),XROW(85000),&
                  RSQ(85143),SSQ(85143)
      IF(NR.EQ.1)GO TO 666
      XNP=NP
      XNQ=NQ
      XNR=NR
  66  FORMAT(6F10.4)
  6   FORMAT(6F15.8)
  96  FORMAT(16X,I6,2X,7F15.8)
  966 FORMAT(I6,2X,7F10.4)
  3   FORMAT(1X,130('*'))
      KPP=0
      KTP=0
  999 CONTINUE
      KPP=KPP+1
      KP=0
      DO 2 I=1,NP
      DO 2 J=1,NR
      XXX(I,J)=0.0
      IF(NQ.EQ.NN)XXX(I,J)=XX(I,J)
      IF(I.GT.NQ)GO TO 2
      ZZZ(I,J)=0.0
      IF(NP.EQ.NN)ZZZ(I,J)=ZZ(I,J)
  2   CONTINUE
  18  CONTINUE
      DO 2007 J=1,NQ
      WRITE(99,966)KP,(ZZ(J,JJ),JJ=1,2),(ZZZ(J,JJ),JJ=1,2)
 2007 CONTINUE
      KP=KP+1
      KTP=KTP+1
      IF(NP.EQ.NN)GO TO 312
      DO 4 I=1,NP
      DO 4 J=1,NQ
      IF(DIST((I-1)*NQ+J).EQ.999.0)GO TO 4
      SUM=0.0
      DO 7 JK=1,NR
  7   SUM=SUM+(XX(I,JK)-ZZ(J,JK))**2
      IF(SUM)50,51,50
  51  XC=1.0
      GO TO 52
  50  XC=(SQRT(DIST((I-1)*NQ+J)))/(SQRT(SUM))
  52  DO 8 K=1,NR
      ZZZ(J,K)=ZZZ(J,K)+(XX(I,K)-XC*(XX(I,K)-ZZ(J,K)&
      ))/(XCOL(J))
  8   CONTINUE
  4   CONTINUE
  213 IF(NQ.EQ.NN)GO TO 210
  312 DO 12 K=1,NR
      DO 99 I=1,NP
      SUMW=0.0
      DO 9 J=1,NQ
      IF(DIST((I-1)*NQ+J).EQ.999.0)GO TO 9
      SUM=0.0
      DO 11 JK=1,NR
      IF(NSYM.EQ.1)SUM=SUM+((XX(I,JK)-ZZ(J,JK)))**2
  11  IF(NSYM.NE.1)SUM=SUM+((XX(I,JK)-ZZZ(J,JK)))**2
      IF(SUM)53,54,53
  54  XC=1.0
      GO TO 55
  53  XC=(SQRT(DIST((I-1)*NQ+J)))/(SQRT(SUM))
  55  IF(NSYM.NE.1)XXX(I,K)=XXX(I,K)+(ZZZ(J,K)-XC*(ZZZ(J,K)-XX(I,K)))
      IF(NSYM.EQ.1)XXX(I,K)=XXX(I,K)+(ZZ(J,K)-XC*(ZZ(J,K)-XX(I,K)))
  9   CONTINUE
  99  XXX(I,K)=XXX(I,K)/XROW(I)
  12  CONTINUE
  210 CONTINUE
      CALL STATKP(NP,NQ,NR,KTP,XNPQT,XCOL,XROW,RSQ,SSQ,DIST,XXX,ZZZ,DAT,&
                                      AMX,BMX,CMX)
      IF(NQ.EQ.NN)GO TO 202
  201 DO 19 I=1,NP
      DO 19 J=1,NR
      XX(I,J)=XXX(I,J)
  19  XXX(I,J)=0.0
      IF(NP.EQ.NN)GO TO 313
  202 DO 21 I=1,NQ
      DO 21 J=1,NR
      ZZ(I,J)=ZZZ(I,J)
  21  ZZZ(I,J)=0.0
  313 IF(KTP.EQ.1)GO TO 211
      XSTOP=(DAT(KTP-1,1)-DAT(KTP,1))/DAT(KTP,1)
      IF(XSTOP.LE.TOL)GO TO 212
  211 IF(KP.LT.KJEP2)GO TO 18
  212 CONTINUE
  666 IF(NR.EQ.1)CALL WHOOPE(NP,NQ,NR,KJEP2,KTP,TOL,XNPQT,XCOL,XROW,&
                     RSQ,SSQ,DIST,XX,ZZ,DAT,AMX,BMX,CMX,NSYM)
  215 WRITE(3,17)
      WRITE(5,177)
  17  FORMAT(16X,'focus',9X,'SSED',11X,'ALPHA',12X,'BETA',12X,'RSQ',&
                  11X,'MAED',10X,'STRESS',8X,'STRESS/ID')
  177 FORMAT(' focus',6X,'SSED',6X,'ALPHA',6X,'BETA',7X,'RSQ',6X,&
                       'MAED',5X,'STRESS',3X,'STRESS/ID')
      DO 25 I=1,KTP
      WRITE(3,96) I,(DAT(I,J),J=1,7)
      WRITE(5,966)I,(DAT(I,J),J=1,7)
  25  CONTINUE
      WRITE(3,3)
      RETURN
      END
!
!
!  *************************************************************************
!   SUBROUTINE WHOOPE---IMPLEMENTS THE CONDITIONAL GLOBAL MINIMUM ALGORITHM
!   DEVELOPED BY POOLE, "LEAST SQUARES METRIC, UNIDIMENSIONAL UNFOLDING,"
!   PSYCHOMETRIKA, 1984.
!  *************************************************************************
!
!
      SUBROUTINE WHOOPE(NP,NQ,NR,KJEP2,KKK,TOL,XNPQT,XCOL,XROW,RSQ,SSQ,&
                           D,X,Z,DAT,AMX,BMX,CMX,NSYM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XCOL(143),XROW(85000),RSQ(85143),SSQ(85143),D(1900999),&
                 X(85000,4),Z(143,4),DAT(100,7),ZBEFORE(143,4)
      KKK=0
      IF(NSYM.EQ.1)GO TO 999
      DO 51 II=1,NQ
      ZBEFORE(II,1)=Z(II,1)
  51  CONTINUE
      DO 9999 IIII=1,KJEP2
      CALL focusX(NP,NQ,D,X,Z)
      DO 52 IJIJ=1,NQ
  52  CONTINUE
      KKK=KKK+1
      CALL STATKP(NP,NQ,NR,KKK,XNPQT,XCOL,XROW,RSQ,SSQ,D,X,Z,DAT,&
                              AMX,BMX,CMX)
      IF(KKK.EQ.1)GO TO 220
      AKKK=(DAT(KKK-1,1)-DAT(KKK,1))/DAT(KKK-1,1)
      IF(AKKK.LE.TOL)GO TO 9998
  220 CALL focusZ(NP,NQ,D,X,Z)
      KKK=KKK+1
      CALL STATKP(NP,NQ,NR,KKK,XNPQT,XCOL,XROW,RSQ,SSQ,D,X,Z,DAT,&
                             AMX,BMX,CMX)
      WRITE(*,1112)KKK,DAT(KKK,1)
 1112 FORMAT(' FOCUSZ',I4,F15.4)
      AKKK=(DAT(KKK-1,1)-DAT(KKK,1))/DAT(KKK-1,1)
      IF(AKKK.LE.TOL)GO TO 9998
 9999 CONTINUE
  999 CONTINUE
      DO 99 II=1,KJEP2
      DO 98 J=1,NQ
      CALL FOCUSW(NP,J,J,NQ,D,X,Z)
  98  X(J,1)=Z(J,1)
      KKK=KKK+1
      CALL STATKP(NP,NQ,NR,KKK,XNPQT,XCOL,XROW,RSQ,SSQ,D,X,Z,DAT,&
                             AMX,BMX,CMX)
      IF(KKK.EQ.1)GO TO 99
      AKKK=(DAT(KKK-1,1)-DAT(KKK,1))/DAT(KKK-1,1)
      IF(AKKK.LE.TOL)GO TO 9998
  99  CONTINUE
 9998 RETURN
      END
!
!
!  *********************************************************************
!   SUBROUTINE FOCUSX--INDIVIDUAL COORDINATES ARE HELD FIXED AND NEW
!   STIMULI COORDINATES ARE ESTIMATED USING THE CGM ALGORITHM
!  *********************************************************************
!
!
      SUBROUTINE FOCUSX(NP,NQ,D,X,Z)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION D(1900999),X(85000,4),Z(143,4),LL(85000),Q(85143),&
                   XX(85000,2)
      DO 303 I=1,NP
      LL(I)=I
  303 Q(I)=X(I,1)
      CALL RSORT(Q,NP,LL)
      DO 11 II=1,NP
  11  CONTINUE
      DO 998 II=1,NQ
      ASUM=0.0
      BSUM=0.0
      WWSUM=0.0
      DO 66 I=1,NP
      IF(D((LL(I)-1)*NQ+II).EQ.999.0)GO TO 66
      XX(I,1)=Q(I)-SQRT(D((LL(I)-1)*NQ+II))
      XX(I,2)=Q(I)+SQRT(D((LL(I)-1)*NQ+II))
      WWSUM=WWSUM+1.0
      ASUM=ASUM+XX(I,1)
      BSUM=BSUM+XX(I,1)**2
  66   CONTINUE
      AA=WWSUM*BSUM-ASUM*ASUM
      KK=1
      DO 77 I=1,NP
      IF(D((LL(I)-1)*NQ+II).EQ.999.0)GO TO 77
      ASUM=ASUM-XX(I,1)+XX(I,2)
      BSUM=BSUM-XX(I,1)**2+XX(I,2)**2
      BB=WWSUM*BSUM-ASUM*ASUM
      CC=DMIN1(AA,BB)
      IF(ABS(CC-AA).LE..00001.AND.KK.GT.1)GO TO 88
      IF(ABS(CC-AA).LE..00001.AND.KK.EQ.1)Z(II,1)= &
                           (ASUM+XX(I,1)-XX(I,2))/WWSUM
      IF(ABS(CC-BB).LE..00001)Z(II,1)=ASUM/WWSUM
!      IF(CC.EQ.AA.AND.KK.GT.1)GO TO 88
!      IF(CC.EQ.AA.AND.KK.EQ.1)Z(II,1)=(ASUM+XX(I,1)-XX(I,2))/WWSUM
!      IF(CC.EQ.BB)Z(II,1)=ASUM/WWSUM
  88   AA=CC
      KK=KK+1
  77   CONTINUE
  998 CONTINUE
      RETURN
      END
!
!
!  **********************************************************************
!    SUBROUTINE FOCUSZ---STIMULI COORDINATES ARE HELD FIXED AND THE
!    INDIVIDUAL COORDINATES ARE ESTIMATED USING THE CGM ALGORITHM
!  **********************************************************************
!
!
      SUBROUTINE FOCUSZ(NP,NQ,D,X,Z)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION D(1900999),X(85000,4),Z(143,4),Q(143),L(143),&
                ZZ(143,4)
      DO 303 I=1,NQ
      L(I)=I
  303 Q(I)=Z(I,1)
      CALL RSORT(Q,NQ,L)
      DO 999 II=1,NP
      ASUM=0.0
      BSUM=0.0
      WWSUM=0.0
      DO 6 I=1,NQ
      IF(D((II-1)*NQ+L(I)).EQ.999.0)GO TO 6
      ZZ(I,1)=Q(I)-SQRT(D((II-1)*NQ+L(I)))
      ZZ(I,2)=Q(I)+SQRT(D((II-1)*NQ+L(I)))
      WWSUM=WWSUM+1.0
      ASUM=ASUM+ZZ(I,1)
      BSUM=BSUM+ZZ(I,1)**2
  6   CONTINUE
      AA=WWSUM*BSUM-ASUM*ASUM
      KK=1
      DO 7 I=1,NQ
      IF(D((II-1)*NQ+L(I)).EQ.999.0)GO TO 7
      ASUM=ASUM-ZZ(I,1)+ZZ(I,2)
      BSUM=BSUM-ZZ(I,1)**2+ZZ(I,2)**2
      BB=WWSUM*BSUM-ASUM*ASUM
      CC=DMIN1(AA,BB)
      IF(ABS(CC-AA).LE..00001.AND.KK.GT.1)GO TO 8
      IF(ABS(CC-AA).LE..00001.AND.KK.EQ.1)X(II,1)= &
                           (ASUM+ZZ(I,1)-ZZ(I,2))/WWSUM
      IF(ABS(CC-BB).LE..00001)X(II,1)=ASUM/WWSUM
!      IF(CC.EQ.AA.AND.KK.GT.1)GO TO 8
!      IF(CC.EQ.AA.AND.KK.EQ.1)X(II,1)=(ASUM+ZZ(I,1)-ZZ(I,2))/WWSUM
!      IF(CC.EQ.BB)X(II,1)=ASUM/WWSUM
  8   AA=CC
      KK=KK+1
  7   CONTINUE
  999 CONTINUE
      RETURN
      END
!
!
!  *********************************************************************
!    SUBROUTINE FOCUSW---PERFORMS LEAST SQUARES METRIC SIMILARITIES
!    ANALYSIS USING THE CONDITIONAL GLOBAL MINIMUM ALGORITHM
!  *********************************************************************
!
!
      SUBROUTINE FOCUSW(NP,NQ1,NQ2,NQ,D,X,Z)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION D(1900999),X(85000,4),Z(143,4),LL(85000),Q(85143),&
                 XX(85000,2)
      DO 303 I=1,NP
      LL(I)=I
  303 Q(I)=X(I,1)
      CALL RSORT(Q,NP,LL)
      DO 998 II=NQ1,NQ2
      ASUM=0.0
      BSUM=0.0
      WWSUM=0.0
      DO 66 I=1,NP
      IF(D((LL(I)-1)*NQ+II).EQ.999.0)GO TO 66
      XX(I,1)=Q(I)-SQRT(D((LL(I)-1)*NQ+II))
      XX(I,2)=Q(I)+SQRT(D((LL(I)-1)*NQ+II))
      WWSUM=WWSUM+1.0
      ASUM=ASUM+XX(I,1)
      BSUM=BSUM+XX(I,1)**2
  66   CONTINUE
      AA=WWSUM*BSUM-ASUM*ASUM
      KK=1
      DO 77 I=1,NP
      IF(D((LL(I)-1)*NQ+II).EQ.999.0)GO TO 77
      ASUM=ASUM-XX(I,1)+XX(I,2)
      BSUM=BSUM-XX(I,1)**2+XX(I,2)**2
      BB=WWSUM*BSUM-ASUM*ASUM
      CC=DMIN1(AA,BB)
      IF(CC.EQ.AA.AND.KK.GT.1)GO TO 88
      IF(CC.EQ.AA.AND.KK.EQ.1)Z(II,1)=(ASUM+XX(I,1)-XX(I,2))/WWSUM
      IF(CC.EQ.BB)Z(II,1)=ASUM/WWSUM
  88   AA=CC
      KK=KK+1
  77   CONTINUE
  998 CONTINUE
      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
