source: LMDZ.3.3/trunk/libf/filtrez/jacobi.F @ 2

Last change on this file since 2 was 2, checked in by lmdz, 25 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 2.5 KB
Line 
1      SUBROUTINE JACOBI(A,N,NP,D,V,NROT)
2      PARAMETER (NMAX=200)
3      DIMENSION A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)
4      IF (n.gt.nmax) THEN
5         print*, 'n, nmax=', n, nmax
6         print*, 'Surdimensionnement insuffisant dans jacobi'
7         CALL abort
8      ENDIF
9      DO 12 IP=1,N
10        DO 11 IQ=1,N
11          V(IP,IQ)=0.
1211      CONTINUE
13        V(IP,IP)=1.
1412    CONTINUE
15      DO 13 IP=1,N
16        B(IP)=A(IP,IP)
17        D(IP)=B(IP)
18        Z(IP)=0.
1913    CONTINUE
20      NROT=0
21      DO 24 I=1,50
22        SM=0.
23        DO 15 IP=1,N-1
24          DO 14 IQ=IP+1,N
25            SM=SM+ABS(A(IP,IQ))
2614        CONTINUE
2715      CONTINUE
28        IF(SM.EQ.0.)RETURN
29        IF(I.LT.4)THEN
30          TRESH=0.2*SM/N**2
31        ELSE
32          TRESH=0.
33        ENDIF
34        DO 22 IP=1,N-1
35          DO 21 IQ=IP+1,N
36            G=100.*ABS(A(IP,IQ))
37            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP)))
38     *         .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
39              A(IP,IQ)=0.
40            ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
41              H=D(IQ)-D(IP)
42              IF(ABS(H)+G.EQ.ABS(H))THEN
43                T=A(IP,IQ)/H
44              ELSE
45                THETA=0.5*H/A(IP,IQ)
46                T=1./(ABS(THETA)+SQRT(1.+THETA**2))
47                IF(THETA.LT.0.)T=-T
48              ENDIF
49              C=1./SQRT(1+T**2)
50              S=T*C
51              TAU=S/(1.+C)
52              H=T*A(IP,IQ)
53              Z(IP)=Z(IP)-H
54              Z(IQ)=Z(IQ)+H
55              D(IP)=D(IP)-H
56              D(IQ)=D(IQ)+H
57              A(IP,IQ)=0.
58              DO 16 J=1,IP-1
59                G=A(J,IP)
60                H=A(J,IQ)
61                A(J,IP)=G-S*(H+G*TAU)
62                A(J,IQ)=H+S*(G-H*TAU)
6316            CONTINUE
64              DO 17 J=IP+1,IQ-1
65                G=A(IP,J)
66                H=A(J,IQ)
67                A(IP,J)=G-S*(H+G*TAU)
68                A(J,IQ)=H+S*(G-H*TAU)
6917            CONTINUE
70              DO 18 J=IQ+1,N
71                G=A(IP,J)
72                H=A(IQ,J)
73                A(IP,J)=G-S*(H+G*TAU)
74                A(IQ,J)=H+S*(G-H*TAU)
7518            CONTINUE
76              DO 19 J=1,N
77                G=V(J,IP)
78                H=V(J,IQ)
79                V(J,IP)=G-S*(H+G*TAU)
80                V(J,IQ)=H+S*(G-H*TAU)
8119            CONTINUE
82              NROT=NROT+1
83            ENDIF
8421        CONTINUE
8522      CONTINUE
86        DO 23 IP=1,N
87          B(IP)=B(IP)+Z(IP)
88          D(IP)=B(IP)
89          Z(IP)=0.
9023      CONTINUE
9124    CONTINUE
92      STOP '50 iterations should never happen'
93      RETURN
94      END
Note: See TracBrowser for help on using the repository browser.