266 SUBROUTINE zlaqr3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
267 $ ihiz, z, ldz, ns, nd, sh, v, ldv, nh, t, ldt,
268 $ nv, wv, ldwv, work, lwork )
276 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
277 $ ldz, lwork, n, nd, nh, ns, nv, nw
281 COMPLEX*16 H( ldh, * ), SH( * ), T( ldt, * ), V( ldv, * ),
282 $ work( * ), wv( ldwv, * ), z( ldz, * )
289 parameter ( zero = ( 0.0d0, 0.0d0 ),
290 $ one = ( 1.0d0, 0.0d0 ) )
291 DOUBLE PRECISION RZERO, RONE
292 parameter ( rzero = 0.0d0, rone = 1.0d0 )
295 COMPLEX*16 BETA, CDUM, S, TAU
296 DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM, ULP
297 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
298 $ knt, krow, kwtop, ltop, lwk1, lwk2, lwk3,
302 DOUBLE PRECISION DLAMCH
304 EXTERNAL dlamch, ilaenv
311 INTRINSIC abs, dble, dcmplx, dconjg, dimag, int, max, min
314 DOUBLE PRECISION CABS1
317 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
323 jw = min( nw, kbot-ktop+1 )
330 CALL zgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
331 lwk1 = int( work( 1 ) )
335 CALL zunmhr(
'R',
'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
337 lwk2 = int( work( 1 ) )
341 CALL zlaqr4( .true., .true., jw, 1, jw, t, ldt, sh, 1, jw, v,
342 $ ldv, work, -1, infqr )
343 lwk3 = int( work( 1 ) )
347 lwkopt = max( jw+max( lwk1, lwk2 ), lwk3 )
352 IF( lwork.EQ.-1 )
THEN
353 work( 1 ) = dcmplx( lwkopt, 0 )
370 safmin = dlamch(
'SAFE MINIMUM' )
371 safmax = rone / safmin
372 CALL dlabad( safmin, safmax )
373 ulp = dlamch(
'PRECISION' )
374 smlnum = safmin*( dble( n ) / ulp )
378 jw = min( nw, kbot-ktop+1 )
379 kwtop = kbot - jw + 1
380 IF( kwtop.EQ.ktop )
THEN
383 s = h( kwtop, kwtop-1 )
386 IF( kbot.EQ.kwtop )
THEN
390 sh( kwtop ) = h( kwtop, kwtop )
393 IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
398 $ h( kwtop, kwtop-1 ) = zero
410 CALL zlacpy(
'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
411 CALL zcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
413 CALL zlaset(
'A', jw, jw, zero, one, v, ldv )
414 nmin = ilaenv( 12,
'ZLAQR3',
'SV', jw, 1, jw, lwork )
415 IF( jw.GT.nmin )
THEN
416 CALL zlaqr4( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
417 $ jw, v, ldv, work, lwork, infqr )
419 CALL zlahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
420 $ jw, v, ldv, infqr )
427 DO 10 knt = infqr + 1, jw
431 foo = cabs1( t( ns, ns ) )
434 IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
446 CALL ztrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
461 DO 30 i = infqr + 1, ns
464 IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
469 $
CALL ztrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
475 DO 40 i = infqr + 1, jw
476 sh( kwtop+i-1 ) = t( i, i )
480 IF( ns.LT.jw .OR. s.EQ.zero )
THEN
481 IF( ns.GT.1 .AND. s.NE.zero )
THEN
485 CALL zcopy( ns, v, ldv, work, 1 )
487 work( i ) = dconjg( work( i ) )
490 CALL zlarfg( ns, beta, work( 2 ), 1, tau )
493 CALL zlaset(
'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
495 CALL zlarf(
'L', ns, jw, work, 1, dconjg( tau ), t, ldt,
497 CALL zlarf(
'R', ns, ns, work, 1, tau, t, ldt,
499 CALL zlarf(
'R', jw, ns, work, 1, tau, v, ldv,
502 CALL zgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
509 $ h( kwtop, kwtop-1 ) = s*dconjg( v( 1, 1 ) )
510 CALL zlacpy(
'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
511 CALL zcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
517 IF( ns.GT.1 .AND. s.NE.zero )
518 $
CALL zunmhr(
'R',
'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
519 $ work( jw+1 ), lwork-jw, info )
528 DO 60 krow = ltop, kwtop - 1, nv
529 kln = min( nv, kwtop-krow )
530 CALL zgemm(
'N',
'N', kln, jw, jw, one, h( krow, kwtop ),
531 $ ldh, v, ldv, zero, wv, ldwv )
532 CALL zlacpy(
'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
538 DO 70 kcol = kbot + 1, n, nh
539 kln = min( nh, n-kcol+1 )
540 CALL zgemm(
'C',
'N', jw, kln, jw, one, v, ldv,
541 $ h( kwtop, kcol ), ldh, zero, t, ldt )
542 CALL zlacpy(
'A', jw, kln, t, ldt, h( kwtop, kcol ),
550 DO 80 krow = iloz, ihiz, nv
551 kln = min( nv, ihiz-krow+1 )
552 CALL zgemm(
'N',
'N', kln, jw, jw, one, z( krow, kwtop ),
553 $ ldz, v, ldv, zero, wv, ldwv )
554 CALL zlacpy(
'A', kln, jw, wv, ldwv, z( krow, kwtop ),
574 work( 1 ) = dcmplx( lwkopt, 0 )
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMHR
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine zlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
subroutine zlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
subroutine ztrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
ZTREXC
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
subroutine zlaqr3(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
ZLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...