LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
dhsein.f
Go to the documentation of this file.
1 *> \brief \b DHSEIN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DHSEIN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhsein.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhsein.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhsein.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
22 * VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL,
23 * IFAILR, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER EIGSRC, INITV, SIDE
27 * INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
28 * ..
29 * .. Array Arguments ..
30 * LOGICAL SELECT( * )
31 * INTEGER IFAILL( * ), IFAILR( * )
32 * DOUBLE PRECISION H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
33 * $ WI( * ), WORK( * ), WR( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> DHSEIN uses inverse iteration to find specified right and/or left
43 *> eigenvectors of a real upper Hessenberg matrix H.
44 *>
45 *> The right eigenvector x and the left eigenvector y of the matrix H
46 *> corresponding to an eigenvalue w are defined by:
47 *>
48 *> H * x = w * x, y**h * H = w * y**h
49 *>
50 *> where y**h denotes the conjugate transpose of the vector y.
51 *> \endverbatim
52 *
53 * Arguments:
54 * ==========
55 *
56 *> \param[in] SIDE
57 *> \verbatim
58 *> SIDE is CHARACTER*1
59 *> = 'R': compute right eigenvectors only;
60 *> = 'L': compute left eigenvectors only;
61 *> = 'B': compute both right and left eigenvectors.
62 *> \endverbatim
63 *>
64 *> \param[in] EIGSRC
65 *> \verbatim
66 *> EIGSRC is CHARACTER*1
67 *> Specifies the source of eigenvalues supplied in (WR,WI):
68 *> = 'Q': the eigenvalues were found using DHSEQR; thus, if
69 *> H has zero subdiagonal elements, and so is
70 *> block-triangular, then the j-th eigenvalue can be
71 *> assumed to be an eigenvalue of the block containing
72 *> the j-th row/column. This property allows DHSEIN to
73 *> perform inverse iteration on just one diagonal block.
74 *> = 'N': no assumptions are made on the correspondence
75 *> between eigenvalues and diagonal blocks. In this
76 *> case, DHSEIN must always perform inverse iteration
77 *> using the whole matrix H.
78 *> \endverbatim
79 *>
80 *> \param[in] INITV
81 *> \verbatim
82 *> INITV is CHARACTER*1
83 *> = 'N': no initial vectors are supplied;
84 *> = 'U': user-supplied initial vectors are stored in the arrays
85 *> VL and/or VR.
86 *> \endverbatim
87 *>
88 *> \param[in,out] SELECT
89 *> \verbatim
90 *> SELECT is LOGICAL array, dimension (N)
91 *> Specifies the eigenvectors to be computed. To select the
92 *> real eigenvector corresponding to a real eigenvalue WR(j),
93 *> SELECT(j) must be set to .TRUE.. To select the complex
94 *> eigenvector corresponding to a complex eigenvalue
95 *> (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)),
96 *> either SELECT(j) or SELECT(j+1) or both must be set to
97 *> .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is
98 *> .FALSE..
99 *> \endverbatim
100 *>
101 *> \param[in] N
102 *> \verbatim
103 *> N is INTEGER
104 *> The order of the matrix H. N >= 0.
105 *> \endverbatim
106 *>
107 *> \param[in] H
108 *> \verbatim
109 *> H is DOUBLE PRECISION array, dimension (LDH,N)
110 *> The upper Hessenberg matrix H.
111 *> If a NaN is detected in H, the routine will return with INFO=-6.
112 *> \endverbatim
113 *>
114 *> \param[in] LDH
115 *> \verbatim
116 *> LDH is INTEGER
117 *> The leading dimension of the array H. LDH >= max(1,N).
118 *> \endverbatim
119 *>
120 *> \param[in,out] WR
121 *> \verbatim
122 *> WR is DOUBLE PRECISION array, dimension (N)
123 *> \endverbatim
124 *>
125 *> \param[in] WI
126 *> \verbatim
127 *> WI is DOUBLE PRECISION array, dimension (N)
128 *>
129 *> On entry, the real and imaginary parts of the eigenvalues of
130 *> H; a complex conjugate pair of eigenvalues must be stored in
131 *> consecutive elements of WR and WI.
132 *> On exit, WR may have been altered since close eigenvalues
133 *> are perturbed slightly in searching for independent
134 *> eigenvectors.
135 *> \endverbatim
136 *>
137 *> \param[in,out] VL
138 *> \verbatim
139 *> VL is DOUBLE PRECISION array, dimension (LDVL,MM)
140 *> On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
141 *> contain starting vectors for the inverse iteration for the
142 *> left eigenvectors; the starting vector for each eigenvector
143 *> must be in the same column(s) in which the eigenvector will
144 *> be stored.
145 *> On exit, if SIDE = 'L' or 'B', the left eigenvectors
146 *> specified by SELECT will be stored consecutively in the
147 *> columns of VL, in the same order as their eigenvalues. A
148 *> complex eigenvector corresponding to a complex eigenvalue is
149 *> stored in two consecutive columns, the first holding the real
150 *> part and the second the imaginary part.
151 *> If SIDE = 'R', VL is not referenced.
152 *> \endverbatim
153 *>
154 *> \param[in] LDVL
155 *> \verbatim
156 *> LDVL is INTEGER
157 *> The leading dimension of the array VL.
158 *> LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
159 *> \endverbatim
160 *>
161 *> \param[in,out] VR
162 *> \verbatim
163 *> VR is DOUBLE PRECISION array, dimension (LDVR,MM)
164 *> On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
165 *> contain starting vectors for the inverse iteration for the
166 *> right eigenvectors; the starting vector for each eigenvector
167 *> must be in the same column(s) in which the eigenvector will
168 *> be stored.
169 *> On exit, if SIDE = 'R' or 'B', the right eigenvectors
170 *> specified by SELECT will be stored consecutively in the
171 *> columns of VR, in the same order as their eigenvalues. A
172 *> complex eigenvector corresponding to a complex eigenvalue is
173 *> stored in two consecutive columns, the first holding the real
174 *> part and the second the imaginary part.
175 *> If SIDE = 'L', VR is not referenced.
176 *> \endverbatim
177 *>
178 *> \param[in] LDVR
179 *> \verbatim
180 *> LDVR is INTEGER
181 *> The leading dimension of the array VR.
182 *> LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
183 *> \endverbatim
184 *>
185 *> \param[in] MM
186 *> \verbatim
187 *> MM is INTEGER
188 *> The number of columns in the arrays VL and/or VR. MM >= M.
189 *> \endverbatim
190 *>
191 *> \param[out] M
192 *> \verbatim
193 *> M is INTEGER
194 *> The number of columns in the arrays VL and/or VR required to
195 *> store the eigenvectors; each selected real eigenvector
196 *> occupies one column and each selected complex eigenvector
197 *> occupies two columns.
198 *> \endverbatim
199 *>
200 *> \param[out] WORK
201 *> \verbatim
202 *> WORK is DOUBLE PRECISION array, dimension ((N+2)*N)
203 *> \endverbatim
204 *>
205 *> \param[out] IFAILL
206 *> \verbatim
207 *> IFAILL is INTEGER array, dimension (MM)
208 *> If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
209 *> eigenvector in the i-th column of VL (corresponding to the
210 *> eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
211 *> eigenvector converged satisfactorily. If the i-th and (i+1)th
212 *> columns of VL hold a complex eigenvector, then IFAILL(i) and
213 *> IFAILL(i+1) are set to the same value.
214 *> If SIDE = 'R', IFAILL is not referenced.
215 *> \endverbatim
216 *>
217 *> \param[out] IFAILR
218 *> \verbatim
219 *> IFAILR is INTEGER array, dimension (MM)
220 *> If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
221 *> eigenvector in the i-th column of VR (corresponding to the
222 *> eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
223 *> eigenvector converged satisfactorily. If the i-th and (i+1)th
224 *> columns of VR hold a complex eigenvector, then IFAILR(i) and
225 *> IFAILR(i+1) are set to the same value.
226 *> If SIDE = 'L', IFAILR is not referenced.
227 *> \endverbatim
228 *>
229 *> \param[out] INFO
230 *> \verbatim
231 *> INFO is INTEGER
232 *> = 0: successful exit
233 *> < 0: if INFO = -i, the i-th argument had an illegal value
234 *> > 0: if INFO = i, i is the number of eigenvectors which
235 *> failed to converge; see IFAILL and IFAILR for further
236 *> details.
237 *> \endverbatim
238 *
239 * Authors:
240 * ========
241 *
242 *> \author Univ. of Tennessee
243 *> \author Univ. of California Berkeley
244 *> \author Univ. of Colorado Denver
245 *> \author NAG Ltd.
246 *
247 *> \ingroup doubleOTHERcomputational
248 *
249 *> \par Further Details:
250 * =====================
251 *>
252 *> \verbatim
253 *>
254 *> Each eigenvector is normalized so that the element of largest
255 *> magnitude has magnitude 1; here the magnitude of a complex number
256 *> (x,y) is taken to be |x|+|y|.
257 *> \endverbatim
258 *>
259 * =====================================================================
260  SUBROUTINE dhsein( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI,
261  $ VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL,
262  $ IFAILR, INFO )
263 *
264 * -- LAPACK computational routine --
265 * -- LAPACK is a software package provided by Univ. of Tennessee, --
266 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
267 *
268 * .. Scalar Arguments ..
269  CHARACTER EIGSRC, INITV, SIDE
270  INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
271 * ..
272 * .. Array Arguments ..
273  LOGICAL SELECT( * )
274  INTEGER IFAILL( * ), IFAILR( * )
275  DOUBLE PRECISION H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
276  $ wi( * ), work( * ), wr( * )
277 * ..
278 *
279 * =====================================================================
280 *
281 * .. Parameters ..
282  DOUBLE PRECISION ZERO, ONE
283  PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
284 * ..
285 * .. Local Scalars ..
286  LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, PAIR, RIGHTV
287  INTEGER I, IINFO, K, KL, KLN, KR, KSI, KSR, LDWORK
288  DOUBLE PRECISION BIGNUM, EPS3, HNORM, SMLNUM, ULP, UNFL, WKI,
289  $ wkr
290 * ..
291 * .. External Functions ..
292  LOGICAL LSAME, DISNAN
293  DOUBLE PRECISION DLAMCH, DLANHS
294  EXTERNAL lsame, dlamch, dlanhs, disnan
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL dlaein, xerbla
298 * ..
299 * .. Intrinsic Functions ..
300  INTRINSIC abs, max
301 * ..
302 * .. Executable Statements ..
303 *
304 * Decode and test the input parameters.
305 *
306  bothv = lsame( side, 'B' )
307  rightv = lsame( side, 'R' ) .OR. bothv
308  leftv = lsame( side, 'L' ) .OR. bothv
309 *
310  fromqr = lsame( eigsrc, 'Q' )
311 *
312  noinit = lsame( initv, 'N' )
313 *
314 * Set M to the number of columns required to store the selected
315 * eigenvectors, and standardize the array SELECT.
316 *
317  m = 0
318  pair = .false.
319  DO 10 k = 1, n
320  IF( pair ) THEN
321  pair = .false.
322  SELECT( k ) = .false.
323  ELSE
324  IF( wi( k ).EQ.zero ) THEN
325  IF( SELECT( k ) )
326  $ m = m + 1
327  ELSE
328  pair = .true.
329  IF( SELECT( k ) .OR. SELECT( k+1 ) ) THEN
330  SELECT( k ) = .true.
331  m = m + 2
332  END IF
333  END IF
334  END IF
335  10 CONTINUE
336 *
337  info = 0
338  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
339  info = -1
340  ELSE IF( .NOT.fromqr .AND. .NOT.lsame( eigsrc, 'N' ) ) THEN
341  info = -2
342  ELSE IF( .NOT.noinit .AND. .NOT.lsame( initv, 'U' ) ) THEN
343  info = -3
344  ELSE IF( n.LT.0 ) THEN
345  info = -5
346  ELSE IF( ldh.LT.max( 1, n ) ) THEN
347  info = -7
348  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
349  info = -11
350  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
351  info = -13
352  ELSE IF( mm.LT.m ) THEN
353  info = -14
354  END IF
355  IF( info.NE.0 ) THEN
356  CALL xerbla( 'DHSEIN', -info )
357  RETURN
358  END IF
359 *
360 * Quick return if possible.
361 *
362  IF( n.EQ.0 )
363  $ RETURN
364 *
365 * Set machine-dependent constants.
366 *
367  unfl = dlamch( 'Safe minimum' )
368  ulp = dlamch( 'Precision' )
369  smlnum = unfl*( n / ulp )
370  bignum = ( one-ulp ) / smlnum
371 *
372  ldwork = n + 1
373 *
374  kl = 1
375  kln = 0
376  IF( fromqr ) THEN
377  kr = 0
378  ELSE
379  kr = n
380  END IF
381  ksr = 1
382 *
383  DO 120 k = 1, n
384  IF( SELECT( k ) ) THEN
385 *
386 * Compute eigenvector(s) corresponding to W(K).
387 *
388  IF( fromqr ) THEN
389 *
390 * If affiliation of eigenvalues is known, check whether
391 * the matrix splits.
392 *
393 * Determine KL and KR such that 1 <= KL <= K <= KR <= N
394 * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
395 * KR = N).
396 *
397 * Then inverse iteration can be performed with the
398 * submatrix H(KL:N,KL:N) for a left eigenvector, and with
399 * the submatrix H(1:KR,1:KR) for a right eigenvector.
400 *
401  DO 20 i = k, kl + 1, -1
402  IF( h( i, i-1 ).EQ.zero )
403  $ GO TO 30
404  20 CONTINUE
405  30 CONTINUE
406  kl = i
407  IF( k.GT.kr ) THEN
408  DO 40 i = k, n - 1
409  IF( h( i+1, i ).EQ.zero )
410  $ GO TO 50
411  40 CONTINUE
412  50 CONTINUE
413  kr = i
414  END IF
415  END IF
416 *
417  IF( kl.NE.kln ) THEN
418  kln = kl
419 *
420 * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
421 * has not ben computed before.
422 *
423  hnorm = dlanhs( 'I', kr-kl+1, h( kl, kl ), ldh, work )
424  IF( disnan( hnorm ) ) THEN
425  info = -6
426  RETURN
427  ELSE IF( hnorm.GT.zero ) THEN
428  eps3 = hnorm*ulp
429  ELSE
430  eps3 = smlnum
431  END IF
432  END IF
433 *
434 * Perturb eigenvalue if it is close to any previous
435 * selected eigenvalues affiliated to the submatrix
436 * H(KL:KR,KL:KR). Close roots are modified by EPS3.
437 *
438  wkr = wr( k )
439  wki = wi( k )
440  60 CONTINUE
441  DO 70 i = k - 1, kl, -1
442  IF( SELECT( i ) .AND. abs( wr( i )-wkr )+
443  $ abs( wi( i )-wki ).LT.eps3 ) THEN
444  wkr = wkr + eps3
445  GO TO 60
446  END IF
447  70 CONTINUE
448  wr( k ) = wkr
449 *
450  pair = wki.NE.zero
451  IF( pair ) THEN
452  ksi = ksr + 1
453  ELSE
454  ksi = ksr
455  END IF
456  IF( leftv ) THEN
457 *
458 * Compute left eigenvector.
459 *
460  CALL dlaein( .false., noinit, n-kl+1, h( kl, kl ), ldh,
461  $ wkr, wki, vl( kl, ksr ), vl( kl, ksi ),
462  $ work, ldwork, work( n*n+n+1 ), eps3, smlnum,
463  $ bignum, iinfo )
464  IF( iinfo.GT.0 ) THEN
465  IF( pair ) THEN
466  info = info + 2
467  ELSE
468  info = info + 1
469  END IF
470  ifaill( ksr ) = k
471  ifaill( ksi ) = k
472  ELSE
473  ifaill( ksr ) = 0
474  ifaill( ksi ) = 0
475  END IF
476  DO 80 i = 1, kl - 1
477  vl( i, ksr ) = zero
478  80 CONTINUE
479  IF( pair ) THEN
480  DO 90 i = 1, kl - 1
481  vl( i, ksi ) = zero
482  90 CONTINUE
483  END IF
484  END IF
485  IF( rightv ) THEN
486 *
487 * Compute right eigenvector.
488 *
489  CALL dlaein( .true., noinit, kr, h, ldh, wkr, wki,
490  $ vr( 1, ksr ), vr( 1, ksi ), work, ldwork,
491  $ work( n*n+n+1 ), eps3, smlnum, bignum,
492  $ iinfo )
493  IF( iinfo.GT.0 ) THEN
494  IF( pair ) THEN
495  info = info + 2
496  ELSE
497  info = info + 1
498  END IF
499  ifailr( ksr ) = k
500  ifailr( ksi ) = k
501  ELSE
502  ifailr( ksr ) = 0
503  ifailr( ksi ) = 0
504  END IF
505  DO 100 i = kr + 1, n
506  vr( i, ksr ) = zero
507  100 CONTINUE
508  IF( pair ) THEN
509  DO 110 i = kr + 1, n
510  vr( i, ksi ) = zero
511  110 CONTINUE
512  END IF
513  END IF
514 *
515  IF( pair ) THEN
516  ksr = ksr + 2
517  ELSE
518  ksr = ksr + 1
519  END IF
520  END IF
521  120 CONTINUE
522 *
523  RETURN
524 *
525 * End of DHSEIN
526 *
527  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlaein(RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO)
DLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
Definition: dlaein.f:172
subroutine dhsein(SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI, VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL, IFAILR, INFO)
DHSEIN
Definition: dhsein.f:263