302 RECURSIVE SUBROUTINE dlaqz0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
303 $ LDA, B, LDB, ALPHAR, ALPHAI, BETA,
304 $ Q, LDQ, Z, LDZ, WORK, LWORK, REC,
309 CHARACTER,
INTENT( IN ) :: wants, wantq, wantz
310 INTEGER,
INTENT( IN ) :: n, ilo, ihi, lda, ldb, ldq, ldz, lwork,
313 INTEGER,
INTENT( OUT ) :: info
315 DOUBLE PRECISION,
INTENT( INOUT ) :: a( lda, * ), b( ldb, * ),
316 $ q( ldq, * ), z( ldz, * ), alphar( * ),
317 $ alphai( * ), beta( * ), work( * )
320 DOUBLE PRECISION :: zero, one, half
321 PARAMETER( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
324 DOUBLE PRECISION :: smlnum, ulp, eshift, safmin, safmax, c1, s1,
326 INTEGER :: istart, istop, iiter, maxit, istart2, k, ld, nshifts,
327 $ nblock, nw, nmin, nibble, n_undeflated, n_deflated,
328 $ ns, sweep_info, shiftpos, lworkreq, k2, istartm,
329 $ istopm, iwants, iwantq, iwantz, norm_info, aed_info,
330 $ nwr, nbr, nsr, itemp1, itemp2, rcost, i
331 LOGICAL :: ilschur, ilq, ilz
332 CHARACTER :: jbcmpz*3
337 DOUBLE PRECISION,
EXTERNAL ::
dlamch
338 LOGICAL,
EXTERNAL ::
lsame
339 INTEGER,
EXTERNAL ::
ilaenv
344 IF(
lsame( wants,
'E' ) )
THEN
347 ELSE IF(
lsame( wants,
'S' ) )
THEN
354 IF(
lsame( wantq,
'N' ) )
THEN
357 ELSE IF(
lsame( wantq,
'V' ) )
THEN
360 ELSE IF(
lsame( wantq,
'I' ) )
THEN
367 IF(
lsame( wantz,
'N' ) )
THEN
370 ELSE IF(
lsame( wantz,
'V' ) )
THEN
373 ELSE IF(
lsame( wantz,
'I' ) )
THEN
383 IF( iwants.EQ.0 )
THEN
385 ELSE IF( iwantq.EQ.0 )
THEN
387 ELSE IF( iwantz.EQ.0 )
THEN
389 ELSE IF( n.LT.0 )
THEN
391 ELSE IF( ilo.LT.1 )
THEN
393 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
395 ELSE IF( lda.LT.n )
THEN
397 ELSE IF( ldb.LT.n )
THEN
399 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) )
THEN
401 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) )
THEN
405 CALL xerbla(
'DLAQZ0', -info )
413 work( 1 ) = dble( 1 )
420 jbcmpz( 1:1 ) = wants
421 jbcmpz( 2:2 ) = wantq
422 jbcmpz( 3:3 ) = wantz
424 nmin =
ilaenv( 12,
'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
426 nwr =
ilaenv( 13,
'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
428 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
430 nibble =
ilaenv( 14,
'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
432 nsr =
ilaenv( 15,
'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
433 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
434 nsr = max( 2, nsr-mod( nsr, 2 ) )
436 rcost =
ilaenv( 17,
'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
437 itemp1 = int( nsr/sqrt( 1+2*nsr/( dble( rcost )/100*n ) ) )
438 itemp1 = ( ( itemp1-1 )/4 )*4+4
441 IF( n .LT. nmin .OR. rec .GE. 2 )
THEN
442 CALL dhgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
443 $ alphar, alphai, beta, q, ldq, z, ldz, work,
453 nw = max( nwr, nmin )
454 CALL dlaqz3( ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b, ldb,
455 $ q, ldq, z, ldz, n_undeflated, n_deflated, alphar,
456 $ alphai, beta, work, nw, work, nw, work, -1, rec,
458 itemp1 = int( work( 1 ) )
460 CALL dlaqz4( ilschur, ilq, ilz, n, ilo, ihi, nsr, nbr, alphar,
461 $ alphai, beta, a, lda, b, ldb, q, ldq, z, ldz, work,
462 $ nbr, work, nbr, work, -1, sweep_info )
463 itemp2 = int( work( 1 ) )
465 lworkreq = max( itemp1+2*nw**2, itemp2+2*nbr**2 )
466 IF ( lwork .EQ.-1 )
THEN
467 work( 1 ) = dble( lworkreq )
469 ELSE IF ( lwork .LT. lworkreq )
THEN
473 CALL xerbla(
'DLAQZ0', info )
479 IF( iwantq.EQ.3 )
CALL dlaset(
'FULL', n, n, zero, one, q, ldq )
480 IF( iwantz.EQ.3 )
CALL dlaset(
'FULL', n, n, zero, one, z, ldz )
483 safmin =
dlamch(
'SAFE MINIMUM' )
485 CALL dlabad( safmin, safmax )
486 ulp =
dlamch(
'PRECISION' )
487 smlnum = safmin*( dble( n )/ulp )
491 maxit = 3*( ihi-ilo+1 )
495 IF( iiter .GE. maxit )
THEN
499 IF ( istart+1 .GE. istop )
THEN
505 IF ( abs( a( istop-1, istop-2 ) ) .LE. max( smlnum,
506 $ ulp*( abs( a( istop-1, istop-1 ) )+abs( a( istop-2,
507 $ istop-2 ) ) ) ) )
THEN
508 a( istop-1, istop-2 ) = zero
512 ELSE IF ( abs( a( istop, istop-1 ) ) .LE. max( smlnum,
513 $ ulp*( abs( a( istop, istop ) )+abs( a( istop-1,
514 $ istop-1 ) ) ) ) )
THEN
515 a( istop, istop-1 ) = zero
521 IF ( abs( a( istart+2, istart+1 ) ) .LE. max( smlnum,
522 $ ulp*( abs( a( istart+1, istart+1 ) )+abs( a( istart+2,
523 $ istart+2 ) ) ) ) )
THEN
524 a( istart+2, istart+1 ) = zero
528 ELSE IF ( abs( a( istart+1, istart ) ) .LE. max( smlnum,
529 $ ulp*( abs( a( istart, istart ) )+abs( a( istart+1,
530 $ istart+1 ) ) ) ) )
THEN
531 a( istart+1, istart ) = zero
537 IF ( istart+1 .GE. istop )
THEN
543 DO k = istop, istart+1, -1
544 IF ( abs( a( k, k-1 ) ) .LE. max( smlnum, ulp*( abs( a( k,
545 $ k ) )+abs( a( k-1, k-1 ) ) ) ) )
THEN
564 DO WHILE ( k.GE.istart2 )
566 IF( k .LT. istop )
THEN
567 temp = temp+abs( b( k, k+1 ) )
569 IF( k .GT. istart2 )
THEN
570 temp = temp+abs( b( k-1, k ) )
573 IF( abs( b( k, k ) ) .LT. max( smlnum, ulp*temp ) )
THEN
577 DO k2 = k, istart2+1, -1
578 CALL dlartg( b( k2-1, k2 ), b( k2-1, k2-1 ), c1, s1,
581 b( k2-1, k2-1 ) = zero
583 CALL drot( k2-2-istartm+1, b( istartm, k2 ), 1,
584 $ b( istartm, k2-1 ), 1, c1, s1 )
585 CALL drot( min( k2+1, istop )-istartm+1, a( istartm,
586 $ k2 ), 1, a( istartm, k2-1 ), 1, c1, s1 )
588 CALL drot( n, z( 1, k2 ), 1, z( 1, k2-1 ), 1, c1,
592 IF( k2.LT.istop )
THEN
593 CALL dlartg( a( k2, k2-1 ), a( k2+1, k2-1 ), c1,
596 a( k2+1, k2-1 ) = zero
598 CALL drot( istopm-k2+1, a( k2, k2 ), lda, a( k2+1,
599 $ k2 ), lda, c1, s1 )
600 CALL drot( istopm-k2+1, b( k2, k2 ), ldb, b( k2+1,
601 $ k2 ), ldb, c1, s1 )
603 CALL drot( n, q( 1, k2 ), 1, q( 1, k2+1 ), 1,
610 IF( istart2.LT.istop )
THEN
611 CALL dlartg( a( istart2, istart2 ), a( istart2+1,
612 $ istart2 ), c1, s1, temp )
613 a( istart2, istart2 ) = temp
614 a( istart2+1, istart2 ) = zero
616 CALL drot( istopm-( istart2+1 )+1, a( istart2,
617 $ istart2+1 ), lda, a( istart2+1,
618 $ istart2+1 ), lda, c1, s1 )
619 CALL drot( istopm-( istart2+1 )+1, b( istart2,
620 $ istart2+1 ), ldb, b( istart2+1,
621 $ istart2+1 ), ldb, c1, s1 )
623 CALL drot( n, q( 1, istart2 ), 1, q( 1,
624 $ istart2+1 ), 1, c1, s1 )
636 IF ( istart2 .GE. istop )
THEN
647 IF ( istop-istart2+1 .LT. nmin )
THEN
651 IF ( istop-istart+1 .LT. nmin )
THEN
662 CALL dlaqz3( ilschur, ilq, ilz, n, istart2, istop, nw, a, lda,
663 $ b, ldb, q, ldq, z, ldz, n_undeflated, n_deflated,
664 $ alphar, alphai, beta, work, nw, work( nw**2+1 ),
665 $ nw, work( 2*nw**2+1 ), lwork-2*nw**2, rec,
668 IF ( n_deflated > 0 )
THEN
669 istop = istop-n_deflated
674 IF ( 100*n_deflated > nibble*( n_deflated+n_undeflated ) .OR.
675 $ istop-istart2+1 .LT. nmin )
THEN
683 ns = min( nshifts, istop-istart2 )
684 ns = min( ns, n_undeflated )
685 shiftpos = istop-n_deflated-n_undeflated+1
690 DO i = shiftpos, shiftpos+n_undeflated-1, 2
691 IF( alphai( i ).NE.-alphai( i+1 ) )
THEN
694 alphar( i ) = alphar( i+1 )
695 alphar( i+1 ) = alphar( i+2 )
699 alphai( i ) = alphai( i+1 )
700 alphai( i+1 ) = alphai( i+2 )
704 beta( i ) = beta( i+1 )
705 beta( i+1 ) = beta( i+2 )
710 IF ( mod( ld, 6 ) .EQ. 0 )
THEN
714 IF( ( dble( maxit )*safmin )*abs( a( istop,
715 $ istop-1 ) ).LT.abs( a( istop-1, istop-1 ) ) )
THEN
716 eshift = a( istop, istop-1 )/b( istop-1, istop-1 )
718 eshift = eshift+one/( safmin*dble( maxit ) )
720 alphar( shiftpos ) = one
721 alphar( shiftpos+1 ) = zero
722 alphai( shiftpos ) = zero
723 alphai( shiftpos+1 ) = zero
724 beta( shiftpos ) = eshift
725 beta( shiftpos+1 ) = eshift
732 CALL dlaqz4( ilschur, ilq, ilz, n, istart2, istop, ns, nblock,
733 $ alphar( shiftpos ), alphai( shiftpos ),
734 $ beta( shiftpos ), a, lda, b, ldb, q, ldq, z, ldz,
735 $ work, nblock, work( nblock**2+1 ), nblock,
736 $ work( 2*nblock**2+1 ), lwork-2*nblock**2,
747 80
CALL dhgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
748 $ alphar, alphai, beta, q, ldq, z, ldz, work, lwork,
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
recursive subroutine dlaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, REC, INFO)
DLAQZ0
subroutine dlaqz4(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS, NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK, INFO)
DLAQZ4
recursive subroutine dlaqz3(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHAR, ALPHAI, BETA, QC, LDQC, ZC, LDZC, WORK, LWORK, REC, INFO)
DLAQZ3