source: LMDZ5/branches/LMDZ6_rc0/libf/filtrez/jacobi.F90 @ 5017

Last change on this file since 5017 was 1910, checked in by Laurent Fairhead, 11 years ago

Merged trunk changes r1860:1909 into testing branch
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • 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 1910 2013-11-29 08:40:25Z dcugnet $
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.