280 RECURSIVE SUBROUTINE claqz0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
281 $ LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z,
282 $ LDZ, WORK, LWORK, RWORK, REC,
287 CHARACTER,
INTENT( IN ) :: wants, wantq, wantz
288 INTEGER,
INTENT( IN ) :: n, ilo, ihi, lda, ldb, ldq, ldz, lwork,
290 INTEGER,
INTENT( OUT ) :: info
291 COMPLEX,
INTENT( INOUT ) :: a( lda, * ), b( ldb, * ), q( ldq, * ),
292 $ z( ldz, * ), alpha( * ), beta( * ), work( * )
293 REAL,
INTENT( OUT ) :: rwork( * )
297 PARAMETER ( czero = ( 0.0, 0.0 ), cone = ( 1.0, 0.0 ) )
298 REAL :: zero, one, half
299 parameter( zero = 0.0, one = 1.0, half = 0.5 )
302 REAL :: smlnum, ulp, safmin, safmax, c1, tempr
303 COMPLEX :: eshift, s1, temp
304 INTEGER :: istart, istop, iiter, maxit, istart2, k, ld, nshifts,
305 $ nblock, nw, nmin, nibble, n_undeflated, n_deflated,
306 $ ns, sweep_info, shiftpos, lworkreq, k2, istartm,
307 $ istopm, iwants, iwantq, iwantz, norm_info, aed_info,
308 $ nwr, nbr, nsr, itemp1, itemp2, rcost
309 LOGICAL :: ilschur, ilq, ilz
310 CHARACTER :: jbcmpz*3
316 LOGICAL,
EXTERNAL ::
lsame
317 INTEGER,
EXTERNAL ::
ilaenv
322 IF(
lsame( wants,
'E' ) )
THEN
325 ELSE IF(
lsame( wants,
'S' ) )
THEN
332 IF(
lsame( wantq,
'N' ) )
THEN
335 ELSE IF(
lsame( wantq,
'V' ) )
THEN
338 ELSE IF(
lsame( wantq,
'I' ) )
THEN
345 IF(
lsame( wantz,
'N' ) )
THEN
348 ELSE IF(
lsame( wantz,
'V' ) )
THEN
351 ELSE IF(
lsame( wantz,
'I' ) )
THEN
361 IF( iwants.EQ.0 )
THEN
363 ELSE IF( iwantq.EQ.0 )
THEN
365 ELSE IF( iwantz.EQ.0 )
THEN
367 ELSE IF( n.LT.0 )
THEN
369 ELSE IF( ilo.LT.1 )
THEN
371 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 )
THEN
373 ELSE IF( lda.LT.n )
THEN
375 ELSE IF( ldb.LT.n )
THEN
377 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) )
THEN
379 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) )
THEN
383 CALL xerbla(
'CLAQZ0', -info )
391 work( 1 ) = real( 1 )
398 jbcmpz( 1:1 ) = wants
399 jbcmpz( 2:2 ) = wantq
400 jbcmpz( 3:3 ) = wantz
402 nmin =
ilaenv( 12,
'CLAQZ0', jbcmpz, n, ilo, ihi, lwork )
404 nwr =
ilaenv( 13,
'CLAQZ0', jbcmpz, n, ilo, ihi, lwork )
406 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
408 nibble =
ilaenv( 14,
'CLAQZ0', jbcmpz, n, ilo, ihi, lwork )
410 nsr =
ilaenv( 15,
'CLAQZ0', jbcmpz, n, ilo, ihi, lwork )
411 nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
412 nsr = max( 2, nsr-mod( nsr, 2 ) )
414 rcost =
ilaenv( 17,
'CLAQZ0', jbcmpz, n, ilo, ihi, lwork )
415 itemp1 = int( nsr/sqrt( 1+2*nsr/( real( rcost )/100*n ) ) )
416 itemp1 = ( ( itemp1-1 )/4 )*4+4
419 IF( n .LT. nmin .OR. rec .GE. 2 )
THEN
420 CALL chgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
421 $ alpha, beta, q, ldq, z, ldz, work, lwork, rwork,
431 nw = max( nwr, nmin )
432 CALL claqz2( ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b, ldb,
433 $ q, ldq, z, ldz, n_undeflated, n_deflated, alpha,
434 $ beta, work, nw, work, nw, work, -1, rwork, rec,
436 itemp1 = int( work( 1 ) )
438 CALL claqz3( ilschur, ilq, ilz, n, ilo, ihi, nsr, nbr, alpha,
439 $ beta, a, lda, b, ldb, q, ldq, z, ldz, work, nbr,
440 $ work, nbr, work, -1, sweep_info )
441 itemp2 = int( work( 1 ) )
443 lworkreq = max( itemp1+2*nw**2, itemp2+2*nbr**2 )
444 IF ( lwork .EQ.-1 )
THEN
445 work( 1 ) = real( lworkreq )
447 ELSE IF ( lwork .LT. lworkreq )
THEN
451 CALL xerbla(
'CLAQZ0', info )
457 IF( iwantq.EQ.3 )
CALL claset(
'FULL', n, n, czero, cone, q,
459 IF( iwantz.EQ.3 )
CALL claset(
'FULL', n, n, czero, cone, z,
463 safmin =
slamch(
'SAFE MINIMUM' )
465 CALL slabad( safmin, safmax )
466 ulp =
slamch(
'PRECISION' )
467 smlnum = safmin*( real( n )/ulp )
471 maxit = 30*( ihi-ilo+1 )
475 IF( iiter .GE. maxit )
THEN
479 IF ( istart+1 .GE. istop )
THEN
485 IF ( abs( a( istop, istop-1 ) ) .LE. max( smlnum,
486 $ ulp*( abs( a( istop, istop ) )+abs( a( istop-1,
487 $ istop-1 ) ) ) ) )
THEN
488 a( istop, istop-1 ) = czero
494 IF ( abs( a( istart+1, istart ) ) .LE. max( smlnum,
495 $ ulp*( abs( a( istart, istart ) )+abs( a( istart+1,
496 $ istart+1 ) ) ) ) )
THEN
497 a( istart+1, istart ) = czero
503 IF ( istart+1 .GE. istop )
THEN
509 DO k = istop, istart+1, -1
510 IF ( abs( a( k, k-1 ) ) .LE. max( smlnum, ulp*( abs( a( k,
511 $ k ) )+abs( a( k-1, k-1 ) ) ) ) )
THEN
530 DO WHILE ( k.GE.istart2 )
532 IF( k .LT. istop )
THEN
533 tempr = tempr+abs( b( k, k+1 ) )
535 IF( k .GT. istart2 )
THEN
536 tempr = tempr+abs( b( k-1, k ) )
539 IF( abs( b( k, k ) ) .LT. max( smlnum, ulp*tempr ) )
THEN
543 DO k2 = k, istart2+1, -1
544 CALL clartg( b( k2-1, k2 ), b( k2-1, k2-1 ), c1, s1,
547 b( k2-1, k2-1 ) = czero
549 CALL crot( k2-2-istartm+1, b( istartm, k2 ), 1,
550 $ b( istartm, k2-1 ), 1, c1, s1 )
551 CALL crot( min( k2+1, istop )-istartm+1, a( istartm,
552 $ k2 ), 1, a( istartm, k2-1 ), 1, c1, s1 )
554 CALL crot( n, z( 1, k2 ), 1, z( 1, k2-1 ), 1, c1,
558 IF( k2.LT.istop )
THEN
559 CALL clartg( a( k2, k2-1 ), a( k2+1, k2-1 ), c1,
562 a( k2+1, k2-1 ) = czero
564 CALL crot( istopm-k2+1, a( k2, k2 ), lda, a( k2+1,
565 $ k2 ), lda, c1, s1 )
566 CALL crot( istopm-k2+1, b( k2, k2 ), ldb, b( k2+1,
567 $ k2 ), ldb, c1, s1 )
569 CALL crot( n, q( 1, k2 ), 1, q( 1, k2+1 ), 1,
576 IF( istart2.LT.istop )
THEN
577 CALL clartg( a( istart2, istart2 ), a( istart2+1,
578 $ istart2 ), c1, s1, temp )
579 a( istart2, istart2 ) = temp
580 a( istart2+1, istart2 ) = czero
582 CALL crot( istopm-( istart2+1 )+1, a( istart2,
583 $ istart2+1 ), lda, a( istart2+1,
584 $ istart2+1 ), lda, c1, s1 )
585 CALL crot( istopm-( istart2+1 )+1, b( istart2,
586 $ istart2+1 ), ldb, b( istart2+1,
587 $ istart2+1 ), ldb, c1, s1 )
589 CALL crot( n, q( 1, istart2 ), 1, q( 1,
590 $ istart2+1 ), 1, c1, conjg( s1 ) )
602 IF ( istart2 .GE. istop )
THEN
613 IF ( istop-istart2+1 .LT. nmin )
THEN
617 IF ( istop-istart+1 .LT. nmin )
THEN
628 CALL claqz2( ilschur, ilq, ilz, n, istart2, istop, nw, a, lda,
629 $ b, ldb, q, ldq, z, ldz, n_undeflated, n_deflated,
630 $ alpha, beta, work, nw, work( nw**2+1 ), nw,
631 $ work( 2*nw**2+1 ), lwork-2*nw**2, rwork, rec,
634 IF ( n_deflated > 0 )
THEN
635 istop = istop-n_deflated
640 IF ( 100*n_deflated > nibble*( n_deflated+n_undeflated ) .OR.
641 $ istop-istart2+1 .LT. nmin )
THEN
649 ns = min( nshifts, istop-istart2 )
650 ns = min( ns, n_undeflated )
651 shiftpos = istop-n_deflated-n_undeflated+1
653 IF ( mod( ld, 6 ) .EQ. 0 )
THEN
657 IF( ( real( maxit )*safmin )*abs( a( istop,
658 $ istop-1 ) ).LT.abs( a( istop-1, istop-1 ) ) )
THEN
659 eshift = a( istop, istop-1 )/b( istop-1, istop-1 )
661 eshift = eshift+cone/( safmin*real( maxit ) )
663 alpha( shiftpos ) = cone
664 beta( shiftpos ) = eshift
671 CALL claqz3( ilschur, ilq, ilz, n, istart2, istop, ns, nblock,
672 $ alpha( shiftpos ), beta( shiftpos ), a, lda, b,
673 $ ldb, q, ldq, z, ldz, work, nblock, work( nblock**
674 $ 2+1 ), nblock, work( 2*nblock**2+1 ),
675 $ lwork-2*nblock**2, sweep_info )
685 80
CALL chgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
686 $ alpha, beta, q, ldq, z, ldz, work, lwork, rwork,
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine chgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
CHGEQZ
subroutine claqz3(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS, NBLOCK_DESIRED, ALPHA, BETA, A, LDA, B, LDB, Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK, INFO)
CLAQZ3
recursive subroutine claqz2(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHA, BETA, QC, LDQC, ZC, LDZC, WORK, LWORK, RWORK, REC, INFO)
CLAQZ2
recursive subroutine claqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, REC, INFO)
CLAQZ0
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.
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
real function slamch(CMACH)
SLAMCH