 LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ slaein()

 subroutine slaein ( logical RIGHTV, logical NOINIT, integer N, real, dimension( ldh, * ) H, integer LDH, real WR, real WI, real, dimension( * ) VR, real, dimension( * ) VI, real, dimension( ldb, * ) B, integer LDB, real, dimension( * ) WORK, real EPS3, real SMLNUM, real BIGNUM, integer INFO )

SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.

Purpose:
``` SLAEIN uses inverse iteration to find a right or left eigenvector
corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
matrix H.```
Parameters
 [in] RIGHTV ``` RIGHTV is LOGICAL = .TRUE. : compute right eigenvector; = .FALSE.: compute left eigenvector.``` [in] NOINIT ``` NOINIT is LOGICAL = .TRUE. : no initial vector supplied in (VR,VI). = .FALSE.: initial vector supplied in (VR,VI).``` [in] N ``` N is INTEGER The order of the matrix H. N >= 0.``` [in] H ``` H is REAL array, dimension (LDH,N) The upper Hessenberg matrix H.``` [in] LDH ``` LDH is INTEGER The leading dimension of the array H. LDH >= max(1,N).``` [in] WR ` WR is REAL` [in] WI ``` WI is REAL The real and imaginary parts of the eigenvalue of H whose corresponding right or left eigenvector is to be computed.``` [in,out] VR ` VR is REAL array, dimension (N)` [in,out] VI ``` VI is REAL array, dimension (N) On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain a real starting vector for inverse iteration using the real eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI must contain the real and imaginary parts of a complex starting vector for inverse iteration using the complex eigenvalue (WR,WI); otherwise VR and VI need not be set. On exit, if WI = 0.0 (real eigenvalue), VR contains the computed real eigenvector; if WI.ne.0.0 (complex eigenvalue), VR and VI contain the real and imaginary parts of the computed complex eigenvector. The eigenvector is normalized so that the component of largest magnitude has magnitude 1; here the magnitude of a complex number (x,y) is taken to be |x| + |y|. VI is not referenced if WI = 0.0.``` [out] B ` B is REAL array, dimension (LDB,N)` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= N+1.``` [out] WORK ` WORK is REAL array, dimension (N)` [in] EPS3 ``` EPS3 is REAL A small machine-dependent value which is used to perturb close eigenvalues, and to replace zero pivots.``` [in] SMLNUM ``` SMLNUM is REAL A machine-dependent value close to the underflow threshold.``` [in] BIGNUM ``` BIGNUM is REAL A machine-dependent value close to the overflow threshold.``` [out] INFO ``` INFO is INTEGER = 0: successful exit = 1: inverse iteration did not converge; VR is set to the last iterate, and so is VI if WI.ne.0.0.```
Date
December 2016

Definition at line 174 of file slaein.f.

174 *
175 * -- LAPACK auxiliary routine (version 3.7.0) --
176 * -- LAPACK is a software package provided by Univ. of Tennessee, --
177 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178 * December 2016
179 *
180 * .. Scalar Arguments ..
181  LOGICAL noinit, rightv
182  INTEGER info, ldb, ldh, n
183  REAL bignum, eps3, smlnum, wi, wr
184 * ..
185 * .. Array Arguments ..
186  REAL b( ldb, * ), h( ldh, * ), vi( * ), vr( * ),
187  \$ work( * )
188 * ..
189 *
190 * =====================================================================
191 *
192 * .. Parameters ..
193  REAL zero, one, tenth
194  parameter( zero = 0.0e+0, one = 1.0e+0, tenth = 1.0e-1 )
195 * ..
196 * .. Local Scalars ..
197  CHARACTER normin, trans
198  INTEGER i, i1, i2, i3, ierr, its, j
199  REAL 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 isamax
205  REAL sasum, slapy2, snrm2
206  EXTERNAL isamax, sasum, slapy2, snrm2
207 * ..
208 * .. External Subroutines ..
210 * ..
211 * .. Intrinsic Functions ..
212  INTRINSIC abs, max, REAL, 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( REAL( 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 = snrm2( n, vr, 1 )
251  CALL sscal( 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 slatrs( '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 = sasum( 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 = isamax( n, vr, 1 )
367  CALL sscal( 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 = slapy2( snrm2( n, vr, 1 ), snrm2( n, vi, 1 ) )
385  rec = ( eps3*rootn ) / max( norm, nrmsml )
386  CALL sscal( n, rec, vr, 1 )
387  CALL sscal( 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 = slapy2( 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 ) = sasum( n-i, b( i, i+1 ), ldb ) +
447  \$ sasum( 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 = slapy2( 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 ) = sasum( j-1, b( 1, j ), 1 ) +
513  \$ sasum( 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 sscal( n, rec, vr, 1 )
538  CALL sscal( 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 sscal( n, rec, vr, 1 )
565  CALL sscal( 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 sladiv( 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 = sasum( n, vr, 1 ) + sasum( 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 sscal( n, one / vnorm, vr, 1 )
624  CALL sscal( n, one / vnorm, vi, 1 )
625 *
626  END IF
627 *
628  RETURN
629 *
630 * End of SLAEIN
631 *
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:73
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:76
subroutine sladiv(A, B, C, D, P, Q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.