Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pycmqr.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pycmqr.f
1 
2 C*********************************************************************
3 
4 C...PYCMQR
5 C...Auxiliary to PYEICG.
6 C
7 C THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
8 C ALGOL PROCEDURE COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN
9 C AND WILKINSON.
10 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(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 OF A COMPLEX
15 C UPPER HESSENBERG MATRIX BY THE QR METHOD.
16 C
17 C ON INPUT
18 C
19 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
20 C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
21 C DIMENSION STATEMENT.
22 C
23 C N IS THE ORDER OF THE MATRIX.
24 C
25 C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
26 C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED,
27 C SET LOW=1, IGH=N.
28 C
29 C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
30 C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
31 C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN
32 C INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN
33 C THE REDUCTION BY CORTH, IF PERFORMED.
34 C
35 C ON OUTPUT
36 C
37 C THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
38 C DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE
39 C CALLING COMQR IF SUBSEQUENT CALCULATION OF
40 C EIGENVECTORS IS TO BE PERFORMED.
41 C
42 C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
43 C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR
44 C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
45 C FOR INDICES IERR+1,...,N.
46 C
47 C IERR IS SET TO
48 C ZERO FOR NORMAL RETURN,
49 C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
50 C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
51 C
52 C CALLS PYCDIV FOR COMPLEX DIVISION.
53 C CALLS PYCSRT FOR COMPLEX SQUARE ROOT.
54 C CALLS PYTHAG FOR DSQRT(A*A + B*B) .
55 C
56 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
57 C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
58 C
59 C THIS VERSION DATED AUGUST 1983.
60 C
61 
62  SUBROUTINE pycmqr(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
63 
64  INTEGER i,j,l,n,en,ll,nm,igh,itn,its,low,lp1,enm1,ierr
65  DOUBLE PRECISION hr(4,4),hi(4,4),wr(4),wi(4)
66  DOUBLE PRECISION si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,norm,tst1,tst2,
67  x pythag
68 
69  ierr = 0
70  IF (low .EQ. igh) goto 130
71 C .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
72  l = low + 1
73 C
74  DO 120 i = l, igh
75  ll = min0(i+1,igh)
76  IF (hi(i,i-1) .EQ. 0.0d0) goto 120
77  norm = pythag(hr(i,i-1),hi(i,i-1))
78  yr = hr(i,i-1) / norm
79  yi = hi(i,i-1) / norm
80  hr(i,i-1) = norm
81  hi(i,i-1) = 0.0d0
82 C
83  DO 100 j = i, igh
84  si = yr * hi(i,j) - yi * hr(i,j)
85  hr(i,j) = yr * hr(i,j) + yi * hi(i,j)
86  hi(i,j) = si
87  100 CONTINUE
88 C
89  DO 110 j = low, ll
90  si = yr * hi(j,i) + yi * hr(j,i)
91  hr(j,i) = yr * hr(j,i) - yi * hi(j,i)
92  hi(j,i) = si
93  110 CONTINUE
94 C
95  120 CONTINUE
96 C .......... STORE ROOTS ISOLATED BY CBAL ..........
97  130 DO 140 i = 1, n
98  IF (i .GE. low .AND. i .LE. igh) goto 140
99  wr(i) = hr(i,i)
100  wi(i) = hi(i,i)
101  140 CONTINUE
102 C
103  en = igh
104  tr = 0.0d0
105  ti = 0.0d0
106  itn = 30*n
107 C .......... SEARCH FOR NEXT EIGENVALUE ..........
108  150 IF (en .LT. low) goto 320
109  its = 0
110  enm1 = en - 1
111 C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
112 C FOR L=EN STEP -1 UNTIL LOW D0 -- ..........
113  160 DO 170 ll = low, en
114  l = en + low - ll
115  IF (l .EQ. low) goto 180
116  tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1))
117  x + dabs(hr(l,l)) + dabs(hi(l,l))
118  tst2 = tst1 + dabs(hr(l,l-1))
119  IF (tst2 .EQ. tst1) goto 180
120  170 CONTINUE
121 C .......... FORM SHIFT ..........
122  180 IF (l .EQ. en) goto 300
123  IF (itn .EQ. 0) goto 310
124  IF (its .EQ. 10 .OR. its .EQ. 20) goto 200
125  sr = hr(en,en)
126  si = hi(en,en)
127  xr = hr(enm1,en) * hr(en,enm1)
128  xi = hi(enm1,en) * hr(en,enm1)
129  IF (xr .EQ. 0.0d0 .AND. xi .EQ. 0.0d0) goto 210
130  yr = (hr(enm1,enm1) - sr) / 2.0d0
131  yi = (hi(enm1,enm1) - si) / 2.0d0
132  CALL pycsrt(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi)
133  IF (yr * zzr + yi * zzi .GE. 0.0d0) goto 190
134  zzr = -zzr
135  zzi = -zzi
136  190 CALL pycdiv(xr,xi,yr+zzr,yi+zzi,xr,xi)
137  sr = sr - xr
138  si = si - xi
139  goto 210
140 C .......... FORM EXCEPTIONAL SHIFT ..........
141  200 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2))
142  si = 0.0d0
143 C
144  210 DO 220 i = low, en
145  hr(i,i) = hr(i,i) - sr
146  hi(i,i) = hi(i,i) - si
147  220 CONTINUE
148 C
149  tr = tr + sr
150  ti = ti + si
151  its = its + 1
152  itn = itn - 1
153 C .......... REDUCE TO TRIANGLE (ROWS) ..........
154  lp1 = l + 1
155 C
156  DO 240 i = lp1, en
157  sr = hr(i,i-1)
158  hr(i,i-1) = 0.0d0
159  norm = pythag(pythag(hr(i-1,i-1),hi(i-1,i-1)),sr)
160  xr = hr(i-1,i-1) / norm
161  wr(i-1) = xr
162  xi = hi(i-1,i-1) / norm
163  wi(i-1) = xi
164  hr(i-1,i-1) = norm
165  hi(i-1,i-1) = 0.0d0
166  hi(i,i-1) = sr / norm
167 C
168  DO 230 j = i, en
169  yr = hr(i-1,j)
170  yi = hi(i-1,j)
171  zzr = hr(i,j)
172  zzi = hi(i,j)
173  hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr
174  hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi
175  hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr
176  hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi
177  230 CONTINUE
178 C
179  240 CONTINUE
180 C
181  si = hi(en,en)
182  IF (si .EQ. 0.0d0) goto 250
183  norm = pythag(hr(en,en),si)
184  sr = hr(en,en) / norm
185  si = si / norm
186  hr(en,en) = norm
187  hi(en,en) = 0.0d0
188 C .......... INVERSE OPERATION (COLUMNS) ..........
189  250 DO 280 j = lp1, en
190  xr = wr(j-1)
191  xi = wi(j-1)
192 C
193  DO 270 i = l, j
194  yr = hr(i,j-1)
195  yi = 0.0d0
196  zzr = hr(i,j)
197  zzi = hi(i,j)
198  IF (i .EQ. j) goto 260
199  yi = hi(i,j-1)
200  hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi
201  260 hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr
202  hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr
203  hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi
204  270 CONTINUE
205 C
206  280 CONTINUE
207 C
208  IF (si .EQ. 0.0d0) goto 160
209 C
210  DO 290 i = l, en
211  yr = hr(i,en)
212  yi = hi(i,en)
213  hr(i,en) = sr * yr - si * yi
214  hi(i,en) = sr * yi + si * yr
215  290 CONTINUE
216 C
217  goto 160
218 C .......... A ROOT FOUND ..........
219  300 wr(en) = hr(en,en) + tr
220  wi(en) = hi(en,en) + ti
221  en = enm1
222  goto 150
223 C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
224 C CONVERGED AFTER 30*N ITERATIONS ..........
225  310 ierr = en
226  320 RETURN
227  END