-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathtred2.f
80 lines (80 loc) · 1.85 KB
/
tred2.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
SUBROUTINE TRED2(A,N,NP,D,E)
DIMENSION A(NP,NP),D(NP),E(NP)
IF(N.GT.1)THEN
DO 18 I=N,2,-1
L=I-1
H=0.
SCALE=0.
IF(L.GT.1)THEN
DO 11 K=1,L
SCALE=SCALE+ABS(A(I,K))
11 CONTINUE
IF(SCALE.EQ.0.)THEN
E(I)=A(I,L)
ELSE
DO 12 K=1,L
A(I,K)=A(I,K)/SCALE
H=H+A(I,K)**2
12 CONTINUE
F=A(I,L)
G=-SIGN(SQRT(H),F)
E(I)=SCALE*G
H=H-F*G
A(I,L)=F-G
F=0.
DO 15 J=1,L
A(J,I)=A(I,J)/H
G=0.
DO 13 K=1,J
G=G+A(J,K)*A(I,K)
13 CONTINUE
IF(L.GT.J)THEN
DO 14 K=J+1,L
G=G+A(K,J)*A(I,K)
14 CONTINUE
ENDIF
E(J)=G/H
F=F+E(J)*A(I,J)
15 CONTINUE
HH=F/(H+H)
DO 17 J=1,L
F=A(I,J)
G=E(J)-HH*F
E(J)=G
DO 16 K=1,J
A(J,K)=A(J,K)-F*E(K)-G*A(I,K)
16 CONTINUE
17 CONTINUE
ENDIF
ELSE
E(I)=A(I,L)
ENDIF
D(I)=H
18 CONTINUE
ENDIF
D(1)=0.
E(1)=0.
DO 23 I=1,N
L=I-1
IF(D(I).NE.0.)THEN
DO 21 J=1,L
G=0.
DO 19 K=1,L
G=G+A(I,K)*A(K,J)
19 CONTINUE
DO 20 K=1,L
A(K,J)=A(K,J)-G*A(K,I)
20 CONTINUE
21 CONTINUE
ENDIF
D(I)=A(I,I)
A(I,I)=1.
IF(L.GE.1)THEN
DO 22 J=1,L
A(I,J)=0.
A(J,I)=0.
22 CONTINUE
ENDIF
23 CONTINUE
RETURN
END