! ! *********************************************************************** ! ! 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