260 INTEGER ihiz, iloz, kacc22, kbot, ktop, ldh, ldu, ldv,
261 $ ldwh, ldwv, ldz, n, nh, nshfts, nv
265 COMPLEX h( ldh, * ), s( * ), u( ldu, * ), v( ldv, * ),
266 $ wh( ldwh, * ), wv( ldwv, * ), z( ldz, * )
272 parameter ( zero = ( 0.0e0, 0.0e0 ),
273 $ one = ( 1.0e0, 0.0e0 ) )
275 parameter ( rzero = 0.0e0, rone = 1.0e0 )
278 COMPLEX alpha, beta, cdum, refsum
279 REAL h11, h12, h21, h22, safmax, safmin, scl,
280 $ smlnum, tst1, tst2, ulp
281 INTEGER i2, i4, incol, j, j2, j4, jbot, jcol, jlen,
282 $ jrow, jtop, k, k1, kdu, kms, knz, krcol, kzs,
283 $ m, m22, mbot, mend, mstart, mtop, nbmps, ndcol,
285 LOGICAL accum, blk22, bmp22
293 INTRINSIC abs, aimag, conjg, max, min, mod, real
306 cabs1( cdum ) = abs(
REAL( CDUM ) ) + abs( aimag( cdum ) )
324 ns = nshfts - mod( nshfts, 2 )
328 safmin =
slamch(
'SAFE MINIMUM' )
329 safmax = rone / safmin
330 CALL slabad( safmin, safmax )
331 ulp =
slamch(
'PRECISION' )
332 smlnum = safmin*(
REAL( N ) / ulp )
337 accum = ( kacc22.EQ.1 ) .OR. ( kacc22.EQ.2 )
341 blk22 = ( ns.GT.2 ) .AND. ( kacc22.EQ.2 )
346 $ h( ktop+2, ktop ) = zero
358 DO 210 incol = 3*( 1-nbmps ) + ktop - 1, kbot - 2, 3*nbmps - 2
361 $
CALL claset(
'ALL', kdu, kdu, zero, one, u, ldu )
375 DO 140 krcol = incol, min( incol+3*nbmps-3, kbot-2 )
384 mtop = max( 1, ( ( ktop-1 )-krcol+2 ) / 3+1 )
385 mbot = min( nbmps, ( kbot-krcol ) / 3 )
387 bmp22 = ( mbot.LT.nbmps ) .AND. ( krcol+3*( m22-1 ) ).EQ.
394 k = krcol + 3*( m-1 )
395 IF( k.EQ.ktop-1 )
THEN
396 CALL claqr1( 3, h( ktop, ktop ), ldh, s( 2*m-1 ),
397 $ s( 2*m ), v( 1, m ) )
399 CALL clarfg( 3, alpha, v( 2, m ), 1, v( 1, m ) )
402 v( 2, m ) = h( k+2, k )
403 v( 3, m ) = h( k+3, k )
404 CALL clarfg( 3, beta, v( 2, m ), 1, v( 1, m ) )
411 IF( h( k+3, k ).NE.zero .OR. h( k+3, k+1 ).NE.
412 $ zero .OR. h( k+3, k+2 ).EQ.zero )
THEN
427 CALL claqr1( 3, h( k+1, k+1 ), ldh, s( 2*m-1 ),
430 CALL clarfg( 3, alpha, vt( 2 ), 1, vt( 1 ) )
431 refsum = conjg( vt( 1 ) )*
432 $ ( h( k+1, k )+conjg( vt( 2 ) )*
435 IF( cabs1( h( k+2, k )-refsum*vt( 2 ) )+
436 $ cabs1( refsum*vt( 3 ) ).GT.ulp*
437 $ ( cabs1( h( k, k ) )+cabs1( h( k+1,
438 $ k+1 ) )+cabs1( h( k+2, k+2 ) ) ) )
THEN
454 h( k+1, k ) = h( k+1, k ) - refsum
467 k = krcol + 3*( m22-1 )
469 IF( k.EQ.ktop-1 )
THEN
470 CALL claqr1( 2, h( k+1, k+1 ), ldh, s( 2*m22-1 ),
471 $ s( 2*m22 ), v( 1, m22 ) )
473 CALL clarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
476 v( 2, m22 ) = h( k+2, k )
477 CALL clarfg( 2, beta, v( 2, m22 ), 1, v( 1, m22 ) )
486 jbot = min( ndcol, kbot )
487 ELSE IF( wantt )
THEN
492 DO 30 j = max( ktop, krcol ), jbot
493 mend = min( mbot, ( j-krcol+2 ) / 3 )
495 k = krcol + 3*( m-1 )
496 refsum = conjg( v( 1, m ) )*
497 $ ( h( k+1, j )+conjg( v( 2, m ) )*h( k+2, j )+
498 $ conjg( v( 3, m ) )*h( k+3, j ) )
499 h( k+1, j ) = h( k+1, j ) - refsum
500 h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m )
501 h( k+3, j ) = h( k+3, j ) - refsum*v( 3, m )
505 k = krcol + 3*( m22-1 )
506 DO 40 j = max( k+1, ktop ), jbot
507 refsum = conjg( v( 1, m22 ) )*
508 $ ( h( k+1, j )+conjg( v( 2, m22 ) )*
510 h( k+1, j ) = h( k+1, j ) - refsum
511 h( k+2, j ) = h( k+2, j ) - refsum*v( 2, m22 )
520 jtop = max( ktop, incol )
521 ELSE IF( wantt )
THEN
527 IF( v( 1, m ).NE.zero )
THEN
528 k = krcol + 3*( m-1 )
529 DO 50 j = jtop, min( kbot, k+3 )
530 refsum = v( 1, m )*( h( j, k+1 )+v( 2, m )*
531 $ h( j, k+2 )+v( 3, m )*h( j, k+3 ) )
532 h( j, k+1 ) = h( j, k+1 ) - refsum
533 h( j, k+2 ) = h( j, k+2 ) -
534 $ refsum*conjg( v( 2, m ) )
535 h( j, k+3 ) = h( j, k+3 ) -
536 $ refsum*conjg( v( 3, m ) )
546 DO 60 j = max( 1, ktop-incol ), kdu
547 refsum = v( 1, m )*( u( j, kms+1 )+v( 2, m )*
548 $ u( j, kms+2 )+v( 3, m )*u( j, kms+3 ) )
549 u( j, kms+1 ) = u( j, kms+1 ) - refsum
550 u( j, kms+2 ) = u( j, kms+2 ) -
551 $ refsum*conjg( v( 2, m ) )
552 u( j, kms+3 ) = u( j, kms+3 ) -
553 $ refsum*conjg( v( 3, m ) )
555 ELSE IF( wantz )
THEN
562 refsum = v( 1, m )*( z( j, k+1 )+v( 2, m )*
563 $ z( j, k+2 )+v( 3, m )*z( j, k+3 ) )
564 z( j, k+1 ) = z( j, k+1 ) - refsum
565 z( j, k+2 ) = z( j, k+2 ) -
566 $ refsum*conjg( v( 2, m ) )
567 z( j, k+3 ) = z( j, k+3 ) -
568 $ refsum*conjg( v( 3, m ) )
576 k = krcol + 3*( m22-1 )
578 IF ( v( 1, m22 ).NE.zero )
THEN
579 DO 90 j = jtop, min( kbot, k+3 )
580 refsum = v( 1, m22 )*( h( j, k+1 )+v( 2, m22 )*
582 h( j, k+1 ) = h( j, k+1 ) - refsum
583 h( j, k+2 ) = h( j, k+2 ) -
584 $ refsum*conjg( v( 2, m22 ) )
589 DO 100 j = max( 1, ktop-incol ), kdu
590 refsum = v( 1, m22 )*( u( j, kms+1 )+
591 $ v( 2, m22 )*u( j, kms+2 ) )
592 u( j, kms+1 ) = u( j, kms+1 ) - refsum
593 u( j, kms+2 ) = u( j, kms+2 ) -
594 $ refsum*conjg( v( 2, m22 ) )
596 ELSE IF( wantz )
THEN
597 DO 110 j = iloz, ihiz
598 refsum = v( 1, m22 )*( z( j, k+1 )+v( 2, m22 )*
600 z( j, k+1 ) = z( j, k+1 ) - refsum
601 z( j, k+2 ) = z( j, k+2 ) -
602 $ refsum*conjg( v( 2, m22 ) )
611 IF( krcol+3*( mstart-1 ).LT.ktop )
612 $ mstart = mstart + 1
616 IF( krcol.EQ.kbot-2 )
618 DO 120 m = mstart, mend
619 k = min( kbot-1, krcol+3*( m-1 ) )
630 IF( h( k+1, k ).NE.zero )
THEN
631 tst1 = cabs1( h( k, k ) ) + cabs1( h( k+1, k+1 ) )
632 IF( tst1.EQ.rzero )
THEN
634 $ tst1 = tst1 + cabs1( h( k, k-1 ) )
636 $ tst1 = tst1 + cabs1( h( k, k-2 ) )
638 $ tst1 = tst1 + cabs1( h( k, k-3 ) )
640 $ tst1 = tst1 + cabs1( h( k+2, k+1 ) )
642 $ tst1 = tst1 + cabs1( h( k+3, k+1 ) )
644 $ tst1 = tst1 + cabs1( h( k+4, k+1 ) )
646 IF( cabs1( h( k+1, k ) ).LE.max( smlnum, ulp*tst1 ) )
648 h12 = max( cabs1( h( k+1, k ) ),
649 $ cabs1( h( k, k+1 ) ) )
650 h21 = min( cabs1( h( k+1, k ) ),
651 $ cabs1( h( k, k+1 ) ) )
652 h11 = max( cabs1( h( k+1, k+1 ) ),
653 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
654 h22 = min( cabs1( h( k+1, k+1 ) ),
655 $ cabs1( h( k, k )-h( k+1, k+1 ) ) )
657 tst2 = h22*( h11 / scl )
659 IF( tst2.EQ.rzero .OR. h21*( h12 / scl ).LE.
660 $ max( smlnum, ulp*tst2 ) )h( k+1, k ) = zero
667 mend = min( nbmps, ( kbot-krcol-1 ) / 3 )
668 DO 130 m = mtop, mend
669 k = krcol + 3*( m-1 )
670 refsum = v( 1, m )*v( 3, m )*h( k+4, k+3 )
671 h( k+4, k+1 ) = -refsum
672 h( k+4, k+2 ) = -refsum*conjg( v( 2, m ) )
673 h( k+4, k+3 ) = h( k+4, k+3 ) - refsum*conjg( v( 3, m ) )
692 IF( ( .NOT.blk22 ) .OR. ( incol.LT.ktop ) .OR.
693 $ ( ndcol.GT.kbot ) .OR. ( ns.LE.2 ) )
THEN
704 k1 = max( 1, ktop-incol )
705 nu = ( kdu-max( 0, ndcol-kbot ) ) - k1 + 1
709 DO 150 jcol = min( ndcol, kbot ) + 1, jbot, nh
710 jlen = min( nh, jbot-jcol+1 )
711 CALL cgemm(
'C',
'N', nu, jlen, nu, one, u( k1, k1 ),
712 $ ldu, h( incol+k1, jcol ), ldh, zero, wh,
714 CALL clacpy(
'ALL', nu, jlen, wh, ldwh,
715 $ h( incol+k1, jcol ), ldh )
720 DO 160 jrow = jtop, max( ktop, incol ) - 1, nv
721 jlen = min( nv, max( ktop, incol )-jrow )
722 CALL cgemm(
'N',
'N', jlen, nu, nu, one,
723 $ h( jrow, incol+k1 ), ldh, u( k1, k1 ),
724 $ ldu, zero, wv, ldwv )
725 CALL clacpy(
'ALL', jlen, nu, wv, ldwv,
726 $ h( jrow, incol+k1 ), ldh )
732 DO 170 jrow = iloz, ihiz, nv
733 jlen = min( nv, ihiz-jrow+1 )
734 CALL cgemm(
'N',
'N', jlen, nu, nu, one,
735 $ z( jrow, incol+k1 ), ldz, u( k1, k1 ),
736 $ ldu, zero, wv, ldwv )
737 CALL clacpy(
'ALL', jlen, nu, wv, ldwv,
738 $ z( jrow, incol+k1 ), ldz )
756 kzs = ( j4-j2 ) - ( ns+1 )
761 DO 180 jcol = min( ndcol, kbot ) + 1, jbot, nh
762 jlen = min( nh, jbot-jcol+1 )
767 CALL clacpy(
'ALL', knz, jlen, h( incol+1+j2, jcol ),
768 $ ldh, wh( kzs+1, 1 ), ldwh )
772 CALL claset(
'ALL', kzs, jlen, zero, zero, wh, ldwh )
773 CALL ctrmm(
'L',
'U',
'C',
'N', knz, jlen, one,
774 $ u( j2+1, 1+kzs ), ldu, wh( kzs+1, 1 ),
779 CALL cgemm(
'C',
'N', i2, jlen, j2, one, u, ldu,
780 $ h( incol+1, jcol ), ldh, one, wh, ldwh )
784 CALL clacpy(
'ALL', j2, jlen, h( incol+1, jcol ), ldh,
785 $ wh( i2+1, 1 ), ldwh )
789 CALL ctrmm(
'L',
'L',
'C',
'N', j2, jlen, one,
790 $ u( 1, i2+1 ), ldu, wh( i2+1, 1 ), ldwh )
794 CALL cgemm(
'C',
'N', i4-i2, jlen, j4-j2, one,
795 $ u( j2+1, i2+1 ), ldu,
796 $ h( incol+1+j2, jcol ), ldh, one,
797 $ wh( i2+1, 1 ), ldwh )
801 CALL clacpy(
'ALL', kdu, jlen, wh, ldwh,
802 $ h( incol+1, jcol ), ldh )
807 DO 190 jrow = jtop, max( incol, ktop ) - 1, nv
808 jlen = min( nv, max( incol, ktop )-jrow )
813 CALL clacpy(
'ALL', jlen, knz, h( jrow, incol+1+j2 ),
814 $ ldh, wv( 1, 1+kzs ), ldwv )
818 CALL claset(
'ALL', jlen, kzs, zero, zero, wv, ldwv )
819 CALL ctrmm(
'R',
'U',
'N',
'N', jlen, knz, one,
820 $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
825 CALL cgemm(
'N',
'N', jlen, i2, j2, one,
826 $ h( jrow, incol+1 ), ldh, u, ldu, one, wv,
831 CALL clacpy(
'ALL', jlen, j2, h( jrow, incol+1 ), ldh,
832 $ wv( 1, 1+i2 ), ldwv )
836 CALL ctrmm(
'R',
'L',
'N',
'N', jlen, i4-i2, one,
837 $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ), ldwv )
841 CALL cgemm(
'N',
'N', jlen, i4-i2, j4-j2, one,
842 $ h( jrow, incol+1+j2 ), ldh,
843 $ u( j2+1, i2+1 ), ldu, one, wv( 1, 1+i2 ),
848 CALL clacpy(
'ALL', jlen, kdu, wv, ldwv,
849 $ h( jrow, incol+1 ), ldh )
855 DO 200 jrow = iloz, ihiz, nv
856 jlen = min( nv, ihiz-jrow+1 )
861 CALL clacpy(
'ALL', jlen, knz,
862 $ z( jrow, incol+1+j2 ), ldz,
863 $ wv( 1, 1+kzs ), ldwv )
867 CALL claset(
'ALL', jlen, kzs, zero, zero, wv,
869 CALL ctrmm(
'R',
'U',
'N',
'N', jlen, knz, one,
870 $ u( j2+1, 1+kzs ), ldu, wv( 1, 1+kzs ),
875 CALL cgemm(
'N',
'N', jlen, i2, j2, one,
876 $ z( jrow, incol+1 ), ldz, u, ldu, one,
881 CALL clacpy(
'ALL', jlen, j2, z( jrow, incol+1 ),
882 $ ldz, wv( 1, 1+i2 ), ldwv )
886 CALL ctrmm(
'R',
'L',
'N',
'N', jlen, i4-i2, one,
887 $ u( 1, i2+1 ), ldu, wv( 1, 1+i2 ),
892 CALL cgemm(
'N',
'N', jlen, i4-i2, j4-j2, one,
893 $ z( jrow, incol+1+j2 ), ldz,
894 $ u( j2+1, i2+1 ), ldu, one,
895 $ wv( 1, 1+i2 ), ldwv )
899 CALL clacpy(
'ALL', jlen, kdu, wv, ldwv,
900 $ z( jrow, incol+1 ), ldz )
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
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...
real function slamch(CMACH)
SLAMCH
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
subroutine claqr1(N, H, LDH, S1, S2, V)
CLAQR1 sets a scalar multiple of the first column of the product of 2-by-2 or 3-by-3 matrix H and spe...
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM