269 SUBROUTINE zlaqr2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
270 $ ihiz, z, ldz, ns, nd, sh, v, ldv, nh, t, ldt,
271 $ nv, wv, ldwv, work, lwork )
279 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
280 $ ldz, lwork, n, nd, nh, ns, nv, nw
284 COMPLEX*16 H( ldh, * ), SH( * ), T( ldt, * ), V( ldv, * ),
285 $ work( * ), wv( ldwv, * ), z( ldz, * )
292 parameter ( zero = ( 0.0d0, 0.0d0 ),
293 $ one = ( 1.0d0, 0.0d0 ) )
294 DOUBLE PRECISION RZERO, RONE
295 parameter ( rzero = 0.0d0, rone = 1.0d0 )
298 COMPLEX*16 BETA, CDUM, S, TAU
299 DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM, ULP
300 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
301 $ knt, krow, kwtop, ltop, lwk1, lwk2, lwkopt
304 DOUBLE PRECISION DLAMCH
312 INTRINSIC abs, dble, dcmplx, dconjg, dimag, int, max, min
315 DOUBLE PRECISION CABS1
318 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
324 jw = min( nw, kbot-ktop+1 )
331 CALL zgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
332 lwk1 = int( work( 1 ) )
336 CALL zunmhr(
'R',
'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
338 lwk2 = int( work( 1 ) )
342 lwkopt = jw + max( lwk1, lwk2 )
347 IF( lwork.EQ.-1 )
THEN
348 work( 1 ) = dcmplx( lwkopt, 0 )
365 safmin = dlamch(
'SAFE MINIMUM' )
366 safmax = rone / safmin
367 CALL dlabad( safmin, safmax )
368 ulp = dlamch(
'PRECISION' )
369 smlnum = safmin*( dble( n ) / ulp )
373 jw = min( nw, kbot-ktop+1 )
374 kwtop = kbot - jw + 1
375 IF( kwtop.EQ.ktop )
THEN
378 s = h( kwtop, kwtop-1 )
381 IF( kbot.EQ.kwtop )
THEN
385 sh( kwtop ) = h( kwtop, kwtop )
388 IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
393 $ h( kwtop, kwtop-1 ) = zero
405 CALL zlacpy(
'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
406 CALL zcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
408 CALL zlaset(
'A', jw, jw, zero, one, v, ldv )
409 CALL zlahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
410 $ jw, v, ldv, infqr )
416 DO 10 knt = infqr + 1, jw
420 foo = cabs1( t( ns, ns ) )
423 IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
435 CALL ztrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
450 DO 30 i = infqr + 1, ns
453 IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
458 $
CALL ztrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
464 DO 40 i = infqr + 1, jw
465 sh( kwtop+i-1 ) = t( i, i )
469 IF( ns.LT.jw .OR. s.EQ.zero )
THEN
470 IF( ns.GT.1 .AND. s.NE.zero )
THEN
474 CALL zcopy( ns, v, ldv, work, 1 )
476 work( i ) = dconjg( work( i ) )
479 CALL zlarfg( ns, beta, work( 2 ), 1, tau )
482 CALL zlaset(
'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
484 CALL zlarf(
'L', ns, jw, work, 1, dconjg( tau ), t, ldt,
486 CALL zlarf(
'R', ns, ns, work, 1, tau, t, ldt,
488 CALL zlarf(
'R', jw, ns, work, 1, tau, v, ldv,
491 CALL zgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
498 $ h( kwtop, kwtop-1 ) = s*dconjg( v( 1, 1 ) )
499 CALL zlacpy(
'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
500 CALL zcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
506 IF( ns.GT.1 .AND. s.NE.zero )
507 $
CALL zunmhr(
'R',
'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
508 $ work( jw+1 ), lwork-jw, info )
517 DO 60 krow = ltop, kwtop - 1, nv
518 kln = min( nv, kwtop-krow )
519 CALL zgemm(
'N',
'N', kln, jw, jw, one, h( krow, kwtop ),
520 $ ldh, v, ldv, zero, wv, ldwv )
521 CALL zlacpy(
'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
527 DO 70 kcol = kbot + 1, n, nh
528 kln = min( nh, n-kcol+1 )
529 CALL zgemm(
'C',
'N', jw, kln, jw, one, v, ldv,
530 $ h( kwtop, kcol ), ldh, zero, t, ldt )
531 CALL zlacpy(
'A', jw, kln, t, ldt, h( kwtop, kcol ),
539 DO 80 krow = iloz, ihiz, nv
540 kln = min( nv, ihiz-krow+1 )
541 CALL zgemm(
'N',
'N', kln, jw, jw, one, z( krow, kwtop ),
542 $ ldz, v, ldv, zero, wv, ldwv )
543 CALL zlacpy(
'A', kln, jw, wv, ldwv, z( krow, kwtop ),
563 work( 1 ) = dcmplx( lwkopt, 0 )
subroutine zunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMHR
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
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 zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
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 zlaqr2(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)
ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
subroutine ztrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
ZTREXC
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).