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