LAPACK 3.12.0
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 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*
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:91
slapy2
real function slapy2(x, y)
SLAPY2 returns sqrt(x2+y2).
Definition slapy2.f:63
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:238
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: