Directory: | ./ |
---|---|
File: | filtrez/jacobi.f90 |
Date: | 2022-01-11 19:19:34 |
Exec | Total | Coverage | |
---|---|---|---|
Lines: | 65 | 66 | 98.5% |
Branches: | 41 | 42 | 97.6% |
Line | Branch | Exec | Source |
---|---|---|---|
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 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 2 times.
|
66 | DO IP=1,N |
21 |
2/2✓ Branch 0 taken 2048 times.
✓ Branch 1 taken 64 times.
|
2112 | DO IQ=1,N |
22 | 2112 | V(IP,IQ)=0. | |
23 | ENDDO | ||
24 | 66 | V(IP,IP)=1. | |
25 | ENDDO | ||
26 |
2/2✓ Branch 0 taken 64 times.
✓ Branch 1 taken 2 times.
|
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 |
1/2✓ Branch 0 taken 24 times.
✗ Branch 1 not taken.
|
24 | DO I=1,50 ! 50? I suspect this should be NP |
33 | ! but convergence is fast enough anyway | ||
34 | SM=0. | ||
35 |
2/2✓ Branch 0 taken 744 times.
✓ Branch 1 taken 24 times.
|
768 | DO IP=1,N-1 |
36 |
2/2✓ Branch 0 taken 744 times.
✓ Branch 1 taken 11904 times.
|
12672 | DO IQ=IP+1,N |
37 | 12648 | SM=SM+ABS(A(IP,IQ)) | |
38 | ENDDO | ||
39 | ENDDO | ||
40 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 22 times.
|
24 | IF(SM.EQ.0.)RETURN |
41 |
2/2✓ Branch 0 taken 6 times.
✓ Branch 1 taken 16 times.
|
22 | IF(I.LT.4)THEN |
42 | 6 | TRESH=0.2*SM/N**2 | |
43 | ELSE | ||
44 | TRESH=0. | ||
45 | ENDIF | ||
46 |
2/2✓ Branch 0 taken 682 times.
✓ Branch 1 taken 22 times.
|
704 | DO IP=1,N-1 |
47 |
2/2✓ Branch 0 taken 682 times.
✓ Branch 1 taken 10912 times.
|
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 |
6/6✓ Branch 0 taken 6944 times.
✓ Branch 1 taken 3968 times.
✓ Branch 2 taken 3424 times.
✓ Branch 3 taken 3520 times.
✓ Branch 4 taken 3402 times.
✓ Branch 5 taken 22 times.
|
11594 | .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN |
51 | 3402 | A(IP,IQ)=0. | |
52 |
2/2✓ Branch 0 taken 7112 times.
✓ Branch 1 taken 398 times.
|
7510 | ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN |
53 | 7112 | H=D(IQ)-D(IP) | |
54 |
2/2✓ Branch 0 taken 112 times.
✓ Branch 1 taken 7000 times.
|
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 |
2/2✓ Branch 0 taken 3706 times.
✓ Branch 1 taken 3294 times.
|
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 |
2/2✓ Branch 0 taken 69492 times.
✓ Branch 1 taken 7112 times.
|
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 |
2/2✓ Branch 0 taken 68938 times.
✓ Branch 1 taken 7112 times.
|
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 |
2/2✓ Branch 0 taken 74930 times.
✓ Branch 1 taken 7112 times.
|
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 |
2/2✓ Branch 0 taken 227584 times.
✓ Branch 1 taken 7112 times.
|
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 |
2/2✓ Branch 0 taken 704 times.
✓ Branch 1 taken 22 times.
|
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 | ||
107 |