 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ 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.