LAPACK  3.10.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.

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.

Definition at line 170 of file slaein.f.

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 *
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:63
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
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 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
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
real function sasum(N, SX, INCX)
SASUM
Definition: sasum.f:72
Here is the call graph for this function:
Here is the caller graph for this function: