LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ slahqr()

 subroutine slahqr ( logical WANTT, logical WANTZ, integer N, integer ILO, integer IHI, real, dimension( ldh, * ) H, integer LDH, real, dimension( * ) WR, real, dimension( * ) WI, integer ILOZ, integer IHIZ, real, dimension( ldz, * ) Z, integer LDZ, integer INFO )

SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.

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

Purpose:
```    SLAHQR is an auxiliary routine called by SHSEQR to update the
eigenvalues and Schur decomposition already computed by SHSEQR, by
dealing with the Hessenberg submatrix in rows and columns ILO to
IHI.```
Parameters
 [in] WANTT ``` WANTT is LOGICAL = .TRUE. : the full Schur form T is required; = .FALSE.: only eigenvalues are required.``` [in] WANTZ ``` WANTZ is LOGICAL = .TRUE. : the matrix of Schur vectors Z is required; = .FALSE.: Schur vectors are not required.``` [in] N ``` N is INTEGER The order of the matrix H. N >= 0.``` [in] ILO ` ILO is INTEGER` [in] IHI ``` IHI is INTEGER It is assumed that H is already upper quasi-triangular in rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless ILO = 1). SLAHQR works primarily with the Hessenberg submatrix in rows and columns ILO to IHI, but applies transformations to all of H if WANTT is .TRUE.. 1 <= ILO <= max(1,IHI); IHI <= N.``` [in,out] H ``` H is REAL array, dimension (LDH,N) On entry, the upper Hessenberg matrix H. On exit, if INFO is zero and if WANTT is .TRUE., H is upper quasi-triangular in rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in standard form. If INFO is zero and WANTT is .FALSE., the contents of H are unspecified on exit. The output state of H if INFO is nonzero is given below under the description of INFO.``` [in] LDH ``` LDH is INTEGER The leading dimension of the array H. LDH >= max(1,N).``` [out] WR ` WR is REAL array, dimension (N)` [out] WI ``` WI is REAL array, dimension (N) The real and imaginary parts, respectively, of the computed eigenvalues ILO to IHI are stored in the corresponding elements of WR and WI. If two eigenvalues are computed as a complex conjugate pair, they are stored in consecutive elements of WR and WI, say the i-th and (i+1)th, with WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the eigenvalues are stored in the same order as on the diagonal of the Schur form returned in H, with WR(i) = H(i,i), and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).``` [in] ILOZ ` ILOZ is INTEGER` [in] IHIZ ``` IHIZ is INTEGER Specify the rows of Z to which transformations must be applied if WANTZ is .TRUE.. 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.``` [in,out] Z ``` Z is REAL array, dimension (LDZ,N) If WANTZ is .TRUE., on entry Z must contain the current matrix Z of transformations accumulated by SHSEQR, and on exit Z has been updated; transformations are applied only to the submatrix Z(ILOZ:IHIZ,ILO:IHI). If WANTZ is .FALSE., Z is not referenced.``` [in] LDZ ``` LDZ is INTEGER The leading dimension of the array Z. LDZ >= max(1,N).``` [out] INFO ``` INFO is INTEGER = 0: successful exit > 0: If INFO = i, SLAHQR failed to compute all the eigenvalues ILO to IHI in a total of 30 iterations per eigenvalue; elements i+1:ihi of WR and WI contain those eigenvalues which have been successfully computed. If INFO > 0 and WANTT is .FALSE., then on exit, the remaining unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix rows and columns ILO through INFO of the final, output value of H. If INFO > 0 and WANTT is .TRUE., then on exit (*) (initial value of H)*U = U*(final value of H) where U is an orthogonal matrix. The final value of H is upper Hessenberg and triangular in rows and columns INFO+1 through IHI. If INFO > 0 and WANTZ is .TRUE., then on exit (final value of Z) = (initial value of Z)*U where U is the orthogonal matrix in (*) (regardless of the value of WANTT.)```
Further Details:
```     02-96 Based on modifications by
David Day, Sandia National Laboratory, USA

12-04 Further modifications by
Ralph Byers, University of Kansas, USA
This is a modified version of SLAHQR from LAPACK version 3.0.
It is (1) more robust against overflow and underflow and
(2) adopts the more conservative Ahues & Tisseur stopping
criterion (LAWN 122, 1997).```

Definition at line 205 of file slahqr.f.

207  IMPLICIT NONE
208 *
209 * -- LAPACK auxiliary routine --
210 * -- LAPACK is a software package provided by Univ. of Tennessee, --
211 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
212 *
213 * .. Scalar Arguments ..
214  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
215  LOGICAL WANTT, WANTZ
216 * ..
217 * .. Array Arguments ..
218  REAL H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
219 * ..
220 *
221 * =========================================================
222 *
223 * .. Parameters ..
224  REAL ZERO, ONE, TWO
225  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
226  REAL DAT1, DAT2
227  parameter( dat1 = 3.0e0 / 4.0e0, dat2 = -0.4375e0 )
228  INTEGER KEXSH
229  parameter( kexsh = 10 )
230 * ..
231 * .. Local Scalars ..
232  REAL AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233  \$ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,
234  \$ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,
235  \$ ULP, V2, V3
236  INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
237  \$ KDEFL
238 * ..
239 * .. Local Arrays ..
240  REAL V( 3 )
241 * ..
242 * .. External Functions ..
243  REAL SLAMCH
244  EXTERNAL slamch
245 * ..
246 * .. External Subroutines ..
247  EXTERNAL scopy, slabad, slanv2, slarfg, srot
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC abs, max, min, real, sqrt
251 * ..
252 * .. Executable Statements ..
253 *
254  info = 0
255 *
256 * Quick return if possible
257 *
258  IF( n.EQ.0 )
259  \$ RETURN
260  IF( ilo.EQ.ihi ) THEN
261  wr( ilo ) = h( ilo, ilo )
262  wi( ilo ) = zero
263  RETURN
264  END IF
265 *
266 * ==== clear out the trash ====
267  DO 10 j = ilo, ihi - 3
268  h( j+2, j ) = zero
269  h( j+3, j ) = zero
270  10 CONTINUE
271  IF( ilo.LE.ihi-2 )
272  \$ h( ihi, ihi-2 ) = zero
273 *
274  nh = ihi - ilo + 1
275  nz = ihiz - iloz + 1
276 *
277 * Set machine-dependent constants for the stopping criterion.
278 *
279  safmin = slamch( 'SAFE MINIMUM' )
280  safmax = one / safmin
281  CALL slabad( safmin, safmax )
282  ulp = slamch( 'PRECISION' )
283  smlnum = safmin*( real( nh ) / ulp )
284 *
285 * I1 and I2 are the indices of the first row and last column of H
286 * to which transformations must be applied. If eigenvalues only are
287 * being computed, I1 and I2 are set inside the main loop.
288 *
289  IF( wantt ) THEN
290  i1 = 1
291  i2 = n
292  END IF
293 *
294 * ITMAX is the total number of QR iterations allowed.
295 *
296  itmax = 30 * max( 10, nh )
297 *
298 * KDEFL counts the number of iterations since a deflation
299 *
300  kdefl = 0
301 *
302 * The main loop begins here. I is the loop index and decreases from
303 * IHI to ILO in steps of 1 or 2. Each iteration of the loop works
304 * with the active submatrix in rows and columns L to I.
305 * Eigenvalues I+1 to IHI have already converged. Either L = ILO or
306 * H(L,L-1) is negligible so that the matrix splits.
307 *
308  i = ihi
309  20 CONTINUE
310  l = ilo
311  IF( i.LT.ilo )
312  \$ GO TO 160
313 *
314 * Perform QR iterations on rows and columns ILO to I until a
315 * submatrix of order 1 or 2 splits off at the bottom because a
316 * subdiagonal element has become negligible.
317 *
318  DO 140 its = 0, itmax
319 *
320 * Look for a single small subdiagonal element.
321 *
322  DO 30 k = i, l + 1, -1
323  IF( abs( h( k, k-1 ) ).LE.smlnum )
324  \$ GO TO 40
325  tst = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
326  IF( tst.EQ.zero ) THEN
327  IF( k-2.GE.ilo )
328  \$ tst = tst + abs( h( k-1, k-2 ) )
329  IF( k+1.LE.ihi )
330  \$ tst = tst + abs( h( k+1, k ) )
331  END IF
332 * ==== The following is a conservative small subdiagonal
333 * . deflation criterion due to Ahues & Tisseur (LAWN 122,
334 * . 1997). It has better mathematical foundation and
335 * . improves accuracy in some cases. ====
336  IF( abs( h( k, k-1 ) ).LE.ulp*tst ) THEN
337  ab = max( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
338  ba = min( abs( h( k, k-1 ) ), abs( h( k-1, k ) ) )
339  aa = max( abs( h( k, k ) ),
340  \$ abs( h( k-1, k-1 )-h( k, k ) ) )
341  bb = min( abs( h( k, k ) ),
342  \$ abs( h( k-1, k-1 )-h( k, k ) ) )
343  s = aa + ab
344  IF( ba*( ab / s ).LE.max( smlnum,
345  \$ ulp*( bb*( aa / s ) ) ) )GO TO 40
346  END IF
347  30 CONTINUE
348  40 CONTINUE
349  l = k
350  IF( l.GT.ilo ) THEN
351 *
352 * H(L,L-1) is negligible
353 *
354  h( l, l-1 ) = zero
355  END IF
356 *
357 * Exit from loop if a submatrix of order 1 or 2 has split off.
358 *
359  IF( l.GE.i-1 )
360  \$ GO TO 150
361  kdefl = kdefl + 1
362 *
363 * Now the active submatrix is in rows and columns L to I. If
364 * eigenvalues only are being computed, only the active submatrix
365 * need be transformed.
366 *
367  IF( .NOT.wantt ) THEN
368  i1 = l
369  i2 = i
370  END IF
371 *
372  IF( mod(kdefl,2*kexsh).EQ.0 ) THEN
373 *
374 * Exceptional shift.
375 *
376  s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
377  h11 = dat1*s + h( i, i )
378  h12 = dat2*s
379  h21 = s
380  h22 = h11
381  ELSE IF( mod(kdefl,kexsh).EQ.0 ) THEN
382 *
383 * Exceptional shift.
384 *
385  s = abs( h( l+1, l ) ) + abs( h( l+2, l+1 ) )
386  h11 = dat1*s + h( l, l )
387  h12 = dat2*s
388  h21 = s
389  h22 = h11
390  ELSE
391 *
392 * Prepare to use Francis' double shift
393 * (i.e. 2nd degree generalized Rayleigh quotient)
394 *
395  h11 = h( i-1, i-1 )
396  h21 = h( i, i-1 )
397  h12 = h( i-1, i )
398  h22 = h( i, i )
399  END IF
400  s = abs( h11 ) + abs( h12 ) + abs( h21 ) + abs( h22 )
401  IF( s.EQ.zero ) THEN
402  rt1r = zero
403  rt1i = zero
404  rt2r = zero
405  rt2i = zero
406  ELSE
407  h11 = h11 / s
408  h21 = h21 / s
409  h12 = h12 / s
410  h22 = h22 / s
411  tr = ( h11+h22 ) / two
412  det = ( h11-tr )*( h22-tr ) - h12*h21
413  rtdisc = sqrt( abs( det ) )
414  IF( det.GE.zero ) THEN
415 *
416 * ==== complex conjugate shifts ====
417 *
418  rt1r = tr*s
419  rt2r = rt1r
420  rt1i = rtdisc*s
421  rt2i = -rt1i
422  ELSE
423 *
424 * ==== real shifts (use only one of them) ====
425 *
426  rt1r = tr + rtdisc
427  rt2r = tr - rtdisc
428  IF( abs( rt1r-h22 ).LE.abs( rt2r-h22 ) ) THEN
429  rt1r = rt1r*s
430  rt2r = rt1r
431  ELSE
432  rt2r = rt2r*s
433  rt1r = rt2r
434  END IF
435  rt1i = zero
436  rt2i = zero
437  END IF
438  END IF
439 *
440 * Look for two consecutive small subdiagonal elements.
441 *
442  DO 50 m = i - 2, l, -1
443 * Determine the effect of starting the double-shift QR
444 * iteration at row M, and see if this would make H(M,M-1)
445 * negligible. (The following uses scaling to avoid
446 * overflows and most underflows.)
447 *
448  h21s = h( m+1, m )
449  s = abs( h( m, m )-rt2r ) + abs( rt2i ) + abs( h21s )
450  h21s = h( m+1, m ) / s
451  v( 1 ) = h21s*h( m, m+1 ) + ( h( m, m )-rt1r )*
452  \$ ( ( h( m, m )-rt2r ) / s ) - rt1i*( rt2i / s )
453  v( 2 ) = h21s*( h( m, m )+h( m+1, m+1 )-rt1r-rt2r )
454  v( 3 ) = h21s*h( m+2, m+1 )
455  s = abs( v( 1 ) ) + abs( v( 2 ) ) + abs( v( 3 ) )
456  v( 1 ) = v( 1 ) / s
457  v( 2 ) = v( 2 ) / s
458  v( 3 ) = v( 3 ) / s
459  IF( m.EQ.l )
460  \$ GO TO 60
461  IF( abs( h( m, m-1 ) )*( abs( v( 2 ) )+abs( v( 3 ) ) ).LE.
462  \$ ulp*abs( v( 1 ) )*( abs( h( m-1, m-1 ) )+abs( h( m,
463  \$ m ) )+abs( h( m+1, m+1 ) ) ) )GO TO 60
464  50 CONTINUE
465  60 CONTINUE
466 *
467 * Double-shift QR step
468 *
469  DO 130 k = m, i - 1
470 *
471 * The first iteration of this loop determines a reflection G
472 * from the vector V and applies it from left and right to H,
473 * thus creating a nonzero bulge below the subdiagonal.
474 *
475 * Each subsequent iteration determines a reflection G to
476 * restore the Hessenberg form in the (K-1)th column, and thus
477 * chases the bulge one step toward the bottom of the active
478 * submatrix. NR is the order of G.
479 *
480  nr = min( 3, i-k+1 )
481  IF( k.GT.m )
482  \$ CALL scopy( nr, h( k, k-1 ), 1, v, 1 )
483  CALL slarfg( nr, v( 1 ), v( 2 ), 1, t1 )
484  IF( k.GT.m ) THEN
485  h( k, k-1 ) = v( 1 )
486  h( k+1, k-1 ) = zero
487  IF( k.LT.i-1 )
488  \$ h( k+2, k-1 ) = zero
489  ELSE IF( m.GT.l ) THEN
490 * ==== Use the following instead of
491 * . H( K, K-1 ) = -H( K, K-1 ) to
492 * . avoid a bug when v(2) and v(3)
493 * . underflow. ====
494  h( k, k-1 ) = h( k, k-1 )*( one-t1 )
495  END IF
496  v2 = v( 2 )
497  t2 = t1*v2
498  IF( nr.EQ.3 ) THEN
499  v3 = v( 3 )
500  t3 = t1*v3
501 *
502 * Apply G from the left to transform the rows of the matrix
503 * in columns K to I2.
504 *
505  DO 70 j = k, i2
506  sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
507  h( k, j ) = h( k, j ) - sum*t1
508  h( k+1, j ) = h( k+1, j ) - sum*t2
509  h( k+2, j ) = h( k+2, j ) - sum*t3
510  70 CONTINUE
511 *
512 * Apply G from the right to transform the columns of the
513 * matrix in rows I1 to min(K+3,I).
514 *
515  DO 80 j = i1, min( k+3, i )
516  sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
517  h( j, k ) = h( j, k ) - sum*t1
518  h( j, k+1 ) = h( j, k+1 ) - sum*t2
519  h( j, k+2 ) = h( j, k+2 ) - sum*t3
520  80 CONTINUE
521 *
522  IF( wantz ) THEN
523 *
524 * Accumulate transformations in the matrix Z
525 *
526  DO 90 j = iloz, ihiz
527  sum = z( j, k ) + v2*z( j, k+1 ) + v3*z( j, k+2 )
528  z( j, k ) = z( j, k ) - sum*t1
529  z( j, k+1 ) = z( j, k+1 ) - sum*t2
530  z( j, k+2 ) = z( j, k+2 ) - sum*t3
531  90 CONTINUE
532  END IF
533  ELSE IF( nr.EQ.2 ) THEN
534 *
535 * Apply G from the left to transform the rows of the matrix
536 * in columns K to I2.
537 *
538  DO 100 j = k, i2
539  sum = h( k, j ) + v2*h( k+1, j )
540  h( k, j ) = h( k, j ) - sum*t1
541  h( k+1, j ) = h( k+1, j ) - sum*t2
542  100 CONTINUE
543 *
544 * Apply G from the right to transform the columns of the
545 * matrix in rows I1 to min(K+3,I).
546 *
547  DO 110 j = i1, i
548  sum = h( j, k ) + v2*h( j, k+1 )
549  h( j, k ) = h( j, k ) - sum*t1
550  h( j, k+1 ) = h( j, k+1 ) - sum*t2
551  110 CONTINUE
552 *
553  IF( wantz ) THEN
554 *
555 * Accumulate transformations in the matrix Z
556 *
557  DO 120 j = iloz, ihiz
558  sum = z( j, k ) + v2*z( j, k+1 )
559  z( j, k ) = z( j, k ) - sum*t1
560  z( j, k+1 ) = z( j, k+1 ) - sum*t2
561  120 CONTINUE
562  END IF
563  END IF
564  130 CONTINUE
565 *
566  140 CONTINUE
567 *
568 * Failure to converge in remaining number of iterations
569 *
570  info = i
571  RETURN
572 *
573  150 CONTINUE
574 *
575  IF( l.EQ.i ) THEN
576 *
577 * H(I,I-1) is negligible: one eigenvalue has converged.
578 *
579  wr( i ) = h( i, i )
580  wi( i ) = zero
581  ELSE IF( l.EQ.i-1 ) THEN
582 *
583 * H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
584 *
585 * Transform the 2-by-2 submatrix to standard Schur form,
586 * and compute and store the eigenvalues.
587 *
588  CALL slanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
589  \$ h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
590  \$ cs, sn )
591 *
592  IF( wantt ) THEN
593 *
594 * Apply the transformation to the rest of H.
595 *
596  IF( i2.GT.i )
597  \$ CALL srot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
598  \$ cs, sn )
599  CALL srot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
600  END IF
601  IF( wantz ) THEN
602 *
603 * Apply the transformation to Z.
604 *
605  CALL srot( nz, z( iloz, i-1 ), 1, z( iloz, i ), 1, cs, sn )
606  END IF
607  END IF
608 * reset deflation counter
609  kdefl = 0
610 *
611 * return to start of the main loop with new value of I.
612 *
613  i = l - 1
614  GO TO 20
615 *
616  160 CONTINUE
617  RETURN
618 *
619 * End of SLAHQR
620 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
Definition: slarfg.f:106
subroutine slanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition: slanv2.f:127
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: