Last change
on this file since 4109 was
576,
checked in by (none), 20 years ago
|
This commit was manufactured by cvs2svn to create tag
'IPSL-CM4_PSMILE'.
|
-
Property svn:eol-style set to
native
-
Property svn:executable set to
*
-
Property svn:keywords set to
Author Date Id Revision
|
File size:
2.5 KB
|
Line | |
---|
1 | SUBROUTINE JACOBI(A,N,NP,D,V,NROT) |
---|
2 | PARAMETER (NMAX=400) |
---|
3 | DIMENSION A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX) |
---|
4 | IF (n.gt.nmax) THEN |
---|
5 | print*, 'n, nmax=', n, nmax |
---|
6 | print*, 'Surdimensionnement insuffisant dans jacobi' |
---|
7 | CALL abort |
---|
8 | ENDIF |
---|
9 | DO 12 IP=1,N |
---|
10 | DO 11 IQ=1,N |
---|
11 | V(IP,IQ)=0. |
---|
12 | 11 CONTINUE |
---|
13 | V(IP,IP)=1. |
---|
14 | 12 CONTINUE |
---|
15 | DO 13 IP=1,N |
---|
16 | B(IP)=A(IP,IP) |
---|
17 | D(IP)=B(IP) |
---|
18 | Z(IP)=0. |
---|
19 | 13 CONTINUE |
---|
20 | NROT=0 |
---|
21 | DO 24 I=1,50 |
---|
22 | SM=0. |
---|
23 | DO 15 IP=1,N-1 |
---|
24 | DO 14 IQ=IP+1,N |
---|
25 | SM=SM+ABS(A(IP,IQ)) |
---|
26 | 14 CONTINUE |
---|
27 | 15 CONTINUE |
---|
28 | IF(SM.EQ.0.)RETURN |
---|
29 | IF(I.LT.4)THEN |
---|
30 | TRESH=0.2*SM/N**2 |
---|
31 | ELSE |
---|
32 | TRESH=0. |
---|
33 | ENDIF |
---|
34 | DO 22 IP=1,N-1 |
---|
35 | DO 21 IQ=IP+1,N |
---|
36 | G=100.*ABS(A(IP,IQ)) |
---|
37 | IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) |
---|
38 | * .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN |
---|
39 | A(IP,IQ)=0. |
---|
40 | ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN |
---|
41 | H=D(IQ)-D(IP) |
---|
42 | IF(ABS(H)+G.EQ.ABS(H))THEN |
---|
43 | T=A(IP,IQ)/H |
---|
44 | ELSE |
---|
45 | THETA=0.5*H/A(IP,IQ) |
---|
46 | T=1./(ABS(THETA)+SQRT(1.+THETA**2)) |
---|
47 | IF(THETA.LT.0.)T=-T |
---|
48 | ENDIF |
---|
49 | C=1./SQRT(1+T**2) |
---|
50 | S=T*C |
---|
51 | TAU=S/(1.+C) |
---|
52 | H=T*A(IP,IQ) |
---|
53 | Z(IP)=Z(IP)-H |
---|
54 | Z(IQ)=Z(IQ)+H |
---|
55 | D(IP)=D(IP)-H |
---|
56 | D(IQ)=D(IQ)+H |
---|
57 | A(IP,IQ)=0. |
---|
58 | DO 16 J=1,IP-1 |
---|
59 | G=A(J,IP) |
---|
60 | H=A(J,IQ) |
---|
61 | A(J,IP)=G-S*(H+G*TAU) |
---|
62 | A(J,IQ)=H+S*(G-H*TAU) |
---|
63 | 16 CONTINUE |
---|
64 | DO 17 J=IP+1,IQ-1 |
---|
65 | G=A(IP,J) |
---|
66 | H=A(J,IQ) |
---|
67 | A(IP,J)=G-S*(H+G*TAU) |
---|
68 | A(J,IQ)=H+S*(G-H*TAU) |
---|
69 | 17 CONTINUE |
---|
70 | DO 18 J=IQ+1,N |
---|
71 | G=A(IP,J) |
---|
72 | H=A(IQ,J) |
---|
73 | A(IP,J)=G-S*(H+G*TAU) |
---|
74 | A(IQ,J)=H+S*(G-H*TAU) |
---|
75 | 18 CONTINUE |
---|
76 | DO 19 J=1,N |
---|
77 | G=V(J,IP) |
---|
78 | H=V(J,IQ) |
---|
79 | V(J,IP)=G-S*(H+G*TAU) |
---|
80 | V(J,IQ)=H+S*(G-H*TAU) |
---|
81 | 19 CONTINUE |
---|
82 | NROT=NROT+1 |
---|
83 | ENDIF |
---|
84 | 21 CONTINUE |
---|
85 | 22 CONTINUE |
---|
86 | DO 23 IP=1,N |
---|
87 | B(IP)=B(IP)+Z(IP) |
---|
88 | D(IP)=B(IP) |
---|
89 | Z(IP)=0. |
---|
90 | 23 CONTINUE |
---|
91 | 24 CONTINUE |
---|
92 | STOP '50 iterations should never happen' |
---|
93 | RETURN |
---|
94 | END |
---|
Note: See
TracBrowser
for help on using the repository browser.