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