250 SUBROUTINE claqr5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S,
251 $ h, ldh, iloz, ihiz, z, ldz, v, ldv, u, ldu, nv,
252 $ wv, ldwv, nh, wh, ldwh )
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 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 claqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
CLAQR5 performs a single small-bulge multi-shift QR sweep.
REAL function slamch(CMACH)
SLAMCH
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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 slabad(SMALL, LARGE)
SLABAD
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM