|
Last change
on this file since 2253 was
1907,
checked in by lguez, 11 years ago
|
|
Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:
$ svn propget copyright libf/phylmd/physiq.F90
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
Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.
|
-
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
|
| Rev | Line | |
|---|
| [524] | 1 | ! |
|---|
| [1289] | 2 | ! $Id: jacobi.F90 1907 2013-11-26 13:10:46Z jyg $ |
|---|
| [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.
