Changeset 1289 for LMDZ4/trunk/libf

Timestamp:

Dec 18, 2009, 3:51:22 PM (15 years ago)

Author:

Ehouarn Millour

Message:

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

File:

• LMDZ4/trunk/libf/filtrez/jacobi.F90 (moved) (moved from LMDZ4/trunk/libf/filtrez/jacobi.F) (3 diffs)

LMDZ4/trunk/libf/filtrez/jacobi.F90

@@ -1 +1 @@
 !
-! $Header$
+! $Id$
 !
       SUBROUTINE JACOBI(A,N,NP,D,V,NROT)
-      PARAMETER (NMAX=400)
-      DIMENSION A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)
-      IF (n.gt.nmax) THEN
-         print*, 'n, nmax=', n, nmax
-         print*, 'Surdimensionnement insuffisant dans jacobi'
-         CALL abort
-      ENDIF
-      DO 12 IP=1,N
-        DO 11 IQ=1,N
+      implicit none
+! Arguments:
+      integer,intent(in) :: N
+      integer,intent(in) :: NP
+      integer,intent(out) :: NROT
+      real,intent(inout) :: A(NP,NP)
+      real,intent(out) :: D(NP)
+      real,intent(out) :: V(NP,NP)
+
+! local variables:
+      integer :: IP,IQ,I,J
+      real :: SM,TRESH,G,H,T,THETA,C,S,TAU
+      real :: B(N)
+      real :: Z(N)
+      
+      DO IP=1,N
+        DO IQ=1,N
           V(IP,IQ)=0.
-11      CONTINUE
+        ENDDO
         V(IP,IP)=1.
-12    CONTINUE
-      DO 13 IP=1,N
+      ENDDO
+      DO IP=1,N
         B(IP)=A(IP,IP)
         D(IP)=B(IP)
         Z(IP)=0.
-13    CONTINUE
+      ENDDO
       NROT=0
-      DO 24 I=1,50
+      DO I=1,50 ! 50? I suspect this should be NP
+                !     but convergence is fast enough anyway
         SM=0.
-        DO 15 IP=1,N-1
-          DO 14 IQ=IP+1,N
+        DO IP=1,N-1
+          DO IQ=IP+1,N
             SM=SM+ABS(A(IP,IQ))
-14        CONTINUE
-15      CONTINUE
+          ENDDO
+        ENDDO
         IF(SM.EQ.0.)RETURN
         IF(I.LT.4)THEN
@@ -35 +44 @@
           TRESH=0.
         ENDIF
-        DO 22 IP=1,N-1
-          DO 21 IQ=IP+1,N
+        DO IP=1,N-1
+          DO IQ=IP+1,N
             G=100.*ABS(A(IP,IQ))
-            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP)))
-     *         .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
+            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) &
+               .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
               A(IP,IQ)=0.
             ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
@@ -59 +68 @@
               D(IQ)=D(IQ)+H
               A(IP,IQ)=0.
-              DO 16 J=1,IP-1
+              DO J=1,IP-1
                 G=A(J,IP)
                 H=A(J,IQ)
                 A(J,IP)=G-S*(H+G*TAU)
                 A(J,IQ)=H+S*(G-H*TAU)
-16            CONTINUE
-              DO 17 J=IP+1,IQ-1
+             ENDDO
+              DO J=IP+1,IQ-1
                 G=A(IP,J)
                 H=A(J,IQ)
                 A(IP,J)=G-S*(H+G*TAU)
                 A(J,IQ)=H+S*(G-H*TAU)
-17            CONTINUE
-              DO 18 J=IQ+1,N
+              ENDDO
+              DO J=IQ+1,N
                 G=A(IP,J)
                 H=A(IQ,J)
                 A(IP,J)=G-S*(H+G*TAU)
                 A(IQ,J)=H+S*(G-H*TAU)
-18            CONTINUE
-              DO 19 J=1,N
+              ENDDO
+              DO J=1,N
                 G=V(J,IP)
                 H=V(J,IQ)
                 V(J,IP)=G-S*(H+G*TAU)
                 V(J,IQ)=H+S*(G-H*TAU)
-19            CONTINUE
+              ENDDO
               NROT=NROT+1
             ENDIF
-21        CONTINUE
-22      CONTINUE
-        DO 23 IP=1,N
+          ENDDO
+        ENDDO
+        DO IP=1,N
           B(IP)=B(IP)+Z(IP)
           D(IP)=B(IP)
           Z(IP)=0.
-23      CONTINUE
-24    CONTINUE
-      STOP '50 iterations should never happen'
+        ENDDO
+      ENDDO ! of DO I=1,50
+      STOP 'Jacobi: 50 iterations should never happen'
       RETURN
-      END
+      END SUBROUTINE JACOBI