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.

Download DLAEIN + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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 ..
206  EXTERNAL dladiv, dlatrs, dscal
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.
Definition: dladiv.f:91
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: