1 |
|
|
! |
2 |
|
|
! $Id: jacobi.F90 1907 2013-11-26 13:10:46Z lguez $ |
3 |
|
|
! |
4 |
|
2 |
SUBROUTINE JACOBI(A,N,NP,D,V,NROT) |
5 |
|
|
implicit none |
6 |
|
|
! Arguments: |
7 |
|
|
integer,intent(in) :: N |
8 |
|
|
integer,intent(in) :: NP |
9 |
|
|
integer,intent(out) :: NROT |
10 |
|
|
real,intent(inout) :: A(NP,NP) |
11 |
|
|
real,intent(out) :: D(NP) |
12 |
|
|
real,intent(out) :: V(NP,NP) |
13 |
|
|
|
14 |
|
|
! local variables: |
15 |
|
|
integer :: IP,IQ,I,J |
16 |
|
|
real :: SM,TRESH,G,H,T,THETA,C,S,TAU |
17 |
|
4 |
real :: B(N) |
18 |
|
2 |
real :: Z(N) |
19 |
|
|
|
20 |
✓✓ |
66 |
DO IP=1,N |
21 |
✓✓ |
2112 |
DO IQ=1,N |
22 |
|
2112 |
V(IP,IQ)=0. |
23 |
|
|
ENDDO |
24 |
|
66 |
V(IP,IP)=1. |
25 |
|
|
ENDDO |
26 |
✓✓ |
66 |
DO IP=1,N |
27 |
|
64 |
B(IP)=A(IP,IP) |
28 |
|
64 |
D(IP)=B(IP) |
29 |
|
66 |
Z(IP)=0. |
30 |
|
|
ENDDO |
31 |
|
2 |
NROT=0 |
32 |
✓✗ |
24 |
DO I=1,50 ! 50? I suspect this should be NP |
33 |
|
|
! but convergence is fast enough anyway |
34 |
|
|
SM=0. |
35 |
✓✓ |
768 |
DO IP=1,N-1 |
36 |
✓✓ |
12672 |
DO IQ=IP+1,N |
37 |
|
12648 |
SM=SM+ABS(A(IP,IQ)) |
38 |
|
|
ENDDO |
39 |
|
|
ENDDO |
40 |
✓✓ |
24 |
IF(SM.EQ.0.)RETURN |
41 |
✓✓ |
22 |
IF(I.LT.4)THEN |
42 |
|
6 |
TRESH=0.2*SM/N**2 |
43 |
|
|
ELSE |
44 |
|
|
TRESH=0. |
45 |
|
|
ENDIF |
46 |
✓✓ |
704 |
DO IP=1,N-1 |
47 |
✓✓ |
11616 |
DO IQ=IP+1,N |
48 |
|
10912 |
G=100.*ABS(A(IP,IQ)) |
49 |
|
|
IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) & |
50 |
✓✓✓✓
✓✓ |
11594 |
.AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN |
51 |
|
3402 |
A(IP,IQ)=0. |
52 |
✓✓ |
7510 |
ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN |
53 |
|
7112 |
H=D(IQ)-D(IP) |
54 |
✓✓ |
7112 |
IF(ABS(H)+G.EQ.ABS(H))THEN |
55 |
|
112 |
T=A(IP,IQ)/H |
56 |
|
|
ELSE |
57 |
|
7000 |
THETA=0.5*H/A(IP,IQ) |
58 |
|
7000 |
T=1./(ABS(THETA)+SQRT(1.+THETA**2)) |
59 |
✓✓ |
7000 |
IF(THETA.LT.0.)T=-T |
60 |
|
|
ENDIF |
61 |
|
7112 |
C=1./SQRT(1+T**2) |
62 |
|
7112 |
S=T*C |
63 |
|
7112 |
TAU=S/(1.+C) |
64 |
|
7112 |
H=T*A(IP,IQ) |
65 |
|
7112 |
Z(IP)=Z(IP)-H |
66 |
|
7112 |
Z(IQ)=Z(IQ)+H |
67 |
|
7112 |
D(IP)=D(IP)-H |
68 |
|
7112 |
D(IQ)=D(IQ)+H |
69 |
|
7112 |
A(IP,IQ)=0. |
70 |
✓✓ |
76604 |
DO J=1,IP-1 |
71 |
|
69492 |
G=A(J,IP) |
72 |
|
69492 |
H=A(J,IQ) |
73 |
|
69492 |
A(J,IP)=G-S*(H+G*TAU) |
74 |
|
76604 |
A(J,IQ)=H+S*(G-H*TAU) |
75 |
|
|
ENDDO |
76 |
✓✓ |
76050 |
DO J=IP+1,IQ-1 |
77 |
|
68938 |
G=A(IP,J) |
78 |
|
68938 |
H=A(J,IQ) |
79 |
|
68938 |
A(IP,J)=G-S*(H+G*TAU) |
80 |
|
76050 |
A(J,IQ)=H+S*(G-H*TAU) |
81 |
|
|
ENDDO |
82 |
✓✓ |
82042 |
DO J=IQ+1,N |
83 |
|
74930 |
G=A(IP,J) |
84 |
|
74930 |
H=A(IQ,J) |
85 |
|
74930 |
A(IP,J)=G-S*(H+G*TAU) |
86 |
|
82042 |
A(IQ,J)=H+S*(G-H*TAU) |
87 |
|
|
ENDDO |
88 |
✓✓ |
234696 |
DO J=1,N |
89 |
|
227584 |
G=V(J,IP) |
90 |
|
227584 |
H=V(J,IQ) |
91 |
|
227584 |
V(J,IP)=G-S*(H+G*TAU) |
92 |
|
234696 |
V(J,IQ)=H+S*(G-H*TAU) |
93 |
|
|
ENDDO |
94 |
|
7112 |
NROT=NROT+1 |
95 |
|
|
ENDIF |
96 |
|
|
ENDDO |
97 |
|
|
ENDDO |
98 |
✓✓ |
726 |
DO IP=1,N |
99 |
|
704 |
B(IP)=B(IP)+Z(IP) |
100 |
|
704 |
D(IP)=B(IP) |
101 |
|
726 |
Z(IP)=0. |
102 |
|
|
ENDDO |
103 |
|
|
ENDDO ! of DO I=1,50 |
104 |
|
|
STOP 'Jacobi: 50 iterations should never happen' |
105 |
|
|
RETURN |
106 |
|
|
END SUBROUTINE JACOBI |