← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/hqr2.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/hqr2.f
 
1
*DECK HQR2
 
2
      SUBROUTINE HQR2 (NM, N, LOW, IGH, H, WR, WI, Z, IERR)
 
3
C***BEGIN PROLOGUE  HQR2
 
4
C***PURPOSE  Compute the eigenvalues and eigenvectors of a real upper
 
5
C            Hessenberg matrix using QR method.
 
6
C***LIBRARY   SLATEC (EISPACK)
 
7
C***CATEGORY  D4C2B
 
8
C***TYPE      SINGLE PRECISION (HQR2-S, COMQR2-C)
 
9
C***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
 
10
C***AUTHOR  Smith, B. T., et al.
 
11
C***DESCRIPTION
 
12
C
 
13
C     This subroutine is a translation of the ALGOL procedure HQR2,
 
14
C     NUM. MATH. 16, 181-204(1970) by Peters and Wilkinson.
 
15
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
 
16
C
 
17
C     This subroutine finds the eigenvalues and eigenvectors
 
18
C     of a REAL UPPER Hessenberg matrix by the QR method.  The
 
19
C     eigenvectors of a REAL GENERAL matrix can also be found
 
20
C     if  ELMHES  and  ELTRAN  or  ORTHES  and  ORTRAN  have
 
21
C     been used to reduce this general matrix to Hessenberg form
 
22
C     and to accumulate the similarity transformations.
 
23
C
 
24
C     On INPUT
 
25
C
 
26
C        NM must be set to the row dimension of the two-dimensional
 
27
C          array parameters, H and Z, as declared in the calling
 
28
C          program dimension statement.  NM is an INTEGER variable.
 
29
C
 
30
C        N is the order of the matrix H.  N is an INTEGER variable.
 
31
C          N must be less than or equal to NM.
 
32
C
 
33
C        LOW and IGH are two INTEGER variables determined by the
 
34
C          balancing subroutine  BALANC.  If  BALANC  has not been
 
35
C          used, set LOW=1 and IGH equal to the order of the matrix, N.
 
36
C
 
37
C        H contains the upper Hessenberg matrix.  H is a two-dimensional
 
38
C          REAL array, dimensioned H(NM,N).
 
39
C
 
40
C        Z contains the transformation matrix produced by  ELTRAN
 
41
C          after the reduction by  ELMHES, or by  ORTRAN  after the
 
42
C          reduction by  ORTHES, if performed.  If the eigenvectors
 
43
C          of the Hessenberg matrix are desired, Z must contain the
 
44
C          identity matrix.  Z is a two-dimensional REAL array,
 
45
C          dimensioned Z(NM,M).
 
46
C
 
47
C     On OUTPUT
 
48
C
 
49
C        H has been destroyed.
 
50
C
 
51
C        WR and WI contain the real and imaginary parts, respectively,
 
52
C          of the eigenvalues.  The eigenvalues are unordered except
 
53
C          that complex conjugate pairs of values appear consecutively
 
54
C          with the eigenvalue having the positive imaginary part first.
 
55
C          If an error exit is made, the eigenvalues should be correct
 
56
C          for indices IERR+1, IERR+2, ..., N.  WR and WI are one-
 
57
C          dimensional REAL arrays, dimensioned WR(N) and WI(N).
 
58
C
 
59
C        Z contains the real and imaginary parts of the eigenvectors.
 
60
C          If the J-th eigenvalue is real, the J-th column of Z
 
61
C          contains its eigenvector.  If the J-th eigenvalue is complex
 
62
C          with positive imaginary part, the J-th and (J+1)-th
 
63
C          columns of Z contain the real and imaginary parts of its
 
64
C          eigenvector.  The eigenvectors are unnormalized.  If an
 
65
C          error exit is made, none of the eigenvectors has been found.
 
66
C
 
67
C        IERR is an INTEGER flag set to
 
68
C          Zero       for normal return,
 
69
C          J          if the J-th eigenvalue has not been
 
70
C                     determined after a total of 30*N iterations.
 
71
C                     The eigenvalues should be correct for indices
 
72
C                     IERR+1, IERR+2, ..., N, but no eigenvectors are
 
73
C                     computed.
 
74
C
 
75
C     Calls CDIV for complex division.
 
76
C
 
77
C     Questions and comments should be directed to B. S. Garbow,
 
78
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
 
79
C     ------------------------------------------------------------------
 
80
C
 
81
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
 
82
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
 
83
C                 system Routines - EISPACK Guide, Springer-Verlag,
 
84
C                 1976.
 
85
C***ROUTINES CALLED  CDIV
 
86
C***REVISION HISTORY  (YYMMDD)
 
87
C   760101  DATE WRITTEN
 
88
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
89
C   890831  Modified array declarations.  (WRB)
 
90
C   890831  REVISION DATE from Version 3.2
 
91
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
92
C   920501  Reformatted the REFERENCES section.  (WRB)
 
93
C***END PROLOGUE  HQR2
 
94
C
 
95
      INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN
 
96
      INTEGER IGH,ITN,ITS,LOW,MP2,ENM2,IERR
 
97
      REAL H(NM,*),WR(*),WI(*),Z(NM,*)
 
98
      REAL P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,S1,S2
 
99
      LOGICAL NOTLAS
 
100
C
 
101
C***FIRST EXECUTABLE STATEMENT  HQR2
 
102
      IERR = 0
 
103
      NORM = 0.0E0
 
104
      K = 1
 
105
C     .......... STORE ROOTS ISOLATED BY BALANC
 
106
C                AND COMPUTE MATRIX NORM ..........
 
107
      DO 50 I = 1, N
 
108
C
 
109
         DO 40 J = K, N
 
110
   40    NORM = NORM + ABS(H(I,J))
 
111
C
 
112
         K = I
 
113
         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
 
114
         WR(I) = H(I,I)
 
115
         WI(I) = 0.0E0
 
116
   50 CONTINUE
 
117
C
 
118
      EN = IGH
 
119
      T = 0.0E0
 
120
      ITN = 30*N
 
121
C     .......... SEARCH FOR NEXT EIGENVALUES ..........
 
122
   60 IF (EN .LT. LOW) GO TO 340
 
123
      ITS = 0
 
124
      NA = EN - 1
 
125
      ENM2 = NA - 1
 
126
C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
 
127
C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
 
128
   70 DO 80 LL = LOW, EN
 
129
         L = EN + LOW - LL
 
130
         IF (L .EQ. LOW) GO TO 100
 
131
         S = ABS(H(L-1,L-1)) + ABS(H(L,L))
 
132
         IF (S .EQ. 0.0E0) S = NORM
 
133
         S2 = S + ABS(H(L,L-1))
 
134
         IF (S2 .EQ. S) GO TO 100
 
135
   80 CONTINUE
 
136
C     .......... FORM SHIFT ..........
 
137
  100 X = H(EN,EN)
 
138
      IF (L .EQ. EN) GO TO 270
 
139
      Y = H(NA,NA)
 
140
      W = H(EN,NA) * H(NA,EN)
 
141
      IF (L .EQ. NA) GO TO 280
 
142
      IF (ITN .EQ. 0) GO TO 1000
 
143
      IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
 
144
C     .......... FORM EXCEPTIONAL SHIFT ..........
 
145
      T = T + X
 
146
C
 
147
      DO 120 I = LOW, EN
 
148
  120 H(I,I) = H(I,I) - X
 
149
C
 
150
      S = ABS(H(EN,NA)) + ABS(H(NA,ENM2))
 
151
      X = 0.75E0 * S
 
152
      Y = X
 
153
      W = -0.4375E0 * S * S
 
154
  130 ITS = ITS + 1
 
155
      ITN = ITN - 1
 
156
C     .......... LOOK FOR TWO CONSECUTIVE SMALL
 
157
C                SUB-DIAGONAL ELEMENTS.
 
158
C                FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
 
159
      DO 140 MM = L, ENM2
 
160
         M = ENM2 + L - MM
 
161
         ZZ = H(M,M)
 
162
         R = X - ZZ
 
163
         S = Y - ZZ
 
164
         P = (R * S - W) / H(M+1,M) + H(M,M+1)
 
165
         Q = H(M+1,M+1) - ZZ - R - S
 
166
         R = H(M+2,M+1)
 
167
         S = ABS(P) + ABS(Q) + ABS(R)
 
168
         P = P / S
 
169
         Q = Q / S
 
170
         R = R / S
 
171
         IF (M .EQ. L) GO TO 150
 
172
         S1 = ABS(P) * (ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1)))
 
173
         S2 = S1 + ABS(H(M,M-1)) * (ABS(Q) + ABS(R))
 
174
         IF (S2 .EQ. S1) GO TO 150
 
175
  140 CONTINUE
 
176
C
 
177
  150 MP2 = M + 2
 
178
C
 
179
      DO 160 I = MP2, EN
 
180
         H(I,I-2) = 0.0E0
 
181
         IF (I .EQ. MP2) GO TO 160
 
182
         H(I,I-3) = 0.0E0
 
183
  160 CONTINUE
 
184
C     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
 
185
C                COLUMNS M TO EN ..........
 
186
      DO 260 K = M, NA
 
187
         NOTLAS = K .NE. NA
 
188
         IF (K .EQ. M) GO TO 170
 
189
         P = H(K,K-1)
 
190
         Q = H(K+1,K-1)
 
191
         R = 0.0E0
 
192
         IF (NOTLAS) R = H(K+2,K-1)
 
193
         X = ABS(P) + ABS(Q) + ABS(R)
 
194
         IF (X .EQ. 0.0E0) GO TO 260
 
195
         P = P / X
 
196
         Q = Q / X
 
197
         R = R / X
 
198
  170    S = SIGN(SQRT(P*P+Q*Q+R*R),P)
 
199
         IF (K .EQ. M) GO TO 180
 
200
         H(K,K-1) = -S * X
 
201
         GO TO 190
 
202
  180    IF (L .NE. M) H(K,K-1) = -H(K,K-1)
 
203
  190    P = P + S
 
204
         X = P / S
 
205
         Y = Q / S
 
206
         ZZ = R / S
 
207
         Q = Q / P
 
208
         R = R / P
 
209
C     .......... ROW MODIFICATION ..........
 
210
         DO 210 J = K, N
 
211
            P = H(K,J) + Q * H(K+1,J)
 
212
            IF (.NOT. NOTLAS) GO TO 200
 
213
            P = P + R * H(K+2,J)
 
214
            H(K+2,J) = H(K+2,J) - P * ZZ
 
215
  200       H(K+1,J) = H(K+1,J) - P * Y
 
216
            H(K,J) = H(K,J) - P * X
 
217
  210    CONTINUE
 
218
C
 
219
         J = MIN(EN,K+3)
 
220
C     .......... COLUMN MODIFICATION ..........
 
221
         DO 230 I = 1, J
 
222
            P = X * H(I,K) + Y * H(I,K+1)
 
223
            IF (.NOT. NOTLAS) GO TO 220
 
224
            P = P + ZZ * H(I,K+2)
 
225
            H(I,K+2) = H(I,K+2) - P * R
 
226
  220       H(I,K+1) = H(I,K+1) - P * Q
 
227
            H(I,K) = H(I,K) - P
 
228
  230    CONTINUE
 
229
C     .......... ACCUMULATE TRANSFORMATIONS ..........
 
230
         DO 250 I = LOW, IGH
 
231
            P = X * Z(I,K) + Y * Z(I,K+1)
 
232
            IF (.NOT. NOTLAS) GO TO 240
 
233
            P = P + ZZ * Z(I,K+2)
 
234
            Z(I,K+2) = Z(I,K+2) - P * R
 
235
  240       Z(I,K+1) = Z(I,K+1) - P * Q
 
236
            Z(I,K) = Z(I,K) - P
 
237
  250    CONTINUE
 
238
C
 
239
  260 CONTINUE
 
240
C
 
241
      GO TO 70
 
242
C     .......... ONE ROOT FOUND ..........
 
243
  270 H(EN,EN) = X + T
 
244
      WR(EN) = H(EN,EN)
 
245
      WI(EN) = 0.0E0
 
246
      EN = NA
 
247
      GO TO 60
 
248
C     .......... TWO ROOTS FOUND ..........
 
249
  280 P = (Y - X) / 2.0E0
 
250
      Q = P * P + W
 
251
      ZZ = SQRT(ABS(Q))
 
252
      H(EN,EN) = X + T
 
253
      X = H(EN,EN)
 
254
      H(NA,NA) = Y + T
 
255
      IF (Q .LT. 0.0E0) GO TO 320
 
256
C     .......... REAL PAIR ..........
 
257
      ZZ = P + SIGN(ZZ,P)
 
258
      WR(NA) = X + ZZ
 
259
      WR(EN) = WR(NA)
 
260
      IF (ZZ .NE. 0.0E0) WR(EN) = X - W / ZZ
 
261
      WI(NA) = 0.0E0
 
262
      WI(EN) = 0.0E0
 
263
      X = H(EN,NA)
 
264
      S = ABS(X) + ABS(ZZ)
 
265
      P = X / S
 
266
      Q = ZZ / S
 
267
      R = SQRT(P*P+Q*Q)
 
268
      P = P / R
 
269
      Q = Q / R
 
270
C     .......... ROW MODIFICATION ..........
 
271
      DO 290 J = NA, N
 
272
         ZZ = H(NA,J)
 
273
         H(NA,J) = Q * ZZ + P * H(EN,J)
 
274
         H(EN,J) = Q * H(EN,J) - P * ZZ
 
275
  290 CONTINUE
 
276
C     .......... COLUMN MODIFICATION ..........
 
277
      DO 300 I = 1, EN
 
278
         ZZ = H(I,NA)
 
279
         H(I,NA) = Q * ZZ + P * H(I,EN)
 
280
         H(I,EN) = Q * H(I,EN) - P * ZZ
 
281
  300 CONTINUE
 
282
C     .......... ACCUMULATE TRANSFORMATIONS ..........
 
283
      DO 310 I = LOW, IGH
 
284
         ZZ = Z(I,NA)
 
285
         Z(I,NA) = Q * ZZ + P * Z(I,EN)
 
286
         Z(I,EN) = Q * Z(I,EN) - P * ZZ
 
287
  310 CONTINUE
 
288
C
 
289
      GO TO 330
 
290
C     .......... COMPLEX PAIR ..........
 
291
  320 WR(NA) = X + P
 
292
      WR(EN) = X + P
 
293
      WI(NA) = ZZ
 
294
      WI(EN) = -ZZ
 
295
  330 EN = ENM2
 
296
      GO TO 60
 
297
C     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND
 
298
C                VECTORS OF UPPER TRIANGULAR FORM ..........
 
299
  340 IF (NORM .EQ. 0.0E0) GO TO 1001
 
300
C     .......... FOR EN=N STEP -1 UNTIL 1 DO -- ..........
 
301
      DO 800 NN = 1, N
 
302
         EN = N + 1 - NN
 
303
         P = WR(EN)
 
304
         Q = WI(EN)
 
305
         NA = EN - 1
 
306
         IF (Q) 710, 600, 800
 
307
C     .......... REAL VECTOR ..........
 
308
  600    M = EN
 
309
         H(EN,EN) = 1.0E0
 
310
         IF (NA .EQ. 0) GO TO 800
 
311
C     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
 
312
         DO 700 II = 1, NA
 
313
            I = EN - II
 
314
            W = H(I,I) - P
 
315
            R = H(I,EN)
 
316
            IF (M .GT. NA) GO TO 620
 
317
C
 
318
            DO 610 J = M, NA
 
319
  610       R = R + H(I,J) * H(J,EN)
 
320
C
 
321
  620       IF (WI(I) .GE. 0.0E0) GO TO 630
 
322
            ZZ = W
 
323
            S = R
 
324
            GO TO 700
 
325
  630       M = I
 
326
            IF (WI(I) .NE. 0.0E0) GO TO 640
 
327
            T = W
 
328
            IF (T .NE. 0.0E0) GO TO 635
 
329
            T = NORM
 
330
  632       T = 0.5E0*T
 
331
            IF (NORM + T .GT. NORM) GO TO 632
 
332
            T = 2.0E0*T
 
333
  635       H(I,EN) = -R / T
 
334
            GO TO 700
 
335
C     .......... SOLVE REAL EQUATIONS ..........
 
336
  640       X = H(I,I+1)
 
337
            Y = H(I+1,I)
 
338
            Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I)
 
339
            T = (X * S - ZZ * R) / Q
 
340
            H(I,EN) = T
 
341
            IF (ABS(X) .LE. ABS(ZZ)) GO TO 650
 
342
            H(I+1,EN) = (-R - W * T) / X
 
343
            GO TO 700
 
344
  650       H(I+1,EN) = (-S - Y * T) / ZZ
 
345
  700    CONTINUE
 
346
C     .......... END REAL VECTOR ..........
 
347
         GO TO 800
 
348
C     .......... COMPLEX VECTOR ..........
 
349
  710    M = NA
 
350
C     .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT
 
351
C                EIGENVECTOR MATRIX IS TRIANGULAR ..........
 
352
         IF (ABS(H(EN,NA)) .LE. ABS(H(NA,EN))) GO TO 720
 
353
         H(NA,NA) = Q / H(EN,NA)
 
354
         H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA)
 
355
         GO TO 730
 
356
  720    CALL CDIV(0.0E0,-H(NA,EN),H(NA,NA)-P,Q,H(NA,NA),H(NA,EN))
 
357
  730    H(EN,NA) = 0.0E0
 
358
         H(EN,EN) = 1.0E0
 
359
         ENM2 = NA - 1
 
360
         IF (ENM2 .EQ. 0) GO TO 800
 
361
C     .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- ..........
 
362
         DO 790 II = 1, ENM2
 
363
            I = NA - II
 
364
            W = H(I,I) - P
 
365
            RA = 0.0E0
 
366
            SA = H(I,EN)
 
367
C
 
368
            DO 760 J = M, NA
 
369
               RA = RA + H(I,J) * H(J,NA)
 
370
               SA = SA + H(I,J) * H(J,EN)
 
371
  760       CONTINUE
 
372
C
 
373
            IF (WI(I) .GE. 0.0E0) GO TO 770
 
374
            ZZ = W
 
375
            R = RA
 
376
            S = SA
 
377
            GO TO 790
 
378
  770       M = I
 
379
            IF (WI(I) .NE. 0.0E0) GO TO 780
 
380
            CALL CDIV(-RA,-SA,W,Q,H(I,NA),H(I,EN))
 
381
            GO TO 790
 
382
C     .......... SOLVE COMPLEX EQUATIONS ..........
 
383
  780       X = H(I,I+1)
 
384
            Y = H(I+1,I)
 
385
            VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q
 
386
            VI = (WR(I) - P) * 2.0E0 * Q
 
387
            IF (VR .NE. 0.0E0 .OR. VI .NE. 0.0E0) GO TO 783
 
388
            S1 = NORM * (ABS(W)+ABS(Q)+ABS(X)+ABS(Y)+ABS(ZZ))
 
389
            VR = S1
 
390
  782       VR = 0.5E0*VR
 
391
            IF (S1 + VR .GT. S1) GO TO 782
 
392
            VR = 2.0E0*VR
 
393
  783       CALL CDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI,
 
394
     1                H(I,NA),H(I,EN))
 
395
            IF (ABS(X) .LE. ABS(ZZ) + ABS(Q)) GO TO 785
 
396
            H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X
 
397
            H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X
 
398
            GO TO 790
 
399
  785       CALL CDIV(-R-Y*H(I,NA),-S-Y*H(I,EN),ZZ,Q,
 
400
     1                H(I+1,NA),H(I+1,EN))
 
401
  790    CONTINUE
 
402
C     .......... END COMPLEX VECTOR ..........
 
403
  800 CONTINUE
 
404
C     .......... END BACK SUBSTITUTION.
 
405
C                VECTORS OF ISOLATED ROOTS ..........
 
406
      DO 840 I = 1, N
 
407
         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
 
408
C
 
409
         DO 820 J = I, N
 
410
  820    Z(I,J) = H(I,J)
 
411
C
 
412
  840 CONTINUE
 
413
C     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
 
414
C                VECTORS OF ORIGINAL FULL MATRIX.
 
415
C                FOR J=N STEP -1 UNTIL LOW DO -- ..........
 
416
      DO 880 JJ = LOW, N
 
417
         J = N + LOW - JJ
 
418
         M = MIN(J,IGH)
 
419
C
 
420
         DO 880 I = LOW, IGH
 
421
            ZZ = 0.0E0
 
422
C
 
423
            DO 860 K = LOW, M
 
424
  860       ZZ = ZZ + Z(I,K) * H(K,J)
 
425
C
 
426
            Z(I,J) = ZZ
 
427
  880 CONTINUE
 
428
C
 
429
      GO TO 1001
 
430
C     .......... SET ERROR -- NO CONVERGENCE TO AN
 
431
C                EIGENVALUE AFTER 30*N ITERATIONS ..........
 
432
 1000 IERR = EN
 
433
 1001 RETURN
 
434
      END