30
31
32
33
34#include "implicit_f.inc"
35
36
37
38 INTEGER N, MV
39
41 . a(*), r(*)
42
43
44
45 INTEGER IQ, J, I, IJ, IA, IND, L, M, MQ, LQ, LM, LL, MM, ILQ, IMQ,
46 . IM, , ILR, IMR, JQ, K
47
49 . range, anorm, anrmx, thr, x, y, sinx, sinx2, cosx, cosx2,
50 . sincs
51 range = 1.0e-7
52 IF (mv/=1)THEN
53 iq = -n
54 DO j = 1,n
55 iq = iq+n
56 DO i = 1,n
57 ij = iq+i
58 r(ij) = zero
59 IF (i==j)r(ij) = one
60
61 ENDDO
62 ENDDO
63 ENDIF
64 anorm = zero
65 DO i = 1,n
66 DO j = 1,n
67 IF (i/=j) THEN
68 ia = i+(j*j-j)/2
69 anorm = anorm+a(ia)*a(ia)
70 ENDIF
71 ENDDO
72 ENDDO
73 IF (anorm<=zero) GOTO 300
74 anorm = onep414*sqrt(anorm)
75 anrmx = anorm*range/float(n)
76 ind = 0
77 thr = anorm
78 thr = thr/float(n)
79 l = 1
80 m = l+1
81
82 100 CONTINUE
83
84 mq = (m*m-m)/2
85 lq = (l*l-l)/2
86 lm = l+mq
87 IF (abs(a(lm))<thr) GOTO 200
88 ind = 1
89 ll = l+lq
90 mm = m+mq
91 x = half*(a(ll)-a(mm))
92 y = -a(lm)/sqrt(a(lm)*a(lm)+x*x)
93 IF (x<zero) THEN
94 y = -y
95 ELSEIF((x==zero).AND.(a(lm)<0)) THEN
96 y = 1
97 ELSEIF((x==zero).AND.(a(lm)>0)) THEN
98 y = -1
99 ENDIF
100 sinx = y/sqrt(two*(one+(sqrt(abs(one-y*y)))))
101 sinx2 = sinx*sinx
102 cosx = sqrt(abs(one-sinx2))
103 cosx2 = cosx*cosx
104 sincs = sinx*cosx
105 ilq = n*(l-1)
106 imq = n*(m-1)
107 DO i = 1,n
108 iq = (i*i-i)/2
109 IF (i/=l.AND.i/=m) THEN
110 IF (i<m) THEN
111 im = i+mq
112 ELSE
113 im = m+iq
114 ENDIF
115 IF (i>=l) THEN
116 il = l+iq
117 ELSE
118 il = i+lq
119 ENDIF
120 x = a(il)*cosx-a(im)*sinx
121 a(im) = a(il)*sinx+a(im)*cosx
122 a(il) = x
123 ENDIF
124 IF (mv/=1) THEN
125 ilr = ilq+i
126 imr = imq+i
127 x = r(ilr)*cosx-r(imr)*sinx
128 r(imr) = r(ilr)*sinx+r(imr)*cosx
129 r(ilr) = x
130 ENDIF
131 ENDDO
132 x = two*a(lm)*sincs
133 y = a(ll)*cosx2+a(mm)*sinx2-x
134 x = a(ll)*sinx2+a(mm)*cosx2+x
135 a(lm) = (a(ll)-a(mm))*sincs+a(lm)*(cosx2-sinx2)
136 a(ll) = y
137 a(mm) = x
138
139 200 CONTINUE
140
141 IF (m/=n) THEN
142 m = m+1
143 GO TO 100
144 ENDIF
145 IF (l/=(n-1)) THEN
146 l = l+1
147 m = l+1
148 GO TO 100
149 ENDIF
150 IF (ind==1) THEN
151 ind = 0
152 l = 1
153 m = l+1
154 GO TO 100
155 ENDIF
156 IF (thr>anrmx) THEN
157 thr = thr/float(n)
158 l = 1
159 m = l+1
160 GO TO 100
161 ENDIF
162
163 300 CONTINUE
164
165 iq = -n
166
167 DO i = 1,n
168 iq = iq+n
169 ll = i+(i*i-i)/2
170 jq = n*(i-2)
171 DO j = i,n
172 jq = jq+n
173 mm = j+(j*j-j)/2
174 IF (a(ll)<a(mm)) THEN
175 x = a(ll)
176 a(ll) = a(mm)
177 a(mm) = x
178 IF (mv/=1) THEN
179 DO k = 1,n
180 ilr = iq+k
181 imr = jq+k
182 x = r(ilr)
183 r(ilr) = r(imr)
184 r(imr) = x
185 ENDDO
186 ENDIF
187 ENDIF
188 ENDDO
189 ENDDO
190
191 anorm=sqrt(r(1)*r(1)+r(2)*r(2)+r(3)*r(3))
192 r(1)=r(1)/anorm
193 r(2)=r(2)/anorm
194 r(3)=r(3)/anorm
195 anorm=sqrt(r(4)*r(4)+r(5)*r(5)+r(6)*r(6))
196 r(4)=r(4)/anorm
197 r(5)=r(5)/anorm
198 r(6)=r(6)/anorm
199 r(7)=r(2)*r(6)-r(3)*r(5)
200 r(8)=r(3)*r(4)-r(1)*r(6)
201 r(9)=r(1)*r(5)-r(2)*r(4)
202 anorm=sqrt(r(7)*r(7)+r(8)*r(8)+r(9)*r(9))
203 r(7)=r(7)/anorm
204 r(8)=r(8)/anorm
205 r(9)=r(9)/anorm
206
207 RETURN