/[lmdze]/trunk/filtrez/jacobi.f
ViewVC logotype

Diff of /trunk/filtrez/jacobi.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/filtrez/jacobi.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/filtrez/jacobi.f90 revision 81 by guez, Wed Mar 5 14:38:41 2014 UTC
# Line 1  Line 1 
1  !  
2  ! $Header: /home/cvsroot/LMDZ4/libf/filtrez/jacobi.F,v 1.1.1.1 2004/05/19 12:53:09 lmdzadmin Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/filtrez/jacobi.F,v 1.1.1.1 2004/05/19
3  !  ! 12:53:09 lmdzadmin Exp $
4        SUBROUTINE JACOBI(A,N,NP,D,V,NROT)  
5        PARAMETER (NMAX=400)  SUBROUTINE jacobi(a, n, np, d, v, nrot)
6        DIMENSION A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)    PARAMETER (nmax=400)
7        IF (n.gt.nmax) THEN    DIMENSION a(np, np), d(np), v(np, np), b(nmax), z(nmax)
8           print*, 'n, nmax=', n, nmax  
9           print*, 'Surdimensionnement insuffisant dans jacobi'    IF (n>nmax) THEN
10           stop 1      PRINT *, 'n, nmax=', n, nmax
11        ENDIF      PRINT *, 'Surdimensionnement insuffisant dans jacobi'
12        DO 12 IP=1,N      STOP 1
13          DO 11 IQ=1,N    END IF
14            V(IP,IQ)=0.    DO ip = 1, n
15  11      CONTINUE      DO iq = 1, n
16          V(IP,IP)=1.        v(ip, iq) = 0.
17  12    CONTINUE      END DO
18        DO 13 IP=1,N      v(ip, ip) = 1.
19          B(IP)=A(IP,IP)    END DO
20          D(IP)=B(IP)    DO ip = 1, n
21          Z(IP)=0.      b(ip) = a(ip, ip)
22  13    CONTINUE      d(ip) = b(ip)
23        NROT=0      z(ip) = 0.
24        DO 24 I=1,50    END DO
25          SM=0.    nrot = 0
26          DO 15 IP=1,N-1    DO i = 1, 50
27            DO 14 IQ=IP+1,N      sm = 0.
28              SM=SM+ABS(A(IP,IQ))      DO ip = 1, n - 1
29  14        CONTINUE        DO iq = ip + 1, n
30  15      CONTINUE          sm = sm + abs(a(ip,iq))
31          IF(SM.EQ.0.)RETURN        END DO
32          IF(I.LT.4)THEN      END DO
33            TRESH=0.2*SM/N**2      IF (sm==0.) RETURN
34          ELSE      IF (i<4) THEN
35            TRESH=0.        tresh = 0.2*sm/n**2
36          ENDIF      ELSE
37          DO 22 IP=1,N-1        tresh = 0.
38            DO 21 IQ=IP+1,N      END IF
39              G=100.*ABS(A(IP,IQ))      DO ip = 1, n - 1
40              IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP)))        DO iq = ip + 1, n
41       *         .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN          g = 100.*abs(a(ip,iq))
42                A(IP,IQ)=0.          IF ((i>4) .AND. (abs(d(ip))+g==abs(d(ip))) .AND. (abs(d(iq))+g==abs( &
43              ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN              d(iq)))) THEN
44                H=D(IQ)-D(IP)            a(ip, iq) = 0.
45                IF(ABS(H)+G.EQ.ABS(H))THEN          ELSE IF (abs(a(ip,iq))>tresh) THEN
46                  T=A(IP,IQ)/H            h = d(iq) - d(ip)
47                ELSE            IF (abs(h)+g==abs(h)) THEN
48                  THETA=0.5*H/A(IP,IQ)              t = a(ip, iq)/h
49                  T=1./(ABS(THETA)+SQRT(1.+THETA**2))            ELSE
50                  IF(THETA.LT.0.)T=-T              theta = 0.5*h/a(ip, iq)
51                ENDIF              t = 1./(abs(theta)+sqrt(1.+theta**2))
52                C=1./SQRT(1+T**2)              IF (theta<0.) t = -t
53                S=T*C            END IF
54                TAU=S/(1.+C)            c = 1./sqrt(1+t**2)
55                H=T*A(IP,IQ)            s = t*c
56                Z(IP)=Z(IP)-H            tau = s/(1.+c)
57                Z(IQ)=Z(IQ)+H            h = t*a(ip, iq)
58                D(IP)=D(IP)-H            z(ip) = z(ip) - h
59                D(IQ)=D(IQ)+H            z(iq) = z(iq) + h
60                A(IP,IQ)=0.            d(ip) = d(ip) - h
61                DO 16 J=1,IP-1            d(iq) = d(iq) + h
62                  G=A(J,IP)            a(ip, iq) = 0.
63                  H=A(J,IQ)            DO j = 1, ip - 1
64                  A(J,IP)=G-S*(H+G*TAU)              g = a(j, ip)
65                  A(J,IQ)=H+S*(G-H*TAU)              h = a(j, iq)
66  16            CONTINUE              a(j, ip) = g - s*(h+g*tau)
67                DO 17 J=IP+1,IQ-1              a(j, iq) = h + s*(g-h*tau)
68                  G=A(IP,J)            END DO
69                  H=A(J,IQ)            DO j = ip + 1, iq - 1
70                  A(IP,J)=G-S*(H+G*TAU)              g = a(ip, j)
71                  A(J,IQ)=H+S*(G-H*TAU)              h = a(j, iq)
72  17            CONTINUE              a(ip, j) = g - s*(h+g*tau)
73                DO 18 J=IQ+1,N              a(j, iq) = h + s*(g-h*tau)
74                  G=A(IP,J)            END DO
75                  H=A(IQ,J)            DO j = iq + 1, n
76                  A(IP,J)=G-S*(H+G*TAU)              g = a(ip, j)
77                  A(IQ,J)=H+S*(G-H*TAU)              h = a(iq, j)
78  18            CONTINUE              a(ip, j) = g - s*(h+g*tau)
79                DO 19 J=1,N              a(iq, j) = h + s*(g-h*tau)
80                  G=V(J,IP)            END DO
81                  H=V(J,IQ)            DO j = 1, n
82                  V(J,IP)=G-S*(H+G*TAU)              g = v(j, ip)
83                  V(J,IQ)=H+S*(G-H*TAU)              h = v(j, iq)
84  19            CONTINUE              v(j, ip) = g - s*(h+g*tau)
85                NROT=NROT+1              v(j, iq) = h + s*(g-h*tau)
86              ENDIF            END DO
87  21        CONTINUE            nrot = nrot + 1
88  22      CONTINUE          END IF
89          DO 23 IP=1,N        END DO
90            B(IP)=B(IP)+Z(IP)      END DO
91            D(IP)=B(IP)      DO ip = 1, n
92            Z(IP)=0.        b(ip) = b(ip) + z(ip)
93  23      CONTINUE        d(ip) = b(ip)
94  24    CONTINUE        z(ip) = 0.
95        STOP '50 iterations should never happen'      END DO
96        RETURN    END DO
97        END    STOP '50 iterations should never happen'
98      RETURN
99    END SUBROUTINE jacobi
Legend:
Removed from v.76  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21