C idioscal.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 Carroll,J.D. & Chang,J.J. (1970), "Analysis of individual differences C in multidimensional scaling via an N-way generalization of C 'Eckart-Young' decomposition" in Psychometrika, 35,283-319, and C Carroll,J.D. (1972). Individual differences and multidimensional C scaling, in R.N.Shepard, A.K.Romney & S.B. Nerlove (Eds.), C Multidimensional Scaling: Theory and Applications in the Behavioral C Sciences (Vol.1, pp.105-155). New York and London: Seminar Press. C Both of these have been reprinted in C P. Davies & A.P.M. Coxon (Eds.), (1984) C "Key Texts in Multidimensional Scaling" Exeter, NH: Heinemann. C First one: pp. 229-252; second: pp. 267-301. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ $ FORTRAN C IDIOSCAL-INDIVIDUAL DIFFERENCES IN ORIENTATION SCALING C 4/11/72, WEIGHTED Q WITH CONSTANT C BY CARROLL AND CHANG DIMENSION TIT(12),BNUM(38),FMT(12),TSUB(4) DIMENSION B(38,38),D(12000),A(780),Q(6,6),GSUB(5),S(21,21) DIMENSION SP(38,38),X(38,6),T(38),COR(38,3),SQR(38,3) DATA GSUB(1)/30HINDSCAL GROUP STIMULUS SPACE / DATA TSUB(1)/24HGROUP STIMULUS SPACE / C N - NUMBER OF POINTS. (MAX.=38) C K - NUMBER OF DIMENSIONS. (MAX.=6) C NSUB - NUMBER OF SUBJECTS. ((N*(N-1))/2)*NSUB .LE. 12000 C AND MAX. NSUB IS 38. C IRDATA- 1, SIMILARITIES C 2, DISSIMILARITIES C 3, EUCLIDEAN DISTANCES C MIDEN- 1, READ IN LABELS FOR PLOTTING. C 0, NO C IPUN - 1, PUNCH TRANSFORMED INDIVIDUAL COORDINATES. C 0, NO C NORM- 0, DO NOT NORMALIZE DATA. C 1, NORMALIZE SQUARED DATA, MAKE SUM SQUARE=1. C 2, NORMALIZE SCALAR PRODUCT, MAKE SUM SQUARE OF SP=1. C IPLOT - 1, PLOT INDIVIDUAL SPACE IN PHASE 1. C 0, NO INDIVIDUAL PLOTS. NK=6 NL=38 NK2=(NK*(NK&1))/2 98 READ100,N,K,NSUB,IRDATA,MIDEN,IPUN,NORM,IPLOT IF(N.EQ.0) GOTO99 READ101,TIT PRINT113,TIT PRINT114,N,K,NSUB,IRDATA,MIDEN,IPUN,NORM,IPLOT IF(MIDEN.EQ.0)GOTO24 READ101,FMT READFMT,(BNUM(I),I=1,N) GOTO25 24 DO26I=1,N 26 CALL PIBCI(I,BNUM(I),4H(I2)) C READ IN DATA 25 READ101,FMT N2=(N*(N-1))/2 K2=(K*(K&1))/2 XN=N XN2=N2 XNS=NSUB DO32I=1,N2 32 A(I)=0. DO5IS=1,NSUB DO41I=2,N II=I-1 41 READFMT,(B(I,J),J=1,II) DO12I=2,N II=I-1 DO12J=1,II IF(IRDATA.NE.1)GOTO11 B(I,J)=-B(I,J) 11 B(J,I)=B(I,J) 12 CONTINUE ADDC=0. IF(IRDATA.EQ.3) GOTO35 C COMPUTE ADDITIVE CONSTANT ADDC=-100. M1=N-1 M2=N-2 DO30I=1,M2 JJ=I&1 DO30J=JJ,M1 KK=J&1 DO30L=KK,N DO30M=1,3 GOTO(21,22,23),M 21 CIJK=B(I,L)-B(J,L)-B(I,J) GOTO34 22 CIJK=B(J,L)-B(I,L)-B(I,J) GOTO34 23 CIJK=B(I,J)-B(J,L)-B(I,L) 34 IF(CIJK.GT.ADDC)ADDC=CIJK 30 CONTINUE 35 CALL SDATA(D,N2,IS ,B,NL,A,N,ADDC,T,SP,XN,NORM) 5 CONTINUE IP=0 C COMPUTE INITIAL GROUP STIMULUS SPACE C COMPUTE AVERAGE DSQ MATRIX 51 SUMA=0. DO45L=1,N2 45 A(L)=A(L)/XNS CALL SCALP(A,N,XN,N2,T,SP,NL) C FACTOR SP MX=1 CALL MTEVV(SP,NL,NL,B,NL,NL,N,MX) C STORE INITIAL STIMULUS SPACE IN X AND PRINT DO9I=1,K DO9J=1,N 9 X(J,I)=B(J,I)*SQRT(SP(I,I)) PRINT103 DO33I=1,K 33 PRINT104,I,(X(J,I),J=1,N) C COMPUTE R FOR AVERAGE SUBJECT TO GET TRANSFORMED GROUP SPACE PRINT 102 CALL COMPR(N,K,1,B,X,B,SP,A,IP,DUMY,Q,BNUM, K2,XN2,TIT,RR, 1NL,N2,IPLOT,IPUN,S,NK2) DO36I=1,K DO36J=1,N 36 X(J,I)=SP(J,I) C COMPUTE CORRELATIONS FOR PHASE 3 56 DO37IS=1,NSUB 37 CALL CORR(N,K,X,NL,D,IS,A,IP,COR(IS,3),SQR(IS,3),XN2,N2) IP=1 PRINT107 PRINT106,IP C ENTER PHASE 1 CALL COMPR(N,K,NSUB,B,X,A,SP,D,IP,COR(1,1),Q,BNUM, K2,XN2, 1TIT,SQR(1,1),NL,N2,IPLOT,IPUN,S,NK2) C ENTER PHASE 2 IP=2 PRINT107 PRINT106,IP CALL COMPX(X,Q,B,N,K,SP,NL) C FACTOR Q, Y*=UY, NORMALIZE Y* AND PUT IN X C PLOT 38 CALL DSUBJ(X,N,K,NL,TIT,BNUM,GSUB,5) IF(IPUN.EQ.0) GOTO58 DO39I=1,K 39 PUNCH105,I,(X(J,I),J=1,N) 105 FORMAT(I3,5E13.5/(3X,5E13.5)) 58 CALL PHAS2(N,K,NSUB,X,SP,B,A,COR(1,2),T,D,XN2,SQR(1,2),NL,N2) C COMPUTE F TESTS PRINT113,TIT CALL FTEST(COR,N,K,K2,NSUB,SQR,NL,B) GOTO98 100 FORMAT(20I4) 101 FORMAT(12A6) 102 FORMAT(32H0OPTIMIZATION ON AVERAGE SUBJECT) 103 FORMAT(23H0INITIAL STIMULUS SPACE) 104 FORMAT(I3,3X,10E12.4/(6X,10E12.4)) 106 FORMAT(6H0PHASEI2) 107 FORMAT(55H0******************************************************) 113 FORMAT(1H112A6) 114 FORMAT(39H0 N K NSUB IRDATA MIDEN IPUN NORM IPLOT/1H 2I2,I4,3I6,2I 15) 1,6I6) 115 FORMAT(3X,10E12.4) 99 STOP END $ FORTRAN $ LIMITS 05,40K,,5000 SUBROUTINE COMPR(N,K,NSUB,B,X,A,SP,D,IP,R,Q,BNUM, K2,XN2, 1TIT,R1,NL,N2,IPLOT,IPUN,S,NK) C USE Q WITH WEIGHTS AND FIRST ROW OF DELTA HAS ALL 1,S DIMENSION S(NK,NK),B(NL,NL),X(NL,6),A(1),SP(NL,NL),D(N2,NSUB) DIMENSION R1(1),TSUB(4) DIMENSION Q(6,6) ,R(1),BNUM(1),FSUB(2),FB(3),T(21),TIT(12),Z(6,6) DATA TSUB(1)/24HAVERAGE STIMULUS SPACE / DATA FB(1)/13H(7HSUBJECTI3)/ C COMPUTE FOR EACH SUBJECT R=S(-1)T KC=K2&1 DO44I=1,KC DO44J=1,I 44 S(I,J)=0. DO45I=1,K DO45J=1,I 45 Q(I,J)=0. C COMPUTE S(-1) L=1 DO11I1=1,K DO11J1=1,I1 L=L&1 M=1 DO10I2=1,K DO10J2=1,I2 M=M&1 DO13M3=2,N MM=M3-1 DO13K3=1,MM 13 S(L,M)=S(L,M)&(X(M3,I1)-X(K3,I1))*(X(M3,J1)-X(K3,J1))*(X(M3,I2)-X( 1K3,I2))*(X(M3,J2)-X(K3,J2)) S(M,L)=S(L,M) IF(M.EQ.L) GOTO11 10 CONTINUE 11 CONTINUE S(1,1)=N2 M=1 DO25I=1,K DO25J=1,I M=M&1 DO23I1=2,N II=I1-1 DO23I2=1,II 23 S(M,1)=S(M,1)&(X(I1,J)-X(I2,J))*(X(I1,I)-X(I2,I)) 25 S(1,M)=S(M,1) MM=1 CALL MTINV(S,NK,NK,B,NK,NK,KC,MM) IF(MM.EQ.1)GOTO15 PRINT102 STOP C COMPUTE T 15 DO50IS=1,NSUB IF(IP.EQ.1) PRINT109,IS 66 DO68I=1,KC 68 T(I)=0. L=1 DO14I1=1,K DO14J1=1,I1 L=L&1 M=0 DO19M3=2,N MM=M3-1 DO19K3=1,MM M=M&1 19 T(L)=T(L)&(X(M3,I1)-X(K3,I1))*(X(M3,J1)-X(K3,J1))*D(M,IS) 14 CONTINUE DO24M=1,N2 24 T(1)=T(1)&D(M,IS) C COMPUTE R=S(-1)T DO16I=1,KC A(I)=0. DO16J=1,KC 16 A(I)=A(I)&S(I,J)*T(J) C UNPACK A, STORE IN SP SUMR=0. L=1 DO33I=1,K DO33J=1,I L=L&1 IF(I.EQ.J) GOTO20 SP(I,J)=A(L)*.5 SP(J,I)=SP(I,J) SUMR=SUMR&SP(J,I)**2 GOTO33 20 SP(I,J)=A(L) 33 CONTINUE SUMR=2.*SUMR PRINT107 DO63I=1,K SUMR=SUMR&SP(I,I)**2 63 PRINT104,I,(SP(I,J),J=1,K) C TEST IF PHASE I, COMPUTE Q IF(IP.NE.1) GOTO61 DO3I=1,K DO3J=1,I Z(I,J)=0. DO3M=1,K 3 Z(I,J)=Z(I,J)&SP(I,M)*SP(J,M) 61 MX=1 CALL MTEVV(SP,NL,NL,B,NL,NL,K,MX) PRINT110 PRINT115,(SP(I,I),I=1,K) PRINT118 DO67I=1,K 67 PRINT104,I,(B(J,I),J=1,K) DO17I=1,K DO17J=1,K IF(SP(I,I).LT.0.) SP(I,I)=0. 17 B(J,I)=SQRT(SP(I,I))*B(J,I) PRINT101 DO12I=1,K 12 PRINT111,I,(B(J,I),J=1,K) C COMPUTE TRANSFORMED INDIVIDUAL SPACE, STORE IN SP DO18M=1,K DO18I=1,N SP(I,M)=0. DO18J=1,K 18 SP(I,M)=SP(I,M)&B(J,M)*X(I,J) IF(IP.EQ.0) GOTO4 PRINT108 GOTO5 4 PRINT105 5 DO65I=1,K IF(IPUN.EQ.0)GOTO64 PUNCH111,I,(SP(J,I),J=1,N) 64 PRINT104,I,(SP(J,I),J=1,N) 65 CONTINUE IF(IP.GT.0) GOTO21 CALL DSUBJ(SP,N,K,NL,TIT,BNUM,TSUB,4) GOTO50 21 IF(IPLOT.EQ.0) GOTO49 CALL FORCNV(IS,1,FSUB,FB,DUM1,DUM2) CALL DSUBJ(SP,N,K,NL,TIT,BNUM,FSUB,2) 49 CALL CORR(N,K,SP,NL,D,IS,WT,IP,R(IS),R1(IS),XN2,N2) IF(IP.NE.1) GOTO50 P=R(IS)**2/SUMR DO22I=1,K DO22J=1,I 22 Q(I,J)=Q(I,J)&P*Z(I,J) 50 CONTINUE 101 FORMAT(32H0COMPOSITE TRANSFORMATION MATRIX) 102 FORMAT(21H S CANNOT BE INVERTED) 103 FORMAT(32H0COMPOSITE TRANSFORMATION MATRIX) 104 FORMAT(I3,3X,10E12.4/(6X,10E12.4)) 105 FORMAT(23H0AVERAGE STIMULUS SPACE) 107 FORMAT(2H0R) 108 FORMAT(29H0TRANSFORMED INDIVIDUAL SPACE) 109 FORMAT(8H0SUBJECTI3) 110 FORMAT(18H0EIGEN VALUES OF R) 111 FORMAT(I3,5E13.5/(3X,5E13.5)) 115 FORMAT(3X,10E12.4) 118 FORMAT(19H0EIGEN VECTORS OF R) 99 RETURN END $ FORTRAN SUBROUTINE COMPX(X,Q,B,N,K,SP,NL) DIMENSION X(NL,6),Q(6,6),B(NL,NL),SP(NL,NL),P(6) DO6I=2,K II=I-1 DO6J=1,II 6 Q(J,I)=Q(I,J) PRINT100 DO10I=1,K 10 PRINT101,I,(Q(I,J),J=1,K) MX=1 CALL MTEVV(Q,6,6,B,NL,NL,K,MX) PRINT102 PRINT101,I,(Q(I,I),I=1,K) PRINT103 DO11J=1,K 11 PRINT101,J,(B(I,J),I=1,K) DO2I=1,K SS=0. DO3J=1,N SUM=0. DO4M=1,K 4 SUM=SUM&B(M,I)*X(J,M) SP(J,I)=SUM 3 SS=SS&SUM**2 2 P(I)=1./SQRT(SS) DO5I=1,K DO5L=1,N 5 X(L,I)=SP(L,I)*P(I) PRINT104 DO12I=1,K 12 PRINT101,I,(X(J,I),J=1,N) 100 FORMAT(2H0Q) 101 FORMAT(I3,2X,10E12.4/(5X,10E12.4)) 102 FORMAT(18H0EIGEN VALUES OF Q) 103 FORMAT(19H0EIGEN VECTORS OF Q) 104 FORMAT(29H0INDSCAL GROUP STIMULUS SPACE) RETURN END $ FORTRAN SUBROUTINE PHAS2(N,K,NSUB,X,SP,B,A,R,W,D,XN2,R1,NL,N2) DIMENSION X(NL,6),SP(NL,NL),B(NL,NL),A(1),R(1),R1(1),W(1),D(N2,NSU 1B) C COMPUTE H STORE IN SP, (K BY K) C WITH CONSTANT VECTOR IN DELTA MATRIX KC=K&1 DO2I=1,K L=I&1 DO2J=1,I M=J&1 SUM=0. DO3I1=2,N II=I1-1 DO3J1=1,II 3 SUM=SUM&(X(I1,I)-X(J1,I))**2*(X(I1,J)-X(J1,J))**2 SP(L,M)=SUM 2 SP(M,L)=SUM SP(1,1)=XN2 DO12I=1,K L=I&1 SP(L,1)=0. DO7I1=2,N II=I1-1 DO7I2=1,II 7 SP(L,1)=SP(L,1)&(X(I1,I)-X(I2,I))**2 12 SP(1,L)=SP(L,1) MM=1 CALL MTINV(SP,NL,NL,B,10,10,KC,MM) IF(MM.EQ.1) GOTO6 PRINT100 100 FORMAT(21H S CANNOT BE INVERTED) STOP 6 PRINT101, (I,I=1,K) C COMPUTE S AND W FOR EACH SUBJ DO10IS=1,NSUB DO4I=1,K M=I&1 A(M)=0. L=0 DO5I1=2,N II=I1-1 DO5J1=1,II L=L&1 5 A(M)=A(M)&(X(I1,I)-X(J1,I))**2*D(L,IS) 4 CONTINUE A(1)=0. DO8I=1,N2 8 A(1)=A(1)&D(I,IS) C COMPUTE W DO9I=1,KC W(I)=0. DO9J=1,KC 9 W(I)=W(I)&A(J)*SP(J,I) PRINT102,IS,(W(I),I=2,KC) CALL CORR(N,K,X,NL,D,IS,W(2),2,R(IS),R1(IS),XN2,N2) 10 CONTINUE 101 FORMAT(23H0SUBJECTS WEIGHT MATRIX/19X,9HDIMENSION/5H0SUBJ,I8,6I15) 102 FORMAT(I4,7E15.5) 103 FORMAT(I2,5E13.5/(2X,5E13.5)) 50 RETURN END $ FORTRAN SUBROUTINE SDATA(D,N2,IS ,B,NL,A,N,ADDC,T,SP,XN,NORM) DIMENSION D(N2,IS), B(NL,NL),A(1),T(1),SP(NL,NL) 35 L=0 SUM=0. DO40I=2,N II=I-1 DO40J=1,II L=L&1 D(L,IS)=(B(I,J)&ADDC)**2 IF(NORM-1)15,16,40 15 A(L)=A(L)&D(L,IS) GOTO40 16 SUM=SUM&D(L,IS)**2 40 CONTINUE IF(NORM-1)50,11,12 11 SUM=1./SQRT(SUM) GOTO13 12 CALL SCALP(D(1,IS),N,XN,N2,T,SP,NL) SUM=0. SUMD=SP(1,1)**2 DO14I=2,N SUMD=SUMD&SP(I,I)**2 II=I-1 DO14J=1,II 14 SUM=SUM&SP(I,J)**2 SUM=1./SQRT(2.*SUM&SUMD) 13 DO43L=1,N2 D(L,IS)=D(L,IS)*SUM 43 A(L)=A(L)&D(L,IS) 50 RETURN END $ FORTRAN SUBROUTINE FTEST( R,N,K,K2,NSUB,SQR,NL,F) DIMENSION NF1(6),NF2(6) DIMENSION RMSR(3),RMSS(3) DIMENSION R(NL,3),DF1(6),DF2(6),DF(6),SQR(NL,3),F(NL,6) XNS=NSUB M=(N*(N-1))/2 XK=K XK2=K2 DF2(1)=M-K2-1 DF2(2)=M-K-1 DF2(3)=M-2 DF2(4)=DF2(1) DF2(5)=DF2(1) DF2(6)=DF2(2) DF1(1)=XK2 DF1(2)=XK DF1(3)=1. DF1(4)=XK2-XK DF1(5)=XK2-1. DF1(6)=XK-1. DO7I=1,6 NF1(I)=DF1(I) 7 NF2(I)=DF2(I) PRINT100 PRINT105 DO10IS=1,NSUB 10 PRINT106,IS,( R(IS,J),J=1,3),(SQR(IS,J),J=1,3) DO2I=1,3 RMSR(I)=0. RMSS(I)=0. DO2IS=1,NSUB RMSR(I)=RMSR(I)&R(IS,I)**2 2 RMSS(I)=RMSS(I)&SQR(IS,I)**2 DO8I=1,3 RMSR(I)=SQRT(RMSR(I)/XNS) 8 RMSS(I)=SQRT(RMSS(I)/XNS) PRINT107 PRINT108,(RMSR(I),I=1,3),(RMSS(I),I=1,3) DO5I=1,6 5 DF(I)=DF2(I)/DF1(I) DO4I=1,NSUB DO3J=1,3 R(I,J)= R(I,J)**2 3 F(I,J)=R(I,J)/(1.-R(I,J))*DF(J) L=3 DO6J1=1,2 II=J1&1 DO6J2=II,3 L=L&1 6 F(I,L)=(R(I,J1)-R(I,J2))/(1.-R(I,J1))*DF(L) 4 CONTINUE PRINT104 PRINT101 PRINT102,(NF1(I),NF2(I),I=1,6) DO11IS=1,NSUB 11 PRINT103,IS,(F(IS,J),J=1,6) 100 FORMAT(82H0TABLE I. CORRELATIONS AND SQUARIANCES BETWEEN DATA(EST 1.)**2 AND OBTAINED DIST**2/26X,12HCORRELATIONS,27X,11HSQUARIANCES) 101 FORMAT(48X,5HPHASE/ 15X,1H1,14X,1H2,14X,1H3,13X,3H1-2,12X 1,3H1-3,12X,3H2-3) 102 FORMAT(4X,2HDF,6X,2I4,5(7X,2I4)/8H SUBJECT) 103 FORMAT(I5,3X,6E15.6) 104 FORMAT(19H0TABLE II. F TESTS) 105 FORMAT(29X,5HPHASE,45X,5HPHASE/8H SUBJECT,7X,1H1,14X,1H2,14X,1H3, 119X,1H1,14X,1H2,14X,1H3) 106 FORMAT(I5,3X,3E15.6,5X,3E15.6) 107 FORMAT(18H0ROOT MEAN SQUARES) 108 FORMAT(8X,3E15.6,5X,3E15.6) RETURN END $ FORTRAN SUBROUTINE SCALP(A,N,XN,N2,T,SP,NL) DIMENSION A(1),T(1),SP(NL,NL) SUMA=0. DO3L=1,N2 3 SUMA=SUMA&A(L) ADIJ=2.*SUMA/XN**2 DO7I=1,N T(I)=0. DO7J=1,N IF(I.EQ.J) GOTO7 M=IFUNF(I,J) T(I)= T(I)&A(M) 7 CONTINUE DO8I=1,N DO8J=1,I QQ=( T(I)& T(J))/XN IF(I.NE.J) GOTO42 SP(I,I)=(QQ-ADIJ)*(.5) GOTO8 42 M=IFUNF(I,J) SP(I,J)=(QQ-A(M)-ADIJ)*(.5) SP(J,I)=SP(I,J) 8 CONTINUE RETURN END $ FORTRAN SUBROUTINE CORR(N,K,X,NL,D,IS,W,IP,R,R1,XN2,N2) DIMENSION X(NL,6),D(N2,IS),W(1) C COMPUTE R BETWEEN DATA AND SOLUTION , USE DSQ. SX=0. SY=0. SQX=0. SQY=0. SXY=0. L=0 DO70I=2,N II=I-1 DO70J=1,II L=L&1 DSQ=0. DO72M=1,K Q=( X(I,M)- X(J,M))**2 IF(IP.EQ.2) GOTO74 DSQ=DSQ&Q GOTO72 74 DSQ=DSQ&Q*W(M) 72 CONTINUE SX=SX&D(L,IS) SY=SY&DSQ SQX=SQX&D(L,IS)**2 SQY=SQY&DSQ**2 SXY=SXY&DSQ*D(L,IS) 70 CONTINUE DEN=(XN2*SQX-SX**2)*(XN2*SQY-SY**2) R =(XN2*SXY-SX*SY)/SQRT(DEN) R1=SXY/SQRT(SQX*SQY) RETURN END $ FORTRAN SUBROUTINE DSUBJ(X,N,NF,NL,TIT,BNUM,FSUB,NW) C PROGRAM TO PLOT N STIMULUS POINTS X IN NF DIMENSIONAL SPACE DIMENSION X(NL,NF),BNUM(1),TIT(1),FSUB(1) DIMENSION FFM(2),FI(2),FJ(2),FDM(4),FFDM(3) DATA FDM(1)/20H(I2,12H DIMENSIONAL)/ DATA FFM(1)/10H(4HDIM.I3)/ C PLOTTING AREA IS 460 BY 460, ORIGIN AT CENTER OF SQUARE NNF=NF-1 CALL FORCNV(NF,1,FFDM,FDM,DUM1,DUM2) 5 XMAX=0. DO10I=1,NF CALL MAXIS(X(1,I),SCL,N) IF(SCL.GT.XMAX) XMAX=SCL 10 CONTINUE XSCALE=460./XMAX 6 DO260I=1,NNF JJ=I&1 DO260J=JJ,NF 34 CALL ROLL CALL REFRNC(-512,512,1024,1024) C PRINT TITLE CALL TEXT(-460,512,TIT,0,12) CALL TEXT(-460,460,FFDM,0,3) C LABEL NF DIMENSIONAL SOLUTION C PRINT SUBTITLE CALL TEXT (-460,500,FSUB,0,NW) C DRAW POINTS DO 22 K=1,N IAX=X(K,I)*XSCALE IF(NF.EQ.1)GOTO23 IAY=X(K,J)*XSCALE GOTO24 23 IAY=0 24 CALL TEXT (IAX,IAY,BNUM(K),0,1) 22 CONTINUE IF(NF.EQ.1)GOTO50 C DRAW AXES CALL DVR1(-460,0,460,0,1) CALL DVR1(0,-460,0,460,1) C LABEL AXES CALL FORCNV(J,1,FJ,FFM,DUM1,DUM2) CALL TEXT(-10,467,FJ,0,2) CALL FORCNV(I,1,FI,FFM,DUM1,DUM2) CALL TEXT(440,7,FI,0,2) 260 CONTINUE 50 RETURN END $ FORTRAN SUBROUTINE MAXI(X,Y,Z,N) C MAXI FINDS THE MAGNITUDE OF ARRAYS X AND Y EACH OF LENGTH N. C ON RETURN ABSOLUTE MAXIMUM IS STORED IN Z. DIMENSION X(1),Y(1) Z=ABS(X(1)) DO10I=1,N AXO=ABS(X(I)) IF(AXO.GT.Z) Z=AXO AYO=ABS(Y(I)) IF(AYO.GT.Z) Z=AYO 10 CONTINUE RETURN ENTRY MAXIS(X,Z,N) Z=ABS(X(1)) DO11I=2,N AXO=ABS(X(I)) IF(AXO.GT.Z) Z=AXO 11 CONTINUE RETURN END $ FORTRAN FUNCTION IFUNF(I,J) IF(I-J)1,2,3 1 IFUNF=((J-1)*(J-2))/2&I GOTO5 2 PRINT100 100 FORMAT(20H0I=J, INDEXING ERROR) STOP 3 IFUNF=((I-1)*(I-2))/2&J 5 RETURN END $ FORTRAN $ INCODE IBMF CMTINV MTINV SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M) MTPK1115 DIMENSION A(I1,I2), B(I3,I4) MTPK1116 C MTPK1117 C THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS- MTPK1118 C JORDAN ELIMINATION SCHEME. ALL CALCULATIONS ARE DONE IN DOUBLE- MTPK1119 C PRECISION ARITHMETIC MTPK1120 C MTPK1121 DOUBLE PRECISION B,BMAX,BMULT MTPK1122 EQUIVALENCE (BMAX,BMULT) MTPK1123 C MTPK1124 C SINGLE TO DOUBLE PRECISION TRANSFER MTPK1125 C MTPK1126 101 DO 103 I = 1,N MTPK1127 102 DO 103 J = 1,N MTPK1128 103 B(I,J) = A(I,J) MTPK1129 C MTPK1130 C FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMATRIX. MTPK1133 C MTPK1134 200 DO 430 K = 1,N MTPK1131 A(K,1) = FLOAT(K) MTPK1135 A(K,2) = A(K,1) 203 BMAX = DABS(B(K,K)) MTPK1137 210 DO 219 I = K,N MTPK1138 211 DO 219 J = K,N MTPK1139 IF (BMAX - DABS(B(I,J)) ) 213,219,219 213 BMAX = DABS ( B(I,J)) MTPK1141 214 A(K,1) = FLOAT(I) MTPK1142 215 A(K,2) = FLOAT(J) MTPK1143 219 CONTINUE MTPK1145 216 IF (BMAX - 1.D-38) 801,301,301 C MTPK1146 C EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL. MTPK1147 C MTPK1148 301 I = IFIX(A(K,1)) MTPK1149 302 J = IFIX(A(K,2)) MTPK1150 313 DO 314 M = 1,N MTPK1151 BMULT = B(I,M) MTPK1152 B(I,M) = B(K,M) MTPK1153 314 B(K,M) = BMULT MTPK1155 303 DO 304 M = 1,N MTPK1156 BMULT = B(M,J) MTPK1157 B(M,J) = B(M,K) MTPK1158 304 B(M,K) = BMULT MTPK1159 C MTPK1160 C ACTUAL INVERSION STEP. MTPK1161 C BEGIN ROW ITERATION MTPK1162 C MTPK1163 401 DO 425 I = 1,N MTPK1164 C IGNORE ROW K. MTPK1165 402 IF (I - K) 405,425,405 C SET ROW MULTIPLIER. MTPK1167 405 BMULT = B(I,K) / B(K,K) C MODIFY ROW ELEMENTS. MTPK1169 410 DO 420 J = 1,N MTPK1170 C IGNORE COLUMN K MTPK1171 411 IF (J - K) 414,412,414 412 B(I,J) = -BMULT MTPK1173 413 GO TO 420 MTPK1174 414 B(I,J) = B(I,J) - B(K,J) * BMULT MTPK1175 420 CONTINUE MTPK1176 425 CONTINUE MTPK1177 C MTPK1178 C DIVIDE PIVOT ROW BY PIVOT ELEMENT MTPK1179 C MTPK1180 426 BMULT = 1.D0/ B(K,K) 427 B(K,K) = 1.D0 C ACTUAL DIVISION MTPK1183 428 DO 429 J = 1,N MTPK1184 429 B(K,J) = B(K,J) * BMULT MTPK1185 430 CONTINUE MTPK1186 C MTPK1187 C REARRANGE AND REASSEMBLE INVERSE MATRIX. MTPK1188 C MTPK1189 501 DO 512 IJK = 1,N MTPK1190 502 I = N-IJK&1 MTPK1191 503 L = IFIX(A(I,2)) MTPK1192 504 J = IFIX(A(I,1)) MTPK1193 505 DO 508 M = 1,N MTPK1194 506 BMULT = B(I,M) MTPK1195 507 B(I,M) = B(L,M) MTPK1196 508 B(L,M) = BMULT MTPK1197 509 DO 512 M = 1,N MTPK1198 510 BMULT = B(M,I) MTPK1199 511 B(M,I) = B(M,J) MTPK1200 512 B(M,J) = BMULT MTPK1201 C MTPK1202 C MOVE A INVERSE BACK INTO A MTPK1203 C MTPK1204 601 DO 603 I = 1,N MTPK1205 602 DO 603 J = 1,N MTPK1206 603 A(I,J) = B(I,J) 701 M = 1 MTPK1208 702 RETURN MTPK1209 C MTPK1210 C ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38 C MTPK1212 801 M = 2 MTPK1213 802 RETURN MTPK1214 999 END MTPK1215 $ FORTRAN $ INCODE IBMF *EIGENJ EIGENVALUE AND EIGENVECTOR SUBROUTINE EIGJ0020 * CD600D4.009 DATE 05/05/65 EIGJ0030 SUBROUTINE EIGENJ (A,V,NN,IV,E,I3,I4) DIMENSION A(5050),V( I3, I4),E(NN) LOGICAL IV EIGJ0060 * EIGENJ FINDS THE EIGENVALUES AND EIGENVECTORS OF EIGJ0070 * A SYMMETRIC MATRIX (A) USING A MODIFIED THRESHOLD JACOBI METHOD EIGJ0080 * ONLY THE UPPER TRIANGLE IS STORED EIGJ0090 * DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER EIGJ0100 * DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER EIGJ0110 N=NN EIGJ0120 IND=1 EIGJ0130 NM2=N-2 EIGJ0140 AMAX=0.0 EIGJ0150 NM=N-1 EIGJ0160 IF(.NOT.IV)GO TO 5 EIGJ0170 * SET UP VECTOR MATRIX EIGJ0180 1 DO 3 I=1,N EIGJ0190 DO 2 J=1,N EIGJ0200 2 V(I,J)=0.0 EIGJ0210 3 V(I,I)=1.0 EIGJ0220 5 K=2 EIGJ0230 DO 28 I=1,NM EIGJ0240 IP=I&1 EIGJ0250 DO 27 J=IP,N EIGJ0260 Y=A(K) EIGJ0270 X=ABS(Y) EIGJ0280 IF(X-AMAX)27,27,6 EIGJ0290 6 GO TO (7,8),IND EIGJ0300 7 AMAX=X EIGJ0310 GO TO 27 EIGJ0320 * TRANSFORMATION SETUP FOLLOWS EIGJ0330 8 II=K&I-J EIGJ0340 JJ=(J*(N&N-J&3))/2-N EIGJ0350 ITEST=1 EIGJ0360 X=A(II)-A(JJ) EIGJ0370 IF(X)10,11,10 EIGJ0380 10 X=Y/X EIGJ0390 Y=.5*X EIGJ0400 IF(ABS(Y)-.414213562)14,11,11 EIGJ0410 11 C=.707106781 EIGJ0420 IF(Y)12,12,13 EIGJ0430 12 S=-C EIGJ0440 GO TO 15 EIGJ0450 13 S=C EIGJ0460 GO TO 15 EIGJ0470 14 X=Y*Y EIGJ0480 C=1.&X EIGJ0490 S=2.*Y/C EIGJ0500 C=(1.-X)/C EIGJ0510 15 X=S*S EIGJ0520 Y=C*C EIGJ0530 XY=S*C EIGJ0540 AXY=2.*A(K)*XY EIGJ0550 Z=A(II)*Y&AXY&A(JJ)*X EIGJ0560 W=A(II)*X-AXY&A(JJ)*Y EIGJ0570 A(K)=A(K)*(Y-X)&XY*(A(JJ)-A(II)) EIGJ0580 A(II)=Z EIGJ0590 A(JJ)=W EIGJ0600 IF(NM2)24,24,16 EIGJ0610 16 IT=I EIGJ0620 JT=J EIGJ0630 * TRANSFORM A EIGJ0640 DO 23 M=1,NM2 EIGJ0650 NP=N-M EIGJ0660 IF(JT-K)20,17,18 EIGJ0670 17 IT=IT&1 EIGJ0680 JT=JT&NP EIGJ0690 ITEST=2 EIGJ0700 18 IF(IT-K)20,19,19 EIGJ0710 19 IT=IT&1 EIGJ0720 JT=JT&1 EIGJ0730 ITEST=3 EIGJ0740 20 X=A(IT)*C&A(JT)*S EIGJ0750 A(JT)=A(JT)*C-A(IT)*S EIGJ0760 A(IT)=X EIGJ0770 GO TO (21,22,23),ITEST EIGJ0780 21 IT=IT&NP EIGJ0790 JT=JT&NP EIGJ0800 GO TO 23 EIGJ0810 22 IT=IT&1 EIGJ0820 JT=JT&NP-1 EIGJ0830 23 CONTINUE EIGJ0840 24 IF(.NOT.IV)GO TO 27 EIGJ0850 * TRANSFORM V EIGJ0860 25 DO 26 M=1,N EIGJ0870 X=V(M,I)*C&V(M,J)*S EIGJ0880 V(M,J)=V(M,J)*C-V(M,I)*S EIGJ0890 26 V(M,I)=X EIGJ0900 27 K=K&1 EIGJ0910 28 K=K&1 EIGJ0920 * SEQUENCE TO ADJUST THRESHOLD FOR EACH SWEEP EIGJ0930 GO TO (29,31),IND EIGJ0940 29 IF(AMAX)30,35,30 EIGJ0950 30 AT=AMAX EIGJ0960 IND=2 EIGJ0970 F=.25 EIGJ0980 AMAX=F*AT EIGJ0990 MAX=6 EIGJ1000 GO TO 5 EIGJ1010 31 MAX=MAX-1 EIGJ1020 IF(MAX)34,33,32 EIGJ1030 32 F=2.0*F*F EIGJ1040 AMAX=AT*F EIGJ1050 33 C=0.0 EIGJ1060 GO TO 5 EIGJ1070 34 IF(C)33,35,33 EIGJ1080 * COPY EIGENVALUES INTO E EIGJ1090 35 MM=1 EIGJ1100 NM=N&1 EIGJ1110 DO 36 M=1,N EIGJ1120 E(M)=A(MM) EIGJ1130 36 MM=MM&NM-M EIGJ1140 RETURN EIGJ1150 END $ FORTRAN $ INCODE IBMF SUBROUTINE MTEVV(X,I1,I2,D,I3,I4,N,M) DIMENSION A(5050),X(I1,I2),D(I3,I4),E(100),NS(100),MS(100) K=1 DO1J=1,N DO1I=J,N A(K)=X(I,J) 1 K=K&1 CALL EIGENJ(A,D,N,M,E,I3,I4) C THE FOLLOWING SYSTEM SORT PACKAGE ROUTINES SORT THE FLOATING POINT C ARRAY E OF N ELEMENTS IN DESCENDING ORDER. CALL SORTFL(-1) IF(M)4,3,4 3 CALL SORT(E(1),E(2),N) DO 50 I=1,N 50 X(I,I)=E(I) RETURN 4 DO 5 I=1,N MS(I)=I 5 NS(I)=I C THE FOLLOWING ROUTINE SORTS THE FLOATING POINT ARRAY E OF N C ELEMENTS IN DESCENDING ORDER AND ARRAY MS IS THE PASSIVE LIST C WHICH WILL BE REARRANGED IN THE SAME ORDER AS E. CALL SORT(E(1),E(2),N,MS(1)) C THE FOLLOWING ROUTINES SORT THE FIXED POINT ARRAY MS OF N ELEMENTS C IN ASCENDING ORDER WITH ARRAY NS AS THE PASSIVE LIST. CALL SORTFX(1) CALL SORT(MS(1),MS(2),N,NS(1)) CALL ARREV(N,NS,D,I3,I4) DO 6 I=1,N 6 X(I,I)=E(I) RETURN END $ FORTRAN $ INCODE IBMF SUBROUTINE ARREV(N,NS,D,I3,I4) DIMENSION NS(N),D(I3,I4) DO 10 I=1,N 1 IF(NS(I)-I)2,10,2 2 CALL XCH(N,I,NS(I),D,NS,I3,I4) GOTO1 10 CONTINUE RETURN END $ FORTRAN $ INCODE IBMF SUBROUTINE XCH(N,L,M,D,NS,I3,I4) DIMENSION NS(N),D(I3,I4) I=L J=M DO 1 K=1,N T=D(K,I) D(K,I)=D(K,J) 1 D(K,J)=T NS(I)=NS(J) NS(J)=J RETURN END