File: DIAGON.FT of Tape: Various/ETH/f2
(Source file text)
SUBROUTINE DIAGON
COMMON /FITPAR/ JIT,FISTOP,AMARQI,UPMARQ,DNMARQ,VARI,
$ UPVAR,DNVAR,VARMIN,MLOOP,LOOPLW,ISTART,ISTOP,ITEST,IDEFIX
COMMON /DIAG/ A(24,24),Z(24,24),D(24),E(24)
COMMON /SIZE/ N,NP,K1,K3
COMMON /VECT/ C(24),D1(6,4),Q(6,4),DIAGEL(24)
C
LOGICAL ITEST
C TEST OUTPUT
600 FORMAT(1H0,16(F7.2,1X))
601 FORMAT(1H )
602 FORMAT(1H ,16(F7.2,1X))
C
DATA EPS/1.E-07/
DATA TOL/1.E-07/
IF (.NOT.ITEST) GO TO 350
WRITE(3,600) (C(J),J=1,N)
WRITE(3,601)
DO 300 K=1,N
300 WRITE(3,602) (A(K,J),J=1,N)
350 Z(1,1)=1.
IF(N.LE.1) GO TO 219
DO 1 I=1,N
DO 1 J=1,N
Z(J,I)=A(I,J)
1 CONTINUE
C
DO 9 II=2,N
L=N-II
I=L+2
F=Z(I,I-1)
G=0.
IF(L.EQ.0) GOTO 100
DO 2 K=1,L
G=G+Z(I,K)**2
2 CONTINUE
100 CONTINUE
H=G+F*F
IF(G.GT.TOL) GOTO 3
E(I)=F
H=0.
GOTO 9
3 L=L+1
G=-SIGN(SQRT(H),F)
E(I)=G
H=H-F*G
Z(I,I-1)=F-G
F=0.
DO 6 J=1,L
Z(J,I)=Z(I,J)/H
G=0.
DO 4 K=1,J
G=G+Z(J,K)*Z(I,K)
4 CONTINUE
J1=J+1
IF(J1.GT.L) GOTO 110
DO 5 K=J1,L
G=G+Z(K,J)*Z(I,K)
5 CONTINUE
110 CONTINUE
E(J)=G/H
F=F+G*Z(J,I)
6 CONTINUE
HH=F/(H+H)
DO 8 J=1,L
F=Z(I,J)
G=E(J)-HH*F
E(J)=G
DO 7 K=1,J
Z(J,K)=Z(J,K)-F*E(K)-G*Z(I,K)
7 CONTINUE
8 CONTINUE
9 D(I)=H
D(1)=0.
E(1)=0.
DO 15 I=1,N
L=I-1
IF(ABS(D(I)).LT.1.E-08)GO TO 13
IF(L.EQ.0) GOTO 13
DO 12 J=1,L
G=0.
DO 10 K=1,L
G=G+Z(I,K)*Z(K,J)
10 CONTINUE
DO 11 K=1,L
Z(K,J)=Z(K,J)-G*Z(K,I)
11 CONTINUE
12 CONTINUE
13 D(I)=Z(I,I)
Z(I,I)=1.
IF(L.EQ.0) GOTO 15
DO 14 J=1,L
Z(I,J)=0.
Z(J,I)=0.
14 CONTINUE
15 CONTINUE
DO 201 I=2,N
E(I-1)=E(I)
201 CONTINUE
E(N)=0.
B=0.
F=0.
J=0
JIT=0
DO 212 L=1,N
J=0
H=EPS*(ABS(D(L))+ABS(E(L)))
IF(B.LT.H) B=H
DO 202 M=L,N
IF(ABS(E(M))-B)203,203,202
202 CONTINUE
203 CONTINUE
IF(M.EQ.L) GOTO 211
204 CONTINUE
IF(J.EQ.30) GOTO 217
J=J+1
P=(D(L+1)-D(L))/(2.*E(L))
R=SQRT(P*P+1.)
H = D(L) - E(L)/(P+SIGN(R,P))
DO 205 I=L,N
D(I)=D(I)-H
205 CONTINUE
F=F+H
P=D(M)
CO=1.
S=0.
M1=M-1
ML=M1+L
DO 210 II=L,M1
I=ML-II
G=CO*E(I)
H=CO*P
IF(ABS(P)-ABS(E(I)))207,206,206
206 CONTINUE
CO=E(I)/P
R=SQRT(CO*CO+1.)
E(I+1)=S*P*R
S=CO/R
CO=1./R
GO TO 208
207 CONTINUE
CO=P/E(I)
R=SQRT(CO*CO+1.)
E(I+1)=S*E(I)*R
S=1./R
CO=CO/R
208 CONTINUE
P=CO*D(I)-S*G
D(I+1)=H+S*(CO*G+S*D(I))
DO 209 K=1,N
H=Z(K,I+1)
Z(K,I+1)=S*Z(K,I)+CO*H
Z(K,I)=CO*Z(K,I)-S*H
209 CONTINUE
210 CONTINUE
E(L)=S*P
D(L)=CO*P
IF(ABS(E(L))-B)211,211,204
211 CONTINUE
D(L)=D(L)+F
JIT=JIT+J
212 CONTINUE
N1=N-1
DO 216 I=1,N1
K=I
P=D(I)
II=I+1
DO 214 J=II,N
IF(D(J)-P) 213,214,214
213 CONTINUE
K=J
P=D(J)
214 CONTINUE
IF(K.EQ.I) GOTO 216
D(K)=D(I)
D(I)=P
DO 215 J=1,N
P=Z(J,I)
Z(J,I)=Z(J,K)
Z(J,K)=P
215 CONTINUE
216 CONTINUE
219 CONTINUE
IF (.NOT.ITEST) GO TO 360
WRITE(3,600) (D(J),J=1,N)
WRITE(3,601)
DO 301 K=1,N
301 WRITE(3,602) (Z(K,J),J=1,N)
360 RETURN
217 CONTINUE
STOP 11111
END