c paramap.f C The authors of this software are J D Carroll and Jih-Jie Chang. C Copyright (c) 1993 by AT&T. C Permission to use, copy, modify, and distribute this software for any C purpose without fee is hereby granted, provided that this entire C notice is included in all copies of any software which is or includes C a copy or modification of this software and in all copies of the C supporting documentation for such software. C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP- C LIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. C This software comes from the FIRST MDS Package of AT&T Bell Laboratories C For explanation of the method this software implements, see C Shepard,R.N. & Carroll,J.D. (1966), "Parametric representation of C nonlinear data structures", In P.R.Krishnaiah (Ed.), Multivariate C Analysis (pp. 561-592). New York: Academic Press C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ REAL KEY,KEYV 00002390 DATA PROB,KEYV,FINISH,YES/'PROB','KEYV','FINI','YES'/ 00002400 DIMENSION FX(20),R(10,10),V(10,10),DIMX(150),DIMY(150),U(11175), 00002410 1 X(150,10),DKX(150,10),TIT(20),COSA(100),RGX(100),RGSA(100), 00002420 2 ALF(100),RLGR(100),XAL(100),AS(100),XL(100),C(100),Y(150), 00002430 3 YY(500),D(11175) 00002440 COMMON U,X,DKX,TIT,KPUNCH,IT,FINI,NS,NF,UEXP,DB,DC,CALF,SA,TK, 00002450 1 DB1,DC1,Z,BC,GV,COSA,DSA,RGX,RGSA,ALF,C,RLGR,XAL,AS,XL,PCUT, 00002460 2 DIMX,DIMY 00002470 EQUIVALENCE (DKX(1),Y(1)), (DKX(151),YY(1)) 00002480 1000 UEXP = 1.0 00002490 DB = 2.0 00002500 DC = -1.0 00002510 GCRIT = 0.0 00002520 CALF = 0.05 00002530 CSA = 5.0 00002540 PCUT = 0.5 00002550 KROW = 4 00002560 KCOL = 4 00002570 NNORM = 1 00002580 READ(5,900)CARD,NS,MR,NFS,NFE,MAXIT,KREADC,KX,KPUNCH,NPO,KEY,IX, 00002590 1 PLTFUN , DSTPRT 00002600 900 FORMAT(A4,2I3,2I2,I3,I2,2I1,I2,A3,I9,2A3) 00002610 IF(CARD .EQ. FINISH) GO TO 926 00002620 IF(CARD .NE. PROB) GO TO 924 00002630 IF(MAXIT .GT. 100) MAXIT = 100 00002640 IF(NPO .EQ. 0) NPO = MAXIT 00002650 WRITE(6,901) NS,MR,NFS,NFE,MAXIT,NPO 00002660 901 FORMAT(1H1,45X,'CARROLL - CHANG PARAMETRIC MAPPING PROGRAM'/55X, 00002670 1'VERSION 2 OF 12/27/67'///10X,'PROBLEM CARD PROCESSED '/35X, 00002680 2'NUMBER OF POINTS = ',I3/35X,'NUMBER OF INPUT DIMENSIONS = ',I3, 00002690 3' (IF = 0, POINTS NOT USED AS INPUT)'/35X,'STARTING NUMBER OF DIME00002700 4NSIONS FOR MAPPING = ',I2/35X,'ENDING NUMBER OF DIMENSIONS FOR MAP00002710 5PING = ',I2/35X,'MAXIMUM NUMBER OF ITERATIONS = ',I3/35X,'X-MATRIX00002720 6 PRINTED AND PLOTTED EVERY',I5,' ITERATIONS') 00002730 WRITE(6,902) 00002740 902 FORMAT(//10X,'INPUT DATA FORM = ') 00002750 IF(KREADC) 903,905,907 00002760 903 WRITE(6,904) 00002770 904 FORMAT(1H+,30X,'COVARIANCE MATRIX') 00002780 GO TO 912 00002790 905 WRITE(6,906) MR 00002800 906 FORMAT(1H+,30X,'DATA POINT COORDINATES IN',I5,' DIMENSIONS, (Y MAT00002810 1RIX)') 00002820 GO TO 912 00002830 907 CONTINUE 00002840 GO TO (908,910,931), KREADC 00002850 908 WRITE(6,909) 00002860 909 FORMAT(1H+,30X,'DISTANCE-SQUARED (U) MATRIX') 00002870 GO TO 912 00002880 931 WRITE(6,932) 00002890 932 FORMAT(1H+,30X,'DISTANCE (DISSIMILARITY) MATRIX - SQRT(U)') 00002900 GO TO 912 00002910 910 WRITE(6,911) 00002920 911 FORMAT(1H+,30X,'CORRELATION MATRIX') 00002930 912 CONTINUE 00002940 IF(KEY .NE. YES) GO TO 915 00002950 READ(5,913)CARD,UEXP,DB,DC,GCRIT,CALF,CSA,PCUT,KROW,KCOL,NNORM 00002960 913 FORMAT(A4,7F8.4,3I3) 00002970 IF(CARD .NE. KEYV) GO TO 924 00002980 WRITE(6,914) 00002990 914 FORMAT(///10X,'KEY VALUES CARD PROCESSED '/) 00003000 GO TO 917 00003010 915 WRITE(6,916) 00003020 916 FORMAT(///10X,'KEY VALUES CARD NOT USED - NORMAL VALUES ARE '/) 00003030 917 WRITE(6,918)UEXP,DB,DC,GCRIT,CALF,CSA,PCUT,KROW,KCOL 00003040 918 FORMAT(15X,'A = ',F8.4/15X,'B = ',F8.4/15X,'C = ',F8.4/15X, 00003050 1'KAPPA CRIT = ',F8.4/15X,'INITIAL STEP SIZE = ',F8.4/15X,'COMPUTED00003060 2 DISTANCE ADDITIVE PARAMETER = ',F8.4/15X,'STEP SIZE CUTBACK COEFF00003070 3ICIENT = ',F8.4/15X,'DATA MATRIX NORMALIZATION OPTIONS -- 00003080 4KROW = ',I4,5X,'KCOL = ',I4//) 00003090 DB = -DB 00003100 NSM1 = NS - 1 00003110 NNS = (NS*(NS-1))/2 00003120 IF(KX .EQ. 0) GO TO 919 00003130 IF(KX .EQ. 1) GO TO 921 00003140 GO TO 924 00003150 919 WRITE(6,920) IX 00003160 920 FORMAT(10X,'INITIAL CONFIGURATION GENERATED BY PROGRAM WITH RANDOM00003170 1LY ASSIGNED POINT NUMBERS'/18X,'RANDOM NUMBER GENERATOR STARTING S00003180 2EED = ',I12//) 00003190 GO TO 923 00003200 921 WRITE(6,922) 00003210 922 FORMAT(10X,'INITIAL CONFIGURATION SUPPLIED AS INPUT DATA -- (X MAT00003220 1RIX)'///) 00003230 923 CONTINUE 00003240 GO TO 928 00003250 924 WRITE(6,925) 00003260 925 FORMAT(//20X,'****CONTROL CARD TERMINATION - ERROR****'/1H1) 00003270 CALL EXIT 00003280 926 WRITE(6,927) 00003290 927 FORMAT(//20X,'****CONTROL CARD TERMINATION - FINISH****'/1H1) 00003300 CALL EXIT 00003310 928 CONTINUE 00003320 READ(5,929) TIT 00003330 929 FORMAT(20A4) 00003340 WRITE(6,930) TIT 00003350 930 FORMAT(10X,'TITLE ',5X,20A4//) 00003360 CALL IDENT 00003370 CALL READY(NS,MR,KREADC,KROW,KCOL) 00003380 NO=NPO 00003390 FINI=1.0 00003400 C PRINT U MATRIX (DISTANCE SQUARED) IF HAS NOT BEEN PRINTED IN READY. 00003410 IF(KREADC .EQ. 1) GO TO 115 00003420 WRITE(6,302) 00003430 WRITE(6,304) (K, K=1,NSM1) 00003440 L1 = 1 00003450 L2 = 1 00003460 DO 66 I = 2,NS 00003470 WRITE(6,303) I,(U(J), J=L1,L2) 00003480 L1 = L2 + 1 00003490 66 L2 = I + L2 00003500 115 CONTINUE 00003510 C CHANGE NEXT STATEMENT IF DESIRE U-MATRIX PUNCHED CONDITIONAL ON 00003520 C THE VARIABLE "KPUNCH". 00003530 GO TO 114 00003540 315 WRITE(7,305) 00003550 L1 = 1 00003560 L2 = 1 00003570 DO 69 I = 2,NS 00003580 WRITE(7,306) I,(U(J), J=L1,L2) 00003590 L1 = L2 + 1 00003600 69 L2 = I + L2 00003610 114 DO 520 NF = NFS,NFE 00003620 WRITE(6,307) NF 00003630 307 FORMAT(1H1,53X,'MAPPING IN',I5,' DIMENSIONS'//) 00003640 CALL CLEAR(X,1500,R,100,V,100) 00003650 SA=CSA 00003660 IF(KX)601,600,601 00003670 601 READ(5,929) FX 00003680 DO 470 I = 1,NS 00003690 470 READ(5,FX) (X(I,J), J=1,NF) 00003700 GOTO113 00003710 600 CALL INITX(IX) 00003720 C BEGIN ITERATION 00003730 113 IT=1 00003740 CALL ITRN(NPO,MAXIT,GCRIT,CMAX) 00003750 IF(NNORM) 121,525,121 00003760 121 CALL NORM 00003770 525 WRITE(6,200) NF 00003780 200 FORMAT(1H1,40X,'HISTORY OF COMPUTATION FOR',I6,' DIMENSION SOLUTI00003790 1ON'///1X,'ITERATION',7X,'KAPPA',6X,'REL GRAD',7X,'RLGX',8X, 00003800 2 'RLGSA',8X,'\X,SA\',9X,'\X\',9X,'COSA',9X,'ALPHA',7X,'SMALL A'//)00003810 DO 201 I = 1,IT 00003820 WRITE(6,202) I,C(I),RLGR(I),RGX(I),RGSA(I),XAL(I),XL(I),COSA(I), 00003830 1 ALF(I),AS(I) 00003840 202 FORMAT(1X,I6,3X,9F13.6) 00003850 RI = I 00003860 201 DIMX(I) = RI 00003870 IF(IT .LT. MAXIT) WRITE(6,206) IT 00003880 206 FORMAT(5X,'*****COMPUTATION TERMINATED AT TIERATION',I5,' DUE TO S00003890 1TEP SIZE UNDERFLOW IN SUBROUTINE "ALPHA"*****') 00003900 WRITE(6,203) CMAX,IT,C(IT) 00003910 203 FORMAT(//5X,'KAPPA MAX = ',F14.6/5X,'FINAL KAPPA AFTER',I4, 00003920 1 ' ITERATIONS = ',F14.6/) 00003930 WRITE(6,204) 00003940 204 FORMAT(1H1,42X,'PLOT OF KAPPA (Y-AXIS) VS. ITERATIONS (X-AXIS)') 00003950 AAA = 0.0 00003960 BBB = 0.0 00003970 CCC = 0.0 00003980 DDD = 0.0 00003990 CALL PLOTX(DIMX,C,IT,-1,+1,0,AAA,BBB,CCC,DDD) 00004000 IF(KREADC .NE. 0) GO TO 210 00004010 IF(PLTFUN .NE. YES) GO TO 210 00004020 WRITE(6,218) MR,NF,C(IT),IT,CMAX 00004030 218 FORMAT(//////////10X,'THE FOLLOWING PLOTS SHOW ALL THE FUNCTIONS R00004040 1EQUIRED TO MAP THE',I4,' DIMENSIONS OF Y INTO THE',I4,' DIMENSIONS00004050 2 OF X'//46X,'KAPPA = ',F10.6,' AFTER',I4,' ITERATIONS'/46X,'KAPPA 00004060 3MAX = ',F10.6) 00004070 DO 216 J = 1,NF 00004080 REWIND 3 00004090 DO 215 I = 1,NS 00004100 215 DIMX(I) = X(I,J) 00004110 DO 217 M = 1,MR 00004120 READ (3) (Y(L), L=1,NS) 00004130 WRITE(6,219) M,J 00004140 219 FORMAT(1H1,34X,'DIMENSION',I4,' COORDINATES OF Y VS. DIMENSION', 00004150 1 I4,' COORDINATES OF X') 00004160 AAA = 0.0 00004170 BBB = 0.0 00004180 CCC = 0.0 00004190 DDD = 0.0 00004200 CALL PLOTX(DIMX,Y,NS,+1,+1,+1,AAA,BBB,CCC,DDD) 00004210 217 CONTINUE 00004220 216 CONTINUE 00004230 210 CALL CLEAR(D,11175) 00004240 REWIND 2 00004250 WRITE (2) (U(M), M=1,NNS) 00004260 L = 0 00004270 DO 230 I = 2,NS 00004280 IM1 = I - 1 00004290 DO 230 J = 1,IM1 00004300 L = L + 1 00004310 DO 231 K = 1,NF 00004320 D(L) = D(L) + (X(I,K) - X(J,K))**2 00004330 231 CONTINUE 00004340 D(L) = SQRT(D(L)) 00004350 U(L) = SQRT(U(L)) 00004360 230 CONTINUE 00004370 IF(DSTPRT .NE. YES) GO TO 320 00004380 WRITE(6,308) 00004390 WRITE(6,304) (K, K=1,NSM1) 00004400 L1 = 1 00004410 L2 = 1 00004420 DO 233 I = 2,NS 00004430 WRITE(6,303) I,(U(J), J=L1,L2) 00004440 L1 = L2 + 1 00004450 233 L2 = I + L2 00004460 WRITE(6,309) 00004470 WRITE(6,304) (K, K=1,NSM1) 00004480 L1 = 1 00004490 L2 = 1 00004500 DO 234 I = 2,NS 00004510 WRITE(6,303) I,(D(J), J=L1,L2) 00004520 L1 = L2 + 1 00004530 234 L2 = I + L2 00004540 320 IF(KPUNCH - 2) 310,311,311 00004550 311 WRITE(7,312) NF 00004560 L1 = 1 00004570 L2 = 1 00004580 DO 313 J = 2,NS 00004590 WRITE(7,314) J,(D(I), I=L1,L2) 00004600 L1 = L2 + 1 00004610 313 L2 = J + L2 00004620 310 WRITE(6,232) NF 00004630 232 FORMAT(1H1,40X,'PLOT OF INPUT DISTANCES-SQRT(U) VS. OUTPUT DISTANC00004640 1ES-(D) FOR',I3,' DIMENSIONAL MAPPING') 00004650 AAA = 0.0 00004660 BBB = 0.0 00004670 CCC = 0.0 00004680 DDD = 0.0 00004690 CALL PLOTX(D,U,NNS,-1,-1,0,AAA,BBB,CCC,DDD) 00004700 REWIND 2 00004710 READ (2) (U(M), M=1,NNS) 00004720 FINI = 1.0 00004730 NPO=NO 00004740 520 CONTINUE 00004750 302 FORMAT(1H1,10X,'U - MATRIX (EQUIVALENT DISTANCE SQUARED)'//) 00004760 303 FORMAT(/I5,10F11.5/(5X,10F11.5)) 00004770 304 FORMAT(10I11) 00004780 305 FORMAT(2X,'U-MATRIX'/2X,'(I3,8F9.5/(3X,8F9.5))') 00004790 306 FORMAT(I3,8F9.5/(3X,8F9.5)) 00004800 308 FORMAT(1H1,10X,'INPUT INTERPOINT DISTANCES (SQRT OF U MATRIX)'//) 00004810 309 FORMAT(1H1,10X,'D - MATRIX, INTERPOINT DISTANCES OF FINAL X CONFIG00004820 1URATION'//) 00004830 312 FORMAT(2X,'D-MATRIX',I4,' DIMENSIONS'/2X,'(I3,9F8.4/(3X,9F8.4))') 00004840 314 FORMAT(I3,9F8.4/(3X,9F8.4)) 00004850 GO TO 1000 00004860 END 00004870 SUBROUTINE READY(NS,MR,KREADC,KROW,KCOL) 00004880 DIMENSION Y(150),U(11175),WR(500),WC(150),CII(150),T(150) 00004890 DIMENSION FMTY(20),FMTC(20),YY(500),DUMMY(1788),DIMX(150) 00004900 DIMENSION DIMY(150),PAD(510) 00004910 COMMON U,WR,WC,CII,T,FMTY,FMTC,PAD,Y,YY,DUMMY,DIMX,DIMY 00004920 NSM1 = NS - 1 00004930 NNS = (NS*(NS-1))/2 00004940 IF(KREADC)83,81,82 00004950 81 XNS = NS 00004960 XM=MR 00004970 C TEST IF MULTIPLYING W 00004980 IF(KROW-10)2,1,1 00004990 1 READ(5,901 ) (WR(I),I=1,MR) 00005000 901 FORMAT(10F10.0) 00005010 2 IF(KCOL-10)7,6,6 00005020 6 READ(5,901 ) (WC(I),I=1,NS) 00005030 7 READ(5,500) FMTY 00005040 C READ IN DATA BY ROW AND COMPUTE DISTANCES 00005050 CALL CLEAR(CII,150,T,150,U,11175) 00005060 WRITE(6,250) NS,MR,(FMTY(I), I=1,10) 00005070 250 FORMAT(1H1,33X,'INPUT DATA (Y-MATRIX) -- COORDINATES OF',I3,' POIN00005080 1TS IN',I3,' DIMENSIONS'/38X,'INPUT FORMAT = ',10A4//) 00005090 WRITE(6,251) (K, K=1,MR) 00005100 251 FORMAT(12I10//) 00005110 REWIND 4 00005120 DO 200 I = 1,NS 00005130 READ(5,FMTY) (YY(J), J=1,MR) 00005140 WRITE(6,252) I,(YY(J), J=1,MR) 00005150 252 FORMAT(I4,12F10.4/(4X,12F10.4)) 00005160 WRITE(4) (YY(J), J=1,MR) 00005170 200 CONTINUE 00005180 REWIND 4 00005190 REWIND 3 00005200 DO 80 I = 1,MR 00005210 DO 201 K = 1,NS 00005220 READ (4) (YY(J), J=1,MR) 00005230 Y(K) = YY(I) 00005240 201 CONTINUE 00005250 REWIND 4 00005260 WRITE (3) (Y(L), L=1,NS) 00005270 IF(KROW-10)11,12,12 00005280 12 KK=KROW-10 00005290 GOTO13 00005300 11 KK=KROW 00005310 13 IF(KK-4)18,24,19 00005320 19 IF(KK-8)18,24,24 00005330 C COMPUTE MEAN 00005340 18 CONTINUE 00005350 SUMY=0 00005360 SSQY=0 00005370 DO20J=1,NS 00005380 SUMY=SUMY+Y(J) 00005390 20 SSQY=SSQY+Y(J)**2 00005400 YMEAN=SUMY/XNS 00005410 SY=SSQY/XNS 00005420 GOTO(21,22,23,24,25,25,26,24),KK 00005430 21 P=0 00005440 GOTO28 00005450 25 P=YMEAN 00005460 28 DO31J=1,NS 00005470 31 Y(J)=(Y(J)-P)/SQRT(SY) 00005480 GOTO24 00005490 22 DO32J=1,NS 00005500 32 Y(J)=Y(J)/SQRT(SY-YMEAN**2) 00005510 GOTO24 00005520 23 P=0 00005530 GOTO29 00005540 26 P=YMEAN 00005550 29 DO36J=1,NS 00005560 36 Y(J)=(Y(J)-P)/YMEAN 00005570 24 IF(KROW-10)34,33,33 00005580 33 DO35J=1,NS 00005590 35 Y(J)=Y(J)*WR(I) 00005600 C COMPUTE CIJ(U) 00005610 34 L = 0 00005620 CII(1)=CII(1)+Y(1)**2 00005630 T(1)=T(1)+Y(1) 00005640 DO40J=2,NS 00005650 II=J-1 00005660 DO41K=1,II 00005670 L=L+1 00005680 41 U(L)=U(L)+Y(J)*Y(K) 00005690 CII(J)=CII(J)+Y(J)**2 00005700 40 T(J)=T(J)+Y(J) 00005710 80 CONTINUE 00005720 REWIND 3 00005730 REWIND 4 00005740 IF(MR .NE. 2) GO TO 270 00005750 READ (3) (Y(J), J=1,NS) 00005760 DO 271 K = 1,NS 00005770 271 DIMX(K) = Y(K) 00005780 READ (3) (Y(J), J=1,NS) 00005790 DO 272 K = 1,NS 00005800 272 DIMY(K) = Y(K) 00005810 WRITE(6,273) 00005820 273 FORMAT(1H1,10X,'PLOT OF 2-DIMENSIONAL Y INPUT MATRIX') 00005830 AAA = 0.0 00005840 BBB = 0.0 00005850 CCC = 0.0 00005860 DDD = 0.0 00005870 CALL PLOTX(DIMX,DIMY,NS,+1,+1,+1,AAA,BBB,CCC,DDD) 00005880 REWIND 3 00005890 270 CONTINUE 00005900 C NORMALIZE COLUMNWISE 00005910 IF(KCOL-10)51,52,52 00005920 52 KK=KCOL-10 00005930 GOTO53 00005940 51 KK=KCOL 00005950 53 IF(KK-4)59,59,56 00005960 56 L=0 00005970 DO60J=2,NS 00005980 JJ=J-1 00005990 DO60K=1,JJ 00006000 L=L+1 00006010 60 U(L)=U(L)-T(J)*T(K)/XM 00006020 C COMPUTE NEW CIJ 00006030 59 L=0 00006040 DO70J=2,NS 00006050 JJ=J-1 00006060 DO70K=1,JJ 00006070 L=L+1 00006080 GOTO(72,73,74,75,72,72,74,75),KK 00006090 72 U(L)=U(L)/SQRT(CII(J)*CII(K)) 00006100 GOTO75 00006110 73 Q=(CII(J)-(T(J)**2)/XNS)*(CII(K)-(T(K)**2)/XNS) 00006120 U(L)=U(L)/SQRT(Q) 00006130 GOTO75 00006140 74 U(L)=U(L)/(T(J)*T(K)) 00006150 75 IF(KCOL-10)70,76,76 00006160 76 U(L)=WC(J)*WC(K)*U(L) 00006170 70 CONTINUE 00006180 C COMPUTE DISTANCES 00006190 L=0 00006200 DO90I=2,NS 00006210 II=I-1 00006220 DO90J=1,II 00006230 L=L+1 00006240 90 U(L)=CII(I)+CII(J)-2.0*U(L) 00006250 GOTO100 00006260 82 READ(5,500) FMTC 00006270 500 FORMAT(20A4) 00006280 GO TO (257,258,259), KREADC 00006290 257 WRITE(6,253) (FMTC(I), I=1,10) 00006300 253 FORMAT(1H1,46X,'INPUT DATA (U-MATRIX) - DISTANCES SQUARED'/ 00006310 1 54X,'INPUT FORMAT = ',10A4//) 00006320 WRITE(6,261) (K, K=1,NSM1) 00006330 GO TO 260 00006340 258 WRITE(6,254) (FMTC(I), I=1,10) 00006350 254 FORMAT(1H1,50X,'INPUT DATA - CORRELATION MATRIX'/54X,'INPUT FORMAT00006360 1 = ',10A4//) 00006370 WRITE(6,261) (K, K=1,NSM1) 00006380 GO TO 260 00006390 259 WRITE(6,255) (FMTC(I), I=1,10) 00006400 255 FORMAT(1H1,38X,'INPUT DATA - DISTANCE OR DISSIMILARITY MATRIX, SQR00006410 1T(U)'/45X,'INPUT FORMAT = ',10A4//) 00006420 WRITE(6,261) (K, K=1,NSM1) 00006430 260 L1 = 1 00006440 L2 = 1 00006450 DO 84 I = 2,NS 00006460 READ(5,FMTC) (U(J), J=L1,L2) 00006470 WRITE(6,262) I,(U(L), L=L1,L2) 00006480 L1 = L2 + 1 00006490 84 L2 = I + L2 00006500 261 FORMAT(10I11) 00006510 262 FORMAT(/I4,10F11.5/(4X,10F11.5)) 00006520 GO TO (100,95,106), KREADC 00006530 106 DO 263 I = 1,NNS 00006540 263 U(I) = U(I)**2 00006550 GO TO 100 00006560 83 READ(5,500) FMTC 00006570 WRITE(6,256) (FMTC(I), I=1,10) 00006580 256 FORMAT(1H1,50X,'INPUT DATA - COVARIANCE MATRIX'/54X,'INPUT FORMAT 00006590 1= ',10A4//) 00006600 WRITE(6,264) (K, K=1,NS) 00006610 264 FORMAT(11I11) 00006620 READ(5,FMTC) CII(1) 00006630 WRITE(6,265) CII(1) 00006640 265 FORMAT(/3X,'1',F11.5) 00006650 L1 = 1 00006660 L2 = 1 00006670 DO 85 I = 2,NS 00006680 READ(5,FMTC) (U(J), J=L1,L2), CII(I) 00006690 WRITE(6,266) I,(U(L), L=L1,L2), CII(I) 00006700 266 FORMAT(/I4,11F11.5/(4X,11F11.5)) 00006710 L1 = L2 + 1 00006720 85 L2 = I + L2 00006730 95 L = 0 00006740 DO 98 I = 2,NS 00006750 II = I - 1 00006760 DO 98 J = 1,II 00006770 L = L + 1 00006780 IF(KREADC) 97,100,99 00006790 97 U(L) = CII(I) + CII(J) - 2.0*U(L) 00006800 GO TO 98 00006810 99 U(L) = 2.0 * (1.0 - U(L)) 00006820 98 CONTINUE 00006830 100 RETURN 00006840 END 00006850 SUBROUTINE INITX(IX) 00006860 C GENERATING INITIAL CONFIGURATION 00006870 DIMENSION U(11175),X(150,10),V(150,10),X1(23),X2(5),NOS(150) 00006880 COMMON U,X,V,X1,NS,NF,X2,NOS 00006890 XN=(NS-1)/NF+1 00006900 IC=0 00006910 DO20I=1,NS 00006920 IC=IC+1 00006930 DO10K=1,NF 00006940 IF(IC-K)5,6,5 00006950 6 V(I,K)=XN 00006960 GOTO10 00006970 5 V(I,K)=0 00006980 10 CONTINUE 00006990 IF(IC-NF)20,8,8 00007000 8 IC=0 00007010 XN=XN-1. 00007020 20 CONTINUE 00007030 DO 40 I = 1,NS 00007040 CALL RANDU(IX,IY,RN) 00007050 IX = IY 00007060 NRN = -100.0 * RN 00007070 40 NOS(I) = NRN 00007080 DO 43 K = 1,NS 00007090 MIN = 999999999 00007100 DO 41 J = 1,NS 00007110 IF(NOS(J)) 44,44,41 00007120 44 IF(NOS(J) .LT. MIN) GO TO 42 00007130 GO TO 41 00007140 42 MIN = NOS(J) 00007150 LOCMIN = J 00007160 41 CONTINUE 00007170 NOS(LOCMIN) = K 00007180 43 CONTINUE 00007190 DO30I=1,NS 00007200 DO30J=1,NF 00007210 N=NOS(I) 00007220 30 X(I,J)=V(N,J) 00007230 RETURN 00007240 END 00007250 SUBROUTINE ITRN(NPO,MAXIT,GCRIT,CMAX) 00007260 DIMENSION U(11175),X(150,10),DKX(150,10),TIT(20),COSA(100), 00007270 1 RGX(100),RGSA(100),ALF(100),RLGR(100),XAL(100),AS(100),XL(100), 00007280 2 C(100) 00007290 COMMON U,X,DKX,TIT,KPUNCH,IT,FINI,NS,NF,UEXP,DB,DC,CALF,SA,TK, 00007300 1 DB1,DC1,Z,BC,GV,COSA,DSA,RGX,RGSA,ALF,C,RLGR,XAL,AS,XL,PCUT 00007310 DB1=DB-1. 00007320 DC1=DC-1. 00007330 BC=-DB/DC 00007340 DV=UEXP+DB 00007350 A1=0 00007360 F1=0 00007370 SUMU=0 00007380 L=0 00007390 NMZ = 0 00007400 DO131I=2,NS 00007410 II=I-1 00007420 DO131J=1,II 00007430 L=L+1 00007440 IF(U(L)) 133,133,132 00007450 133 NMZ = NMZ + 1 00007460 GO TO 131 00007470 132 A1=A1+U(L)**DV 00007480 F1=F1+U(L)**DC 00007490 U(L)=U(L)**UEXP 00007500 SUMU=SUMU+U(L) 00007510 131 CONTINUE 00007520 ZD=(2.*F1)**BC 00007530 Z=1./(2.*A1*ZD) 00007540 XNS=NS*(NS-1) 00007550 CMAX=Z*2.*SUMU*XNS**BC 00007560 WRITE(6,321) CMAX 00007570 IF(NMZ .GT. 0) WRITE(6,134) NMZ 00007580 134 FORMAT(//5X,'THERE WERE',I6,' ZERO DISTANCES'/) 00007590 IF(NMZ .EQ. 0) WRITE(6,135) 00007600 135 FORMAT(//5X,'THERE WERE NO ZERO DISTANCES'/) 00007610 XNS=NS 00007620 XM=XNS-1. 00007630 TK=2./XM 00007640 NO=NPO 00007650 IT=0 00007660 CALL NRMX(NF,NS,X,IT,SA) 00007670 IT=1 00007680 500 CALL MINK 00007690 81 IF(RGX(IT)-GCRIT)281,281,282 00007700 281 IF(RGSA(IT)-GCRIT)284,284,282 00007710 284 FINI=0 00007720 PRINT316,GCRIT 00007730 GOTO94 00007740 C TEST END OF ITERATION 00007750 282 IF(IT-MAXIT)91,283,283 00007760 283 FINI=0 00007770 WRITE(6,315) MAXIT 00007780 GOTO94 00007790 91 IF(IT-1)90,94,95 00007800 95 IF(IT-NPO)97,90,90 00007810 90 NPO=NPO+NO 00007820 94 CALL PPX 00007830 97 CALL ALPHA 00007840 IF(FINI)89,89,15 00007850 15 IT=IT+1 00007860 GOTO500 00007870 315 FORMAT(///5X,'REACHED MAXIMUM OF',I5,' ITERATIONS'/) 00007880 316 FORMAT('0REACHED CRITERION GCRIT='E14.7) 00007890 321 FORMAT(///5X,'KAPPA MAX = ',F14.6/) 00007900 89 RETURN 00007910 END 00007920 SUBROUTINE MINK 00007930 DIMENSION RR(150),U(11175),X(150,10),DKX(150,10),TIT(20), 00007940 1 COSA(100),RGX(100),RGSA(100),ALF(100),RLGR(100),XAL(100),AS(100),00007950 2 XL(100),C(100) 00007960 COMMON U,X,DKX,TIT,KPUNCH,IT,FINI,NS,NF,UEXP,DB,DC,CALF,SA,TK, 00007970 1 DB1,DC1,Z,BC,GV,COSA,DSA,RGX,RGSA,ALF,C,RLGR,XAL,AS,XL,PCUT 00007980 XFUNF(I,J)=((I-1)*(I-2))/2+J 00007990 21 A=0 00008000 B=0 00008010 SDXSQ=0 00008020 SXSQ=0 00008030 GV=0 00008040 L=0 00008050 TKSA=TK*SA 00008060 SASQ=TKSA*SA 00008070 DO40I=2,NS 00008080 JJ=I-1 00008090 DO40J=1,JJ 00008100 L=L+1 00008110 SUML=0 00008120 DO35M=1,NF 00008130 35 SUML=SUML+(X(I,M)-X(J,M))**2 00008140 DEL=SUML+SASQ 00008150 B=B+DEL**DC 00008160 40 A=A+U(L)*DEL**(DB1+1.) 00008170 A=2.*A 00008180 B=2.*B 00008190 C(IT)=Z*A*B**BC 00008200 ZB=4.*Z*DB*B**(BC-1.) 00008210 C COMPUTE S 00008220 L=0 00008230 SUMS=0 00008240 DO45J=1,NS 00008250 RR(J)=0 00008260 DO44I=1,NS 00008270 IF(I-J)41,44,43 00008280 43 L=XFUNF(I,J) 00008290 GOTO46 00008300 41 L=XFUNF(J,I) 00008310 46 SUML=0 00008320 DO54M=1,NF 00008330 54 SUML=SUML+(X(I,M)-X(J,M))**2 00008340 DEL=SUML+SASQ 00008350 RR(I)=+ZB*(B*U(L)*DEL**DB1-A*DEL**DC1) 00008360 132 RR(J)=RR(J)-RR(I) 00008370 44 CONTINUE 00008380 SUMS=SUMS+RR(J) 00008390 C COMPUTE DELTA X 00008400 DO110K=1,NF 00008410 DX=0 00008420 DO111I=1,NS 00008430 111 DX=DX+X(I,K)*RR(I) 00008440 IF(IT-1)61,62,61 00008450 61 GV=GV+DX*DKX(J,K) 00008460 62 DKX(J,K)=DX 00008470 SXSQ=SXSQ+X(J,K)**2 00008480 110 SDXSQ=SDXSQ+DKX(J,K)**2 00008490 45 CONTINUE 00008500 IF(IT .EQ. 1) GO TO 74 00008510 GV=GV/XGI 00008520 C COMPUTE DELTA SA 00008530 74 DSA = TKSA/2.*SUMS 00008540 XSAL=SXSQ+SA**2 00008550 GI=SQRT (SDXSQ+DSA**2) 00008560 XAL(IT)=SQRT (XSAL) 00008570 XGI=XAL(IT)/GI 00008580 RLGR(IT)=GI*XAL(IT) 00008590 XL(IT)=SQRT (SXSQ) 00008600 IF(IT-1)64,64,19 00008610 19 GV=GV+DSA*DSA1 00008620 COSA(IT)=GV/(GI*GI1) 00008630 64 DSA1=DSA 00008640 GI1=GI 00008650 RGX(IT)=SQRT (SDXSQ)*XL(IT) 00008660 RGSA(IT)=DSA*ABS (SA) 00008670 AS(IT)=SA 00008680 DO120K=1,NF 00008690 DO120J=1,NS 00008700 120 DKX(J,K)=DKX(J,K)*XGI 00008710 DSA=DSA*XGI 00008720 RETURN 00008730 END 00008740 SUBROUTINE NRMX(NF,NS,X,IT,SA) 00008750 DIMENSION X(150,10) 00008760 C NORMALIZE X 00008770 SSQX=0 00008780 DO535I=1,NF 00008790 SUMX=0 00008800 DO540J=1,NS 00008810 540 SUMX=SUMX+X(J,I) 00008820 XMEAN=SUMX/FLOAT (NS) 00008830 DO541J=1,NS 00008840 X(J,I)=X(J,I)-XMEAN 00008850 541 SSQX=SSQX+X(J,I)**2 00008860 535 CONTINUE 00008870 XS=SSQX+SA**2 00008880 IF(IT)100,534,100 00008890 534 CX=SQRT (FLOAT (NF)/XS) 00008900 GOTO20 00008910 100 CX=1./SQRT (XS) 00008920 SA=SA*CX 00008930 20 DO538I=1,NF 00008940 DO538J=1,NS 00008950 538 X(J,I)=X(J,I)*CX 00008960 RETURN 00008970 END 00008980 SUBROUTINE PPX 00008990 DIMENSION DIMA(150),DIMB(150),U(11175),X(150,10),DKX(150,10), 00009000 1 TIT(20),COSA(100),RGX(100),RGSA(100),ALF(100),RLGR(100),XAL(100),00009010 2 AS(100),XL(100),C(100) 00009020 COMMON U,X,DKX,TIT,KPUNCH,IT,FINI,NS,NF,UEXP,DB,DC,CALF,SA,TK, 00009030 1 DB1,DC1,Z,BC,GV,COSA,DSA,RGX,RGSA,ALF,C,RLGR,XAL,AS,XL,PCUT, 00009040 2 DIMA,DIMB 00009050 100 IF(FINI) 70,16,20 00009060 20 IF(IT - 1) 5,5,16 00009070 5 WRITE(6,315) 00009080 315 FORMAT(//5X,'INITIAL X MATRIX') 00009090 GO TO 58 00009100 70 WRITE(6,308) 00009110 308 FORMAT(//5X,'NORMALIZED X MATRIX') 00009120 IF(NF .GT. 1) WRITE(6,324) 00009130 324 FORMAT(1H+,25X,' -- ROTATED TO PRINCIPAL AXES') 00009140 GO TO 58 00009150 16 WRITE(6,311) IT 00009160 311 FORMAT(//5X,'X MATRIX, ITERATION',I5) 00009170 58 WRITE(6,320) (M,M=1,NF) 00009180 320 FORMAT(//10I12) 00009190 DO 53 I = 1,NS 00009200 53 WRITE(6,306) I,(X(I,J), J = 1,NF) 00009210 306 FORMAT(I4,10F12.5) 00009220 IF(KPUNCH .EQ. 0) GO TO 55 00009230 GO TO (61,61,62), KPUNCH 00009240 62 IF(FINI .NE. -1.0) GO TO 55 00009250 61 WRITE(7,313) IT 00009260 313 FORMAT(2X,'X - MATRIX, ITERATION',I4/2X,'(I3,9F8.4/(3X,9F8.4))') 00009270 DO 57 J = 1,NS 00009280 57 WRITE(7,307) J,(X(J,K), K=1,NF) 00009290 307 FORMAT(I3,9F8.4/(3X,9F8.4)) 00009300 55 CONTINUE 00009310 IF(NF .EQ. 1) GO TO 73 00009320 NFM1 = NF - 1 00009330 DO 80 I = 1,NFM1 00009340 IP1 = I + 1 00009350 DO 71 J = IP1,NF 00009360 IF(FINI) 330,332,334 00009370 330 WRITE(6,331) 00009380 331 FORMAT(1H1,5X,'NORMALIZED-ROTATED X-MATRIX') 00009390 GO TO 340 00009400 332 WRITE(6,333) 00009410 333 FORMAT(1H1,5X,'FINAL ITERATION X-MATRIX') 00009420 GO TO 340 00009430 334 IF(IT .EQ. 1) GO TO 335 00009440 GO TO 337 00009450 335 WRITE(6,336) 00009460 336 FORMAT(1H1,5X,'PLOT OF INITIAL X-MATRIX') 00009470 GO TO 340 00009480 337 WRITE(6,338) 00009490 338 FORMAT(1H1,5X,'PLOT OF X-MATRIX') 00009500 340 WRITE(6,322) I,J,IT,NF 00009510 322 FORMAT(1H+,35X,'DIM',I3,' (X-AXIS) VS. DIM',I3,' (Y-AXIS) -- ITER'00009520 1,I3,' OF',I3,' DIM SOLUTION') 00009530 DO 72 K = 1,NS 00009540 DIMA(K) = X(K,I) 00009550 DIMB(K) = X(K,J) 00009560 72 CONTINUE 00009570 AAA = 0.0 00009580 BBB = 0.0 00009590 CCC = 0.0 00009600 DDD = 0.0 00009610 CALL PLOTX(DIMA,DIMB,NS,1,1,1,AAA,BBB,CCC,DDD) 00009620 71 CONTINUE 00009630 80 CONTINUE 00009640 GO TO 75 00009650 73 WRITE(6,323) IT 00009660 323 FORMAT(1H1,40X,'X PLOT FOR ONE-DIMENSIONAL SOLUTION -- ITERATION',00009670 1 I3) 00009680 CALL CLEAR(DIMB,150) 00009690 DO 74 K = 1,NS 00009700 74 DIMA(K) = X(K,1) 00009710 AAA = 0.0 00009720 BBB = 0.0 00009730 CCC = +1.0 00009740 DDD = -1.0 00009750 CALL PLOTX(DIMA,DIMB,NS,1,1,0,AAA,BBB,CCC,DDD) 00009760 75 CONTINUE 00009770 IF(FINI .GT. 0.0) WRITE(6,321) NF 00009780 321 FORMAT(////30X,'CONTINUE ITERATIONS IN',I5,' DIMENSIONS'////) 00009790 RETURN 00009800 END 00009810 SUBROUTINE ALPHA 00009820 DIMENSION U(11175),X(150,10),DKX(150,10),TIT(20),COSA(100), 00009830 1 RGX(100),RGSA(100),ALF(100),RLGR(100),XAL(100),AS(100),XL(100), 00009840 2 C(100) 00009850 COMMON U,X,DKX,TIT,KPUNCH,IT,FINI,NS,NF,UEXP,DB,DC,CALF,SA,TK, 00009860 1 DB1,DC1,Z,BC,GV,COSA,DSA,RGX,RGSA,ALF,C,RLGR,XAL,AS,XL,PCUT 00009870 BOARG = RLGR(IT) * (C(IT) - 1.0) 00009880 IF(BOARG) 200,201,201 00009890 200 BOARG = -BOARG 00009900 201 BO = SQRT(BOARG) 00009910 IF(IT-1)2,1,2 00009920 1 ALF(IT)=CALF 00009930 QP=1.-PCUT 00009940 GO TO 91 00009950 2 Q=8.**(COSA(IT)**5) 00009960 ALF(IT)=ALF(IT-1)*Q 00009970 91 AALF=BO*ALF(IT) 00009980 41 IF(FINI)96,100,96 00009990 96 KQ=0 00010000 98 DO13I=1,NF 00010010 DO13J=1,NS 00010020 13 X(J,I)=X(J,I)+AALF*DKX(J,I) 00010030 SA=SA+AALF*DSA 00010040 C COMPUTE NEW K 00010050 21 A=0 00010060 B=0 00010070 L=0 00010080 TKSA=TK*SA 00010090 SASQ=TKSA*SA 00010100 DO40I=2,NS 00010110 JJ=I-1 00010120 DO40J=1,JJ 00010130 L=L+1 00010140 SUML=0 00010150 DO35M=1,NF 00010160 35 SUML=SUML+(X(I,M)-X(J,M))**2 00010170 DEL=SUML+SASQ 00010180 B=B+DEL**DC 00010190 40 A=A+U(L)*DEL**(DB1+1.) 00010200 A=2.*A 00010210 B=2.*B 00010220 CA=Z*A*B**BC 00010230 IF(CA-C(IT))95,95,99 00010240 99 IF(KQ)105,97,105 00010250 97 AALF=-PCUT*AALF 00010260 GOTO85 00010270 105 AALF=QP*AALF 00010280 85 KQ=KQ+1 00010290 IF(KQ .GE. 10) GO TO 300 00010300 GOTO98 00010310 95 ALF(IT)=ALF(IT)*QP**KQ 00010320 WRITE(6,101) IT,KQ,BO 00010330 101 FORMAT(30X,'ITERATION',I5,' ',5X,'KQ = ',I3,5X,'BO = ',F12.6) 00010340 45 CALL NRMX(NF,NS,X,IT,SA) 00010350 100 RETURN 00010360 300 FINI = 0.0 00010370 CALL NRMX(NF,NS,X,IT,SA) 00010380 WRITE(6,102) IT 00010390 102 FORMAT(//5X,'*****COMPUTATION TERMINATED AT ITERATION',I5,' DUE TO00010400 1 STEP SIZE UNDERFLOW IN SUBROUTINE "ALPHA"*****') 00010410 RETURN 00010420 END 00010430 SUBROUTINE NORM 00010440 DIMENSION XM(10),XN(150,10),R(10,10),V(10,10),U(11175),X(150,10), 00010450 1 SYMR(55), XV(100),TIT(20) 00010460 COMMON U,X,XN,TIT,KPUNCH,IT,FINI,NS,NF 00010470 EQUIVALENCE (XN(1),SYMR(1)), (XN(56),XM(1)),(XN(66),R(1)), 00010480 1 (XN(166),XV(1)) 00010490 XNS=NS 00010500 NN=(NS*(NS-1))/2 00010510 CALL CLEAR(R,100,XN,1500) 00010520 DO6J=1,NF 00010530 SMEAN=0 00010540 DO5I=1,NS 00010550 5 SMEAN=SMEAN+X(I,J) 00010560 6 XM(J)=SMEAN/XNS 00010570 WRITE(6,209) 00010580 WRITE(6,213) (K, K=1,NF) 00010590 WRITE(6,206) (XM(J), J=1,NF) 00010600 DO 7 I = 1,NF 00010610 DO 7 J = 1,NS 00010620 7 X(J,I)=X(J,I)-XM(I) 00010630 WRITE(6,214) 00010640 WRITE(6,213) (K, K=1,NF) 00010650 DO 14 I = 1,NS 00010660 14 WRITE(6,212) I,(X(I,J), J=1,NF) 00010670 IF(NF .EQ. 1) GO TO 300 00010680 DO 9 I = 1,NF 00010690 DO 9 J = 1,NF 00010700 DO 9 K = 1,NS 00010710 9 R(I,J)=R(I,J)+X(K,I)*X(K,J) 00010720 WRITE(6,210) 00010730 DO 15 I = 1,NF 00010740 15 WRITE(6,212) I,(R(I,J), J = 1,NF) 00010750 L = 0 00010760 DO 8 J = 1,NF 00010770 DO 8 I = 1,J 00010780 L = L + 1 00010790 8 SYMR(L) = R(I,J) 00010800 CALL EIGEN(SYMR,XV,NF,0) 00010810 L = 0 00010820 DO 13 J = 1,NF 00010830 DO 13 I = 1,NF 00010840 L = L + 1 00010850 13 V(I,J) = XV(L) 00010860 L = 0 00010870 DO 11 I = 1,NF 00010880 L = L + I 00010890 11 R(I,I) = SYMR(L) 00010900 WRITE(6,205) 00010910 WRITE(6,213) (K, K=1,NF) 00010920 WRITE(6,206) (R(I,I), I = 1,NF) 00010930 WRITE(6,211) 00010940 WRITE(6,213) (K, K=1,NF) 00010950 DO 18 I = 1,NF 00010960 18 WRITE(6,212) I,(V(I,J), J = 1,NF) 00010970 CALL CLEAR(XN,1500) 00010980 DO 10 I = 1,NS 00010990 DO 10 J = 1,NF 00011000 DO 10 K = 1,NF 00011010 10 XN(I,J)=XN(I,J)+X(I,K)*V(K,J) 00011020 DO 12 J = 1,NF 00011030 DO 12 I = 1,NS 00011040 12 X(I,J)=XN(I,J) 00011050 300 FINI = -1.0 00011060 CALL PPX 00011070 205 FORMAT(//10X,'EIGEN ROOTS OF R - MATRIX (FACTOR IMPORTANCE)'/) 00011080 206 FORMAT(5X,10F12.5) 00011090 209 FORMAT(1H1,30X,'***** NORMALIZATION OF X - MATRIX, AND ROTATION TO00011100 1 PRINCIPAL AXES *****'////10X,'PRE-NORMALIZED X-MATRIX FACTOR MEAN00011110 2S'//) 00011120 210 FORMAT(///10X,'R - MATRIX, SUMS OF CROSS-PRODUCTS OF DEVIATIONS FR00011130 1OM MEANS OF THE X - MATRIX'/) 00011140 211 FORMAT(///10X,'V - MATRIX, EIGEN VECTORS OF R IN COLUMNS'/) 00011150 212 FORMAT(I5,10F12.5) 00011160 213 FORMAT(10I12) 00011170 214 FORMAT(/10X,'X-MATRIX, DEVIATIONS FROM MEANS'/) 00011180 RETURN 00011190 END 00011200 SUBROUTINE EIGEN(A,R,N,MV) 00011210 DIMENSION A(1),R(1) 00011220 C 00011230 C SUBROUTINE EIGEN 00011240 C 00011250 C PURPOSE 00011260 C COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC 00011270 C MATRIX 00011280 C 00011290 C USAGE 00011300 C CALL EIGEN(A,R,N,MV) 00011310 C 00011320 C DESCRIPTION OF PARAMETERS 00011330 C A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION. 00011340 C RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF 00011350 C MATRIX A IN DECENDING ORDER. 00011360 C R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, 00011370 C IN SAME SEQUENCE AS EIGENVALUES) 00011380 C N - ORDER OF MATRICES A AND R 00011390 C MV- INPUT CODE 00011400 C 0 COMPUTE EIGENVALUES AND EIGENVECTORS 00011410 C 1 COMPUTE EIGENVALUES ONLY (R NEED NOT BE 00011420 C DIMENSIONED BUT MUST STILL APPEAR IN CALLING 00011430 C SEQUENCE) 00011440 C 00011450 C REMARKS 00011460 C ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1) 00011470 C MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R 00011480 C 00011490 C SUBROUTINES USED BY THIS SUBROUTINE 00011500 C NONE 00011510 C 00011520 C METHOD 00011530 C DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED 00011540 C BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN MATHEMATICAL 00011550 C METHODS FOR DIGITAL COMPUTERS, EDITED BY A. RALSTON AND 00011560 C H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7 00011570 C 00011580 C ..................................................................00011590 C 00011600 C 00011610 C GENERATE IDENTITY MATRIX 00011620 C 00011630 IF(MV-1) 22,14,22 00011640 22 IQ=-N 00011650 DO 1 J=1,N 00011660 IQ=IQ+N 00011670 DO 1 I=1,N 00011680 IJ=IQ+I 00011690 R(IJ)=0.0 00011700 IF(I-J) 1,31,1 00011710 31 R(IJ)=1.0 00011720 1 CONTINUE 00011730 C 00011740 C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) 00011750 C 00011760 14 ANORM=0.0 00011770 DO 2 I=1,N 00011780 DO 2 J=I,N 00011790 IF(I-J) 15,2,15 00011800 15 IA=I+(J*J-J)/2 00011810 ANORM=ANORM+A(IA)*A(IA) 00011820 2 CONTINUE 00011830 IF(ANORM) 39, 39, 50 00011840 50 ANORM=2.0*SQRT(ANORM) 00011850 ANRMX =ANORM*1.0E-8 00011860 C 00011870 C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR 00011880 C 00011890 IND=0 00011900 THR=ANORM 00011910 3 THR=THR/FLOAT(N) 00011920 13 L=1 00011930 4 M=L+1 00011940 C 00011950 C COMPUTE SIN AND COS 00011960 C 00011970 5 MQ=(M*M-M)/2 00011980 LQ=(L*L-L)/2 00011990 LM=L+MQ 00012000 IF(ABS(A(LM))-THR) 7,32,32 00012010 32 IND=1 00012020 LL=L+LQ 00012030 MM=M+MQ 00012040 X=0.5*(A(LL)-A(MM)) 00012050 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) 00012060 IF(X) 33,34,34 00012070 33 Y=-Y 00012080 34 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y)))) 00012090 SINX2=SINX*SINX 00012100 COSX=SQRT(1.0-SINX2) 00012110 COSX2=COSX*COSX 00012120 SINCS =SINX*COSX 00012130 C 00012140 C ROTATE L AND M COLUMNS 00012150 C 00012160 ILQ=N*(L-1) 00012170 IMQ=N*(M-1) 00012180 DO 6 I=1,N 00012190 IQ=(I*I-I)/2 00012200 IF(I-L) 16,11,16 00012210 16 IF(I-M) 17,11,18 00012220 17 IM=I+MQ 00012230 GO TO 19 00012240 18 IM=M+IQ 00012250 19 IF(I-L) 37,23,23 00012260 37 IL=I+LQ 00012270 GO TO 24 00012280 23 IL=L+IQ 00012290 24 X=A(IL)*COSX-A(IM)*SINX 00012300 A(IM)=A(IL)*SINX+A(IM)*COSX 00012310 A(IL)=X 00012320 11 IF(MV-1) 26,6,26 00012330 26 ILR=ILQ+I 00012340 IMR=IMQ+I 00012350 X=R(ILR)*COSX-R(IMR)*SINX 00012360 R(IMR)=R(ILR)*SINX+R(IMR)*COSX 00012370 R(ILR)=X 00012380 6 CONTINUE 00012390 X=2.0*A(LM)*SINCS 00012400 Y=A(LL)*COSX2+A(MM)*SINX2-X 00012410 X=A(LL)*SINX2+A(MM)*COSX2+X 00012420 A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) 00012430 A(LL)=Y 00012440 A(MM)=X 00012450 C 00012460 C TESTS FOR COMPLETION 00012470 C 00012480 C TEST FOR M = LAST COLUMN 00012490 7 IF(M-N) 35,8,35 00012500 35 M=M+1 00012510 GO TO 5 00012520 C TEST FOR L = SECOND FROM LAST COLUMN 00012530 8 IF(L-(N-1)) 36,9,36 00012540 36 L=L+1 00012550 GO TO 4 00012560 9 IF(IND-1) 10,38,10 00012570 38 IND=0 00012580 GO TO 13 00012590 C COMPARE THRESHOLD WITH FINAL NORM 00012600 10 IF(THR-ANRMX ) 39,39,3 00012610 C 00012620 C SORT EIGENVALUES AND EIGENVECTORS 00012630 C 00012640 39 IQ=-N 00012650 DO 30 I=1,N 00012660 IQ=IQ+N 00012670 LL=I+(I*I-I)/2 00012680 JQ=N*(I-2) 00012690 DO 30 J=I,N 00012700 JQ=JQ+N 00012710 MM=J+(J*J-J)/2 00012720 IF(A(LL)-A(MM)) 40,30,30 00012730 40 X=A(LL) 00012740 A(LL)=A(MM) 00012750 A(MM)=X 00012760 IF(MV-1) 28,30,28 00012770 28 DO 20 K=1,N 00012780 ILR=IQ+K 00012790 IMR=JQ+K 00012800 X=R(ILR) 00012810 R(ILR)=R(IMR) 00012820 20 R(IMR)=X 00012830 30 CONTINUE 00012840 RETURN 00012850 END 00012860 SUBROUTINE PLOTX(Z,Y,NPOI,AXES,IDENT,IP,XMAX,XMIN,YMAX,YMIN) 00012870 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC00012880 C 00012890 C ROUTINE TO GENERATE A SQUARE PLOT OF ARRAY X VS. Y. 00012900 C 00012910 C THE INPUT VECTORS HAVE NPOI ENTRIES (A NUMBER LESS THAN OR EQUAL TO 00012920 C THE LENGTH OF THE ARGUMENT ARRAYS IN THE CALLING PROGRAM. 00012930 C 00012940 C THE PLOTTING IS DONE ON TAPE -OUT-. 00012950 C 00012960 C IF -AXES- IS POSITIVE, AXES WILL APPEAR ON THE PLOT. 00012970 C IF -AXES- IS NEGATIVE, NO AXES WILL APPEAR. 00012980 C 00012990 C 00013000 C IF -XMAX- = -XMIN- THE PROGRAM WILL DETERMINE ITS OWN SCALE FOR 00013010 C THE X AXIS, AND SIMILARLY FOR THE Y AXIS IF -YMAX- = -YMIN-. 00013020 C IF -XMAX- IS NOT EQUAL TO -XMIN-, THE PROGRAM WILL USE THESE VALUES 00013030 C FOR THE UPPER AND LOWER BOUNDS OF THE X AXIS, AND SIMILARLY FOR THE 00013040 C Y AXIS. 00013050 C IF -IP- = ZERO, AXIS SCALES WILL NOT BE MADE EQUAL. IF -IP- = 1, 00013060 C THE AXIS SCALES WILL BE MADE EQUAL. 00013070 C 00013080 C IF -IDENT- IS POSITIVE, THE PLOTTED POINTS WILL BE IDENTIFIED BY 00013090 C NUMBER. THERE ARE FIFTY (50) IDENTIFICATION CHARACTERS AVAILABLE.00013100 C IF -IDENT- IS NEGATIVE, THE PLOTTED POINTS WILL APPEAR AS "X", AND 00013110 C MULTIPLE POINTS WILL BE IDENTIFIED BY THE NUMBER OF SUPERIMPOSED 00013120 C POINTS AT THAT LOCATION. ***NOTE*** THIS OPTION SHOULD BE USED WHEN 00013130 C NPOI SIGNIFICANTLY EXCEEDS 50. 00013140 C 00013150 C WRITTEN BY FORREST W. YOUNG, NOVEMBER 1965. 00013160 C VERSION 3, APRIL 1967 00013170 C 00013180 C----------ADAPTED BY F.J. CARMONE, JR., 5/1/67. 00013190 C----------MODIFIED BY R.F. MCCRACKEN 11/15/67. 00013200 C 00013210 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC00013220 C 00013230 DIMENSION Z(1),Y(1),SMALL(13),CHAR(50),HOLL(11) 00013240 REAL ITEM(55,121) 00013250 INTEGER OUT,AXES 00013260 LOGICAL FLAG/.FALSE./ 00013270 DATA CHAR/'1','2','3','4','5','6','7','8','9','A','B','C','D', 00013280 1'E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T', 00013290 2'U','V','W','X','Y','Z','+','/','=','*','&','$','@','%','?','<', 00013300 3'(',')','"','#',' '/ 00013310 DATA HOLL/' ','X','2','3','4','5','6','7','8','9','M'/ 00013320 DATA BLANK/' '/,DD/'.'/,PLUS/'+'/,AIE/'\'/,AMINUS/'-'/,ORIG/'0'/, 00013330 1 AMULT/' '/,GT50/'>'/ 00013340 OUT = 6 00013350 IF(FLAG) GO TO 400 00013360 DO 115 I = 1,55 00013370 DO 115 J=1,121 00013380 115 ITEM(I,J) = BLANK 00013390 FLAG = .TRUE. 00013400 400 CONTINUE 00013410 IF(XMAX .NE. XMIN) GO TO 222 00013420 XMAX = -1.0E7 00013430 XMIN = +1.0E7 00013440 DO 221 I = 1,NPOI 00013450 IF(Z(I) .GT. XMAX) XMAX = Z(I) 00013460 IF(Z(I) .LT. XMIN) XMIN = Z(I) 00013470 221 CONTINUE 00013480 XEXP = 0.1 * (XMAX - XMIN) 00013490 XMAX = XMAX + XEXP 00013500 XMIN = XMIN - XEXP 00013510 222 IF(YMAX .NE. YMIN) GO TO 224 00013520 YMAX = -1.0E7 00013530 YMIN = +1.0E7 00013540 DO 223 I = 1,NPOI 00013550 IF(Y(I) .GT. YMAX) YMAX = Y(I) 00013560 IF(Y(I) .LT. YMIN) YMIN = Y(I) 00013570 223 CONTINUE 00013580 YEXP = 0.1 * (YMAX - YMIN) 00013590 YMAX = YMAX + YEXP 00013600 YMIN = YMIN - YEXP 00013610 224 IF(IP .NE. 1) GO TO 226 00013620 IF(YMAX - XMAX) 230,232,231 00013630 230 YMAX = XMAX 00013640 GO TO 232 00013650 231 XMAX = YMAX 00013660 232 CONTINUE 00013670 IF(YMIN - XMIN) 233,225,234 00013680 233 XMIN = YMIN 00013690 GO TO 225 00013700 234 YMIN = XMIN 00013710 225 CONTINUE 00013720 EXP = 0.1667 * (XMAX - XMIN) 00013730 XMAX = XMAX + EXP 00013740 XMIN = XMIN - EXP 00013750 226 CONTINUE 00013760 SPANX = XMAX - XMIN 00013770 SPANY = YMAX - YMIN 00013780 DELX = SPANX/120.0 00013790 DELY = SPANY/54.0 00013800 IF(AXES .LE. 0) GO TO 119 00013810 K = 0 00013820 M = 0 00013830 IF(XMIN .LT. 0 .AND. XMAX .GT. 0) GO TO 200 00013840 GO TO 202 00013850 200 K = (-XMIN/DELX) + 1.5 00013860 DO 201 I = 1,55 00013870 201 ITEM(I,K) = AIE 00013880 202 CONTINUE 00013890 IF(YMIN .LT. 0 .AND. YMAX .GT. 0) GO TO 203 00013900 GO TO 119 00013910 203 M = (-YMIN/DELY) + 1.5 00013920 M = 56 - M 00013930 DO 204 J = 1,121 00013940 204 ITEM(M,J) = AMINUS 00013950 IF(K .EQ. 0 .OR. M .EQ. 0) GO TO 119 00013960 ITEM(M,K) = ORIG 00013970 119 CONTINUE 00013980 VALUE = YMAX + DELY 00013990 SMALL(1) = XMIN 00014000 DO 120 I = 1,12 00014010 120 SMALL(I+1) = SMALL(I) + 10.0*DELX 00014020 IF(IDENT .GT. 0) GO TO 140 00014030 DO 135 II=1,NPOI 00014040 I = (YMAX - Y(II))/DELY + 1.5 00014050 J = (Z(II) - XMIN)/DELX + 1.5 00014060 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 135 00014070 DO 134 JJ=1,10 00014080 IF(ITEM(I,J) .EQ. HOLL(JJ)) GO TO 136 00014090 134 CONTINUE 00014100 IF(ITEM(I,J) .EQ. AIE .OR. ITEM(I,J) .EQ. AMINUS .OR. ITEM(I,J) 00014110 1 .EQ. ORIG) ITEM(I,J) = HOLL(2) 00014120 GO TO 135 00014130 136 ITEM(I,J) = HOLL(JJ+1) 00014140 135 CONTINUE 00014150 GO TO 141 00014160 140 CONTINUE 00014170 DO 137 II = 1,NPOI 00014180 I = (YMAX - Y(II))/DELY + 1.5 00014190 J = (Z(II) - XMIN)/DELX + 1.5 00014200 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 137 00014210 IF(ITEM(I,J) .EQ. BLANK .OR. ITEM(I,J) .EQ. AIE .OR. ITEM(I,J) 00014220 1 .EQ. AMINUS .OR. ITEM(I,J) .EQ. ORIG) GO TO 138 00014230 ITEM(I,J) = AMULT 00014240 GO TO 137 00014250 138 CONTINUE 00014260 IF(II .GT. 50) GO TO 139 00014270 ITEM(I,J) = CHAR(II) 00014280 GO TO 137 00014290 139 ITEM(I,J) = GT50 00014300 137 CONTINUE 00014310 141 CONTINUE 00014320 WRITE(OUT,3001) 00014330 3001 FORMAT(9X,' +.........+.........+.........+.........+.........+...00014340 1......+.........+.........+.........+.........+.........+.........00014350 2+') 00014360 DO 330 I = 1,55 00014370 VALUE = VALUE - DELY 00014380 II = I - 1 00014390 L = II + 6 00014400 IF(L/6*6 - L) 310,320,310 00014410 320 B = PLUS 00014420 WRITE(OUT,3300) VALUE,B, (ITEM(I,J), J=1,121), B 00014430 3300 FORMAT(1X,F8.3,123A1) 00014440 GO TO 330 00014450 310 B = DD 00014460 WRITE(OUT,3303) B, (ITEM(I,J), J=1,121), B 00014470 3303 FORMAT(9X,123A1) 00014480 330 CONTINUE 00014490 WRITE(OUT,3001) 00014500 WRITE(OUT,3302) (SMALL(K), K=1,12) 00014510 3302 FORMAT(4X,12F10.3) 00014520 WRITE(OUT,3301) SMALL(13) 00014530 3301 FORMAT(1H+,123X,F8.2/) 00014540 IF(AXES .LE. 0) GO TO 302 00014550 K = (-XMIN/DELX) + 1.5 00014560 DO 300 I = 1,55 00014570 300 ITEM(I,K) = BLANK 00014580 M = (-YMIN/DELY) + 1.5 00014590 M = 56 - M 00014600 DO 301 J = 1,121 00014610 301 ITEM(M,J) = BLANK 00014620 302 CONTINUE 00014630 DO 303 II = 1,NPOI 00014640 I = (YMAX - Y(II))/DELY + 1.5 00014650 J = (Z(II) - XMIN)/DELX + 1.5 00014660 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 303 00014670 ITEM(I,J) = BLANK 00014680 303 CONTINUE 00014690 RETURN 00014700 END 00014710 SUBROUTINE RANDU(IX,IY,YFL) 00014720 IY = IX * 65539 00014730 IF(IY) 5,6,6 00014740 5 IY = IY + 2147483647 + 1 00014750 6 YFL = IY 00014760 YFL = YFL * .4656613E-9 00014770 RETURN 00014780 END 00014790 SUBROUTINE IDENT 00014800 DIMENSION CHAR(50) 00014810 DATA CHAR/'1','2','3','4','5','6','7','8','9','A','B','C','D', 00014820 1'E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T', 00014830 2'U','V','W','X','Y','Z','+','/','=','*','&','$','@','%','?','<', 00014840 3'(',')','"','#',' '/ 00014850 WRITE(6,151) (J, J=1,25) 00014860 151 FORMAT(/45X,'POINT IDENTIFICATION KEY FOR PLOTS'//' PT #',25I5/) 00014870 WRITE(6,152) (CHAR(I), I=1,25) 00014880 152 FORMAT(' CHAR',25(4X,A1)/) 00014890 WRITE(6,153) (J, J=26,50) 00014900 153 FORMAT(1X,'PT #',25I5/) 00014910 WRITE(6,154) (CHAR(I), I=26,50) 00014920 154 FORMAT(1X,'CHAR',25(4X,A1)) 00014930 WRITE(6,155) 00014940 155 FORMAT(////1X,'ALL POINT NUMBERS ABOVE 50 ARE DESIGNATED BY THE SY00014950 1MBOL >'//1X,'MULTIPLE POINTS AT THE SAME LOCATION ARE DESIGNATED00014960 2 BY THE SYMBOL ;') 00014970 RETURN 00014980 END 00014990