LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
slaein.f
Go to the documentation of this file.
1 *> \brief \b SLAEIN 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 SLAEIN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaein.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaein.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaein.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAEIN( 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 * REAL BIGNUM, EPS3, SMLNUM, WI, WR
28 * ..
29 * .. Array Arguments ..
30 * REAL B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
31 * $ WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> SLAEIN 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 REAL 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 REAL
83 *> \endverbatim
84 *>
85 *> \param[in] WI
86 *> \verbatim
87 *> WI is REAL
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 REAL array, dimension (N)
95 *> \endverbatim
96 *>
97 *> \param[in,out] VI
98 *> \verbatim
99 *> VI is REAL 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 REAL 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 REAL array, dimension (N)
130 *> \endverbatim
131 *>
132 *> \param[in] EPS3
133 *> \verbatim
134 *> EPS3 is REAL
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 REAL
142 *> A machine-dependent value close to the underflow threshold.
143 *> \endverbatim
144 *>
145 *> \param[in] BIGNUM
146 *> \verbatim
147 *> BIGNUM is REAL
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 *> \ingroup realOTHERauxiliary
168 *
169 * =====================================================================
170  SUBROUTINE slaein( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
171  $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
172 *
173 * -- LAPACK auxiliary routine --
174 * -- LAPACK is a software package provided by Univ. of Tennessee, --
175 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176 *
177 * .. Scalar Arguments ..
178  LOGICAL NOINIT, RIGHTV
179  INTEGER INFO, LDB, LDH, N
180  REAL BIGNUM, EPS3, SMLNUM, WI, WR
181 * ..
182 * .. Array Arguments ..
183  REAL B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
184  $ work( * )
185 * ..
186 *
187 * =====================================================================
188 *
189 * .. Parameters ..
190  REAL ZERO, ONE, TENTH
191  parameter( zero = 0.0e+0, one = 1.0e+0, tenth = 1.0e-1 )
192 * ..
193 * .. Local Scalars ..
194  CHARACTER NORMIN, TRANS
195  INTEGER I, I1, I2, I3, IERR, ITS, J
196  REAL ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
197  $ rec, rootn, scale, temp, vcrit, vmax, vnorm, w,
198  $ w1, x, xi, xr, y
199 * ..
200 * .. External Functions ..
201  INTEGER ISAMAX
202  REAL SASUM, SLAPY2, SNRM2
203  EXTERNAL isamax, sasum, slapy2, snrm2
204 * ..
205 * .. External Subroutines ..
206  EXTERNAL sladiv, slatrs, sscal
207 * ..
208 * .. Intrinsic Functions ..
209  INTRINSIC abs, max, real, sqrt
210 * ..
211 * .. Executable Statements ..
212 *
213  info = 0
214 *
215 * GROWTO is the threshold used in the acceptance test for an
216 * eigenvector.
217 *
218  rootn = sqrt( real( n ) )
219  growto = tenth / rootn
220  nrmsml = max( one, eps3*rootn )*smlnum
221 *
222 * Form B = H - (WR,WI)*I (except that the subdiagonal elements and
223 * the imaginary parts of the diagonal elements are not stored).
224 *
225  DO 20 j = 1, n
226  DO 10 i = 1, j - 1
227  b( i, j ) = h( i, j )
228  10 CONTINUE
229  b( j, j ) = h( j, j ) - wr
230  20 CONTINUE
231 *
232  IF( wi.EQ.zero ) THEN
233 *
234 * Real eigenvalue.
235 *
236  IF( noinit ) THEN
237 *
238 * Set initial vector.
239 *
240  DO 30 i = 1, n
241  vr( i ) = eps3
242  30 CONTINUE
243  ELSE
244 *
245 * Scale supplied initial vector.
246 *
247  vnorm = snrm2( n, vr, 1 )
248  CALL sscal( n, ( eps3*rootn ) / max( vnorm, nrmsml ), vr,
249  $ 1 )
250  END IF
251 *
252  IF( rightv ) THEN
253 *
254 * LU decomposition with partial pivoting of B, replacing zero
255 * pivots by EPS3.
256 *
257  DO 60 i = 1, n - 1
258  ei = h( i+1, i )
259  IF( abs( b( i, i ) ).LT.abs( ei ) ) THEN
260 *
261 * Interchange rows and eliminate.
262 *
263  x = b( i, i ) / ei
264  b( i, i ) = ei
265  DO 40 j = i + 1, n
266  temp = b( i+1, j )
267  b( i+1, j ) = b( i, j ) - x*temp
268  b( i, j ) = temp
269  40 CONTINUE
270  ELSE
271 *
272 * Eliminate without interchange.
273 *
274  IF( b( i, i ).EQ.zero )
275  $ b( i, i ) = eps3
276  x = ei / b( i, i )
277  IF( x.NE.zero ) THEN
278  DO 50 j = i + 1, n
279  b( i+1, j ) = b( i+1, j ) - x*b( i, j )
280  50 CONTINUE
281  END IF
282  END IF
283  60 CONTINUE
284  IF( b( n, n ).EQ.zero )
285  $ b( n, n ) = eps3
286 *
287  trans = 'N'
288 *
289  ELSE
290 *
291 * UL decomposition with partial pivoting of B, replacing zero
292 * pivots by EPS3.
293 *
294  DO 90 j = n, 2, -1
295  ej = h( j, j-1 )
296  IF( abs( b( j, j ) ).LT.abs( ej ) ) THEN
297 *
298 * Interchange columns and eliminate.
299 *
300  x = b( j, j ) / ej
301  b( j, j ) = ej
302  DO 70 i = 1, j - 1
303  temp = b( i, j-1 )
304  b( i, j-1 ) = b( i, j ) - x*temp
305  b( i, j ) = temp
306  70 CONTINUE
307  ELSE
308 *
309 * Eliminate without interchange.
310 *
311  IF( b( j, j ).EQ.zero )
312  $ b( j, j ) = eps3
313  x = ej / b( j, j )
314  IF( x.NE.zero ) THEN
315  DO 80 i = 1, j - 1
316  b( i, j-1 ) = b( i, j-1 ) - x*b( i, j )
317  80 CONTINUE
318  END IF
319  END IF
320  90 CONTINUE
321  IF( b( 1, 1 ).EQ.zero )
322  $ b( 1, 1 ) = eps3
323 *
324  trans = 'T'
325 *
326  END IF
327 *
328  normin = 'N'
329  DO 110 its = 1, n
330 *
331 * Solve U*x = scale*v for a right eigenvector
332 * or U**T*x = scale*v for a left eigenvector,
333 * overwriting x on v.
334 *
335  CALL slatrs( 'Upper', trans, 'Nonunit', normin, n, b, ldb,
336  $ vr, scale, work, ierr )
337  normin = 'Y'
338 *
339 * Test for sufficient growth in the norm of v.
340 *
341  vnorm = sasum( n, vr, 1 )
342  IF( vnorm.GE.growto*scale )
343  $ GO TO 120
344 *
345 * Choose new orthogonal starting vector and try again.
346 *
347  temp = eps3 / ( rootn+one )
348  vr( 1 ) = eps3
349  DO 100 i = 2, n
350  vr( i ) = temp
351  100 CONTINUE
352  vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
353  110 CONTINUE
354 *
355 * Failure to find eigenvector in N iterations.
356 *
357  info = 1
358 *
359  120 CONTINUE
360 *
361 * Normalize eigenvector.
362 *
363  i = isamax( n, vr, 1 )
364  CALL sscal( n, one / abs( vr( i ) ), vr, 1 )
365  ELSE
366 *
367 * Complex eigenvalue.
368 *
369  IF( noinit ) THEN
370 *
371 * Set initial vector.
372 *
373  DO 130 i = 1, n
374  vr( i ) = eps3
375  vi( i ) = zero
376  130 CONTINUE
377  ELSE
378 *
379 * Scale supplied initial vector.
380 *
381  norm = slapy2( snrm2( n, vr, 1 ), snrm2( n, vi, 1 ) )
382  rec = ( eps3*rootn ) / max( norm, nrmsml )
383  CALL sscal( n, rec, vr, 1 )
384  CALL sscal( n, rec, vi, 1 )
385  END IF
386 *
387  IF( rightv ) THEN
388 *
389 * LU decomposition with partial pivoting of B, replacing zero
390 * pivots by EPS3.
391 *
392 * The imaginary part of the (i,j)-th element of U is stored in
393 * B(j+1,i).
394 *
395  b( 2, 1 ) = -wi
396  DO 140 i = 2, n
397  b( i+1, 1 ) = zero
398  140 CONTINUE
399 *
400  DO 170 i = 1, n - 1
401  absbii = slapy2( b( i, i ), b( i+1, i ) )
402  ei = h( i+1, i )
403  IF( absbii.LT.abs( ei ) ) THEN
404 *
405 * Interchange rows and eliminate.
406 *
407  xr = b( i, i ) / ei
408  xi = b( i+1, i ) / ei
409  b( i, i ) = ei
410  b( i+1, i ) = zero
411  DO 150 j = i + 1, n
412  temp = b( i+1, j )
413  b( i+1, j ) = b( i, j ) - xr*temp
414  b( j+1, i+1 ) = b( j+1, i ) - xi*temp
415  b( i, j ) = temp
416  b( j+1, i ) = zero
417  150 CONTINUE
418  b( i+2, i ) = -wi
419  b( i+1, i+1 ) = b( i+1, i+1 ) - xi*wi
420  b( i+2, i+1 ) = b( i+2, i+1 ) + xr*wi
421  ELSE
422 *
423 * Eliminate without interchanging rows.
424 *
425  IF( absbii.EQ.zero ) THEN
426  b( i, i ) = eps3
427  b( i+1, i ) = zero
428  absbii = eps3
429  END IF
430  ei = ( ei / absbii ) / absbii
431  xr = b( i, i )*ei
432  xi = -b( i+1, i )*ei
433  DO 160 j = i + 1, n
434  b( i+1, j ) = b( i+1, j ) - xr*b( i, j ) +
435  $ xi*b( j+1, i )
436  b( j+1, i+1 ) = -xr*b( j+1, i ) - xi*b( i, j )
437  160 CONTINUE
438  b( i+2, i+1 ) = b( i+2, i+1 ) - wi
439  END IF
440 *
441 * Compute 1-norm of offdiagonal elements of i-th row.
442 *
443  work( i ) = sasum( n-i, b( i, i+1 ), ldb ) +
444  $ sasum( n-i, b( i+2, i ), 1 )
445  170 CONTINUE
446  IF( b( n, n ).EQ.zero .AND. b( n+1, n ).EQ.zero )
447  $ b( n, n ) = eps3
448  work( n ) = zero
449 *
450  i1 = n
451  i2 = 1
452  i3 = -1
453  ELSE
454 *
455 * UL decomposition with partial pivoting of conjg(B),
456 * replacing zero pivots by EPS3.
457 *
458 * The imaginary part of the (i,j)-th element of U is stored in
459 * B(j+1,i).
460 *
461  b( n+1, n ) = wi
462  DO 180 j = 1, n - 1
463  b( n+1, j ) = zero
464  180 CONTINUE
465 *
466  DO 210 j = n, 2, -1
467  ej = h( j, j-1 )
468  absbjj = slapy2( b( j, j ), b( j+1, j ) )
469  IF( absbjj.LT.abs( ej ) ) THEN
470 *
471 * Interchange columns and eliminate
472 *
473  xr = b( j, j ) / ej
474  xi = b( j+1, j ) / ej
475  b( j, j ) = ej
476  b( j+1, j ) = zero
477  DO 190 i = 1, j - 1
478  temp = b( i, j-1 )
479  b( i, j-1 ) = b( i, j ) - xr*temp
480  b( j, i ) = b( j+1, i ) - xi*temp
481  b( i, j ) = temp
482  b( j+1, i ) = zero
483  190 CONTINUE
484  b( j+1, j-1 ) = wi
485  b( j-1, j-1 ) = b( j-1, j-1 ) + xi*wi
486  b( j, j-1 ) = b( j, j-1 ) - xr*wi
487  ELSE
488 *
489 * Eliminate without interchange.
490 *
491  IF( absbjj.EQ.zero ) THEN
492  b( j, j ) = eps3
493  b( j+1, j ) = zero
494  absbjj = eps3
495  END IF
496  ej = ( ej / absbjj ) / absbjj
497  xr = b( j, j )*ej
498  xi = -b( j+1, j )*ej
499  DO 200 i = 1, j - 1
500  b( i, j-1 ) = b( i, j-1 ) - xr*b( i, j ) +
501  $ xi*b( j+1, i )
502  b( j, i ) = -xr*b( j+1, i ) - xi*b( i, j )
503  200 CONTINUE
504  b( j, j-1 ) = b( j, j-1 ) + wi
505  END IF
506 *
507 * Compute 1-norm of offdiagonal elements of j-th column.
508 *
509  work( j ) = sasum( j-1, b( 1, j ), 1 ) +
510  $ sasum( j-1, b( j+1, 1 ), ldb )
511  210 CONTINUE
512  IF( b( 1, 1 ).EQ.zero .AND. b( 2, 1 ).EQ.zero )
513  $ b( 1, 1 ) = eps3
514  work( 1 ) = zero
515 *
516  i1 = 1
517  i2 = n
518  i3 = 1
519  END IF
520 *
521  DO 270 its = 1, n
522  scale = one
523  vmax = one
524  vcrit = bignum
525 *
526 * Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
527 * or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
528 * overwriting (xr,xi) on (vr,vi).
529 *
530  DO 250 i = i1, i2, i3
531 *
532  IF( work( i ).GT.vcrit ) THEN
533  rec = one / vmax
534  CALL sscal( n, rec, vr, 1 )
535  CALL sscal( n, rec, vi, 1 )
536  scale = scale*rec
537  vmax = one
538  vcrit = bignum
539  END IF
540 *
541  xr = vr( i )
542  xi = vi( i )
543  IF( rightv ) THEN
544  DO 220 j = i + 1, n
545  xr = xr - b( i, j )*vr( j ) + b( j+1, i )*vi( j )
546  xi = xi - b( i, j )*vi( j ) - b( j+1, i )*vr( j )
547  220 CONTINUE
548  ELSE
549  DO 230 j = 1, i - 1
550  xr = xr - b( j, i )*vr( j ) + b( i+1, j )*vi( j )
551  xi = xi - b( j, i )*vi( j ) - b( i+1, j )*vr( j )
552  230 CONTINUE
553  END IF
554 *
555  w = abs( b( i, i ) ) + abs( b( i+1, i ) )
556  IF( w.GT.smlnum ) THEN
557  IF( w.LT.one ) THEN
558  w1 = abs( xr ) + abs( xi )
559  IF( w1.GT.w*bignum ) THEN
560  rec = one / w1
561  CALL sscal( n, rec, vr, 1 )
562  CALL sscal( n, rec, vi, 1 )
563  xr = vr( i )
564  xi = vi( i )
565  scale = scale*rec
566  vmax = vmax*rec
567  END IF
568  END IF
569 *
570 * Divide by diagonal element of B.
571 *
572  CALL sladiv( xr, xi, b( i, i ), b( i+1, i ), vr( i ),
573  $ vi( i ) )
574  vmax = max( abs( vr( i ) )+abs( vi( i ) ), vmax )
575  vcrit = bignum / vmax
576  ELSE
577  DO 240 j = 1, n
578  vr( j ) = zero
579  vi( j ) = zero
580  240 CONTINUE
581  vr( i ) = one
582  vi( i ) = one
583  scale = zero
584  vmax = one
585  vcrit = bignum
586  END IF
587  250 CONTINUE
588 *
589 * Test for sufficient growth in the norm of (VR,VI).
590 *
591  vnorm = sasum( n, vr, 1 ) + sasum( n, vi, 1 )
592  IF( vnorm.GE.growto*scale )
593  $ GO TO 280
594 *
595 * Choose a new orthogonal starting vector and try again.
596 *
597  y = eps3 / ( rootn+one )
598  vr( 1 ) = eps3
599  vi( 1 ) = zero
600 *
601  DO 260 i = 2, n
602  vr( i ) = y
603  vi( i ) = zero
604  260 CONTINUE
605  vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
606  270 CONTINUE
607 *
608 * Failure to find eigenvector in N iterations
609 *
610  info = 1
611 *
612  280 CONTINUE
613 *
614 * Normalize eigenvector.
615 *
616  vnorm = zero
617  DO 290 i = 1, n
618  vnorm = max( vnorm, abs( vr( i ) )+abs( vi( i ) ) )
619  290 CONTINUE
620  CALL sscal( n, one / vnorm, vr, 1 )
621  CALL sscal( n, one / vnorm, vi, 1 )
622 *
623  END IF
624 *
625  RETURN
626 *
627 * End of SLAEIN
628 *
629  END
subroutine slatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
SLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition: slatrs.f:238
subroutine slaein(RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO)
SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
Definition: slaein.f:172
subroutine sladiv(A, B, C, D, P, Q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: sladiv.f:91
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79