source: LMDZ5/branches/LMDZ5_AR5/libf/filtrez/jacobi.F90 @ 4085

Last change on this file since 4085 was 1289, checked in by Ehouarn Millour, 15 years ago

Minor improvement: Bring Jacobi into the XXIst century.
Routine can now handle matix of any size (tested on local Linux machines and IBM Power6 that this doesn't change results compared to old version).

EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.7 KB
Line 
1!
2! $Id: jacobi.F90 1289 2009-12-18 14:51:22Z fhourdin $
3!
4      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      real :: B(N)
18      real :: Z(N)
19     
20      DO IP=1,N
21        DO IQ=1,N
22          V(IP,IQ)=0.
23        ENDDO
24        V(IP,IP)=1.
25      ENDDO
26      DO IP=1,N
27        B(IP)=A(IP,IP)
28        D(IP)=B(IP)
29        Z(IP)=0.
30      ENDDO
31      NROT=0
32      DO I=1,50 ! 50? I suspect this should be NP
33                !     but convergence is fast enough anyway
34        SM=0.
35        DO IP=1,N-1
36          DO IQ=IP+1,N
37            SM=SM+ABS(A(IP,IQ))
38          ENDDO
39        ENDDO
40        IF(SM.EQ.0.)RETURN
41        IF(I.LT.4)THEN
42          TRESH=0.2*SM/N**2
43        ELSE
44          TRESH=0.
45        ENDIF
46        DO IP=1,N-1
47          DO IQ=IP+1,N
48            G=100.*ABS(A(IP,IQ))
49            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) &
50               .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
51              A(IP,IQ)=0.
52            ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
53              H=D(IQ)-D(IP)
54              IF(ABS(H)+G.EQ.ABS(H))THEN
55                T=A(IP,IQ)/H
56              ELSE
57                THETA=0.5*H/A(IP,IQ)
58                T=1./(ABS(THETA)+SQRT(1.+THETA**2))
59                IF(THETA.LT.0.)T=-T
60              ENDIF
61              C=1./SQRT(1+T**2)
62              S=T*C
63              TAU=S/(1.+C)
64              H=T*A(IP,IQ)
65              Z(IP)=Z(IP)-H
66              Z(IQ)=Z(IQ)+H
67              D(IP)=D(IP)-H
68              D(IQ)=D(IQ)+H
69              A(IP,IQ)=0.
70              DO J=1,IP-1
71                G=A(J,IP)
72                H=A(J,IQ)
73                A(J,IP)=G-S*(H+G*TAU)
74                A(J,IQ)=H+S*(G-H*TAU)
75             ENDDO
76              DO J=IP+1,IQ-1
77                G=A(IP,J)
78                H=A(J,IQ)
79                A(IP,J)=G-S*(H+G*TAU)
80                A(J,IQ)=H+S*(G-H*TAU)
81              ENDDO
82              DO J=IQ+1,N
83                G=A(IP,J)
84                H=A(IQ,J)
85                A(IP,J)=G-S*(H+G*TAU)
86                A(IQ,J)=H+S*(G-H*TAU)
87              ENDDO
88              DO J=1,N
89                G=V(J,IP)
90                H=V(J,IQ)
91                V(J,IP)=G-S*(H+G*TAU)
92                V(J,IQ)=H+S*(G-H*TAU)
93              ENDDO
94              NROT=NROT+1
95            ENDIF
96          ENDDO
97        ENDDO
98        DO IP=1,N
99          B(IP)=B(IP)+Z(IP)
100          D(IP)=B(IP)
101          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
Note: See TracBrowser for help on using the repository browser.