source: LMDZ5/branches/testing/libf/filtrez/jacobi.F90 @ 4328

Last change on this file since 4328 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
RevLine 
[524]1!
[1289]2! $Id: jacobi.F90 1910 2013-11-29 08:40:25Z dcugnet $
[524]3!
4      SUBROUTINE JACOBI(A,N,NP,D,V,NROT)
[1289]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
[524]22          V(IP,IQ)=0.
[1289]23        ENDDO
[524]24        V(IP,IP)=1.
[1289]25      ENDDO
26      DO IP=1,N
[524]27        B(IP)=A(IP,IP)
28        D(IP)=B(IP)
29        Z(IP)=0.
[1289]30      ENDDO
[524]31      NROT=0
[1289]32      DO I=1,50 ! 50? I suspect this should be NP
33                !     but convergence is fast enough anyway
[524]34        SM=0.
[1289]35        DO IP=1,N-1
36          DO IQ=IP+1,N
[524]37            SM=SM+ABS(A(IP,IQ))
[1289]38          ENDDO
39        ENDDO
[524]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
[1289]46        DO IP=1,N-1
47          DO IQ=IP+1,N
[524]48            G=100.*ABS(A(IP,IQ))
[1289]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
[524]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.
[1289]70              DO J=1,IP-1
[524]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)
[1289]75             ENDDO
76              DO J=IP+1,IQ-1
[524]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)
[1289]81              ENDDO
82              DO J=IQ+1,N
[524]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)
[1289]87              ENDDO
88              DO J=1,N
[524]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)
[1289]93              ENDDO
[524]94              NROT=NROT+1
95            ENDIF
[1289]96          ENDDO
97        ENDDO
98        DO IP=1,N
[524]99          B(IP)=B(IP)+Z(IP)
100          D(IP)=B(IP)
101          Z(IP)=0.
[1289]102        ENDDO
103      ENDDO ! of DO I=1,50
104      STOP 'Jacobi: 50 iterations should never happen'
[524]105      RETURN
[1289]106      END SUBROUTINE JACOBI
Note: See TracBrowser for help on using the repository browser.