trunk/filtrez/jacobi.f


! $Header: /home/cvsroot/LMDZ4/libf/filtrez/jacobi.F,v 1.1.1.1 2004/05/19
! 12:53:09 lmdzadmin Exp $

SUBROUTINE jacobi(a, n, np, d, v, nrot)
  implicit none
  integer, PARAMETER:: nmax=400
  integer n, np
  real a(np, np), d(np), v(np, np), b(nmax), z(nmax)
  integer ip, iq, nrot, i, j
  real sm, g, h, t, theta, c, s, tau, tresh

  IF (n>nmax) THEN
    PRINT *, 'n, nmax=', n, nmax
    PRINT *, 'Surdimensionnement insuffisant dans jacobi'
    STOP 1
  END IF
  DO ip = 1, n
    DO iq = 1, n
      v(ip, iq) = 0.
    END DO
    v(ip, ip) = 1.
  END DO
  DO ip = 1, n
    b(ip) = a(ip, ip)
    d(ip) = b(ip)
    z(ip) = 0.
  END DO
  nrot = 0
  DO i = 1, 50
    sm = 0.
    DO ip = 1, n - 1
      DO iq = ip + 1, n
        sm = sm + abs(a(ip,iq))
      END DO
    END DO
    IF (sm==0.) RETURN
    IF (i<4) THEN
      tresh = 0.2*sm/n**2
    ELSE
      tresh = 0.
    END IF
    DO ip = 1, n - 1
      DO iq = ip + 1, n
        g = 100.*abs(a(ip,iq))
        IF ((i>4) .AND. (abs(d(ip))+g==abs(d(ip))) .AND. (abs(d(iq))+g==abs( &
            d(iq)))) THEN
          a(ip, iq) = 0.
        ELSE IF (abs(a(ip,iq))>tresh) THEN
          h = d(iq) - d(ip)
          IF (abs(h)+g==abs(h)) THEN
            t = a(ip, iq)/h
          ELSE
            theta = 0.5*h/a(ip, iq)
            t = 1./(abs(theta)+sqrt(1.+theta**2))
            IF (theta<0.) t = -t
          END IF
          c = 1./sqrt(1+t**2)
          s = t*c
          tau = s/(1.+c)
          h = t*a(ip, iq)
          z(ip) = z(ip) - h
          z(iq) = z(iq) + h
          d(ip) = d(ip) - h
          d(iq) = d(iq) + h
          a(ip, iq) = 0.
          DO j = 1, ip - 1
            g = a(j, ip)
            h = a(j, iq)
            a(j, ip) = g - s*(h+g*tau)
            a(j, iq) = h + s*(g-h*tau)
          END DO
          DO j = ip + 1, iq - 1
            g = a(ip, j)
            h = a(j, iq)
            a(ip, j) = g - s*(h+g*tau)
            a(j, iq) = h + s*(g-h*tau)
          END DO
          DO j = iq + 1, n
            g = a(ip, j)
            h = a(iq, j)
            a(ip, j) = g - s*(h+g*tau)
            a(iq, j) = h + s*(g-h*tau)
          END DO
          DO j = 1, n
            g = v(j, ip)
            h = v(j, iq)
            v(j, ip) = g - s*(h+g*tau)
            v(j, iq) = h + s*(g-h*tau)
          END DO
          nrot = nrot + 1
        END IF
      END DO
    END DO
    DO ip = 1, n
      b(ip) = b(ip) + z(ip)
      d(ip) = b(ip)
      z(ip) = 0.
    END DO
  END DO
  STOP '50 iterations should never happen'
  RETURN
END SUBROUTINE jacobi
1
2	! $Header: /home/cvsroot/LMDZ4/libf/filtrez/jacobi.F,v 1.1.1.1 2004/05/19
3	! 12:53:09 lmdzadmin Exp $
4
5	SUBROUTINE jacobi(a, n, np, d, v, nrot)
6	implicit none
7	integer, PARAMETER:: nmax=400
8	integer n, np
9	real a(np, np), d(np), v(np, np), b(nmax), z(nmax)
10	integer ip, iq, nrot, i, j
11	real sm, g, h, t, theta, c, s, tau, tresh
12
13	IF (n>nmax) THEN
14	PRINT *, 'n, nmax=', n, nmax
15	PRINT *, 'Surdimensionnement insuffisant dans jacobi'
16	STOP 1
17	END IF
18	DO ip = 1, n
19	DO iq = 1, n
20	v(ip, iq) = 0.
21	END DO
22	v(ip, ip) = 1.
23	END DO
24	DO ip = 1, n
25	b(ip) = a(ip, ip)
26	d(ip) = b(ip)
27	z(ip) = 0.
28	END DO
29	nrot = 0
30	DO i = 1, 50
31	sm = 0.
32	DO ip = 1, n - 1
33	DO iq = ip + 1, n
34	sm = sm + abs(a(ip,iq))
35	END DO
36	END DO
37	IF (sm==0.) RETURN
38	IF (i<4) THEN
39	tresh = 0.2sm/n*2
40	ELSE
41	tresh = 0.
42	END IF
43	DO ip = 1, n - 1
44	DO iq = ip + 1, n
45	g = 100.*abs(a(ip,iq))
46	IF ((i>4) .AND. (abs(d(ip))+g==abs(d(ip))) .AND. (abs(d(iq))+g==abs( &
47	d(iq)))) THEN
48	a(ip, iq) = 0.
49	ELSE IF (abs(a(ip,iq))>tresh) THEN
50	h = d(iq) - d(ip)
51	IF (abs(h)+g==abs(h)) THEN
52	t = a(ip, iq)/h
53	ELSE
54	theta = 0.5*h/a(ip, iq)
55	t = 1./(abs(theta)+sqrt(1.+theta**2))
56	IF (theta<0.) t = -t
57	END IF
58	c = 1./sqrt(1+t**2)
59	s = t*c
60	tau = s/(1.+c)
61	h = t*a(ip, iq)
62	z(ip) = z(ip) - h
63	z(iq) = z(iq) + h
64	d(ip) = d(ip) - h
65	d(iq) = d(iq) + h
66	a(ip, iq) = 0.
67	DO j = 1, ip - 1
68	g = a(j, ip)
69	h = a(j, iq)
70	a(j, ip) = g - s(h+gtau)
71	a(j, iq) = h + s(g-htau)
72	END DO
73	DO j = ip + 1, iq - 1
74	g = a(ip, j)
75	h = a(j, iq)
76	a(ip, j) = g - s(h+gtau)
77	a(j, iq) = h + s(g-htau)
78	END DO
79	DO j = iq + 1, n
80	g = a(ip, j)
81	h = a(iq, j)
82	a(ip, j) = g - s(h+gtau)
83	a(iq, j) = h + s(g-htau)
84	END DO
85	DO j = 1, n
86	g = v(j, ip)
87	h = v(j, iq)
88	v(j, ip) = g - s(h+gtau)
89	v(j, iq) = h + s(g-htau)
90	END DO
91	nrot = nrot + 1
92	END IF
93	END DO
94	END DO
95	DO ip = 1, n
96	b(ip) = b(ip) + z(ip)
97	d(ip) = b(ip)
98	z(ip) = 0.
99	END DO
100	END DO
101	STOP '50 iterations should never happen'
102	RETURN
103	END SUBROUTINE jacobi