|
Last change
on this file since 3663 was
1289,
checked in by Ehouarn Millour, 15 years ago
|
|
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
|
-
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 1289 2009-12-18 14:51:22Z oboucher $ |
|---|
| [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.
