LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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 168 of file slaein.f.

171*
172* -- LAPACK auxiliary routine --
173* -- LAPACK is a software package provided by Univ. of Tennessee, --
174* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
175*
176* .. Scalar Arguments ..
177 LOGICAL NOINIT, RIGHTV
178 INTEGER INFO, LDB, LDH, N
179 REAL BIGNUM, EPS3, SMLNUM, WI, WR
180* ..
181* .. Array Arguments ..
182 REAL B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
183 $ WORK( * )
184* ..
185*
186* =====================================================================
187*
188* .. Parameters ..
189 REAL ZERO, ONE, TENTH
190 parameter( zero = 0.0e+0, one = 1.0e+0, tenth = 1.0e-1 )
191* ..
192* .. Local Scalars ..
193 CHARACTER NORMIN, TRANS
194 INTEGER I, I1, I2, I3, IERR, ITS, J
195 REAL ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
196 $ REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
197 $ W1, X, XI, XR, Y
198* ..
199* .. External Functions ..
200 INTEGER ISAMAX
201 REAL SASUM, SLAPY2, SNRM2
202 EXTERNAL isamax, sasum, slapy2, snrm2
203* ..
204* .. External Subroutines ..
205 EXTERNAL sladiv, slatrs, sscal
206* ..
207* .. Intrinsic Functions ..
208 INTRINSIC abs, max, real, sqrt
209* ..
210* .. Executable Statements ..
211*
212 info = 0
213*
214* GROWTO is the threshold used in the acceptance test for an
215* eigenvector.
216*
217 rootn = sqrt( real( n ) )
218 growto = tenth / rootn
219 nrmsml = max( one, eps3*rootn )*smlnum
220*
221* Form B = H - (WR,WI)*I (except that the subdiagonal elements and
222* the imaginary parts of the diagonal elements are not stored).
223*
224 DO 20 j = 1, n
225 DO 10 i = 1, j - 1
226 b( i, j ) = h( i, j )
227 10 CONTINUE
228 b( j, j ) = h( j, j ) - wr
229 20 CONTINUE
230*
231 IF( wi.EQ.zero ) THEN
232*
233* Real eigenvalue.
234*
235 IF( noinit ) THEN
236*
237* Set initial vector.
238*
239 DO 30 i = 1, n
240 vr( i ) = eps3
241 30 CONTINUE
242 ELSE
243*
244* Scale supplied initial vector.
245*
246 vnorm = snrm2( n, vr, 1 )
247 CALL sscal( n, ( eps3*rootn ) / max( vnorm, nrmsml ), vr,
248 $ 1 )
249 END IF
250*
251 IF( rightv ) THEN
252*
253* LU decomposition with partial pivoting of B, replacing zero
254* pivots by EPS3.
255*
256 DO 60 i = 1, n - 1
257 ei = h( i+1, i )
258 IF( abs( b( i, i ) ).LT.abs( ei ) ) THEN
259*
260* Interchange rows and eliminate.
261*
262 x = b( i, i ) / ei
263 b( i, i ) = ei
264 DO 40 j = i + 1, n
265 temp = b( i+1, j )
266 b( i+1, j ) = b( i, j ) - x*temp
267 b( i, j ) = temp
268 40 CONTINUE
269 ELSE
270*
271* Eliminate without interchange.
272*
273 IF( b( i, i ).EQ.zero )
274 $ b( i, i ) = eps3
275 x = ei / b( i, i )
276 IF( x.NE.zero ) THEN
277 DO 50 j = i + 1, n
278 b( i+1, j ) = b( i+1, j ) - x*b( i, j )
279 50 CONTINUE
280 END IF
281 END IF
282 60 CONTINUE
283 IF( b( n, n ).EQ.zero )
284 $ b( n, n ) = eps3
285*
286 trans = 'N'
287*
288 ELSE
289*
290* UL decomposition with partial pivoting of B, replacing zero
291* pivots by EPS3.
292*
293 DO 90 j = n, 2, -1
294 ej = h( j, j-1 )
295 IF( abs( b( j, j ) ).LT.abs( ej ) ) THEN
296*
297* Interchange columns and eliminate.
298*
299 x = b( j, j ) / ej
300 b( j, j ) = ej
301 DO 70 i = 1, j - 1
302 temp = b( i, j-1 )
303 b( i, j-1 ) = b( i, j ) - x*temp
304 b( i, j ) = temp
305 70 CONTINUE
306 ELSE
307*
308* Eliminate without interchange.
309*
310 IF( b( j, j ).EQ.zero )
311 $ b( j, j ) = eps3
312 x = ej / b( j, j )
313 IF( x.NE.zero ) THEN
314 DO 80 i = 1, j - 1
315 b( i, j-1 ) = b( i, j-1 ) - x*b( i, j )
316 80 CONTINUE
317 END IF
318 END IF
319 90 CONTINUE
320 IF( b( 1, 1 ).EQ.zero )
321 $ b( 1, 1 ) = eps3
322*
323 trans = 'T'
324*
325 END IF
326*
327 normin = 'N'
328 DO 110 its = 1, n
329*
330* Solve U*x = scale*v for a right eigenvector
331* or U**T*x = scale*v for a left eigenvector,
332* overwriting x on v.
333*
334 CALL slatrs( 'Upper', trans, 'Nonunit', normin, n, b,
335 $ 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 ),
382 $ snrm2( n, vi, 1 ) )
383 rec = ( eps3*rootn ) / max( norm, nrmsml )
384 CALL sscal( n, rec, vr, 1 )
385 CALL sscal( n, rec, vi, 1 )
386 END IF
387*
388 IF( rightv ) THEN
389*
390* LU decomposition with partial pivoting of B, replacing zero
391* pivots by EPS3.
392*
393* The imaginary part of the (i,j)-th element of U is stored in
394* B(j+1,i).
395*
396 b( 2, 1 ) = -wi
397 DO 140 i = 2, n
398 b( i+1, 1 ) = zero
399 140 CONTINUE
400*
401 DO 170 i = 1, n - 1
402 absbii = slapy2( b( i, i ), b( i+1, i ) )
403 ei = h( i+1, i )
404 IF( absbii.LT.abs( ei ) ) THEN
405*
406* Interchange rows and eliminate.
407*
408 xr = b( i, i ) / ei
409 xi = b( i+1, i ) / ei
410 b( i, i ) = ei
411 b( i+1, i ) = zero
412 DO 150 j = i + 1, n
413 temp = b( i+1, j )
414 b( i+1, j ) = b( i, j ) - xr*temp
415 b( j+1, i+1 ) = b( j+1, i ) - xi*temp
416 b( i, j ) = temp
417 b( j+1, i ) = zero
418 150 CONTINUE
419 b( i+2, i ) = -wi
420 b( i+1, i+1 ) = b( i+1, i+1 ) - xi*wi
421 b( i+2, i+1 ) = b( i+2, i+1 ) + xr*wi
422 ELSE
423*
424* Eliminate without interchanging rows.
425*
426 IF( absbii.EQ.zero ) THEN
427 b( i, i ) = eps3
428 b( i+1, i ) = zero
429 absbii = eps3
430 END IF
431 ei = ( ei / absbii ) / absbii
432 xr = b( i, i )*ei
433 xi = -b( i+1, i )*ei
434 DO 160 j = i + 1, n
435 b( i+1, j ) = b( i+1, j ) - xr*b( i, j ) +
436 $ xi*b( j+1, i )
437 b( j+1, i+1 ) = -xr*b( j+1, i ) - xi*b( i, j )
438 160 CONTINUE
439 b( i+2, i+1 ) = b( i+2, i+1 ) - wi
440 END IF
441*
442* Compute 1-norm of offdiagonal elements of i-th row.
443*
444 work( i ) = sasum( n-i, b( i, i+1 ), ldb ) +
445 $ sasum( n-i, b( i+2, i ), 1 )
446 170 CONTINUE
447 IF( b( n, n ).EQ.zero .AND. b( n+1, n ).EQ.zero )
448 $ b( n, n ) = eps3
449 work( n ) = zero
450*
451 i1 = n
452 i2 = 1
453 i3 = -1
454 ELSE
455*
456* UL decomposition with partial pivoting of conjg(B),
457* replacing zero pivots by EPS3.
458*
459* The imaginary part of the (i,j)-th element of U is stored in
460* B(j+1,i).
461*
462 b( n+1, n ) = wi
463 DO 180 j = 1, n - 1
464 b( n+1, j ) = zero
465 180 CONTINUE
466*
467 DO 210 j = n, 2, -1
468 ej = h( j, j-1 )
469 absbjj = slapy2( b( j, j ), b( j+1, j ) )
470 IF( absbjj.LT.abs( ej ) ) THEN
471*
472* Interchange columns and eliminate
473*
474 xr = b( j, j ) / ej
475 xi = b( j+1, j ) / ej
476 b( j, j ) = ej
477 b( j+1, j ) = zero
478 DO 190 i = 1, j - 1
479 temp = b( i, j-1 )
480 b( i, j-1 ) = b( i, j ) - xr*temp
481 b( j, i ) = b( j+1, i ) - xi*temp
482 b( i, j ) = temp
483 b( j+1, i ) = zero
484 190 CONTINUE
485 b( j+1, j-1 ) = wi
486 b( j-1, j-1 ) = b( j-1, j-1 ) + xi*wi
487 b( j, j-1 ) = b( j, j-1 ) - xr*wi
488 ELSE
489*
490* Eliminate without interchange.
491*
492 IF( absbjj.EQ.zero ) THEN
493 b( j, j ) = eps3
494 b( j+1, j ) = zero
495 absbjj = eps3
496 END IF
497 ej = ( ej / absbjj ) / absbjj
498 xr = b( j, j )*ej
499 xi = -b( j+1, j )*ej
500 DO 200 i = 1, j - 1
501 b( i, j-1 ) = b( i, j-1 ) - xr*b( i, j ) +
502 $ xi*b( j+1, i )
503 b( j, i ) = -xr*b( j+1, i ) - xi*b( i, j )
504 200 CONTINUE
505 b( j, j-1 ) = b( j, j-1 ) + wi
506 END IF
507*
508* Compute 1-norm of offdiagonal elements of j-th column.
509*
510 work( j ) = sasum( j-1, b( 1, j ), 1 ) +
511 $ sasum( j-1, b( j+1, 1 ), ldb )
512 210 CONTINUE
513 IF( b( 1, 1 ).EQ.zero .AND. b( 2, 1 ).EQ.zero )
514 $ b( 1, 1 ) = eps3
515 work( 1 ) = zero
516*
517 i1 = 1
518 i2 = n
519 i3 = 1
520 END IF
521*
522 DO 270 its = 1, n
523 scale = one
524 vmax = one
525 vcrit = bignum
526*
527* Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
528* or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
529* overwriting (xr,xi) on (vr,vi).
530*
531 DO 250 i = i1, i2, i3
532*
533 IF( work( i ).GT.vcrit ) THEN
534 rec = one / vmax
535 CALL sscal( n, rec, vr, 1 )
536 CALL sscal( n, rec, vi, 1 )
537 scale = scale*rec
538 vmax = one
539 vcrit = bignum
540 END IF
541*
542 xr = vr( i )
543 xi = vi( i )
544 IF( rightv ) THEN
545 DO 220 j = i + 1, n
546 xr = xr - b( i, j )*vr( j ) + b( j+1, i )*vi( j )
547 xi = xi - b( i, j )*vi( j ) - b( j+1, i )*vr( j )
548 220 CONTINUE
549 ELSE
550 DO 230 j = 1, i - 1
551 xr = xr - b( j, i )*vr( j ) + b( i+1, j )*vi( j )
552 xi = xi - b( j, i )*vi( j ) - b( i+1, j )*vr( j )
553 230 CONTINUE
554 END IF
555*
556 w = abs( b( i, i ) ) + abs( b( i+1, i ) )
557 IF( w.GT.smlnum ) THEN
558 IF( w.LT.one ) THEN
559 w1 = abs( xr ) + abs( xi )
560 IF( w1.GT.w*bignum ) THEN
561 rec = one / w1
562 CALL sscal( n, rec, vr, 1 )
563 CALL sscal( n, rec, vi, 1 )
564 xr = vr( i )
565 xi = vi( i )
566 scale = scale*rec
567 vmax = vmax*rec
568 END IF
569 END IF
570*
571* Divide by diagonal element of B.
572*
573 CALL sladiv( xr, xi, b( i, i ), b( i+1, i ),
574 $ vr( i ),
575 $ vi( i ) )
576 vmax = max( abs( vr( i ) )+abs( vi( i ) ), vmax )
577 vcrit = bignum / vmax
578 ELSE
579 DO 240 j = 1, n
580 vr( j ) = zero
581 vi( j ) = zero
582 240 CONTINUE
583 vr( i ) = one
584 vi( i ) = one
585 scale = zero
586 vmax = one
587 vcrit = bignum
588 END IF
589 250 CONTINUE
590*
591* Test for sufficient growth in the norm of (VR,VI).
592*
593 vnorm = sasum( n, vr, 1 ) + sasum( n, vi, 1 )
594 IF( vnorm.GE.growto*scale )
595 $ GO TO 280
596*
597* Choose a new orthogonal starting vector and try again.
598*
599 y = eps3 / ( rootn+one )
600 vr( 1 ) = eps3
601 vi( 1 ) = zero
602*
603 DO 260 i = 2, n
604 vr( i ) = y
605 vi( i ) = zero
606 260 CONTINUE
607 vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
608 270 CONTINUE
609*
610* Failure to find eigenvector in N iterations
611*
612 info = 1
613*
614 280 CONTINUE
615*
616* Normalize eigenvector.
617*
618 vnorm = zero
619 DO 290 i = 1, n
620 vnorm = max( vnorm, abs( vr( i ) )+abs( vi( i ) ) )
621 290 CONTINUE
622 CALL sscal( n, one / vnorm, vr, 1 )
623 CALL sscal( n, one / vnorm, vi, 1 )
624*
625 END IF
626*
627 RETURN
628*
629* End of SLAEIN
630*
sasum
real function sasum(n, sx, incx)
SASUM
Definition sasum.f:72
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
sladiv
subroutine sladiv(a, b, c, d, p, q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
Definition sladiv.f:89
slapy2
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
Definition slapy2.f:61
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:237
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: