source: LMDZ4/branches/pre_V3/libf/filtrez/jacobi.F

Last change on this file was 524, checked in by lmdzadmin, 21 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.5 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE JACOBI(A,N,NP,D,V,NROT)
5      PARAMETER (NMAX=400)
6      DIMENSION A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)
7      IF (n.gt.nmax) THEN
8         print*, 'n, nmax=', n, nmax
9         print*, 'Surdimensionnement insuffisant dans jacobi'
10         CALL abort
11      ENDIF
12      DO 12 IP=1,N
13        DO 11 IQ=1,N
14          V(IP,IQ)=0.
1511      CONTINUE
16        V(IP,IP)=1.
1712    CONTINUE
18      DO 13 IP=1,N
19        B(IP)=A(IP,IP)
20        D(IP)=B(IP)
21        Z(IP)=0.
2213    CONTINUE
23      NROT=0
24      DO 24 I=1,50
25        SM=0.
26        DO 15 IP=1,N-1
27          DO 14 IQ=IP+1,N
28            SM=SM+ABS(A(IP,IQ))
2914        CONTINUE
3015      CONTINUE
31        IF(SM.EQ.0.)RETURN
32        IF(I.LT.4)THEN
33          TRESH=0.2*SM/N**2
34        ELSE
35          TRESH=0.
36        ENDIF
37        DO 22 IP=1,N-1
38          DO 21 IQ=IP+1,N
39            G=100.*ABS(A(IP,IQ))
40            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP)))
41     *         .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
42              A(IP,IQ)=0.
43            ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
44              H=D(IQ)-D(IP)
45              IF(ABS(H)+G.EQ.ABS(H))THEN
46                T=A(IP,IQ)/H
47              ELSE
48                THETA=0.5*H/A(IP,IQ)
49                T=1./(ABS(THETA)+SQRT(1.+THETA**2))
50                IF(THETA.LT.0.)T=-T
51              ENDIF
52              C=1./SQRT(1+T**2)
53              S=T*C
54              TAU=S/(1.+C)
55              H=T*A(IP,IQ)
56              Z(IP)=Z(IP)-H
57              Z(IQ)=Z(IQ)+H
58              D(IP)=D(IP)-H
59              D(IQ)=D(IQ)+H
60              A(IP,IQ)=0.
61              DO 16 J=1,IP-1
62                G=A(J,IP)
63                H=A(J,IQ)
64                A(J,IP)=G-S*(H+G*TAU)
65                A(J,IQ)=H+S*(G-H*TAU)
6616            CONTINUE
67              DO 17 J=IP+1,IQ-1
68                G=A(IP,J)
69                H=A(J,IQ)
70                A(IP,J)=G-S*(H+G*TAU)
71                A(J,IQ)=H+S*(G-H*TAU)
7217            CONTINUE
73              DO 18 J=IQ+1,N
74                G=A(IP,J)
75                H=A(IQ,J)
76                A(IP,J)=G-S*(H+G*TAU)
77                A(IQ,J)=H+S*(G-H*TAU)
7818            CONTINUE
79              DO 19 J=1,N
80                G=V(J,IP)
81                H=V(J,IQ)
82                V(J,IP)=G-S*(H+G*TAU)
83                V(J,IQ)=H+S*(G-H*TAU)
8419            CONTINUE
85              NROT=NROT+1
86            ENDIF
8721        CONTINUE
8822      CONTINUE
89        DO 23 IP=1,N
90          B(IP)=B(IP)+Z(IP)
91          D(IP)=B(IP)
92          Z(IP)=0.
9323      CONTINUE
9424    CONTINUE
95      STOP '50 iterations should never happen'
96      RETURN
97      END
Note: See TracBrowser for help on using the repository browser.