← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/cinvit.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/cinvit.f
 
1
*DECK CINVIT
 
2
      SUBROUTINE CINVIT (NM, N, AR, AI, WR, WI, SELECT, MM, M, ZR, ZI,
 
3
     +   IERR, RM1, RM2, RV1, RV2)
 
4
C***BEGIN PROLOGUE  CINVIT
 
5
C***PURPOSE  Compute the eigenvectors of a complex upper Hessenberg
 
6
C            associated with specified eigenvalues using inverse
 
7
C            iteration.
 
8
C***LIBRARY   SLATEC (EISPACK)
 
9
C***CATEGORY  D4C2B
 
10
C***TYPE      COMPLEX (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 CXINVIT
 
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 COMPLEX 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, AR, AI, ZR and ZI, as declared in the
 
27
C          calling program dimension statement.  NM is an INTEGER
 
28
C          variable.
 
29
C
 
30
C        N is the order of the matrix A=(AR,AI).  N is an INTEGER
 
31
C          variable.  N must be less than or equal to NM.
 
32
C
 
33
C        AR and AI contain the real and imaginary parts, respectively,
 
34
C          of the complex upper Hessenberg matrix.  AR and AI are
 
35
C          two-dimensional REAL arrays, dimensioned AR(NM,N)
 
36
C          and AI(NM,N).
 
37
C
 
38
C        WR and WI contain the real and imaginary parts, respectively,
 
39
C          of the eigenvalues of the matrix.  The eigenvalues must be
 
40
C          stored in a manner identical to that of subroutine  COMLR,
 
41
C          which recognizes possible splitting of the matrix.  WR and
 
42
C          WI are one-dimensional REAL arrays, dimensioned WR(N) and
 
43
C          WI(N).
 
44
C
 
45
C        SELECT specifies the eigenvectors to be found.  The
 
46
C          eigenvector corresponding to the J-th eigenvalue is
 
47
C          specified by setting SELECT(J) to .TRUE.  SELECT is a
 
48
C          one-dimensional LOGICAL array, dimensioned SELECT(N).
 
49
C
 
50
C        MM should be set to an upper bound for the number of
 
51
C          eigenvectors to be found.  MM is an INTEGER variable.
 
52
C
 
53
C     On OUTPUT
 
54
C
 
55
C        AR, AI, WI, and SELECT are unaltered.
 
56
C
 
57
C        WR may have been altered since close eigenvalues are perturbed
 
58
C          slightly in searching for independent eigenvectors.
 
59
C
 
60
C        M is the number of eigenvectors actually found.  M is an
 
61
C          INTEGER variable.
 
62
C
 
63
C        ZR and ZI contain the real and imaginary parts, respectively,
 
64
C          of the eigenvectors corresponding to the flagged eigenvalues.
 
65
C          The eigenvectors are normalized so that the component of
 
66
C          largest magnitude is 1.  Any vector which fails the
 
67
C          acceptance test is set to zero.  ZR and ZI are
 
68
C          two-dimensional REAL arrays, dimensioned ZR(NM,MM) and
 
69
C          ZI(NM,MM).
 
70
C
 
71
C        IERR is an INTEGER flag set to
 
72
C          Zero       for normal return,
 
73
C          -(2*N+1)   if more than MM eigenvectors have been requested
 
74
C                     (the MM eigenvectors calculated to this point are
 
75
C                     in ZR and ZI),
 
76
C          -K         if the iteration corresponding to the K-th
 
77
C                     value fails (if this occurs more than once, K
 
78
C                     is the index of the last occurrence); the
 
79
C                     corresponding columns of ZR and ZI are set to
 
80
C                     zero vectors,
 
81
C          -(N+K)     if both error situations occur.
 
82
C
 
83
C        RV1 and RV2 are one-dimensional REAL arrays used for
 
84
C          temporary storage, dimensioned RV1(N) and RV2(N).
 
85
C          They hold the approximate eigenvectors during the inverse
 
86
C          iteration process.
 
87
C
 
88
C        RM1 and RM2 are two-dimensional REAL arrays used for
 
89
C          temporary storage, dimensioned RM1(N,N) and RM2(N,N).
 
90
C          These arrays hold the triangularized form of the upper
 
91
C          Hessenberg matrix used in the inverse iteration process.
 
92
C
 
93
C     The ALGOL procedure GUESSVEC appears in CINVIT in-line.
 
94
C
 
95
C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
 
96
C     Calls CDIV for complex division.
 
97
C
 
98
C     Questions and comments should be directed to B. S. Garbow,
 
99
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
 
100
C     ------------------------------------------------------------------
 
101
C
 
102
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
 
103
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
 
104
C                 system Routines - EISPACK Guide, Springer-Verlag,
 
105
C                 1976.
 
106
C***ROUTINES CALLED  CDIV, PYTHAG
 
107
C***REVISION HISTORY  (YYMMDD)
 
108
C   760101  DATE WRITTEN
 
109
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
110
C   890831  Modified array declarations.  (WRB)
 
111
C   890831  REVISION DATE from Version 3.2
 
112
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
113
C   920501  Reformatted the REFERENCES section.  (WRB)
 
114
C***END PROLOGUE  CINVIT
 
115
C
 
116
      INTEGER I,J,K,M,N,S,II,MM,MP,NM,UK,IP1,ITS,KM1,IERR
 
117
      REAL AR(NM,*),AI(NM,*),WR(*),WI(*),ZR(NM,*),ZI(NM,*)
 
118
      REAL RM1(N,*),RM2(N,*),RV1(*),RV2(*)
 
119
      REAL X,Y,EPS3,NORM,NORMV,GROWTO,ILAMBD,RLAMBD,UKROOT
 
120
      REAL PYTHAG
 
121
      LOGICAL SELECT(N)
 
122
C
 
123
C***FIRST EXECUTABLE STATEMENT  CINVIT
 
124
      IERR = 0
 
125
      UK = 0
 
126
      S = 1
 
127
C
 
128
      DO 980 K = 1, N
 
129
         IF (.NOT. SELECT(K)) GO TO 980
 
130
         IF (S .GT. MM) GO TO 1000
 
131
         IF (UK .GE. K) GO TO 200
 
132
C     .......... CHECK FOR POSSIBLE SPLITTING ..........
 
133
         DO 120 UK = K, N
 
134
            IF (UK .EQ. N) GO TO 140
 
135
            IF (AR(UK+1,UK) .EQ. 0.0E0 .AND. AI(UK+1,UK) .EQ. 0.0E0)
 
136
     1         GO TO 140
 
137
  120    CONTINUE
 
138
C     .......... COMPUTE INFINITY NORM OF LEADING UK BY UK
 
139
C                (HESSENBERG) MATRIX ..........
 
140
  140    NORM = 0.0E0
 
141
         MP = 1
 
142
C
 
143
         DO 180 I = 1, UK
 
144
            X = 0.0E0
 
145
C
 
146
            DO 160 J = MP, UK
 
147
  160       X = X + PYTHAG(AR(I,J),AI(I,J))
 
148
C
 
149
            IF (X .GT. NORM) NORM = X
 
150
            MP = I
 
151
  180    CONTINUE
 
152
C     .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION
 
153
C                AND CLOSE ROOTS ARE MODIFIED BY EPS3 ..........
 
154
         IF (NORM .EQ. 0.0E0) NORM = 1.0E0
 
155
         EPS3 = NORM
 
156
  190    EPS3 = 0.5E0*EPS3
 
157
         IF (NORM + EPS3 .GT. NORM) GO TO 190
 
158
         EPS3 = 2.0E0*EPS3
 
159
C     .......... GROWTO IS THE CRITERION FOR GROWTH ..........
 
160
         UKROOT = SQRT(REAL(UK))
 
161
         GROWTO = 0.1E0 / UKROOT
 
162
  200    RLAMBD = WR(K)
 
163
         ILAMBD = WI(K)
 
164
         IF (K .EQ. 1) GO TO 280
 
165
         KM1 = K - 1
 
166
         GO TO 240
 
167
C     .......... PERTURB EIGENVALUE IF IT IS CLOSE
 
168
C                TO ANY PREVIOUS EIGENVALUE ..........
 
169
  220    RLAMBD = RLAMBD + EPS3
 
170
C     .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- ..........
 
171
  240    DO 260 II = 1, KM1
 
172
            I = K - II
 
173
            IF (SELECT(I) .AND. ABS(WR(I)-RLAMBD) .LT. EPS3 .AND.
 
174
     1         ABS(WI(I)-ILAMBD) .LT. EPS3) GO TO 220
 
175
  260    CONTINUE
 
176
C
 
177
         WR(K) = RLAMBD
 
178
C     .......... FORM UPPER HESSENBERG (AR,AI)-(RLAMBD,ILAMBD)*I
 
179
C                AND INITIAL COMPLEX VECTOR ..........
 
180
  280    MP = 1
 
181
C
 
182
         DO 320 I = 1, UK
 
183
C
 
184
            DO 300 J = MP, UK
 
185
               RM1(I,J) = AR(I,J)
 
186
               RM2(I,J) = AI(I,J)
 
187
  300       CONTINUE
 
188
C
 
189
            RM1(I,I) = RM1(I,I) - RLAMBD
 
190
            RM2(I,I) = RM2(I,I) - ILAMBD
 
191
            MP = I
 
192
            RV1(I) = EPS3
 
193
  320    CONTINUE
 
194
C     .......... TRIANGULAR DECOMPOSITION WITH INTERCHANGES,
 
195
C                REPLACING ZERO PIVOTS BY EPS3 ..........
 
196
         IF (UK .EQ. 1) GO TO 420
 
197
C
 
198
         DO 400 I = 2, UK
 
199
            MP = I - 1
 
200
            IF (PYTHAG(RM1(I,MP),RM2(I,MP)) .LE.
 
201
     1         PYTHAG(RM1(MP,MP),RM2(MP,MP))) GO TO 360
 
202
C
 
203
            DO 340 J = MP, UK
 
204
               Y = RM1(I,J)
 
205
               RM1(I,J) = RM1(MP,J)
 
206
               RM1(MP,J) = Y
 
207
               Y = RM2(I,J)
 
208
               RM2(I,J) = RM2(MP,J)
 
209
               RM2(MP,J) = Y
 
210
  340       CONTINUE
 
211
C
 
212
  360       IF (RM1(MP,MP) .EQ. 0.0E0 .AND. RM2(MP,MP) .EQ. 0.0E0)
 
213
     1         RM1(MP,MP) = EPS3
 
214
            CALL CDIV(RM1(I,MP),RM2(I,MP),RM1(MP,MP),RM2(MP,MP),X,Y)
 
215
            IF (X .EQ. 0.0E0 .AND. Y .EQ. 0.0E0) GO TO 400
 
216
C
 
217
            DO 380 J = I, UK
 
218
               RM1(I,J) = RM1(I,J) - X * RM1(MP,J) + Y * RM2(MP,J)
 
219
               RM2(I,J) = RM2(I,J) - X * RM2(MP,J) - Y * RM1(MP,J)
 
220
  380       CONTINUE
 
221
C
 
222
  400    CONTINUE
 
223
C
 
224
  420    IF (RM1(UK,UK) .EQ. 0.0E0 .AND. RM2(UK,UK) .EQ. 0.0E0)
 
225
     1      RM1(UK,UK) = EPS3
 
226
         ITS = 0
 
227
C     .......... BACK SUBSTITUTION
 
228
C                FOR I=UK STEP -1 UNTIL 1 DO -- ..........
 
229
  660    DO 720 II = 1, UK
 
230
            I = UK + 1 - II
 
231
            X = RV1(I)
 
232
            Y = 0.0E0
 
233
            IF (I .EQ. UK) GO TO 700
 
234
            IP1 = I + 1
 
235
C
 
236
            DO 680 J = IP1, UK
 
237
               X = X - RM1(I,J) * RV1(J) + RM2(I,J) * RV2(J)
 
238
               Y = Y - RM1(I,J) * RV2(J) - RM2(I,J) * RV1(J)
 
239
  680       CONTINUE
 
240
C
 
241
  700       CALL CDIV(X,Y,RM1(I,I),RM2(I,I),RV1(I),RV2(I))
 
242
  720    CONTINUE
 
243
C     .......... ACCEPTANCE TEST FOR EIGENVECTOR
 
244
C                AND NORMALIZATION ..........
 
245
         ITS = ITS + 1
 
246
         NORM = 0.0E0
 
247
         NORMV = 0.0E0
 
248
C
 
249
         DO 780 I = 1, UK
 
250
            X = PYTHAG(RV1(I),RV2(I))
 
251
            IF (NORMV .GE. X) GO TO 760
 
252
            NORMV = X
 
253
            J = I
 
254
  760       NORM = NORM + X
 
255
  780    CONTINUE
 
256
C
 
257
         IF (NORM .LT. GROWTO) GO TO 840
 
258
C     .......... ACCEPT VECTOR ..........
 
259
         X = RV1(J)
 
260
         Y = RV2(J)
 
261
C
 
262
         DO 820 I = 1, UK
 
263
            CALL CDIV(RV1(I),RV2(I),X,Y,ZR(I,S),ZI(I,S))
 
264
  820    CONTINUE
 
265
C
 
266
         IF (UK .EQ. N) GO TO 940
 
267
         J = UK + 1
 
268
         GO TO 900
 
269
C     .......... IN-LINE PROCEDURE FOR CHOOSING
 
270
C                A NEW STARTING VECTOR ..........
 
271
  840    IF (ITS .GE. UK) GO TO 880
 
272
         X = UKROOT
 
273
         Y = EPS3 / (X + 1.0E0)
 
274
         RV1(1) = EPS3
 
275
C
 
276
         DO 860 I = 2, UK
 
277
  860    RV1(I) = Y
 
278
C
 
279
         J = UK - ITS + 1
 
280
         RV1(J) = RV1(J) - EPS3 * X
 
281
         GO TO 660
 
282
C     .......... SET ERROR -- UNACCEPTED EIGENVECTOR ..........
 
283
  880    J = 1
 
284
         IERR = -K
 
285
C     .......... SET REMAINING VECTOR COMPONENTS TO ZERO ..........
 
286
  900    DO 920 I = J, N
 
287
            ZR(I,S) = 0.0E0
 
288
            ZI(I,S) = 0.0E0
 
289
  920    CONTINUE
 
290
C
 
291
  940    S = S + 1
 
292
  980 CONTINUE
 
293
C
 
294
      GO TO 1001
 
295
C     .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR
 
296
C                SPACE REQUIRED ..........
 
297
 1000 IF (IERR .NE. 0) IERR = IERR - N
 
298
      IF (IERR .EQ. 0) IERR = -(2 * N + 1)
 
299
 1001 M = S - 1
 
300
      RETURN
 
301
      END