Trailing-Edge
-
PDP-10 Archives
-
-
There are no other files named in the archive.
C
C***********************************************************************
C
C PROJECTED REVISIONS IN THE FORTRAN IV (7094) VERSION OF BMD72X
C TO ALLOW ITS EXECUTION ON THE IBM SYSTEM/360. AUGUST 9, 1966.
C
C***********************************************************************
C
C FACTOR ANALYSIS - MAIN PROGRAM MARCH 28, 1966
C
DIMENSION DATE(5)
DOUBLE PRECISION FINISH
DOUBLE PRECISION P,PC,PROBLM
DATA DATE/'APRI','L 25',', 19','69 ',' '/
DOUBLE PRECISION CUNST,GAMO
DATA MAXST/6700/
DIMENSION A(6700)
DATA BINARY,YES,PROBLM,FINISH,ONO/4HBNRY,3HYES,6HPROBLM,6HFINISH,
.2HNO/
C EXTERNAL SIGN STATEMENT WAS REMOVED FROM HERE
LOGICAL FSC
COMMON PC,NC,NV,MAT,ICOM,IROT,GAMA,NROT,CONST,COR,VEC,PFS,WFS,MT,N
XF,RE,MFT,NFO,REO,AABBCC,FSC,MXC,MXR,AMIS,NVV
C
C
AMIS=ONO
C THIS VARIABLE (AMIS) IS A DUMMY VARIABLE AT THIS TIME. IT WAS
C INTENDED THAT AT SOME FUTURE TIME IT WOULD BE USED IN
C ALLOWING FOR MISSING VARIABLES. THIS HAS NOT YET BEEN DONE.
C
C
CALL USAGEB('BMDX72')
55 REWIND 1
REWIND 2
REWIND 3
81 FORMAT(37H1BMDX72 - FACTOR ANALYSIS - REVISED ,5A4
X/41H0HEALTH SCIENCES COMPUTING FACILITY, UCLA)
READ(5,2) P,PC,NC,NV,MAT,ICOM,MXC,IROT,MXR,GAMO,NROT,CUNST,AKN,
1COR,VEC,PFS,WFS,MT,NF,RE,MFT,NFO,REO,UPL
*,NLF,KL,KRL
IF(P.EQ.FINISH) STOP
IF(P.NE.PROBLM) GO TO 121
C PRIOR TO MAY 1,1968 THIS PROGRAM
C PRINTED THE PROBLM CARD AS IT WAS READ IN
C UNDER A-FORMAT, THEN USED 'USEBUF' TO
C REREAD THIS CARD. THIS A-FORMAT INFORMATION
C WAS STORED IN A REAL * 8 ' FF(14) ' ARRAY
C USING ' FORMAT (13A6,A2) ' TYPE READ. THIS
C SEEMED UNNECESSARY AND WAS DROPPED.
2 FORMAT(2A6,I6,I3,2I1,I2,I1,I2,A6,I3,A6,A2,4A3,2(2I2,A2),F3.3,3I2)
KN=1
IF(AKN.EQ.ONO) KN=0
WRITE(6,81) DATE
IF(UPL.EQ.0.0) UPL=.95
IF(ICOM.EQ.0) ICOM=1
MAT=4-MAT
IF(MAT.LE.0) MAT=5-MAT
IF (NROT.EQ.0) NROT=NV*0.5
IF(MT.EQ.0) MT=5
IF(NFO.EQ.0 .AND. WFS.EQ.YES) NFO=1
IF(NF.EQ.0) NF=1
CONST=1.0
CALL ATOF (CUNST,6,CONST)
GAMA=100.00
CALL ATOF (GAMO,6,GAMA)
IF(MXC.EQ.0) MXC=1
IF(MXR.EQ.0) MXR=50
NF=18*NF
NFO=18*NFO
IF(NF.LT.0) NF=-1
IF(NFO.LT.0) NFO=-1
L00=1+IABS(NF)
L0=L00+IABS(NFO)
L8=L0+NV
L9=L8+NV
M3=L9+NV
L1=M3+NV
L2=L1+NV
L3=L2+NV
L4=L3+NV
L5=L4+NV
L6=L5+NV
L7=MAX0(L6+NV,L4+255)
LL=L7+NV*NROT
NNR=0
IF(IROT.GT.1) NNR=NROT*NROT
NVV=(NV*(NV+1))/2
NLF=NLF*20
K8=L7+NVV
K9=K8+NLF
K91=K9-1
IF(K9.GT.MAXST) GO TO 791
LA=1
A(1)=BINARY
A(L00)=BINARY
FSC=.FALSE.
IF(PFS.EQ.YES .OR. WFS.EQ.YES) FSC=.TRUE.
CALL DOIT1(A,A(L00),A(L0),A(L8),A(L9),A(L7),A(L4),UPL)
IF(NLF.NE.0) READ (5,20) (A(I),I=K8,K91)
20 FORMAT(20A4)
17 CALL DOIT2(A,A(L4),A(L00),A(M3),A(L0),A(L8),A(L9),A(L7),A(L7),NV,
XLA,LAL)
GO TO (11,11,12,13,57),LA
11 CALL MULTR(A(L7),A(M3),A(L0),A(L1),NV,A(L4))
GO TO 17
12 A(L0)=AABBCC
CALL RAV(A(L7),A(M3),NV,A(L0),A(L1),A(L2),A(L3),A(L4),A(L5),A(L6),
X.00001,CONST,NROT,MXC)
LL=L7+NV*NROT
NNR=0
IF(IROT.GT.1) NNR=NROT*NROT
IF(LL.GT.MAXST) GO TO 371
GO TO 17
13 IF(KL.NE.0) CALL LOUT(A(L7),NV,NROT,KL,A(K8))
IF(IROT.EQ.0) GO TO 57
IF(LL+NNR.GT.MAXST) GO TO 371
GO TO (8,9,10),IROT
8 CALL ROTAT1(A(L7),A(LL),NV,NROT,GAMA,.00001,MXR,1,KN,A(L0),A(L1)
X,A(L2),A(L4),A(L3),A(L5),A(L6))
IF(KRL.NE.0) CALL LOUT(A(L7),NV,NROT,KRL,A(K8))
GO TO 17
9 CALL ROTAT2(A(L7),A(LL),NV,NROT,GAMA,.00001,MXR,1,KN,A(L0),A(L1)
X,A(L2),A(L4),A(L3),A(L5),A(L6),UPL)
IF(KRL.NE.0) CALL LOUT(A(L7),NV,NROT,KRL,A(K8))
GO TO 17
10 CALL ROTAT3(A(L7),A(LL),NV,NROT,GAMA,.00001,MXR,1,KN,A(L0),A(L1)
X,A(L2),A(L4),A(L3),A(L5),A(L6))
IF(KRL.NE.0) CALL LOUT(A(L7),NV,NROT,KRL,A(K8))
GO TO 17
121 WRITE (6,122)
122 FORMAT(28H0UNRECOGNIZABLE PROBLEM CARD)
STOP
371 WRITE (6,372)
372 FORMAT(47H0INSUFFICIENT STORAGE IS AVAILABLE FOR ROTATION)
STOP
791 WRITE (6,792)
792 FORMAT(36H0THIS PROBLEM HAS TOO MANY VARIABLES)
STOP
57 IF(.NOT.FSC) GO TO 55
502 FORMAT(1H0,'HENCE FACTOR SCORES ARE NOT COMPUTED OR PRINTED, ALTHO
XUGH REQUESTED')
IF(MAT.EQ.5) WRITE(6,500)
500 FORMAT(1H0,'A CORRELATION OR COVARIANCE MATRIX IS READ AS DATA')
IF(MAT.EQ.6) WRITE(6,501)
501 FORMAT(1H0,'A FACTOR LOADING MATRIX IS READ AS DATA ')
IF(MAT.GE.5) WRITE(6,502)
LL7=LL+NNR
IF(MAT.GE.5) GO TO 55
IF(LL7+NV*NROT.GT.MAXST) GO TO 356
CALL FACSCO(A(L00),A(L7),A(LL),A(LL7),A(L1),A(L2),A(L3),A(L8),A(L9
X),A(L4),NV,NROT)
GO TO 55
356 WRITE (6,558)
558 FORMAT(67H0INSUFFICIENT STORAGE IS AVAILABLE FOR COMPUTATION OF FA
XCTOR SCORES)
GO TO 55
END
SUBROUTINE LOUT(A,NV,NR,KL,F)
DIMENSION A(NV,NR),F(16)
DO 1 I=1,NV
1 WRITE(KL,F) I,(A(I,J),J=1,NR)
RETURN
END
C FUNCTION NBL FOR BMDX72 MARCH 28, 1966
LOGICAL FUNCTION NBL(X)
EXTERNAL SIGN
NBL=.TRUE.
IF (X.EQ.0.0.AND.SIGN(1.0,X).LT.0.0) NBL=.FALSE.
RETURN
END
C SUBROUTINE FACSCO FOR BMDX72 MARCH 28, 1966
SUBROUTINE FACSCO(FO,A,T,H,X,Y,Z,S1,D1,C,NV,NROT)
DATA ST,BL,YES/1H*,1H ,3HYES/
DOUBLE PRECISION PC
DIMENSION A(NV,NROT),T(NROT,NROT),H(NV,NROT),S1(2),D1(2),X(2),Y(2)
X,Z(2),C(255),FO(180)
COMMON PC,NC,DM,MAT,ICOM,IROT,GAMA,DUM2,CONST,COR,VEC,PFS,WFS,MT,N
XF,RE,MFT,NFO,REO,AABBCC,FSC,MXC,MXR,AMIS,NVV
LOGICAL NBL
IF(IROT.LT.2) GO TO 37
DO 1 I=1,NV
DO 2 K=1,NROT
X(K)=A(I,K)
2 A(I,K)=0.
DO 1 J=1,NROT
DO 1 K=1,NROT
1 A(I,J)=A(I,J)+T(J,K)*X(K)
37 M=255
DO 3 I=1,NV
DO 21 K=1,NV
IF(M.LT.255) GO TO 20
M=0
READ (1) C
200 FORMAT (20A4)
20 M=M+1
21 X(K)=-C(M)
DO 3 J=1,NROT
H(I,J)=0.
DO 3 K=1,NV
3 H(I,J)=H(I,J)+X(K)*A(K,J)
REWIND 1
WRITE (6,300)
300 FORMAT(26H1FACTOR SCORE COEFFICIENTS)
L1=0
301 L0=L1+1
L1=MIN0(L1+10,NROT)
WRITE (6,302) (L,L=L0,L1)
302 FORMAT(//20X,6HFACTOR,/2X,10I12)
WRITE (6,303)
303 FORMAT(9H VARIABLE)
DO 304 I=1,NV
304 WRITE (6,305) I,(H(I,L),L=L0,L1)
305 FORMAT(I5,10F12.5)
IF(L1.NE.NROT) GO TO 301
M=255
IF(PFS.NE.YES) GO TO 100
IF(AMIS.EQ.YES) WRITE (6,101)
IF(AMIS.NE.YES) WRITE (6,102)
101 FORMAT(1H1,10X,47HFACTOR SCORES (* COMPUTED FROM INCOMPLETE DATA)/
X5H CASE)
102 FORMAT(1H1,10X,13HFACTOR SCORES/5H CASE)
100 DO 71 L=1,NC
HH=BL
DO 5 I=1,NV
IF(M.LT.255) GO TO 6
M=0
READ (2) C
6 M=M+1
Z(I)=C(M)
IF (MAT.EQ.4) X(I)=(C(M)-S1(I))/D1(I)
IF (MAT.EQ.3) X(I)=C(M)/D1(I)
IF (MAT.EQ.2) X(I)=C(M)-S1(I)
IF (MAT.EQ.1) X(I)=C(M)
4 IF(AMIS.NE.YES .OR. NBL(C(M))) GO TO 5
X(I)=0.
HH=ST
5 CONTINUE
DO 17 I=1,NROT
Y(I)=0.
DO 17 J=1,NV
17 Y(I)=Y(I)+H(J,I)*X(J)
IF(WFS.NE.YES) GO TO 71
IF(NFO.GT.0) WRITE (MFT,FO) (Z(I),I=1,NV),(Y(I),I=1,NROT)
GO TO 1000
1000 IF(NFO.LT.0) WRITE (MFT) (Z(I),I=1,NV),(Y(I),I=1,NROT)
GO TO 71
71 IF(PFS.EQ.YES) WRITE (6,75) L,HH,(Y(J),J=1,NROT)
75 FORMAT(1X,I4,A1,10F12.5/(6X,10F12.5))
RETURN
END
C SUBROUTINE DOIT1 FOR BMDX72 MARCH 28, 1966
SUBROUTINE DOIT1(F,FO,X,S1,D1,A,C,UPL)
DIMENSION C(255)
DIMENSION F(180),FO(180),X(2),S1(2),D1(2),A(2)
DOUBLE PRECISION PC
COMMON PC,NC,NV,MAT,ICOM,IROT,GAMA,NROT,CONST,COR,VEC,PFS,WFS,MT,N
XF,RE,MFT,NFO,REO,AABBCC,FSC,MXC,MXR,AMIS,NVV
LOGICAL NBL
DATA ONO,YES/2HNO,3HYES/
LOGICAL FSC
IF(GAMA.NE.100.00) GO TO 70
GAMA=0.0
IF(IROT.EQ.2) GAMA=0.5
IF(IROT.EQ.1) GAMA=1.0
70 MOT=MOD(MAT,2)
IF(RE.NE.ONO .AND. MT.NE.5) REWIND MT
IF(MFT.EQ.0) REO=ONO
IF(REO.NE.ONO .AND. MFT.NE.6) REWIND MFT
IF(NF.GT.0) READ (5,87) (F(I),I=1,NF)
IF(NFO.GT.0) READ (5,87) (FO(I),I=1,NFO)
87 FORMAT (18A4)
WRITE (6,50) PC,NC,NV,MT,MXC,MXR,NROT,CONST
50 FORMAT(40H0PROBLEM CODE
XA6/40H0NUMBER OF CASES
XI6/40H0NUMBER OF VARIABLES
XI6/40H0INPUT TAPE
XI6,/,'0MAX. ITERATIONS FOR COMMUNALITIES',8X,I3,/,'0MAX. ITERATION
XS FOR ROTATION',13X,I3,/,'0NUMBER OF FACTORS TO BE ROTATED',10X,I4
X,/,'0CONSTANT',31X,F10.6)
WRITE(6,2000) UPL
2000 FORMAT(40H0UPPER LIMIT ON CORRELATION COEFFICIENT ,2X,F6.5)
WRITE (6,51) (F(I),I=1,NF)
51 FORMAT(40H0INPUT FORMAT
X18A4/(40X,18A4))
IF(AMIS.EQ.YES) WRITE (6,500)
500 FORMAT(35H0BLANKS ARE TREATED AS MISSING DATA)
IF(WFS.EQ.YES) WRITE (6,52) MFT,(FO(I),I=1,NFO)
52 FORMAT(40H0OUTPUT TAPE
XI6/40H0OUTPUT FORMAT
X18A4/(40X,18A4))
IF(MAT.EQ.1) WRITE (6,601)
IF(MAT.EQ.2) WRITE (6,602)
IF(MAT.EQ.3) WRITE (6,603)
IF(MAT.EQ.4) WRITE (6,604)
IF(MAT.EQ.5) WRITE (6,605)
IF(MAT.EQ.6) WRITE (6,606)
IF(MAT.GE.6) GO TO 666
601 FORMAT(49H0THE COVARIANCE MATRIX ABOUT THE ORIGIN IS FORMED)
602 FORMAT(32H0THE COVARIANCE MATRIX IS FORMED)
603 FORMAT(50H0THE CORRELATION MATRIX ABOUT THE ORIGIN IS FORMED)
604 FORMAT(33H0THE CORRELATION MATRIX IS FORMED)
605 FORMAT(37H0THE MATRIX TO BE FACTORED IS READ IN)
606 FORMAT(27H0A FACTOR MATRIX IS READ IN)
IF(ICOM.EQ.1) WRITE (6,701)
IF(ICOM.EQ.2) WRITE (6,702)
IF(ICOM.EQ.3) WRITE (6,703)
IF(ICOM.EQ.4) WRITE (6,704)
701 FORMAT(32H0DIAGONAL ELEMENTS ARE UNALTERED)
702 FORMAT(64H0INITIAL COMMUNALITY ESTIMATES ARE SQUARED MULTIPLE CORR
XELATIONS)
703 FORMAT(63H0INITIAL COMMUNALITY ESTIMATES ARE LARGEST ABSOLUTE ROW
X VALUES)
704 FORMAT(42H0INITIAL COMMUNALITY ESTIMATES ARE READ IN)
666 IF(IROT.EQ.0) WRITE (6,800)
IF(IROT.EQ.1 .AND. GAMA.EQ.0.0) WRITE (6,811)
IF(IROT.EQ.1 .AND. GAMA.EQ.1.0) WRITE (6,821)
IF(IROT.EQ.1 .AND.GAMA.NE.1.0 .AND. GAMA.NE.0.0) WRITE (6,831)
XGAMA
IF(IROT.EQ.2) WRITE (6,802) GAMA
IF(IROT.EQ.3) WRITE (6,803) GAMA
800 FORMAT(25H0NO ROTATION IS PERFORMED)
811 FORMAT(32H0QUARTIMAX ROTATION IS PERFORMED)
821 FORMAT(30H0VARIMAX ROTATION IS PERFORMED)
831 FORMAT(53H0ORTHOGONAL ROTATION IS PERFORMED WITH GAMMA EQUAL TOF7.
X4)
801 FORMAT(33H0ORTHOGONAL ROTATION IS PERFORMED)
802 FORMAT(50H0OBLIMIN ROTATION IS PERFORMED WITH GAMMA EQUAL TOF7.4)
803 FORMAT(70H0OBLIQUE ROTATION FOR SIMPLE LOADINGS IS PERFORMED WITH
XGAMMA EQUAL TOF7.4)
IF(MAT.NE.5) GO TO 20
L1=0
DO 21 J=1,NV
L0=L1+1
L1=L0+NV-J
IF(NF.GT.0) READ (MT,F) (X(I),I=1,NV)
1000 IF(NF.LE.0) READ (MT) (X(I),I=1,NV)
1001 I=J
DO 21 L=L0,L1
A(L)=X(I)
21 I=I+1
RETURN
20 IF(MAT.EQ.6) RETURN
22 L=0
DO 1 I=1,NV
D1(I)=0.
S1(I)=0.
DO 1 J=I,NV
L=L+1
1 A(L)=0.
ANC=NC
IF(AMIS.EQ.YES) RETURN
30 M=0
DO 4 LL=1,NC
IF(NF.GT.0) READ (MT,F) (X(I),I=1,NV)
GO TO 1004
1004 IF(NF.LE.0) READ (MT) (X(I),I=1,NV)
GO TO 1005
1005 IF(.NOT.FSC) GO TO 35
DO 36 I=1,NV
IF(M.NE.255) GO TO 37
M=0
WRITE(2) C
200 FORMAT (20A4)
37 M=M+1
36 C(M)=X(I)
35 H=LL
DO 3 I=1,NV
D1(I)=(X(I)-S1(I))/H
S1(I)=S1(I)+D1(I)
3 IF(MOT.EQ.1) D1(I)=X(I)
H=H*(H-1.)
IF(MOT.EQ.1) H=1.
L=0
DO 4 I=1,NV
DD=D1(I)*H
DO 4 J=I,NV
L=L+1
4 A(L)=A(L)+DD*D1(J)
IF(.NOT.FSC) GO TO 10
WRITE (2) C
ENDFILE 2
REWIND 2
10 L=1
H=ANC
M=255
DO 60 I=1,NV
AL=A(L)
IF(MOT.EQ.1) AL=A(L)-H*S1(I)*S1(I)
D12=SQRT(AL/(H-1.))
LL=L
DO 220 J=I,NV
7 GO TO (400,6,400,6),MAT
400 A(L)=A(L)/H
GO TO 220
6 A(L)=A(L)/(H-1.)
220 L=L+1
X(I)=A(LL)
60 D1(I)=D12
WRITE (6,80) (I,S1(I),D1(I),I=1,NV)
80 FORMAT('1VARIABLE MEAN ST.DEV.'//(I6,2F15.6))
DO 61 I=1,NV
61 D1(I)=SQRT(X(I))
IF(MAT.LT.3) RETURN
L=0
DO 11 I=1,NV
DO 11 J=I,NV
L=L+1
11 A(L)=A(L)/(D1(I)*D1(J))
RETURN
END
C SUBROUTINE DOIT2 FOR BMDX72 MARCH 28, 1966
SUBROUTINE DOIT2(F,H,FO,R,X,S1,D1,A,AB,NV,LA,LAL)
DOUBLE PRECISION COMM,COM
DIMENSION FO(180),R(180),X(2),S1(2),D1(2),A(2),AB(NV,2)
X,F(2),H(2)
DOUBLE PRECISION PC
COMMON PC,NC,DUMMY,MAT,ICOM,IROT,GAMA,NROT,CONST,COR,VEC,PFS,WFS,
XMT,NF,RE,MFT,NFO,REO,AABBCC,FSC,MXC,MXR,AMIS,NVV
LOGICAL FSC
DATA COMM,YES/6HCOMMUN,3HYES/
GO TO (571,31,573,574),LA
571 IF(MAT.NE.6) GO TO 771
DO 723 I=1,NV
IF(NF.GT.0) READ (MT,F) (AB(I,J),J=1,NROT)
723 IF(NF.LE.0) READ (MT) (AB(I,J),J=1,NROT)
GO TO 100
771 L=1
K=NV
AABBCC=0.
DO 50 I=1,NV
AABBCC=AABBCC+A(L)
L=L+K
50 K=K-1
IF(.NOT.(FSC.OR.ICOM.EQ.2.OR.MXC.GT.1)) GO TO 310
C
C
C
WRITE (1) (A(L),L=1,NVV)
200 FORMAT (20A4)
LA=2
IF(FSC.OR.ICOM.EQ.2) RETURN
31 ENDFILE 1
REWIND 1
IF (FSC.OR.ICOM.EQ.2) GO TO 1798
310 GO TO (36,42,32,33),ICOM
32 DO 40 I=1,NV
40 R(I)=0.
L=0
NV1=NV-1
DO 41 I=1,NV1
L=L+1
I1=I+1
DO 41 J=I1,NV
L=L+1
R(I)=AMAX1(ABS(A(L)),R(I))
41 R(J)=AMAX1(ABS(A(L)),R(J))
GO TO 42
1798 CONTINUE
C
C
READ (1) (A(L),L=1,NVV)
C
GO TO 310
33 L1=0
38 L0=L1+1
L1=MIN0(L1+11,NV)
37 FORMAT(A6,11F6.6)
READ (5,37) COM,(R(L),L=L0,L1)
IF(COM.NE.COMM) GO TO 89
IF(L1.LT.NV) GO TO 38
GO TO 42
89 WRITE (6,98) COM,(R(L),L=L0,L1)
98 FORMAT(42H0UNRECOGNIZABLE COMMUNALITY ESTIMATES CARD/1X,A6,11F7.6)
STOP
42 L=1
K=NV
DO 34 I=1,NV
A(L)=R(I)
L=L+K
34 K=K-1
36 IF(COR.NE.YES) GO TO 85
IF(MAT.EQ.1) WRITE (6,831)
IF(MAT.EQ.2) WRITE (6,832)
IF(MAT.EQ.3) WRITE (6,833)
IF(MAT.EQ.4) WRITE (6,834)
831 FORMAT(35H1COVARIANCE MATRIX ABOUT THE ORIGIN)
832 FORMAT(18H1COVARIANCE MATRIX)
833 FORMAT(36H1CORRELATION MATRIX ABOUT THE ORIGIN)
834 FORMAT(19H1CORRELATION MATRIX)
L1=0
81 L0=L1+1
L1=MIN0(NV,L1+10)
L11=L0
K1=((L0-1)*(2*NV-L0+2))/2+1
WRITE (6,82) (L,L=L0,L1)
82 FORMAT(1H-,10I12)
DO 84 I=L0,NV
K=K1
K1=K1+1
K0=NV-L0
DO 90 L=L0,L11
X(L)=A(K)
K=K+K0
90 K0=K0-1
WRITE (6,86) I,(X(L),L=L0,L11)
86 FORMAT(I4,10F12.5)
84 L11=MIN0(L11+1,L1)
IF(L1.NE.NV) GO TO 81
85 IF(MXC.GT.1) WRITE (6,928)
928 FORMAT(28H-ITERATION FOR COMMUNALITIES/1H0,5X,14HMEAN ABS. DEV. 4X
X,14HMAX. ABS. DEV.)
LA=3
RETURN
573 ENDFILE 3
REWIND 3
IF (NROT.LE.0) GO TO 123
DO 12 I=1,NROT
H(I)=SQRT(H(I))
12 READ (3) (AB(J,I),J=1,NV)
REWIND 3
123 DO 812 I=1,NV
X(I)=0.
DO 812 J=1,NROT
AB(I,J)=AB(I,J)*H(J)
812 X(I)=X(I)+AB(I,J)*AB(I,J)
WRITE (6,813) (I,R(I),X(I),I=1,NV)
813 FORMAT(37H-VARIABLE ESTIMATED FINAL/12X,28HCOMMUNALITY
X COMMUNALITY/(I5,2F17.6))
IF(VEC.NE.YES) GO TO 100
LLL=1
WRITE (6,61)
61 FORMAT(30H1FACTOR MATRIX BEFORE ROTATION)
68 L1=0
62 L0=L1+1
L1=MIN0(L1+10,NROT)
WRITE (6,63) (L,L=L0,L1)
63 FORMAT(//20X,6HFACTOR/2X,10I12)
WRITE (6,633)
633 FORMAT(1X,8HVARIABLE)
DO 64 I=1,NV
64 WRITE (6,65) I,(AB(I,J),J=L0,L1)
IF(L1.NE.NROT) GO TO 62
65 FORMAT(I5,10F12.5)
LA=5
IF(LLL.NE.1) RETURN
100 LA=4
IF(IROT.EQ.0) LA=5
RETURN
574 WRITE (6,69)
69 FORMAT(22H1ROTATED FACTOR MATRIX)
LLL=2
GO TO 68
END
C SUBROUTINE MULTR FOR BMDX72 MARCH 28, 1966
SUBROUTINE MULTR(A,U,V,W,N,C)
DIMENSION A(2),U(2),V(2),W(2),C(255)
LOGICAL W
L=N
M=1
K=0
DO 1 I=1,N
V(I)=A(M)
M=M+L
L=L-1
IF(V(I).NE.0.) K=I
1 W(I)=.FALSE.
IF (K.EQ.0) GO TO 50
6 M=K-N
L=N
DO 2 I=1,N
IF(I.GT.K) L=1
M=M+L
L=L-1
U(I)=A(M)
2 A(M)=0.
P=U(K)
W(K)=.TRUE.
U(K)=-1.
M=1
T=0.
K=0
DO 5 I=1,N
Y=-U(I)/P
MM=M
DO 4 J=I,N
A(M)=A(M)+Y*U(J)
4 M=M+1
IF(W(I)) GO TO 5
IF (V(I).EQ.0.0) GO TO 5
H=A(MM)/V(I)
IF(H.LT.T) GO TO 5
T=H
K=I
5 CONTINUE
IF(T.GT.1.E-5) GO TO 6
50 L=N
M=1
DO 9 K=1,N
IF(W(K)) GO TO 7
U(K)=V(K)
GO TO 11
7 L1=N
U(K)=V(K)+1./A(M)
M1=1
L2=N
M2=K-N
DO 8 I=1,N
IF(I.GT.K) L2=1
M2=M2+L2
L2=L2-1
IF(W(I)) GO TO 10
IF(A(M1)-A(M2)*A(M2)/A(M) .GT. 1.E-5*V(I)) U(K)=V(K)
10 M1=M1+L1
L1=L1-1
8 CONTINUE
11 M=M+L
9 L=L-1
MM=0
DO 12 K=1,N
M=K-N
L=N
DO 12 I=1,N
IF(I.GT.K) L=1
M=M+L
L=L-1
IF(MM.NE.255) GO TO 20
MM=0
WRITE ( 1) C
200 FORMAT (20A4)
20 MM=MM+1
C(MM)=A(M)
12 IF(.NOT.(W(I).AND.W(K))) C(MM)=0.
WRITE (1) C
RETURN
END
C SUBROUTINE RAV FOR BMDX72 MARCH 28, 1966
SUBROUTINE RAV(A,C,N,U,V,W,G,H,R,P,ACC,CONST,NNN,MXC)
LOGICAL ITR
DIMENSION A(N),C(N)
DIMENSION U(N),V(N),W(N),G(N),H(N),R(N),P(N)
DATA CON/O201400000400/
EXTERNAL SIGN
ITR=.TRUE.
MX=0
L=1
M=N
DO 167 I=1,N
C(I)=A(L)
L=L+M
167 M=M-1
AG=U(1)
NNNN=NNN
158 DO 98 I=1,N
V(I)=0.
98 U(I)=0.
JJ=0
MX=MX+1
IF(MX.GE.MXC) ITR=.FALSE.
IF(ITR) REWIND 1
NM2=N-2
DO 91 KK=1,NM2
KP1=KK+1
AN=0.
II=JJ+1
DO 92 I=KP1,N
G(I)=0.
II=II+1
R(I)=A(II)+U(KK)*V(I)+U(I)*V(KK)
92 AN=AN+R(I)*R(I)
AN=SQRT(AN)
A2=R(KP1)
V(KK)=A(JJ+1)+2.*U(KK)*V(KK)
X=SQRT(1.+ABS(A2)/AN)
IF(X.GT.1.1) GO TO 45
A2=-A2
X=SQRT(1.-ABS(A2)/AN)
45 U(KK)=-SIGN(AN,A2)
Y=SIGN(1./(X*AN), A2)
W(KP1)=X
JJ=JJ+2
A(JJ)=X
KP2=KK+2
DO 93 I=KP2,N
JJ=JJ+1
W(I)=R(I)*Y
93 A(JJ)=W(I)
DO 95 K=KP1,N
II=II+1
A(II)=A(II)+U(K)*V(K)*2.
G(K)=G(K)+A(II)*W(K)
IF(K.EQ.N) GO TO 99
KPP1=K+1
DO 95 I=KPP1,N
II=II+1
A(II)=A(II)+U(I)*V(K)+V(I)*U(K)
G(I)=G(I)+A(II)*W(K)
95 G(K)=G(K)+A(II)*W(I)
99 UAU=0.
DO 96 I=KP1,N
UAU=UAU+G(I)*W(I)
96 U(I)=W(I)
DO 97 I=KP1,N
97 V(I)=UAU*U(I)/2.-G(I)
91 CONTINUE
X=U(N)*V(N-1)+U(N-1)*V(N)
Y=2.*U(N-1)*V(N-1)
Z=2.*U(N)*V(N)
U(N-1)=A(JJ+2)+X
V(N-1)=A(JJ+1)+Y
V(N)=A(JJ+3)+Z
JJJ=JJ-2
NM1=N-1
DO 15 I=1,NM1
U(I)=SIGN(AMAX1(1.E-9*ABS(V(I)),ABS(U(I))),U(I))
15 W(I)=U(I)*U(I)
D=AMAX1(ABS(V(1))+ABS(U(1)),ABS(V(N))+ABS(U(NM1)))
DO 10 I=2,NM1
10 D=AMAX1(ABS(U(I-1))+ABS(V(I))+ABS(U(I)),D)
DO 20 I=1,N
G(I)=-D
20 H(I)=D
DO 30 I=1,N
70 X=(G(I)+H(I))/2.
IF (X.EQ.G(I).OR.X.EQ.H(I)) GO TO 30
K=0
GG=V(1)-X
IF(GG.GT.0.) K=K+1
DO 40 J=2,N
IF (GG.NE.0.0) GG=V(J)-X-W(J-1)/GG
IF (GG.EQ.0.0) GG=V(J)-X
40 IF(GG.GT.0.) K=K+1
DO 50 J=I,N
IF(K.LT.J) H(J)=AMIN1(H(J),X)
50 IF(K.GE.J) G(J)=AMAX1(G(J),X)
GO TO 70
30 CONTINUE
GO TO 207
2077 WRITE (6,200) (H(I),I=1,N)
200 FORMAT(13H1EIGENVALUES //(10F12.5))
GG=0.
DO 201 I=1,N
GG=GG+H(I)
IF(H(I).GT.0.) N1N=I
IF (AG.NE.0.0) G(I)=GG/AG
IF (AG.EQ.0.0) G(I)=0.0
201 CONTINUE
WRITE (6,202) (G(I),I=1,N1N)
202 FORMAT(40H0CUMULATIVE PROPORTION OF TOTAL VARIANCE//(10F12.5))
RETURN
207 T=1.
DO 203 I=1,NM1
IF(H(I).NE.H(I+1) .OR. H(I).EQ.0.) GO TO 203
T=0.
WRITE(6,1) I
1 FORMAT(' EIGEN VALUE OF I= 'I4,' AND I+1 EQUAL, CHECK ACCURACY ')
H(I)=H(I)*CON
203 CONTINUE
IF(T.EQ.0.) GO TO 207
L=1
M=N
DO 102 I=1,N
A(L)=0.
L=L+M
102 M=M-1
510 DO 107 II=1,NNNN
IF(H(II).LE.CONST) GO TO 108
DO 57 I=1,N
57 R(I)=1.E-8
LIP=0
56 LIP=LIP+1
DO 51 I=1,N
G(I)=V(I)-H(II)
P(I)=U(I-1)
IF(ABS(R(I)).LT.1.E-17) R(I)=SIGN(1.E-17,R(I))
51 W(I)=U(I)
P(1)=0.
W(N)=0.
DO 52 I=1,NM1
IF(ABS(G(I)).GT.ABS(P(I+1))) GO TO 53
T=P(I+1)
P(I+1)=G(I)
G(I)=T
T=G(I+1)
G(I+1)=W(I)
W(I)=T
P(I)=W(I+1)
W(I+1)=0.
T=R(I+1)
R(I+1)=R(I)
R(I)=T
53 IF(G(I).NE.0.0) D=P(I+1)/G(I)
IF (G(I).EQ.0.0) D=0.0
G(I+1)=G(I+1)-W(I)*D
W(I+1)=W(I+1)-P(I)*D
R(I+1)=R(I+1)-D*R(I)
52 P(I+1)=0.
L=N
IF(G(L).EQ.0.0) G(L)=1.E-30
DO 54 I=1,N
IF (G(L).NE.0.0) R(L)=(R(L)-W(L)*R(L+1)-P(L)*R(L+2))/G(L)
IF (G(L).EQ.0.0) R(L)=0.0
IF(ABS(R(L)).LT.1.E8) GO TO 54
DO 85 J=1,N
IF(R(L).NE.0.0) R(J)=R(J)/R(L)
85 IF(R(L).EQ.0.0) R(J)=0.0
54 L=L-1
T=0.
DO 86 I=1,N
86 T=T+R(I)*R(I)
T=SQRT(T)
IF(LIP.EQ.1) T=T*1.E8
DO 83 I=1,N
IF (T.NE.0.0) R(I)=R(I)/T
IF (T.EQ.0.0) R(I)=0.0
83 CONTINUE
GO TO (56,55),LIP
55 JJ=JJJ
100 KQK=6
DO 64 K=1,NM2
NK=N-K
T=0.
DO 63 I=NK,N
JJ=JJ+1
G(I)=A(JJ)
63 T=T+G(I)*R(I)
JJ=JJ-KQK
KQK=KQK+2
DO 64 I=NK,N
64 R(I)=R(I)-T*G(I)
T=0.
TM=0
DO 65 I=1,N
IF(ABS(R(I)).GT.ABS(TM)) TM=R(I)
65 T=T+R(I)*R(I)
T=SIGN(SQRT(T),TM)
DO 66 I=1,N
IF (T.NE.0.0) R(I)=R(I)/T
IF (T.EQ.0.0) R(I)=0.0
66 CONTINUE
L=1
M=N
DO 101 I=1,N
A(L)=A(L)+R(I)*R(I)*H(II)
L=L+M
101 M=M-1
IF(.NOT.ITR) WRITE (3) (R(I),I=1,N)
107 CONTINUE
GO TO 109
108 NNN=II-1
109 IF (NNN.NE.0) GO TO 1099
WRITE (6,1098) H(I)
STOP
1098 FORMAT (' PROGRAM TERMINATED BECAUSE THE MAX. EIGENVALUE IS ',
1 F10.6/' OR THE FUNCTION SPECIFIED IS ZERO' )
1099 CONTINUE
X=0.
Y=0.
L=1
M=N
DO 110 I=1,N
Z=ABS(A(L)-C(I))
IF(ITR) C(I)=A(L)
L=L+M
M=M-1
X=X+Z
110 Y=AMAX1(Y,Z)
IF (N.NE.0.0) X=X/FLOAT(N)
IF (N.EQ.0.0) X=0.0
IF(.NOT.ITR) GO TO 2077
WRITE (6,916) MX,X,Y
IF(Y.LT.0.001) ITR=.FALSE.
916 FORMAT(I4,F11.4,F18.4)
NVV=(N*(N+1))/2
C
C
C
READ (1) (A(L),L=1,NVV)
300 FORMAT (20A4)
C
L=1
M=N
DO 152 I=1,N
A(L)=C(I)
L=L+M
152 M=M-1
GO TO 158
END
C SUBROUTINE ROTAT1 FOR BMDX72 MARCH 28, 1966
SUBROUTINE ROTAT1(A,TT,LV,LR,GG,ACC,MR,INI,NOR,S,C,V,TL,FL,XL,YL)
DIMENSION A(LV,LR),TT(LR,LR),C(2),S(2),V(2),TL(2),FL(2),XL(2),YL(2
X)
WRITE (6,882)
882 FORMAT(30H1ORTHOGONAL ROTATION
X/23H0ITERATION SIMPLICITY/13X,9HCRITERION)
NV=LV
NR=LR
GA=GG/FLOAT(NV)
IF(NOR.NE.1) GO TO 420
DO 40 I=1,NV
S(I)=0.
DO 4 J=1,NR
4 S(I)=S(I)+A(I,J)*A(I,J)
FL(I)=SQRT(S(I))
DO 2 J=1,NR
2 A(I,J)=A(I,J)/FL(I)
40 CONTINUE
420 DO 43 I=1,NR
41 C(I)=0.
DO 43 J=1,NV
43 C(I)=C(I)+A(J,I)*A(J,I)
D=0.
DO 44 J=1,NR
44 D=D+C(J)
F0=1.E20
50 DO 32 L=1,MR
G=0.
H=0.
DO 100 J=1,NR
100 V(J)=0.
DO 36 I=1,NV
S(I)=0.
DO 101 J=1,NR
T=A(I,J)*A(I,J)
S(I)=S(I)+T
101 V(J)=V(J)+T*T
36 H=H+S(I)*S(I)
DO 102 J=1,NR
102 G=G+V(J)-GA*C(J)*C(J)
F=H-GA*D*D-G
L1=L-1
WRITE (6,78) L1,F
78 FORMAT(I6,F16.6)
IF(L1.EQ.0) FCC=ABS(F*ACC)
IF(ABS(F-F0).LE.FCC) GO TO 38
F0=F
DO 32 IP=1,NR
DO 32 IQ=1,IP
IF(IP-IQ) 33,32,33
33 Y=0.
T=0.
R=0.
Z=0.
DO 34 I=1,NV
PQ=(A(I,IQ)+A(I,IP))*(A(I,IQ)-A(I,IP))
AB=A(I,IP)*A(I,IQ)
Z=Z+AB
T=T+AB*AB
Y=Y+AB*PQ
34 R=R+PQ*PQ
H=C(IQ)-C(IP)
X=T-GA*Z*Z-.25*(R-GA*H*H)
Y=Y-GA*Z*H
PHI=ATAN2(-Y,-X)/4.
S1=SIN(PHI)
C1=COS(PHI)
DO 35 I=1,NV
X=C1*A(I,IP)+S1*A(I,IQ)
A(I,IQ)=S1*A(I,IP)-C1*A(I,IQ)
35 A(I,IP)=X
Z=2.*S1*C1*Z
S1=S1*S1
C1=C1*C1
T=C(IP)
C(IP)=C1*T+S1*C(IQ)+Z
C(IQ)=S1*T+C1*C(IQ)-Z
32 CONTINUE
38 IF(NOR) 88,88,711
711 DO 46 I=1,NR
DO 46 J=1,NV
46 A(J,I)=A(J,I)*FL(J)
88 RETURN
END
C SUBROUTINE PFC FOR BMDX72 MARCH 28, 1966
SUBROUTINE PFC(A,N)
DIMENSION A(N,N)
WRITE (6,1)
1 FORMAT(26H-FACTOR CORRELATION MATRIX)
L1=0
2 L0=L1+1
L1=MIN0(N,L1+10)
WRITE (6,3) (L,L=L0,L1)
3 FORMAT(2H0 10I12)
DO 4 I=1,N
4 WRITE (6,5) I,(A(I,J),J=L0,L1)
5 FORMAT(I5,10F12.5)
IF(L1.NE.N) GO TO 2
RETURN
END
SUBROUTINE ROTAT2 (A,TT,LV,LR,GG,ACC,MR,INI,NOR,S,C,V,TL,FL,XL,YL,
XUPL)
DIMENSION A(LV,LR),TT(LR,LR),C(2),S(2),V(2),TL(2),FL(2),XL(2),
XYL(2)
WRITE (6,882)
882 FORMAT(30H1OBLIMIN ROTATION
X/23H0ITERATION SIMPLICITY/13X,9HCRITERION)
NR=LR
NV=LV
GA=GG/FLOAT(NV)
DO 30 I=1,NV
S(I)=0.0
DO 4 J=1,NR
4 S(I)=S(I)+A(I,J)*A(I,J)
IF(NOR) 30,30,32
32 FL(I)=SQRT(S(I))
S(I)=1.
DO 2 J=1,NR
2 A(I,J)=A(I,J)/FL(I)
30 CONTINUE
D=0.
DO 34 I=1,NR
IF(INI) 31,31,3
3 DO 6 J=1,NR
6 TT(I,J)=0.
TT(I,I)=1.
31 C(I)=0.
DO 33 J=1,NV
33 C(I)=C(I)+A(J,I)*A(J,I)
34 D=D+C(I)
60 IF(INI) 38,38,39
38 CALL MATINV(TT,NR,TOL,TL,XL,YL)
DO 40 I=1,NR
TL(I)=SQRT(TT(I,I))
DO 40 J=1,I
TT(I,J)=TT(I,J)/(TL(I)*TL(J))
40 TT(J,I)=TT(I,J)
F0=1.E20
39 DO 11 L=1,MR
F=0.
DO 200 J=1,NR
F=F-GA*C(J)*C(J)
DO 200 I=1,NV
T=A(I,J)*A(I,J)
200 F=F+T*T
F=-F-GA*D*D
DO 201 I=1,NV
201 F=F+S(I)*S(I)
L1=L-1
IF(L1.EQ.0) FCC=ABS(F*ACC)
WRITE (6,52) L1,F
52 FORMAT(I6,F16.6)
IF(ABS(F-F0).LE.FCC) GO TO 17
F0=F
DO 11 IP=1,NR
DO 11 IQ=1,NR
IF(IP-IQ) 18,11,18
45 A1=1.
A2=0.
GO TO 46
18 X=0.
Y=0.
Z=0.
D=D-C(IP)
DO 12 I=1,NV
AA=A(I,IP)*A(I,IP)
S(I)=S(I)-AA
W=S(I)-GA*D
X=X+AA*W
Y=Y+A(I,IP)*A(I,IQ)*W
12 Z=Z+A(I,IQ)*A(I,IQ)*W
O=TT(IP,IQ)
A1=1.-O*O
B1=Y*O-.5*(X+Z)
C1=X*Z-Y*Y
AL=B1*B1-A1*C1
IF(AL) 45,49,49
49 AL=-(B1+SQRT(AL))/A1
A1=AL-Z
A2=Y-AL*O
IF(A1) 23,24,23
24 A1=A2
A2=AL-X
23 W=A1*A1+A2*(A2+2.*A1*O)
IF(ABS(A1).LE.ABS(.001*A2) .OR. W.LE.0.) GO TO 45
W=SQRT(W)
A1=A1/W
A2=A2/W
DO 14 I=1,NR
IF(I.EQ.IP) GO TO 14
V(I)=A1*TT(I,IP)+A2*TT(I,IQ)
IF(ABS(V(I)).GT.UPL) GO TO 45
14 CONTINUE
DO 210 I=1,NR
TT(I,IP)=V(I)
210 TT(IP,I)=V(I)
46 C(IP)=0.
DO 13 I=1,NV
A(I,IP)=A1*A(I,IP)+A2*A(I,IQ)
AA=A(I,IP)*A(I,IP)
S(I)=S(I)+AA
13 C(IP)=C(IP)+AA
D=D+C(IP)
TT(IP,IP)=1.
11 CONTINUE
17 CALL MATINV(TT,NR,TOL,TL,XL,YL)
DO 19 I=1,NR
S(I)=SQRT(TT(I,I))
DO 20 J=1,I
TT(I,J)=TT(I,J)/(S(I)*S(J))
20 TT(J,I)=TT(I,J)
DO 19 J=1,NV
19 A(J,I)=A(J,I)*S(I)
89 CALL PFC(TT,NR)
IF(NOR) 88,88,71
71 DO 36 I=1,NR
DO 36 J=1,NV
36 A(J,I)=A(J,I)*FL(J)
88 RETURN
END
SUBROUTINE MATINV(A,N,T,IN,V,U)
DIMENSION A(N,N),IN(2),V(2),U(2)
DO 1 I=1,N
V(I)=A(I,I)
1 IN(I)=0
K=1
DO 7 L=1,N
DO 2 I=1,N
U(I)=A(K,I)
IF(I.GT.K) U(I)=A(I,K)
A(I,K)=0.
2 A(K,I)=0.
P=U(K)
T=H
H=-1.E20
IN(K)=1
U(K)=-1.
DO 7 I=1,N
Y=U(I)/P
DO 4 J=1,I
4 A(I,J)=A(I,J)-Y*U(J)
IF(IN(I)) 5,5,7
5 IF(H-A(I,I)/V(I)) 6,7,7
6 H=A(I,I)/V(I)
K=I
7 CONTINUE
DO 8 I=1,N
DO 8 J=1,I
A(I,J)=-A(I,J)
8 A(J,I)=A(I,J)
RETURN
END
C SUBROUTINE ROTAT3 FOR BMDX72 MARCH 28, 1966
SUBROUTINE ROTAT3(A,TT,LV,LR,GG,ACC,MR,INI,NOR,S,C,V,TL,FL,XL,YL)
DIMENSION A(LV,LR),TT(LR,LR),C(2),S(2),V(2),TL(2),FL(2),XL(2),YL(2
X)
WRITE (6,882)
882 FORMAT(40H1ROTATION FOR SIMPLE LOADINGS
X/23H0ITERATION SIMPLICITY/13X,9HCRITERION)
NV=LV
NR=LR
GA=GG/FLOAT(NV)
DO 30 I=1,NV
S(I)=0.
DO 4 J=1,NR
4 S(I)=S(I)+A(I,J)*A(I,J)
IF(NOR) 30,30,32
32 FL(I)=SQRT(S(I))
S(I)=1.
DO 2 J=1,NR
2 A(I,J)=A(I,J)/FL(I)
30 CONTINUE
DO 33 I=1,NR
IF(INI) 31,31,3
3 DO 6 J=1,NR
6 TT(I,J)=0.
TT(I,I)=1.
31 C(I)=0.
V(I)=0.
DO 33 J=1,NV
AA=A(J,I)*A(J,I)
C(I)=C(I)+AA
33 V(I)=V(I)+AA*AA
G=0.
D=0.
DO 34 J=1,NR
D=D+C(J)
V(J)=V(J)-GA*C(J)*C(J)
34 G=G+V(J)
H=0.
DO 5 I=1,NV
5 H=H+S(I)*S(I)
H=H-GA*D*D
F0=H-G
FCC=F0*ACC
L=0
WRITE (6,52) L,F0
52 FORMAT(I6,F16.6)
70 DO 80 L=1,MR
DO 81 IP=1,NR
DO 81 IQ=1,NR
IF(IP-IQ) 82,81,82
82 D=D-C(IP)-C(IQ)
G=G-V(IP)-V(IQ)
P=0.
R=0.
T=0.
O=0.
Y=0.
Z=0.
DO 83 I=1,NV
A1=A(I,IP)
A2=A(I,IQ)
AA=A1*A1
BB=A2*A2
AB=A1*A2
S(I)=S(I)-AA-BB
Z=Z+AB
R=R+AA*BB
P=P+AB*AA
T=T+AA*S(I)
83 O=O+AB*S(I)
X=TT(IP,IQ)
GAC=GA*C(IP)
R=R-GAC*C(IQ)
P=P-GAC*Z
O=O-GA*Z*D
U=U-GA*C(IQ)*D
T=T-GAC*D
100 P1=1.5*(X-P/V(IP))
Q1=.5*(V(IP)-4.*X*P+R+2.*T)/V(IP)
R1=.5*(X*(T+R)-P-O)/V(IP)
CALL ROOT(P1,Q1,R1,A2)
A22=A2*A2
P=1.+2.*X*A2+A22
IF(P)89,90,90
89 WRITE (6,91) P
91 FORMAT(2H $F20.7)
P=-P
90 A1=SQRT(P)
A3=A2/A1
A11=P*P
C(IP)=P*C(IP)
V(IP)=A11*V(IP)
Z=0.
Y=0.
DO 84 I=1,NV
A(I,IQ)=A(I,IQ)-A2*A(I,IP)
A(I,IP)=A1*A(I,IP)
BB=A(I,IQ)*A(I,IQ)
Z=Z+BB*BB
Y=Y+BB
84 S(I)=S(I)+A(I,IP)*A(I,IP)+BB
V(IQ)=Z-GA*Y*Y
C(IQ)=Y
D=D+Y+C(IP)
G=G+V(IP)+V(IQ)
A1=1./A1
DO 85 I=1,NR
TT(I,IP)=A1*TT(I,IP)+A3*TT(I,IQ)
85 TT(IP,I)=TT(I,IP)
TT(IP,IP)=1.
81 CONTINUE
H=0.
DO 86 I=1,NV
86 H=H+S(I)*S(I)
F=H-GA*D*D-G
WRITE (6,52) L,F
157 IF(ABS(F-F0)-FCC) 38,38,80
80 F0=F
38 CONTINUE
CALL PFC(TT,NR)
IF(NOR) 88,88,71
71 DO 36 I=1,NR
DO 36 J=1,NV
36 A(J,I)=A(J,I)*FL(J)
88 RETURN
END
C SUBROUTINE ROOT FOR BMDX72 MARCH 28, 1966
SUBROUTINE ROOT(P,Q,R,X)
H=(ABS(P)+ABS(Q)+ABS(R))*1.E-5
P2=2.*P
A=P*P-3.*Q
IF(A) 1,1,2
1 X=0.
GO TO 5
2 A=SQRT(A)
X=(A-P)/3.
X1=-(P+A)/3.
IF(X*(Q+X*(P+X))+X1*(Q+X1*(P+X1))+R+R)3,4,4
3 X=X+1.
GO TO5
4 X=X1-1.
5 DO 7 I=1,50
F=R+X*(Q+X*(P+X))
FP=Q+X*(P2+3.*X)
DX=F/FP
X=X-DX
IF(ABS(F)-H) 6,6,7
7 CONTINUE
WRITE (6,8) P,Q,R,X
8 FORMAT(2H *4E20.8)
6 RETURN
END
SUBROUTINE ATOF(A,N,F)
DIMENSION A(1)
LOGICAL BLANK
BLANK=.TRUE.
S=1.0
NUMB=0
TEN=1.0
DIV=1.0
DO 10 I=1,N
L=INTCHR(A,I)
IF(L.EQ.36) GO TO 10
BLANK=.FALSE.
IF(L.NE.38) GO TO 2
S=-1.0
GO TO 10
2 IF(L.NE.44) GO TO 4
TEN=10.0
GO TO 10
4 IF(L.GT.9) GO TO 9
NUMB=NUMB*10+L
DIV=DIV*TEN
9 CONTINUE
10 CONTINUE
IF(BLANK)RETURN
F=S*FLOAT(NUMB)/DIV
RETURN
END
FUNCTION INTCHR(STRING,N)
DIMENSION SEQ(50),STRING(1),EBCD(5)
DATA SEQ/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,
X 1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,
X 1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,
X 1HU,1HV,1HW,1HX,1HY,1HZ,1H ,1H+,1H-,1H*,
X 1H/,1H(,1H),1H,,1H.,1H',1H=,1H$,1H ,1H /
DATA EBCD/1H+,1H(,1H),1H',1H=/
CALL GETCHR(STRING,N,CHR)
IF (CHR.NE.SEQ(37)) GO TO 2
INTCHR = 36
GO TO 10
2 DO 1 I=1,48
IF(SEQ(I).EQ.CHR) GO TO 9
1 CONTINUE
I=51
IF(EBCD(1).EQ.CHR) I=38
IF(EBCD(2).EQ.CHR) I=42
IF(EBCD(3).EQ.CHR) I=43
IF(EBCD(4).EQ.CHR) I=46
IF(EBCD(5).EQ.CHR) I=47
9 INTCHR=I-1
10 RETURN
END