Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pycmq2.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pycmq2.f
1 
2 C*********************************************************************
3 
4 C...PYCMQ2
5 C...Auxiliary to PYEICG.
6 C
7 C THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
8 C ALGOL PROCEDURE COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS
9 C AND WILKINSON.
10 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
11 C THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
12 C (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
13 C
14 C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
15 C OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR
16 C METHOD. THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
17 C CAN ALSO BE FOUND IF CORTH HAS BEEN USED TO REDUCE
18 C THIS GENERAL MATRIX TO HESSENBERG FORM.
19 C
20 C ON INPUT
21 C
22 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
23 C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
24 C DIMENSION STATEMENT.
25 C
26 C N IS THE ORDER OF THE MATRIX.
27 C
28 C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
29 C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED,
30 C SET LOW=1, IGH=N.
31 C
32 C ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
33 C FORMATIONS USED IN THE REDUCTION BY CORTH, IF PERFORMED.
34 C ONLY ELEMENTS LOW THROUGH IGH ARE USED. IF THE EIGENVECTORS
35 C OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND
36 C ORTI(J) TO 0.0D0 FOR THESE ELEMENTS.
37 C
38 C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
39 C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
40 C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER
41 C INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE
42 C REDUCTION BY CORTH, IF PERFORMED. IF THE EIGENVECTORS OF
43 C THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE
44 C ARBITRARY.
45 C
46 C ON OUTPUT
47 C
48 C ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI
49 C HAVE BEEN DESTROYED.
50 C
51 C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
52 C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR
53 C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
54 C FOR INDICES IERR+1,...,N.
55 C
56 C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
57 C RESPECTIVELY, OF THE EIGENVECTORS. THE EIGENVECTORS
58 C ARE UNNORMALIZED. IF AN ERROR EXIT IS MADE, NONE OF
59 C THE EIGENVECTORS HAS BEEN FOUND.
60 C
61 C IERR IS SET TO
62 C ZERO FOR NORMAL RETURN,
63 C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
64 C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
65 C
66 C CALLS PYCDIV FOR COMPLEX DIVISION.
67 C CALLS PYCSRT FOR COMPLEX SQUARE ROOT.
68 C CALLS PYTHAG FOR DSQRT(A*A + B*B) .
69 C
70 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
71 C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
72 C
73 C THIS VERSION DATED OCTOBER 1989.
74 C
75 C MESHED OVERFLOW CONTROL WITH VECTORS OF ISOLATED ROOTS (10/19/89 BSG)
76 C MESHED OVERFLOW CONTROL WITH TRIANGULAR MULTIPLY (10/30/89 BSG)
77 C
78 
79  SUBROUTINE pycmq2(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR)
80 
81  INTEGER i,j,k,l,m,n,en,ii,jj,ll,nm,nn,igh,ip1,
82  x itn,its,low,lp1,enm1,iend,ierr
83  DOUBLE PRECISION hr(4,4),hi(4,4),wr(4),wi(4),zr(4,4),zi(4,4),
84  x ortr(4),orti(4)
85  DOUBLE PRECISION si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,norm,tst1,tst2,
86  x pythag
87 
88  ierr = 0
89 C .......... INITIALIZE EIGENVECTOR MATRIX ..........
90  DO 110 j = 1, n
91 C
92  DO 100 i = 1, n
93  zr(i,j) = 0.0d0
94  zi(i,j) = 0.0d0
95  100 CONTINUE
96  zr(j,j) = 1.0d0
97  110 CONTINUE
98 C .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
99 C FROM THE INFORMATION LEFT BY CORTH ..........
100  iend = igh - low - 1
101  IF (iend.LT.0) goto 220
102  IF (iend.EQ.0) goto 170
103 C .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
104  DO 160 ii = 1, iend
105  i = igh - ii
106  IF (ortr(i) .EQ. 0.0d0 .AND. orti(i) .EQ. 0.0d0) goto 160
107  IF (hr(i,i-1) .EQ. 0.0d0 .AND. hi(i,i-1) .EQ. 0.0d0) goto 160
108 C .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
109  norm = hr(i,i-1) * ortr(i) + hi(i,i-1) * orti(i)
110  ip1 = i + 1
111 C
112  DO 120 k = ip1, igh
113  ortr(k) = hr(k,i-1)
114  orti(k) = hi(k,i-1)
115  120 CONTINUE
116 C
117  DO 150 j = i, igh
118  sr = 0.0d0
119  si = 0.0d0
120 C
121  DO 130 k = i, igh
122  sr = sr + ortr(k) * zr(k,j) + orti(k) * zi(k,j)
123  si = si + ortr(k) * zi(k,j) - orti(k) * zr(k,j)
124  130 CONTINUE
125 C
126  sr = sr / norm
127  si = si / norm
128 C
129  DO 140 k = i, igh
130  zr(k,j) = zr(k,j) + sr * ortr(k) - si * orti(k)
131  zi(k,j) = zi(k,j) + sr * orti(k) + si * ortr(k)
132  140 CONTINUE
133 C
134  150 CONTINUE
135 C
136  160 CONTINUE
137 C .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
138  170 l = low + 1
139 C
140  DO 210 i = l, igh
141  ll = min0(i+1,igh)
142  IF (hi(i,i-1) .EQ. 0.0d0) goto 210
143  norm = pythag(hr(i,i-1),hi(i,i-1))
144  yr = hr(i,i-1) / norm
145  yi = hi(i,i-1) / norm
146  hr(i,i-1) = norm
147  hi(i,i-1) = 0.0d0
148 C
149  DO 180 j = i, n
150  si = yr * hi(i,j) - yi * hr(i,j)
151  hr(i,j) = yr * hr(i,j) + yi * hi(i,j)
152  hi(i,j) = si
153  180 CONTINUE
154 C
155  DO 190 j = 1, ll
156  si = yr * hi(j,i) + yi * hr(j,i)
157  hr(j,i) = yr * hr(j,i) - yi * hi(j,i)
158  hi(j,i) = si
159  190 CONTINUE
160 C
161  DO 200 j = low, igh
162  si = yr * zi(j,i) + yi * zr(j,i)
163  zr(j,i) = yr * zr(j,i) - yi * zi(j,i)
164  zi(j,i) = si
165  200 CONTINUE
166 C
167  210 CONTINUE
168 C .......... STORE ROOTS ISOLATED BY CBAL ..........
169  220 DO 230 i = 1, n
170  IF (i .GE. low .AND. i .LE. igh) goto 230
171  wr(i) = hr(i,i)
172  wi(i) = hi(i,i)
173  230 CONTINUE
174 C
175  en = igh
176  tr = 0.0d0
177  ti = 0.0d0
178  itn = 30*n
179 C .......... SEARCH FOR NEXT EIGENVALUE ..........
180  240 IF (en .LT. low) goto 430
181  its = 0
182  enm1 = en - 1
183 C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
184 C FOR L=EN STEP -1 UNTIL LOW DO -- ..........
185  250 DO 260 ll = low, en
186  l = en + low - ll
187  IF (l .EQ. low) goto 270
188  tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1))
189  x + dabs(hr(l,l)) + dabs(hi(l,l))
190  tst2 = tst1 + dabs(hr(l,l-1))
191  IF (tst2 .EQ. tst1) goto 270
192  260 CONTINUE
193 C .......... FORM SHIFT ..........
194  270 IF (l .EQ. en) goto 420
195  IF (itn .EQ. 0) goto 550
196  IF (its .EQ. 10 .OR. its .EQ. 20) goto 290
197  sr = hr(en,en)
198  si = hi(en,en)
199  xr = hr(enm1,en) * hr(en,enm1)
200  xi = hi(enm1,en) * hr(en,enm1)
201  IF (xr .EQ. 0.0d0 .AND. xi .EQ. 0.0d0) goto 300
202  yr = (hr(enm1,enm1) - sr) / 2.0d0
203  yi = (hi(enm1,enm1) - si) / 2.0d0
204  CALL pycsrt(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi)
205  IF (yr * zzr + yi * zzi .GE. 0.0d0) goto 280
206  zzr = -zzr
207  zzi = -zzi
208  280 CALL pycdiv(xr,xi,yr+zzr,yi+zzi,xr,xi)
209  sr = sr - xr
210  si = si - xi
211  goto 300
212 C .......... FORM EXCEPTIONAL SHIFT ..........
213  290 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2))
214  si = 0.0d0
215 C
216  300 DO 310 i = low, en
217  hr(i,i) = hr(i,i) - sr
218  hi(i,i) = hi(i,i) - si
219  310 CONTINUE
220 C
221  tr = tr + sr
222  ti = ti + si
223  its = its + 1
224  itn = itn - 1
225 C .......... REDUCE TO TRIANGLE (ROWS) ..........
226  lp1 = l + 1
227 C
228  DO 330 i = lp1, en
229  sr = hr(i,i-1)
230  hr(i,i-1) = 0.0d0
231  norm = pythag(pythag(hr(i-1,i-1),hi(i-1,i-1)),sr)
232  xr = hr(i-1,i-1) / norm
233  wr(i-1) = xr
234  xi = hi(i-1,i-1) / norm
235  wi(i-1) = xi
236  hr(i-1,i-1) = norm
237  hi(i-1,i-1) = 0.0d0
238  hi(i,i-1) = sr / norm
239 C
240  DO 320 j = i, n
241  yr = hr(i-1,j)
242  yi = hi(i-1,j)
243  zzr = hr(i,j)
244  zzi = hi(i,j)
245  hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr
246  hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi
247  hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr
248  hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi
249  320 CONTINUE
250 C
251  330 CONTINUE
252 C
253  si = hi(en,en)
254  IF (si .EQ. 0.0d0) goto 350
255  norm = pythag(hr(en,en),si)
256  sr = hr(en,en) / norm
257  si = si / norm
258  hr(en,en) = norm
259  hi(en,en) = 0.0d0
260  IF (en .EQ. n) goto 350
261  ip1 = en + 1
262 C
263  DO 340 j = ip1, n
264  yr = hr(en,j)
265  yi = hi(en,j)
266  hr(en,j) = sr * yr + si * yi
267  hi(en,j) = sr * yi - si * yr
268  340 CONTINUE
269 C .......... INVERSE OPERATION (COLUMNS) ..........
270  350 DO 390 j = lp1, en
271  xr = wr(j-1)
272  xi = wi(j-1)
273 C
274  DO 370 i = 1, j
275  yr = hr(i,j-1)
276  yi = 0.0d0
277  zzr = hr(i,j)
278  zzi = hi(i,j)
279  IF (i .EQ. j) goto 360
280  yi = hi(i,j-1)
281  hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi
282  360 hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr
283  hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr
284  hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi
285  370 CONTINUE
286 C
287  DO 380 i = low, igh
288  yr = zr(i,j-1)
289  yi = zi(i,j-1)
290  zzr = zr(i,j)
291  zzi = zi(i,j)
292  zr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr
293  zi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi
294  zr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr
295  zi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi
296  380 CONTINUE
297 C
298  390 CONTINUE
299 C
300  IF (si .EQ. 0.0d0) goto 250
301 C
302  DO 400 i = 1, en
303  yr = hr(i,en)
304  yi = hi(i,en)
305  hr(i,en) = sr * yr - si * yi
306  hi(i,en) = sr * yi + si * yr
307  400 CONTINUE
308 C
309  DO 410 i = low, igh
310  yr = zr(i,en)
311  yi = zi(i,en)
312  zr(i,en) = sr * yr - si * yi
313  zi(i,en) = sr * yi + si * yr
314  410 CONTINUE
315 C
316  goto 250
317 C .......... A ROOT FOUND ..........
318  420 hr(en,en) = hr(en,en) + tr
319  wr(en) = hr(en,en)
320  hi(en,en) = hi(en,en) + ti
321  wi(en) = hi(en,en)
322  en = enm1
323  goto 240
324 C .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND
325 C VECTORS OF UPPER TRIANGULAR FORM ..........
326  430 norm = 0.0d0
327 C
328  DO 440 i = 1, n
329 C
330  DO 440 j = i, n
331  tr = dabs(hr(i,j)) + dabs(hi(i,j))
332  IF (tr .GT. norm) norm = tr
333  440 CONTINUE
334 C
335  IF (n .EQ. 1 .OR. norm .EQ. 0.0d0) goto 560
336 C .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........
337  DO 500 nn = 2, n
338  en = n + 2 - nn
339  xr = wr(en)
340  xi = wi(en)
341  hr(en,en) = 1.0d0
342  hi(en,en) = 0.0d0
343  enm1 = en - 1
344 C .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
345  DO 490 ii = 1, enm1
346  i = en - ii
347  zzr = 0.0d0
348  zzi = 0.0d0
349  ip1 = i + 1
350 C
351  DO 450 j = ip1, en
352  zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en)
353  zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en)
354  450 CONTINUE
355 C
356  yr = xr - wr(i)
357  yi = xi - wi(i)
358  IF (yr .NE. 0.0d0 .OR. yi .NE. 0.0d0) goto 470
359  tst1 = norm
360  yr = tst1
361  460 yr = 0.01d0 * yr
362  tst2 = norm + yr
363  IF (tst2 .GT. tst1) goto 460
364  470 CONTINUE
365  CALL pycdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en))
366 C .......... OVERFLOW CONTROL ..........
367  tr = dabs(hr(i,en)) + dabs(hi(i,en))
368  IF (tr .EQ. 0.0d0) goto 490
369  tst1 = tr
370  tst2 = tst1 + 1.0d0/tst1
371  IF (tst2 .GT. tst1) goto 490
372  DO 480 j = i, en
373  hr(j,en) = hr(j,en)/tr
374  hi(j,en) = hi(j,en)/tr
375  480 CONTINUE
376 C
377  490 CONTINUE
378 C
379  500 CONTINUE
380 C .......... END BACKSUBSTITUTION ..........
381 C .......... VECTORS OF ISOLATED ROOTS ..........
382  DO 520 i = 1, n
383  IF (i .GE. low .AND. i .LE. igh) goto 520
384 C
385  DO 510 j = i, n
386  zr(i,j) = hr(i,j)
387  zi(i,j) = hi(i,j)
388  510 CONTINUE
389 C
390  520 CONTINUE
391 C .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
392 C VECTORS OF ORIGINAL FULL MATRIX.
393 C FOR J=N STEP -1 UNTIL LOW DO -- ..........
394  DO 540 jj = low, n
395  j = n + low - jj
396  m = min0(j,igh)
397 C
398  DO 540 i = low, igh
399  zzr = 0.0d0
400  zzi = 0.0d0
401 C
402  DO 530 k = low, m
403  zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j)
404  zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j)
405  530 CONTINUE
406 C
407  zr(i,j) = zzr
408  zi(i,j) = zzi
409  540 CONTINUE
410 C
411  goto 560
412 C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
413 C CONVERGED AFTER 30*N ITERATIONS ..........
414  550 ierr = en
415  560 RETURN
416  END