LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.

Download SLAEIN + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 174 of file slaein.f.

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  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 ..
209  EXTERNAL sladiv, slatrs, sscal
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 *
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
sladiv
subroutine sladiv(A, B, C, D, P, Q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition: sladiv.f:93
slapy2
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
snrm2
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
sasum
real function sasum(N, SX, INCX)
SASUM
Definition: sasum.f:54
slatrs
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:240
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: