LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlaein.f
Go to the documentation of this file.
1 *> \brief \b DLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLAEIN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaein.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaein.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaein.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
22 * LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
23 *
24 * .. Scalar Arguments ..
25 * LOGICAL NOINIT, RIGHTV
26 * INTEGER INFO, LDB, LDH, N
27 * DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR
28 * ..
29 * .. Array Arguments ..
30 * DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
31 * $ WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> DLAEIN uses inverse iteration to find a right or left eigenvector
41 *> corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
42 *> matrix H.
43 *> \endverbatim
44 *
45 * Arguments:
46 * ==========
47 *
48 *> \param[in] RIGHTV
49 *> \verbatim
50 *> RIGHTV is LOGICAL
51 *> = .TRUE. : compute right eigenvector;
52 *> = .FALSE.: compute left eigenvector.
53 *> \endverbatim
54 *>
55 *> \param[in] NOINIT
56 *> \verbatim
57 *> NOINIT is LOGICAL
58 *> = .TRUE. : no initial vector supplied in (VR,VI).
59 *> = .FALSE.: initial vector supplied in (VR,VI).
60 *> \endverbatim
61 *>
62 *> \param[in] N
63 *> \verbatim
64 *> N is INTEGER
65 *> The order of the matrix H. N >= 0.
66 *> \endverbatim
67 *>
68 *> \param[in] H
69 *> \verbatim
70 *> H is DOUBLE PRECISION array, dimension (LDH,N)
71 *> The upper Hessenberg matrix H.
72 *> \endverbatim
73 *>
74 *> \param[in] LDH
75 *> \verbatim
76 *> LDH is INTEGER
77 *> The leading dimension of the array H. LDH >= max(1,N).
78 *> \endverbatim
79 *>
80 *> \param[in] WR
81 *> \verbatim
82 *> WR is DOUBLE PRECISION
83 *> \endverbatim
84 *>
85 *> \param[in] WI
86 *> \verbatim
87 *> WI is DOUBLE PRECISION
88 *> The real and imaginary parts of the eigenvalue of H whose
89 *> corresponding right or left eigenvector is to be computed.
90 *> \endverbatim
91 *>
92 *> \param[in,out] VR
93 *> \verbatim
94 *> VR is DOUBLE PRECISION array, dimension (N)
95 *> \endverbatim
96 *>
97 *> \param[in,out] VI
98 *> \verbatim
99 *> VI is DOUBLE PRECISION array, dimension (N)
100 *> On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
101 *> a real starting vector for inverse iteration using the real
102 *> eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
103 *> must contain the real and imaginary parts of a complex
104 *> starting vector for inverse iteration using the complex
105 *> eigenvalue (WR,WI); otherwise VR and VI need not be set.
106 *> On exit, if WI = 0.0 (real eigenvalue), VR contains the
107 *> computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
108 *> VR and VI contain the real and imaginary parts of the
109 *> computed complex eigenvector. The eigenvector is normalized
110 *> so that the component of largest magnitude has magnitude 1;
111 *> here the magnitude of a complex number (x,y) is taken to be
112 *> |x| + |y|.
113 *> VI is not referenced if WI = 0.0.
114 *> \endverbatim
115 *>
116 *> \param[out] B
117 *> \verbatim
118 *> B is DOUBLE PRECISION array, dimension (LDB,N)
119 *> \endverbatim
120 *>
121 *> \param[in] LDB
122 *> \verbatim
123 *> LDB is INTEGER
124 *> The leading dimension of the array B. LDB >= N+1.
125 *> \endverbatim
126 *>
127 *> \param[out] WORK
128 *> \verbatim
129 *> WORK is DOUBLE PRECISION array, dimension (N)
130 *> \endverbatim
131 *>
132 *> \param[in] EPS3
133 *> \verbatim
134 *> EPS3 is DOUBLE PRECISION
135 *> A small machine-dependent value which is used to perturb
136 *> close eigenvalues, and to replace zero pivots.
137 *> \endverbatim
138 *>
139 *> \param[in] SMLNUM
140 *> \verbatim
141 *> SMLNUM is DOUBLE PRECISION
142 *> A machine-dependent value close to the underflow threshold.
143 *> \endverbatim
144 *>
145 *> \param[in] BIGNUM
146 *> \verbatim
147 *> BIGNUM is DOUBLE PRECISION
148 *> A machine-dependent value close to the overflow threshold.
149 *> \endverbatim
150 *>
151 *> \param[out] INFO
152 *> \verbatim
153 *> INFO is INTEGER
154 *> = 0: successful exit
155 *> = 1: inverse iteration did not converge; VR is set to the
156 *> last iterate, and so is VI if WI.ne.0.0.
157 *> \endverbatim
158 *
159 * Authors:
160 * ========
161 *
162 *> \author Univ. of Tennessee
163 *> \author Univ. of California Berkeley
164 *> \author Univ. of Colorado Denver
165 *> \author NAG Ltd.
166 *
167 *> \date September 2012
168 *
169 *> \ingroup doubleOTHERauxiliary
170 *
171 * =====================================================================
172  SUBROUTINE dlaein( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
173  $ ldb, work, eps3, smlnum, bignum, info )
174 *
175 * -- LAPACK auxiliary routine (version 3.4.2) --
176 * -- LAPACK is a software package provided by Univ. of Tennessee, --
177 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178 * September 2012
179 *
180 * .. Scalar Arguments ..
181  LOGICAL noinit, rightv
182  INTEGER info, ldb, ldh, n
183  DOUBLE PRECISION bignum, eps3, smlnum, wi, wr
184 * ..
185 * .. Array Arguments ..
186  DOUBLE PRECISION b( ldb, * ), h( ldh, * ), vi( * ), vr( * ),
187  $ work( * )
188 * ..
189 *
190 * =====================================================================
191 *
192 * .. Parameters ..
193  DOUBLE PRECISION zero, one, tenth
194  parameter( zero = 0.0d+0, one = 1.0d+0, tenth = 1.0d-1 )
195 * ..
196 * .. Local Scalars ..
197  CHARACTER normin, trans
198  INTEGER i, i1, i2, i3, ierr, its, j
199  DOUBLE PRECISION absbii, absbjj, ei, ej, growto, norm, nrmsml,
200  $ rec, rootn, scale, temp, vcrit, vmax, vnorm, w,
201  $ w1, x, xi, xr, y
202 * ..
203 * .. External Functions ..
204  INTEGER idamax
205  DOUBLE PRECISION dasum, dlapy2, dnrm2
206  EXTERNAL idamax, dasum, dlapy2, dnrm2
207 * ..
208 * .. External Subroutines ..
209  EXTERNAL dladiv, dlatrs, dscal
210 * ..
211 * .. Intrinsic Functions ..
212  INTRINSIC abs, dble, max, sqrt
213 * ..
214 * .. Executable Statements ..
215 *
216  info = 0
217 *
218 * GROWTO is the threshold used in the acceptance test for an
219 * eigenvector.
220 *
221  rootn = sqrt( dble( n ) )
222  growto = tenth / rootn
223  nrmsml = max( one, eps3*rootn )*smlnum
224 *
225 * Form B = H - (WR,WI)*I (except that the subdiagonal elements and
226 * the imaginary parts of the diagonal elements are not stored).
227 *
228  DO 20 j = 1, n
229  DO 10 i = 1, j - 1
230  b( i, j ) = h( i, j )
231  10 continue
232  b( j, j ) = h( j, j ) - wr
233  20 continue
234 *
235  IF( wi.EQ.zero ) THEN
236 *
237 * Real eigenvalue.
238 *
239  IF( noinit ) THEN
240 *
241 * Set initial vector.
242 *
243  DO 30 i = 1, n
244  vr( i ) = eps3
245  30 continue
246  ELSE
247 *
248 * Scale supplied initial vector.
249 *
250  vnorm = dnrm2( n, vr, 1 )
251  CALL dscal( n, ( eps3*rootn ) / max( vnorm, nrmsml ), vr,
252  $ 1 )
253  END IF
254 *
255  IF( rightv ) THEN
256 *
257 * LU decomposition with partial pivoting of B, replacing zero
258 * pivots by EPS3.
259 *
260  DO 60 i = 1, n - 1
261  ei = h( i+1, i )
262  IF( abs( b( i, i ) ).LT.abs( ei ) ) THEN
263 *
264 * Interchange rows and eliminate.
265 *
266  x = b( i, i ) / ei
267  b( i, i ) = ei
268  DO 40 j = i + 1, n
269  temp = b( i+1, j )
270  b( i+1, j ) = b( i, j ) - x*temp
271  b( i, j ) = temp
272  40 continue
273  ELSE
274 *
275 * Eliminate without interchange.
276 *
277  IF( b( i, i ).EQ.zero )
278  $ b( i, i ) = eps3
279  x = ei / b( i, i )
280  IF( x.NE.zero ) THEN
281  DO 50 j = i + 1, n
282  b( i+1, j ) = b( i+1, j ) - x*b( i, j )
283  50 continue
284  END IF
285  END IF
286  60 continue
287  IF( b( n, n ).EQ.zero )
288  $ b( n, n ) = eps3
289 *
290  trans = 'N'
291 *
292  ELSE
293 *
294 * UL decomposition with partial pivoting of B, replacing zero
295 * pivots by EPS3.
296 *
297  DO 90 j = n, 2, -1
298  ej = h( j, j-1 )
299  IF( abs( b( j, j ) ).LT.abs( ej ) ) THEN
300 *
301 * Interchange columns and eliminate.
302 *
303  x = b( j, j ) / ej
304  b( j, j ) = ej
305  DO 70 i = 1, j - 1
306  temp = b( i, j-1 )
307  b( i, j-1 ) = b( i, j ) - x*temp
308  b( i, j ) = temp
309  70 continue
310  ELSE
311 *
312 * Eliminate without interchange.
313 *
314  IF( b( j, j ).EQ.zero )
315  $ b( j, j ) = eps3
316  x = ej / b( j, j )
317  IF( x.NE.zero ) THEN
318  DO 80 i = 1, j - 1
319  b( i, j-1 ) = b( i, j-1 ) - x*b( i, j )
320  80 continue
321  END IF
322  END IF
323  90 continue
324  IF( b( 1, 1 ).EQ.zero )
325  $ b( 1, 1 ) = eps3
326 *
327  trans = 'T'
328 *
329  END IF
330 *
331  normin = 'N'
332  DO 110 its = 1, n
333 *
334 * Solve U*x = scale*v for a right eigenvector
335 * or U**T*x = scale*v for a left eigenvector,
336 * overwriting x on v.
337 *
338  CALL dlatrs( 'Upper', trans, 'Nonunit', normin, n, b, ldb,
339  $ vr, scale, work, ierr )
340  normin = 'Y'
341 *
342 * Test for sufficient growth in the norm of v.
343 *
344  vnorm = dasum( n, vr, 1 )
345  IF( vnorm.GE.growto*scale )
346  $ go to 120
347 *
348 * Choose new orthogonal starting vector and try again.
349 *
350  temp = eps3 / ( rootn+one )
351  vr( 1 ) = eps3
352  DO 100 i = 2, n
353  vr( i ) = temp
354  100 continue
355  vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
356  110 continue
357 *
358 * Failure to find eigenvector in N iterations.
359 *
360  info = 1
361 *
362  120 continue
363 *
364 * Normalize eigenvector.
365 *
366  i = idamax( n, vr, 1 )
367  CALL dscal( n, one / abs( vr( i ) ), vr, 1 )
368  ELSE
369 *
370 * Complex eigenvalue.
371 *
372  IF( noinit ) THEN
373 *
374 * Set initial vector.
375 *
376  DO 130 i = 1, n
377  vr( i ) = eps3
378  vi( i ) = zero
379  130 continue
380  ELSE
381 *
382 * Scale supplied initial vector.
383 *
384  norm = dlapy2( dnrm2( n, vr, 1 ), dnrm2( n, vi, 1 ) )
385  rec = ( eps3*rootn ) / max( norm, nrmsml )
386  CALL dscal( n, rec, vr, 1 )
387  CALL dscal( n, rec, vi, 1 )
388  END IF
389 *
390  IF( rightv ) THEN
391 *
392 * LU decomposition with partial pivoting of B, replacing zero
393 * pivots by EPS3.
394 *
395 * The imaginary part of the (i,j)-th element of U is stored in
396 * B(j+1,i).
397 *
398  b( 2, 1 ) = -wi
399  DO 140 i = 2, n
400  b( i+1, 1 ) = zero
401  140 continue
402 *
403  DO 170 i = 1, n - 1
404  absbii = dlapy2( b( i, i ), b( i+1, i ) )
405  ei = h( i+1, i )
406  IF( absbii.LT.abs( ei ) ) THEN
407 *
408 * Interchange rows and eliminate.
409 *
410  xr = b( i, i ) / ei
411  xi = b( i+1, i ) / ei
412  b( i, i ) = ei
413  b( i+1, i ) = zero
414  DO 150 j = i + 1, n
415  temp = b( i+1, j )
416  b( i+1, j ) = b( i, j ) - xr*temp
417  b( j+1, i+1 ) = b( j+1, i ) - xi*temp
418  b( i, j ) = temp
419  b( j+1, i ) = zero
420  150 continue
421  b( i+2, i ) = -wi
422  b( i+1, i+1 ) = b( i+1, i+1 ) - xi*wi
423  b( i+2, i+1 ) = b( i+2, i+1 ) + xr*wi
424  ELSE
425 *
426 * Eliminate without interchanging rows.
427 *
428  IF( absbii.EQ.zero ) THEN
429  b( i, i ) = eps3
430  b( i+1, i ) = zero
431  absbii = eps3
432  END IF
433  ei = ( ei / absbii ) / absbii
434  xr = b( i, i )*ei
435  xi = -b( i+1, i )*ei
436  DO 160 j = i + 1, n
437  b( i+1, j ) = b( i+1, j ) - xr*b( i, j ) +
438  $ xi*b( j+1, i )
439  b( j+1, i+1 ) = -xr*b( j+1, i ) - xi*b( i, j )
440  160 continue
441  b( i+2, i+1 ) = b( i+2, i+1 ) - wi
442  END IF
443 *
444 * Compute 1-norm of offdiagonal elements of i-th row.
445 *
446  work( i ) = dasum( n-i, b( i, i+1 ), ldb ) +
447  $ dasum( n-i, b( i+2, i ), 1 )
448  170 continue
449  IF( b( n, n ).EQ.zero .AND. b( n+1, n ).EQ.zero )
450  $ b( n, n ) = eps3
451  work( n ) = zero
452 *
453  i1 = n
454  i2 = 1
455  i3 = -1
456  ELSE
457 *
458 * UL decomposition with partial pivoting of conjg(B),
459 * replacing zero pivots by EPS3.
460 *
461 * The imaginary part of the (i,j)-th element of U is stored in
462 * B(j+1,i).
463 *
464  b( n+1, n ) = wi
465  DO 180 j = 1, n - 1
466  b( n+1, j ) = zero
467  180 continue
468 *
469  DO 210 j = n, 2, -1
470  ej = h( j, j-1 )
471  absbjj = dlapy2( b( j, j ), b( j+1, j ) )
472  IF( absbjj.LT.abs( ej ) ) THEN
473 *
474 * Interchange columns and eliminate
475 *
476  xr = b( j, j ) / ej
477  xi = b( j+1, j ) / ej
478  b( j, j ) = ej
479  b( j+1, j ) = zero
480  DO 190 i = 1, j - 1
481  temp = b( i, j-1 )
482  b( i, j-1 ) = b( i, j ) - xr*temp
483  b( j, i ) = b( j+1, i ) - xi*temp
484  b( i, j ) = temp
485  b( j+1, i ) = zero
486  190 continue
487  b( j+1, j-1 ) = wi
488  b( j-1, j-1 ) = b( j-1, j-1 ) + xi*wi
489  b( j, j-1 ) = b( j, j-1 ) - xr*wi
490  ELSE
491 *
492 * Eliminate without interchange.
493 *
494  IF( absbjj.EQ.zero ) THEN
495  b( j, j ) = eps3
496  b( j+1, j ) = zero
497  absbjj = eps3
498  END IF
499  ej = ( ej / absbjj ) / absbjj
500  xr = b( j, j )*ej
501  xi = -b( j+1, j )*ej
502  DO 200 i = 1, j - 1
503  b( i, j-1 ) = b( i, j-1 ) - xr*b( i, j ) +
504  $ xi*b( j+1, i )
505  b( j, i ) = -xr*b( j+1, i ) - xi*b( i, j )
506  200 continue
507  b( j, j-1 ) = b( j, j-1 ) + wi
508  END IF
509 *
510 * Compute 1-norm of offdiagonal elements of j-th column.
511 *
512  work( j ) = dasum( j-1, b( 1, j ), 1 ) +
513  $ dasum( j-1, b( j+1, 1 ), ldb )
514  210 continue
515  IF( b( 1, 1 ).EQ.zero .AND. b( 2, 1 ).EQ.zero )
516  $ b( 1, 1 ) = eps3
517  work( 1 ) = zero
518 *
519  i1 = 1
520  i2 = n
521  i3 = 1
522  END IF
523 *
524  DO 270 its = 1, n
525  scale = one
526  vmax = one
527  vcrit = bignum
528 *
529 * Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
530 * or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
531 * overwriting (xr,xi) on (vr,vi).
532 *
533  DO 250 i = i1, i2, i3
534 *
535  IF( work( i ).GT.vcrit ) THEN
536  rec = one / vmax
537  CALL dscal( n, rec, vr, 1 )
538  CALL dscal( n, rec, vi, 1 )
539  scale = scale*rec
540  vmax = one
541  vcrit = bignum
542  END IF
543 *
544  xr = vr( i )
545  xi = vi( i )
546  IF( rightv ) THEN
547  DO 220 j = i + 1, n
548  xr = xr - b( i, j )*vr( j ) + b( j+1, i )*vi( j )
549  xi = xi - b( i, j )*vi( j ) - b( j+1, i )*vr( j )
550  220 continue
551  ELSE
552  DO 230 j = 1, i - 1
553  xr = xr - b( j, i )*vr( j ) + b( i+1, j )*vi( j )
554  xi = xi - b( j, i )*vi( j ) - b( i+1, j )*vr( j )
555  230 continue
556  END IF
557 *
558  w = abs( b( i, i ) ) + abs( b( i+1, i ) )
559  IF( w.GT.smlnum ) THEN
560  IF( w.LT.one ) THEN
561  w1 = abs( xr ) + abs( xi )
562  IF( w1.GT.w*bignum ) THEN
563  rec = one / w1
564  CALL dscal( n, rec, vr, 1 )
565  CALL dscal( n, rec, vi, 1 )
566  xr = vr( i )
567  xi = vi( i )
568  scale = scale*rec
569  vmax = vmax*rec
570  END IF
571  END IF
572 *
573 * Divide by diagonal element of B.
574 *
575  CALL dladiv( xr, xi, b( i, i ), b( i+1, i ), vr( i ),
576  $ vi( i ) )
577  vmax = max( abs( vr( i ) )+abs( vi( i ) ), vmax )
578  vcrit = bignum / vmax
579  ELSE
580  DO 240 j = 1, n
581  vr( j ) = zero
582  vi( j ) = zero
583  240 continue
584  vr( i ) = one
585  vi( i ) = one
586  scale = zero
587  vmax = one
588  vcrit = bignum
589  END IF
590  250 continue
591 *
592 * Test for sufficient growth in the norm of (VR,VI).
593 *
594  vnorm = dasum( n, vr, 1 ) + dasum( n, vi, 1 )
595  IF( vnorm.GE.growto*scale )
596  $ go to 280
597 *
598 * Choose a new orthogonal starting vector and try again.
599 *
600  y = eps3 / ( rootn+one )
601  vr( 1 ) = eps3
602  vi( 1 ) = zero
603 *
604  DO 260 i = 2, n
605  vr( i ) = y
606  vi( i ) = zero
607  260 continue
608  vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
609  270 continue
610 *
611 * Failure to find eigenvector in N iterations
612 *
613  info = 1
614 *
615  280 continue
616 *
617 * Normalize eigenvector.
618 *
619  vnorm = zero
620  DO 290 i = 1, n
621  vnorm = max( vnorm, abs( vr( i ) )+abs( vi( i ) ) )
622  290 continue
623  CALL dscal( n, one / vnorm, vr, 1 )
624  CALL dscal( n, one / vnorm, vi, 1 )
625 *
626  END IF
627 *
628  return
629 *
630 * End of DLAEIN
631 *
632  END