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)
  PARAMETER (nmax=400)
  DIMENSION a(np, np), d(np), v(np, np), b(nmax), z(nmax)

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