LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dlaein()

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

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

Purpose:
``` DLAEIN 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 DOUBLE PRECISION 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 DOUBLE PRECISION` [in] WI ``` WI is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)` [in,out] VI ``` VI is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDB,N)` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= N+1.``` [out] WORK ` WORK is DOUBLE PRECISION array, dimension (N)` [in] EPS3 ``` EPS3 is DOUBLE PRECISION A small machine-dependent value which is used to perturb close eigenvalues, and to replace zero pivots.``` [in] SMLNUM ``` SMLNUM is DOUBLE PRECISION A machine-dependent value close to the underflow threshold.``` [in] BIGNUM ``` BIGNUM is DOUBLE PRECISION 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.```

Definition at line 170 of file dlaein.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 DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR
181* ..
182* .. Array Arguments ..
183 DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
184 \$ WORK( * )
185* ..
186*
187* =====================================================================
188*
189* .. Parameters ..
190 DOUBLE PRECISION ZERO, ONE, TENTH
191 parameter( zero = 0.0d+0, one = 1.0d+0, tenth = 1.0d-1 )
192* ..
193* .. Local Scalars ..
194 CHARACTER NORMIN, TRANS
195 INTEGER I, I1, I2, I3, IERR, ITS, J
196 DOUBLE PRECISION 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 IDAMAX
202 DOUBLE PRECISION DASUM, DLAPY2, DNRM2
203 EXTERNAL idamax, dasum, dlapy2, dnrm2
204* ..
205* .. External Subroutines ..
207* ..
208* .. Intrinsic Functions ..
209 INTRINSIC abs, dble, max, 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( dble( 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 = dnrm2( n, vr, 1 )
248 CALL dscal( 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 dlatrs( '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 = dasum( 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 = idamax( n, vr, 1 )
364 CALL dscal( 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 = dlapy2( dnrm2( n, vr, 1 ), dnrm2( n, vi, 1 ) )
382 rec = ( eps3*rootn ) / max( norm, nrmsml )
383 CALL dscal( n, rec, vr, 1 )
384 CALL dscal( 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 = dlapy2( 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 ) = dasum( n-i, b( i, i+1 ), ldb ) +
444 \$ dasum( 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 = dlapy2( 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 ) = dasum( j-1, b( 1, j ), 1 ) +
510 \$ dasum( 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 dscal( n, rec, vr, 1 )
535 CALL dscal( 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 dscal( n, rec, vr, 1 )
562 CALL dscal( 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 dladiv( 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 = dasum( n, vr, 1 ) + dasum( 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 dscal( n, one / vnorm, vr, 1 )
621 CALL dscal( n, one / vnorm, vi, 1 )
622*
623 END IF
624*
625 RETURN
626*
627* End of DLAEIN
628*
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:63
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
double precision function dasum(N, DX, INCX)
DASUM
Definition: dasum.f:71
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
DLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition: dlatrs.f:238
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.