LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine claqr3 ( logical  WANTT,
logical  WANTZ,
integer  N,
integer  KTOP,
integer  KBOT,
integer  NW,
complex, dimension( ldh, * )  H,
integer  LDH,
integer  ILOZ,
integer  IHIZ,
complex, dimension( ldz, * )  Z,
integer  LDZ,
integer  NS,
integer  ND,
complex, dimension( * )  SH,
complex, dimension( ldv, * )  V,
integer  LDV,
integer  NH,
complex, dimension( ldt, * )  T,
integer  LDT,
integer  NV,
complex, dimension( ldwv, * )  WV,
integer  LDWV,
complex, dimension( * )  WORK,
integer  LWORK 
)

CLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).

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

Purpose:
    Aggressive early deflation:

    CLAQR3 accepts as input an upper Hessenberg matrix
    H and performs an unitary similarity transformation
    designed to detect and deflate fully converged eigenvalues from
    a trailing principal submatrix.  On output H has been over-
    written by a new Hessenberg matrix that is a perturbation of
    an unitary similarity transformation of H.  It is to be
    hoped that the final version of H has many zero subdiagonal
    entries.
Parameters
[in]WANTT
          WANTT is LOGICAL
          If .TRUE., then the Hessenberg matrix H is fully updated
          so that the triangular Schur factor may be
          computed (in cooperation with the calling subroutine).
          If .FALSE., then only enough of H is updated to preserve
          the eigenvalues.
[in]WANTZ
          WANTZ is LOGICAL
          If .TRUE., then the unitary matrix Z is updated so
          so that the unitary Schur factor may be computed
          (in cooperation with the calling subroutine).
          If .FALSE., then Z is not referenced.
[in]N
          N is INTEGER
          The order of the matrix H and (if WANTZ is .TRUE.) the
          order of the unitary matrix Z.
[in]KTOP
          KTOP is INTEGER
          It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
          KBOT and KTOP together determine an isolated block
          along the diagonal of the Hessenberg matrix.
[in]KBOT
          KBOT is INTEGER
          It is assumed without a check that either
          KBOT = N or H(KBOT+1,KBOT)=0.  KBOT and KTOP together
          determine an isolated block along the diagonal of the
          Hessenberg matrix.
[in]NW
          NW is INTEGER
          Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).
[in,out]H
          H is COMPLEX array, dimension (LDH,N)
          On input the initial N-by-N section of H stores the
          Hessenberg matrix undergoing aggressive early deflation.
          On output H has been transformed by a unitary
          similarity transformation, perturbed, and the returned
          to Hessenberg form that (it is to be hoped) has some
          zero subdiagonal entries.
[in]LDH
          LDH is integer
          Leading dimension of H just as declared in the calling
          subroutine.  N .LE. LDH
[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 .LE. ILOZ .LE. IHIZ .LE. N.
[in,out]Z
          Z is COMPLEX array, dimension (LDZ,N)
          IF WANTZ is .TRUE., then on output, the unitary
          similarity transformation mentioned above has been
          accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
          If WANTZ is .FALSE., then Z is unreferenced.
[in]LDZ
          LDZ is integer
          The leading dimension of Z just as declared in the
          calling subroutine.  1 .LE. LDZ.
[out]NS
          NS is integer
          The number of unconverged (ie approximate) eigenvalues
          returned in SR and SI that may be used as shifts by the
          calling subroutine.
[out]ND
          ND is integer
          The number of converged eigenvalues uncovered by this
          subroutine.
[out]SH
          SH is COMPLEX array, dimension KBOT
          On output, approximate eigenvalues that may
          be used for shifts are stored in SH(KBOT-ND-NS+1)
          through SR(KBOT-ND).  Converged eigenvalues are
          stored in SH(KBOT-ND+1) through SH(KBOT).
[out]V
          V is COMPLEX array, dimension (LDV,NW)
          An NW-by-NW work array.
[in]LDV
          LDV is integer scalar
          The leading dimension of V just as declared in the
          calling subroutine.  NW .LE. LDV
[in]NH
          NH is integer scalar
          The number of columns of T.  NH.GE.NW.
[out]T
          T is COMPLEX array, dimension (LDT,NW)
[in]LDT
          LDT is integer
          The leading dimension of T just as declared in the
          calling subroutine.  NW .LE. LDT
[in]NV
          NV is integer
          The number of rows of work array WV available for
          workspace.  NV.GE.NW.
[out]WV
          WV is COMPLEX array, dimension (LDWV,NW)
[in]LDWV
          LDWV is integer
          The leading dimension of W just as declared in the
          calling subroutine.  NW .LE. LDV
[out]WORK
          WORK is COMPLEX array, dimension LWORK.
          On exit, WORK(1) is set to an estimate of the optimal value
          of LWORK for the given values of N, NW, KTOP and KBOT.
[in]LWORK
          LWORK is integer
          The dimension of the work array WORK.  LWORK = 2*NW
          suffices, but greater efficiency may result from larger
          values of LWORK.

          If LWORK = -1, then a workspace query is assumed; CLAQR3
          only estimates the optimal workspace size for the given
          values of N, NW, KTOP and KBOT.  The estimate is returned
          in WORK(1).  No error message related to LWORK is issued
          by XERBLA.  Neither H nor Z are accessed.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016
Contributors:
Karen Braman and Ralph Byers, Department of Mathematics, University of Kansas, USA

Definition at line 268 of file claqr3.f.

268 *
269 * -- LAPACK auxiliary routine (version 3.6.1) --
270 * -- LAPACK is a software package provided by Univ. of Tennessee, --
271 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
272 * June 2016
273 *
274 * .. Scalar Arguments ..
275  INTEGER ihiz, iloz, kbot, ktop, ldh, ldt, ldv, ldwv,
276  $ ldz, lwork, n, nd, nh, ns, nv, nw
277  LOGICAL wantt, wantz
278 * ..
279 * .. Array Arguments ..
280  COMPLEX h( ldh, * ), sh( * ), t( ldt, * ), v( ldv, * ),
281  $ work( * ), wv( ldwv, * ), z( ldz, * )
282 * ..
283 *
284 * ================================================================
285 *
286 * .. Parameters ..
287  COMPLEX zero, one
288  parameter ( zero = ( 0.0e0, 0.0e0 ),
289  $ one = ( 1.0e0, 0.0e0 ) )
290  REAL rzero, rone
291  parameter ( rzero = 0.0e0, rone = 1.0e0 )
292 * ..
293 * .. Local Scalars ..
294  COMPLEX beta, cdum, s, tau
295  REAL foo, safmax, safmin, smlnum, ulp
296  INTEGER i, ifst, ilst, info, infqr, j, jw, kcol, kln,
297  $ knt, krow, kwtop, ltop, lwk1, lwk2, lwk3,
298  $ lwkopt, nmin
299 * ..
300 * .. External Functions ..
301  REAL slamch
302  INTEGER ilaenv
303  EXTERNAL slamch, ilaenv
304 * ..
305 * .. External Subroutines ..
306  EXTERNAL ccopy, cgehrd, cgemm, clacpy, clahqr, claqr4,
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC abs, aimag, cmplx, conjg, int, max, min, real
311 * ..
312 * .. Statement Functions ..
313  REAL cabs1
314 * ..
315 * .. Statement Function definitions ..
316  cabs1( cdum ) = abs( REAL( CDUM ) ) + abs( aimag( cdum ) )
317 * ..
318 * .. Executable Statements ..
319 *
320 * ==== Estimate optimal workspace. ====
321 *
322  jw = min( nw, kbot-ktop+1 )
323  IF( jw.LE.2 ) THEN
324  lwkopt = 1
325  ELSE
326 *
327 * ==== Workspace query call to CGEHRD ====
328 *
329  CALL cgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
330  lwk1 = int( work( 1 ) )
331 *
332 * ==== Workspace query call to CUNMHR ====
333 *
334  CALL cunmhr( 'R', 'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
335  $ work, -1, info )
336  lwk2 = int( work( 1 ) )
337 *
338 * ==== Workspace query call to CLAQR4 ====
339 *
340  CALL claqr4( .true., .true., jw, 1, jw, t, ldt, sh, 1, jw, v,
341  $ ldv, work, -1, infqr )
342  lwk3 = int( work( 1 ) )
343 *
344 * ==== Optimal workspace ====
345 *
346  lwkopt = max( jw+max( lwk1, lwk2 ), lwk3 )
347  END IF
348 *
349 * ==== Quick return in case of workspace query. ====
350 *
351  IF( lwork.EQ.-1 ) THEN
352  work( 1 ) = cmplx( lwkopt, 0 )
353  RETURN
354  END IF
355 *
356 * ==== Nothing to do ...
357 * ... for an empty active block ... ====
358  ns = 0
359  nd = 0
360  work( 1 ) = one
361  IF( ktop.GT.kbot )
362  $ RETURN
363 * ... nor for an empty deflation window. ====
364  IF( nw.LT.1 )
365  $ RETURN
366 *
367 * ==== Machine constants ====
368 *
369  safmin = slamch( 'SAFE MINIMUM' )
370  safmax = rone / safmin
371  CALL slabad( safmin, safmax )
372  ulp = slamch( 'PRECISION' )
373  smlnum = safmin*( REAL( N ) / ulp )
374 *
375 * ==== Setup deflation window ====
376 *
377  jw = min( nw, kbot-ktop+1 )
378  kwtop = kbot - jw + 1
379  IF( kwtop.EQ.ktop ) THEN
380  s = zero
381  ELSE
382  s = h( kwtop, kwtop-1 )
383  END IF
384 *
385  IF( kbot.EQ.kwtop ) THEN
386 *
387 * ==== 1-by-1 deflation window: not much to do ====
388 *
389  sh( kwtop ) = h( kwtop, kwtop )
390  ns = 1
391  nd = 0
392  IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
393  $ kwtop ) ) ) ) THEN
394  ns = 0
395  nd = 1
396  IF( kwtop.GT.ktop )
397  $ h( kwtop, kwtop-1 ) = zero
398  END IF
399  work( 1 ) = one
400  RETURN
401  END IF
402 *
403 * ==== Convert to spike-triangular form. (In case of a
404 * . rare QR failure, this routine continues to do
405 * . aggressive early deflation using that part of
406 * . the deflation window that converged using INFQR
407 * . here and there to keep track.) ====
408 *
409  CALL clacpy( 'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
410  CALL ccopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
411 *
412  CALL claset( 'A', jw, jw, zero, one, v, ldv )
413  nmin = ilaenv( 12, 'CLAQR3', 'SV', jw, 1, jw, lwork )
414  IF( jw.GT.nmin ) THEN
415  CALL claqr4( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
416  $ jw, v, ldv, work, lwork, infqr )
417  ELSE
418  CALL clahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
419  $ jw, v, ldv, infqr )
420  END IF
421 *
422 * ==== Deflation detection loop ====
423 *
424  ns = jw
425  ilst = infqr + 1
426  DO 10 knt = infqr + 1, jw
427 *
428 * ==== Small spike tip deflation test ====
429 *
430  foo = cabs1( t( ns, ns ) )
431  IF( foo.EQ.rzero )
432  $ foo = cabs1( s )
433  IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
434  $ THEN
435 *
436 * ==== One more converged eigenvalue ====
437 *
438  ns = ns - 1
439  ELSE
440 *
441 * ==== One undeflatable eigenvalue. Move it up out of the
442 * . way. (CTREXC can not fail in this case.) ====
443 *
444  ifst = ns
445  CALL ctrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
446  ilst = ilst + 1
447  END IF
448  10 CONTINUE
449 *
450 * ==== Return to Hessenberg form ====
451 *
452  IF( ns.EQ.0 )
453  $ s = zero
454 *
455  IF( ns.LT.jw ) THEN
456 *
457 * ==== sorting the diagonal of T improves accuracy for
458 * . graded matrices. ====
459 *
460  DO 30 i = infqr + 1, ns
461  ifst = i
462  DO 20 j = i + 1, ns
463  IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
464  $ ifst = j
465  20 CONTINUE
466  ilst = i
467  IF( ifst.NE.ilst )
468  $ CALL ctrexc( 'V', jw, t, ldt, v, ldv, ifst, ilst, info )
469  30 CONTINUE
470  END IF
471 *
472 * ==== Restore shift/eigenvalue array from T ====
473 *
474  DO 40 i = infqr + 1, jw
475  sh( kwtop+i-1 ) = t( i, i )
476  40 CONTINUE
477 *
478 *
479  IF( ns.LT.jw .OR. s.EQ.zero ) THEN
480  IF( ns.GT.1 .AND. s.NE.zero ) THEN
481 *
482 * ==== Reflect spike back into lower triangle ====
483 *
484  CALL ccopy( ns, v, ldv, work, 1 )
485  DO 50 i = 1, ns
486  work( i ) = conjg( work( i ) )
487  50 CONTINUE
488  beta = work( 1 )
489  CALL clarfg( ns, beta, work( 2 ), 1, tau )
490  work( 1 ) = one
491 *
492  CALL claset( 'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
493 *
494  CALL clarf( 'L', ns, jw, work, 1, conjg( tau ), t, ldt,
495  $ work( jw+1 ) )
496  CALL clarf( 'R', ns, ns, work, 1, tau, t, ldt,
497  $ work( jw+1 ) )
498  CALL clarf( 'R', jw, ns, work, 1, tau, v, ldv,
499  $ work( jw+1 ) )
500 *
501  CALL cgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
502  $ lwork-jw, info )
503  END IF
504 *
505 * ==== Copy updated reduced window into place ====
506 *
507  IF( kwtop.GT.1 )
508  $ h( kwtop, kwtop-1 ) = s*conjg( v( 1, 1 ) )
509  CALL clacpy( 'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
510  CALL ccopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
511  $ ldh+1 )
512 *
513 * ==== Accumulate orthogonal matrix in order update
514 * . H and Z, if requested. ====
515 *
516  IF( ns.GT.1 .AND. s.NE.zero )
517  $ CALL cunmhr( 'R', 'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
518  $ work( jw+1 ), lwork-jw, info )
519 *
520 * ==== Update vertical slab in H ====
521 *
522  IF( wantt ) THEN
523  ltop = 1
524  ELSE
525  ltop = ktop
526  END IF
527  DO 60 krow = ltop, kwtop - 1, nv
528  kln = min( nv, kwtop-krow )
529  CALL cgemm( 'N', 'N', kln, jw, jw, one, h( krow, kwtop ),
530  $ ldh, v, ldv, zero, wv, ldwv )
531  CALL clacpy( 'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
532  60 CONTINUE
533 *
534 * ==== Update horizontal slab in H ====
535 *
536  IF( wantt ) THEN
537  DO 70 kcol = kbot + 1, n, nh
538  kln = min( nh, n-kcol+1 )
539  CALL cgemm( 'C', 'N', jw, kln, jw, one, v, ldv,
540  $ h( kwtop, kcol ), ldh, zero, t, ldt )
541  CALL clacpy( 'A', jw, kln, t, ldt, h( kwtop, kcol ),
542  $ ldh )
543  70 CONTINUE
544  END IF
545 *
546 * ==== Update vertical slab in Z ====
547 *
548  IF( wantz ) THEN
549  DO 80 krow = iloz, ihiz, nv
550  kln = min( nv, ihiz-krow+1 )
551  CALL cgemm( 'N', 'N', kln, jw, jw, one, z( krow, kwtop ),
552  $ ldz, v, ldv, zero, wv, ldwv )
553  CALL clacpy( 'A', kln, jw, wv, ldwv, z( krow, kwtop ),
554  $ ldz )
555  80 CONTINUE
556  END IF
557  END IF
558 *
559 * ==== Return the number of deflations ... ====
560 *
561  nd = jw - ns
562 *
563 * ==== ... and the number of shifts. (Subtracting
564 * . INFQR from the spike length takes care
565 * . of the case of a rare QR failure while
566 * . calculating eigenvalues of the deflation
567 * . window.) ====
568 *
569  ns = ns - infqr
570 *
571 * ==== Return optimal workspace. ====
572 *
573  work( 1 ) = cmplx( lwkopt, 0 )
574 *
575 * ==== End of CLAQR3 ====
576 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine clahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
Definition: clahqr.f:197
subroutine claqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
CLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: claqr4.f:251
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:108
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130
subroutine cunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMHR
Definition: cunmhr.f:181
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
Definition: cgehrd.f:169
subroutine ctrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
CTREXC
Definition: ctrexc.f:126
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.f:108

Here is the call graph for this function:

Here is the caller graph for this function: