My Project
 All Classes Files Functions Variables Macros
jacobi.F90
Go to the documentation of this file.
1 !
2 ! $Id: jacobi.F90 1289 2009-12-18 14:51:22Z emillour $
3 !
4  SUBROUTINE jacobi(A,N,NP,D,V,NROT)
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
22  v(ip,iq)=0.
23  ENDDO
24  v(ip,ip)=1.
25  ENDDO
26  DO ip=1,n
27  b(ip)=a(ip,ip)
28  d(ip)=b(ip)
29  z(ip)=0.
30  ENDDO
31  nrot=0
32  DO i=1,50 ! 50? I suspect this should be NP
33  ! but convergence is fast enough anyway
34  sm=0.
35  DO ip=1,n-1
36  DO iq=ip+1,n
37  sm=sm+abs(a(ip,iq))
38  ENDDO
39  ENDDO
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
46  DO ip=1,n-1
47  DO iq=ip+1,n
48  g=100.*abs(a(ip,iq))
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
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.
70  DO j=1,ip-1
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)
75  ENDDO
76  DO j=ip+1,iq-1
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)
81  ENDDO
82  DO j=iq+1,n
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)
87  ENDDO
88  DO j=1,n
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)
93  ENDDO
94  nrot=nrot+1
95  ENDIF
96  ENDDO
97  ENDDO
98  DO ip=1,n
99  b(ip)=b(ip)+z(ip)
100  d(ip)=b(ip)
101  z(ip)=0.
102  ENDDO
103  ENDDO ! of DO I=1,50
104  stop 'Jacobi: 50 iterations should never happen'
105  RETURN
106  END SUBROUTINE jacobi