← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/invit.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/invit.f
 
1
*DECK INVIT
 
2
      SUBROUTINE INVIT (NM, N, A, WR, WI, SELECT, MM, M, Z, IERR, RM1,
 
3
     +   RV1, RV2)
 
4
C***BEGIN PROLOGUE  INVIT
 
5
C***PURPOSE  Compute the eigenvectors of a real upper Hessenberg
 
6
C            matrix associated with specified eigenvalues by inverse
 
7
C            iteration.
 
8
C***LIBRARY   SLATEC (EISPACK)
 
9
C***CATEGORY  D4C2B
 
10
C***TYPE      SINGLE PRECISION (INVIT-S, CINVIT-C)
 
11
C***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
 
12
C***AUTHOR  Smith, B. T., et al.
 
13
C***DESCRIPTION
 
14
C
 
15
C     This subroutine is a translation of the ALGOL procedure INVIT
 
16
C     by Peters and Wilkinson.
 
17
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
 
18
C
 
19
C     This subroutine finds those eigenvectors of a REAL UPPER
 
20
C     Hessenberg matrix corresponding to specified eigenvalues,
 
21
C     using inverse iteration.
 
22
C
 
23
C     On INPUT
 
24
C
 
25
C        NM must be set to the row dimension of the two-dimensional
 
26
C          array parameters, A and Z, as declared in the calling
 
27
C          program dimension statement.  NM is an INTEGER variable.
 
28
C
 
29
C        N is the order of the matrix A.  N is an INTEGER variable.
 
30
C          N must be less than or equal to NM.
 
31
C
 
32
C        A contains the upper Hessenberg matrix.  A is a two-dimensional
 
33
C          REAL array, dimensioned A(NM,N).
 
34
C
 
35
C        WR and WI contain the real and imaginary parts, respectively,
 
36
C          of the eigenvalues of the Hessenberg matrix.  The eigenvalues
 
37
C          must be stored in a manner identical to that output by
 
38
C          subroutine  HQR,  which recognizes possible splitting of the
 
39
C          matrix.  WR and WI are one-dimensional REAL arrays,
 
40
C          dimensioned WR(N) and WI(N).
 
41
C
 
42
C        SELECT specifies the eigenvectors to be found. The
 
43
C          eigenvector corresponding to the J-th eigenvalue is
 
44
C          specified by setting SELECT(J) to .TRUE.  SELECT is a
 
45
C          one-dimensional LOGICAL array, dimensioned SELECT(N).
 
46
C
 
47
C        MM should be set to an upper bound for the number of
 
48
C          columns required to store the eigenvectors to be found.
 
49
C          NOTE that two columns are required to store the
 
50
C          eigenvector corresponding to a complex eigenvalue.  One
 
51
C          column is required to store the eigenvector corresponding
 
52
C          to a real eigenvalue.  MM is an INTEGER variable.
 
53
C
 
54
C     On OUTPUT
 
55
C
 
56
C        A and WI are unaltered.
 
57
C
 
58
C        WR may have been altered since close eigenvalues are perturbed
 
59
C          slightly in searching for independent eigenvectors.
 
60
C
 
61
C        SELECT may have been altered.  If the elements corresponding
 
62
C          to a pair of conjugate complex eigenvalues were each
 
63
C          initially set to .TRUE., the program resets the second of
 
64
C          the two elements to .FALSE.
 
65
C
 
66
C        M is the number of columns actually used to store the
 
67
C          eigenvectors.  M is an INTEGER variable.
 
68
C
 
69
C        Z contains the real and imaginary parts of the eigenvectors.
 
70
C          The eigenvectors are packed into the columns of Z starting
 
71
C          at the first column.  If the next selected eigenvalue is
 
72
C          real, the next column of Z contains its eigenvector.  If the
 
73
C          eigenvalue is complex, the next two columns of Z contain the
 
74
C          real and imaginary parts of its eigenvector, with the real
 
75
C          part first.  The eigenvectors are normalized so that the
 
76
C          component of largest magnitude is 1. Any vector which fails
 
77
C          the acceptance test is set to zero.  Z is a two-dimensional
 
78
C          REAL array, dimensioned Z(NM,MM).
 
79
C
 
80
C        IERR is an INTEGER flag set to
 
81
C          Zero       for normal return,
 
82
C          -(2*N+1)   if more than MM columns of Z are necessary
 
83
C                     to store the eigenvectors corresponding to
 
84
C                     the specified eigenvalues (in this case, M is
 
85
C                     equal to the number of columns of Z containing
 
86
C                     eigenvectors already computed),
 
87
C          -K         if the iteration corresponding to the K-th
 
88
C                     value fails (if this occurs more than once, K
 
89
C                     is the index of the last occurrence); the
 
90
C                     corresponding columns of Z are set to zero
 
91
C                     vectors,
 
92
C          -(N+K)     if both error situations occur.
 
93
C
 
94
C        RM1 is a two-dimensional REAL array used for temporary storage.
 
95
C          This array holds the triangularized form of the upper
 
96
C          Hessenberg matrix used in the inverse iteration process.
 
97
C          RM1 is dimensioned RM1(N,N).
 
98
C
 
99
C        RV1 and RV2 are one-dimensional REAL arrays used for temporary
 
100
C          storage.  They hold the approximate eigenvectors during the
 
101
C          inverse iteration process.  RV1 and RV2 are dimensioned
 
102
C          RV1(N) and RV2(N).
 
103
C
 
104
C     The ALGOL procedure GUESSVEC appears in INVIT in-line.
 
105
C
 
106
C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
 
107
C     Calls CDIV for complex division.
 
108
C
 
109
C     Questions and comments should be directed to B. S. Garbow,
 
110
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
 
111
C     ------------------------------------------------------------------
 
112
C
 
113
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
 
114
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
 
115
C                 system Routines - EISPACK Guide, Springer-Verlag,
 
116
C                 1976.
 
117
C***ROUTINES CALLED  CDIV, PYTHAG
 
118
C***REVISION HISTORY  (YYMMDD)
 
119
C   760101  DATE WRITTEN
 
120
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
121
C   890831  Modified array declarations.  (WRB)
 
122
C   890831  REVISION DATE from Version 3.2
 
123
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
124
C   920501  Reformatted the REFERENCES section.  (WRB)
 
125
C***END PROLOGUE  INVIT
 
126
C
 
127
      INTEGER I,J,K,L,M,N,S,II,IP,MM,MP,NM,NS,N1,UK,IP1,ITS,KM1,IERR
 
128
      REAL A(NM,*),WR(*),WI(*),Z(NM,*)
 
129
      REAL RM1(N,*),RV1(*),RV2(*)
 
130
      REAL T,W,X,Y,EPS3
 
131
      REAL NORM,NORMV,GROWTO,ILAMBD,RLAMBD,UKROOT
 
132
      REAL PYTHAG
 
133
      LOGICAL SELECT(N)
 
134
C
 
135
C***FIRST EXECUTABLE STATEMENT  INVIT
 
136
      IERR = 0
 
137
      UK = 0
 
138
      S = 1
 
139
C     .......... IP = 0, REAL EIGENVALUE
 
140
C                     1, FIRST OF CONJUGATE COMPLEX PAIR
 
141
C                    -1, SECOND OF CONJUGATE COMPLEX PAIR ..........
 
142
      IP = 0
 
143
      N1 = N - 1
 
144
C
 
145
      DO 980 K = 1, N
 
146
         IF (WI(K) .EQ. 0.0E0 .OR. IP .LT. 0) GO TO 100
 
147
         IP = 1
 
148
         IF (SELECT(K) .AND. SELECT(K+1)) SELECT(K+1) = .FALSE.
 
149
  100    IF (.NOT. SELECT(K)) GO TO 960
 
150
         IF (WI(K) .NE. 0.0E0) S = S + 1
 
151
         IF (S .GT. MM) GO TO 1000
 
152
         IF (UK .GE. K) GO TO 200
 
153
C     .......... CHECK FOR POSSIBLE SPLITTING ..........
 
154
         DO 120 UK = K, N
 
155
            IF (UK .EQ. N) GO TO 140
 
156
            IF (A(UK+1,UK) .EQ. 0.0E0) GO TO 140
 
157
  120    CONTINUE
 
158
C     .......... COMPUTE INFINITY NORM OF LEADING UK BY UK
 
159
C                (HESSENBERG) MATRIX ..........
 
160
  140    NORM = 0.0E0
 
161
         MP = 1
 
162
C
 
163
         DO 180 I = 1, UK
 
164
            X = 0.0E0
 
165
C
 
166
            DO 160 J = MP, UK
 
167
  160       X = X + ABS(A(I,J))
 
168
C
 
169
            IF (X .GT. NORM) NORM = X
 
170
            MP = I
 
171
  180    CONTINUE
 
172
C     .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION
 
173
C                AND CLOSE ROOTS ARE MODIFIED BY EPS3 ..........
 
174
         IF (NORM .EQ. 0.0E0) NORM = 1.0E0
 
175
         EPS3 = NORM
 
176
  190    EPS3 = 0.5E0*EPS3
 
177
         IF (NORM + EPS3 .GT. NORM) GO TO 190
 
178
         EPS3 = 2.0E0*EPS3
 
179
C     .......... GROWTO IS THE CRITERION FOR THE GROWTH ..........
 
180
         UKROOT = SQRT(REAL(UK))
 
181
         GROWTO = 0.1E0 / UKROOT
 
182
  200    RLAMBD = WR(K)
 
183
         ILAMBD = WI(K)
 
184
         IF (K .EQ. 1) GO TO 280
 
185
         KM1 = K - 1
 
186
         GO TO 240
 
187
C     .......... PERTURB EIGENVALUE IF IT IS CLOSE
 
188
C                TO ANY PREVIOUS EIGENVALUE ..........
 
189
  220    RLAMBD = RLAMBD + EPS3
 
190
C     .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- ..........
 
191
  240    DO 260 II = 1, KM1
 
192
            I = K - II
 
193
            IF (SELECT(I) .AND. ABS(WR(I)-RLAMBD) .LT. EPS3 .AND.
 
194
     1         ABS(WI(I)-ILAMBD) .LT. EPS3) GO TO 220
 
195
  260    CONTINUE
 
196
C
 
197
         WR(K) = RLAMBD
 
198
C     .......... PERTURB CONJUGATE EIGENVALUE TO MATCH ..........
 
199
         IP1 = K + IP
 
200
         WR(IP1) = RLAMBD
 
201
C     .......... FORM UPPER HESSENBERG A-RLAMBD*I (TRANSPOSED)
 
202
C                AND INITIAL REAL VECTOR ..........
 
203
  280    MP = 1
 
204
C
 
205
         DO 320 I = 1, UK
 
206
C
 
207
            DO 300 J = MP, UK
 
208
  300       RM1(J,I) = A(I,J)
 
209
C
 
210
            RM1(I,I) = RM1(I,I) - RLAMBD
 
211
            MP = I
 
212
            RV1(I) = EPS3
 
213
  320    CONTINUE
 
214
C
 
215
         ITS = 0
 
216
         IF (ILAMBD .NE. 0.0E0) GO TO 520
 
217
C     .......... REAL EIGENVALUE.
 
218
C                TRIANGULAR DECOMPOSITION WITH INTERCHANGES,
 
219
C                REPLACING ZERO PIVOTS BY EPS3 ..........
 
220
         IF (UK .EQ. 1) GO TO 420
 
221
C
 
222
         DO 400 I = 2, UK
 
223
            MP = I - 1
 
224
            IF (ABS(RM1(MP,I)) .LE. ABS(RM1(MP,MP))) GO TO 360
 
225
C
 
226
            DO 340 J = MP, UK
 
227
               Y = RM1(J,I)
 
228
               RM1(J,I) = RM1(J,MP)
 
229
               RM1(J,MP) = Y
 
230
  340       CONTINUE
 
231
C
 
232
  360       IF (RM1(MP,MP) .EQ. 0.0E0) RM1(MP,MP) = EPS3
 
233
            X = RM1(MP,I) / RM1(MP,MP)
 
234
            IF (X .EQ. 0.0E0) GO TO 400
 
235
C
 
236
            DO 380 J = I, UK
 
237
  380       RM1(J,I) = RM1(J,I) - X * RM1(J,MP)
 
238
C
 
239
  400    CONTINUE
 
240
C
 
241
  420    IF (RM1(UK,UK) .EQ. 0.0E0) RM1(UK,UK) = EPS3
 
242
C     .......... BACK SUBSTITUTION FOR REAL VECTOR
 
243
C                FOR I=UK STEP -1 UNTIL 1 DO -- ..........
 
244
  440    DO 500 II = 1, UK
 
245
            I = UK + 1 - II
 
246
            Y = RV1(I)
 
247
            IF (I .EQ. UK) GO TO 480
 
248
            IP1 = I + 1
 
249
C
 
250
            DO 460 J = IP1, UK
 
251
  460       Y = Y - RM1(J,I) * RV1(J)
 
252
C
 
253
  480       RV1(I) = Y / RM1(I,I)
 
254
  500    CONTINUE
 
255
C
 
256
         GO TO 740
 
257
C     .......... COMPLEX EIGENVALUE.
 
258
C                TRIANGULAR DECOMPOSITION WITH INTERCHANGES,
 
259
C                REPLACING ZERO PIVOTS BY EPS3.  STORE IMAGINARY
 
260
C                PARTS IN UPPER TRIANGLE STARTING AT (1,3) ..........
 
261
  520    NS = N - S
 
262
         Z(1,S-1) = -ILAMBD
 
263
         Z(1,S) = 0.0E0
 
264
         IF (N .EQ. 2) GO TO 550
 
265
         RM1(1,3) = -ILAMBD
 
266
         Z(1,S-1) = 0.0E0
 
267
         IF (N .EQ. 3) GO TO 550
 
268
C
 
269
         DO 540 I = 4, N
 
270
  540    RM1(1,I) = 0.0E0
 
271
C
 
272
  550    DO 640 I = 2, UK
 
273
            MP = I - 1
 
274
            W = RM1(MP,I)
 
275
            IF (I .LT. N) T = RM1(MP,I+1)
 
276
            IF (I .EQ. N) T = Z(MP,S-1)
 
277
            X = RM1(MP,MP) * RM1(MP,MP) + T * T
 
278
            IF (W * W .LE. X) GO TO 580
 
279
            X = RM1(MP,MP) / W
 
280
            Y = T / W
 
281
            RM1(MP,MP) = W
 
282
            IF (I .LT. N) RM1(MP,I+1) = 0.0E0
 
283
            IF (I .EQ. N) Z(MP,S-1) = 0.0E0
 
284
C
 
285
            DO 560 J = I, UK
 
286
               W = RM1(J,I)
 
287
               RM1(J,I) = RM1(J,MP) - X * W
 
288
               RM1(J,MP) = W
 
289
               IF (J .LT. N1) GO TO 555
 
290
               L = J - NS
 
291
               Z(I,L) = Z(MP,L) - Y * W
 
292
               Z(MP,L) = 0.0E0
 
293
               GO TO 560
 
294
  555          RM1(I,J+2) = RM1(MP,J+2) - Y * W
 
295
               RM1(MP,J+2) = 0.0E0
 
296
  560       CONTINUE
 
297
C
 
298
            RM1(I,I) = RM1(I,I) - Y * ILAMBD
 
299
            IF (I .LT. N1) GO TO 570
 
300
            L = I - NS
 
301
            Z(MP,L) = -ILAMBD
 
302
            Z(I,L) = Z(I,L) + X * ILAMBD
 
303
            GO TO 640
 
304
  570       RM1(MP,I+2) = -ILAMBD
 
305
            RM1(I,I+2) = RM1(I,I+2) + X * ILAMBD
 
306
            GO TO 640
 
307
  580       IF (X .NE. 0.0E0) GO TO 600
 
308
            RM1(MP,MP) = EPS3
 
309
            IF (I .LT. N) RM1(MP,I+1) = 0.0E0
 
310
            IF (I .EQ. N) Z(MP,S-1) = 0.0E0
 
311
            T = 0.0E0
 
312
            X = EPS3 * EPS3
 
313
  600       W = W / X
 
314
            X = RM1(MP,MP) * W
 
315
            Y = -T * W
 
316
C
 
317
            DO 620 J = I, UK
 
318
               IF (J .LT. N1) GO TO 610
 
319
               L = J - NS
 
320
               T = Z(MP,L)
 
321
               Z(I,L) = -X * T - Y * RM1(J,MP)
 
322
               GO TO 615
 
323
  610          T = RM1(MP,J+2)
 
324
               RM1(I,J+2) = -X * T - Y * RM1(J,MP)
 
325
  615          RM1(J,I) = RM1(J,I) - X * RM1(J,MP) + Y * T
 
326
  620       CONTINUE
 
327
C
 
328
            IF (I .LT. N1) GO TO 630
 
329
            L = I - NS
 
330
            Z(I,L) = Z(I,L) - ILAMBD
 
331
            GO TO 640
 
332
  630       RM1(I,I+2) = RM1(I,I+2) - ILAMBD
 
333
  640    CONTINUE
 
334
C
 
335
         IF (UK .LT. N1) GO TO 650
 
336
         L = UK - NS
 
337
         T = Z(UK,L)
 
338
         GO TO 655
 
339
  650    T = RM1(UK,UK+2)
 
340
  655    IF (RM1(UK,UK) .EQ. 0.0E0 .AND. T .EQ. 0.0E0) RM1(UK,UK) = EPS3
 
341
C     .......... BACK SUBSTITUTION FOR COMPLEX VECTOR
 
342
C                FOR I=UK STEP -1 UNTIL 1 DO -- ..........
 
343
  660    DO 720 II = 1, UK
 
344
            I = UK + 1 - II
 
345
            X = RV1(I)
 
346
            Y = 0.0E0
 
347
            IF (I .EQ. UK) GO TO 700
 
348
            IP1 = I + 1
 
349
C
 
350
            DO 680 J = IP1, UK
 
351
               IF (J .LT. N1) GO TO 670
 
352
               L = J - NS
 
353
               T = Z(I,L)
 
354
               GO TO 675
 
355
  670          T = RM1(I,J+2)
 
356
  675          X = X - RM1(J,I) * RV1(J) + T * RV2(J)
 
357
               Y = Y - RM1(J,I) * RV2(J) - T * RV1(J)
 
358
  680       CONTINUE
 
359
C
 
360
  700       IF (I .LT. N1) GO TO 710
 
361
            L = I - NS
 
362
            T = Z(I,L)
 
363
            GO TO 715
 
364
  710       T = RM1(I,I+2)
 
365
  715       CALL CDIV(X,Y,RM1(I,I),T,RV1(I),RV2(I))
 
366
  720    CONTINUE
 
367
C     .......... ACCEPTANCE TEST FOR REAL OR COMPLEX
 
368
C                EIGENVECTOR AND NORMALIZATION ..........
 
369
  740    ITS = ITS + 1
 
370
         NORM = 0.0E0
 
371
         NORMV = 0.0E0
 
372
C
 
373
         DO 780 I = 1, UK
 
374
            IF (ILAMBD .EQ. 0.0E0) X = ABS(RV1(I))
 
375
            IF (ILAMBD .NE. 0.0E0) X = PYTHAG(RV1(I),RV2(I))
 
376
            IF (NORMV .GE. X) GO TO 760
 
377
            NORMV = X
 
378
            J = I
 
379
  760       NORM = NORM + X
 
380
  780    CONTINUE
 
381
C
 
382
         IF (NORM .LT. GROWTO) GO TO 840
 
383
C     .......... ACCEPT VECTOR ..........
 
384
         X = RV1(J)
 
385
         IF (ILAMBD .EQ. 0.0E0) X = 1.0E0 / X
 
386
         IF (ILAMBD .NE. 0.0E0) Y = RV2(J)
 
387
C
 
388
         DO 820 I = 1, UK
 
389
            IF (ILAMBD .NE. 0.0E0) GO TO 800
 
390
            Z(I,S) = RV1(I) * X
 
391
            GO TO 820
 
392
  800       CALL CDIV(RV1(I),RV2(I),X,Y,Z(I,S-1),Z(I,S))
 
393
  820    CONTINUE
 
394
C
 
395
         IF (UK .EQ. N) GO TO 940
 
396
         J = UK + 1
 
397
         GO TO 900
 
398
C     .......... IN-LINE PROCEDURE FOR CHOOSING
 
399
C                A NEW STARTING VECTOR ..........
 
400
  840    IF (ITS .GE. UK) GO TO 880
 
401
         X = UKROOT
 
402
         Y = EPS3 / (X + 1.0E0)
 
403
         RV1(1) = EPS3
 
404
C
 
405
         DO 860 I = 2, UK
 
406
  860    RV1(I) = Y
 
407
C
 
408
         J = UK - ITS + 1
 
409
         RV1(J) = RV1(J) - EPS3 * X
 
410
         IF (ILAMBD .EQ. 0.0E0) GO TO 440
 
411
         GO TO 660
 
412
C     .......... SET ERROR -- UNACCEPTED EIGENVECTOR ..........
 
413
  880    J = 1
 
414
         IERR = -K
 
415
C     .......... SET REMAINING VECTOR COMPONENTS TO ZERO ..........
 
416
  900    DO 920 I = J, N
 
417
            Z(I,S) = 0.0E0
 
418
            IF (ILAMBD .NE. 0.0E0) Z(I,S-1) = 0.0E0
 
419
  920    CONTINUE
 
420
C
 
421
  940    S = S + 1
 
422
  960    IF (IP .EQ. (-1)) IP = 0
 
423
         IF (IP .EQ. 1) IP = -1
 
424
  980 CONTINUE
 
425
C
 
426
      GO TO 1001
 
427
C     .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR
 
428
C                SPACE REQUIRED ..........
 
429
 1000 IF (IERR .NE. 0) IERR = IERR - N
 
430
      IF (IERR .EQ. 0) IERR = -(2 * N + 1)
 
431
 1001 M = S - 1 - ABS(IP)
 
432
      RETURN
 
433
      END