165 SUBROUTINE zhbgst( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
166 $ ldx, work, rwork, info )
175 INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
178 DOUBLE PRECISION RWORK( * )
179 COMPLEX*16 AB( ldab, * ), BB( ldbb, * ), WORK( * ),
186 COMPLEX*16 CZERO, CONE
188 parameter ( czero = ( 0.0d+0, 0.0d+0 ),
189 $ cone = ( 1.0d+0, 0.0d+0 ), one = 1.0d+0 )
192 LOGICAL UPDATE, UPPER, WANTX
193 INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
194 $ ka1, kb1, kbt, l, m, nr, nrt, nx
196 COMPLEX*16 RA, RA1, T
207 INTRINSIC dble, dconjg, max, min
213 wantx = lsame( vect,
'V' )
214 upper = lsame( uplo,
'U' )
218 IF( .NOT.wantx .AND. .NOT.lsame( vect,
'N' ) )
THEN
220 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo,
'L' ) )
THEN
222 ELSE IF( n.LT.0 )
THEN
224 ELSE IF( ka.LT.0 )
THEN
226 ELSE IF( kb.LT.0 .OR. kb.GT.ka )
THEN
228 ELSE IF( ldab.LT.ka+1 )
THEN
230 ELSE IF( ldbb.LT.kb+1 )
THEN
232 ELSE IF( ldx.LT.1 .OR. wantx .AND. ldx.LT.max( 1, n ) )
THEN
236 CALL xerbla(
'ZHBGST', -info )
250 $
CALL zlaset(
'Full', n, n, czero, cone, x, ldx )
348 bii = dble( bb( kb1, i ) )
349 ab( ka1, i ) = ( dble( ab( ka1, i ) ) / bii ) / bii
351 ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
353 DO 30 j = max( 1, i-ka ), i - 1
354 ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
356 DO 60 k = i - kbt, i - 1
358 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
360 $ dconjg( ab( k-i+ka1, i ) ) -
361 $ dconjg( bb( k-i+kb1, i ) )*
363 $ dble( ab( ka1, i ) )*
365 $ dconjg( bb( k-i+kb1, i ) )
367 DO 50 j = max( 1, i-ka ), i - kbt - 1
368 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
369 $ dconjg( bb( k-i+kb1, i ) )*
374 DO 70 k = max( j-ka, i-kbt ), i - 1
375 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
376 $ bb( k-i+kb1, i )*ab( i-j+ka1, j )
384 CALL zdscal( n-m, one / bii, x( m+1, i ), 1 )
386 $
CALL zgerc( n-m, kbt, -cone, x( m+1, i ), 1,
387 $ bb( kb1-kbt, i ), 1, x( m+1, i-kbt ),
393 ra1 = ab( i-i1+ka1, i1 )
406 IF( i-k+ka.LT.n .AND. i-k.GT.1 )
THEN
410 CALL zlartg( ab( k+1, i-k+ka ), ra1,
411 $ rwork( i-k+ka-m ), work( i-k+ka-m ), ra )
416 t = -bb( kb1-k, i )*ra1
417 work( i-k ) = rwork( i-k+ka-m )*t -
418 $ dconjg( work( i-k+ka-m ) )*
420 ab( 1, i-k+ka ) = work( i-k+ka-m )*t +
421 $ rwork( i-k+ka-m )*ab( 1, i-k+ka )
425 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
426 nr = ( n-j2+ka ) / ka1
427 j1 = j2 + ( nr-1 )*ka1
429 j2t = max( j2, i+2*ka-k+1 )
433 nrt = ( n-j2t+ka ) / ka1
434 DO 90 j = j2t, j1, ka1
439 work( j-m ) = work( j-m )*ab( 1, j+1 )
440 ab( 1, j+1 ) = rwork( j-m )*ab( 1, j+1 )
447 $
CALL zlargv( nrt, ab( 1, j2t ), inca, work( j2t-m ), ka1,
448 $ rwork( j2t-m ), ka1 )
454 CALL zlartv( nr, ab( ka1-l, j2 ), inca,
455 $ ab( ka-l, j2+1 ), inca, rwork( j2-m ),
456 $ work( j2-m ), ka1 )
462 CALL zlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
463 $ ab( ka, j2+1 ), inca, rwork( j2-m ),
464 $ work( j2-m ), ka1 )
466 CALL zlacgv( nr, work( j2-m ), ka1 )
471 DO 110 l = ka - 1, kb - k + 1, -1
472 nrt = ( n-j2+l ) / ka1
474 $
CALL zlartv( nrt, ab( l, j2+ka1-l ), inca,
475 $ ab( l+1, j2+ka1-l ), inca, rwork( j2-m ),
476 $ work( j2-m ), ka1 )
483 DO 120 j = j2, j1, ka1
484 CALL zrot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
485 $ rwork( j-m ), dconjg( work( j-m ) ) )
491 IF( i2.LE.n .AND. kbt.GT.0 )
THEN
496 work( i-kbt ) = -bb( kb1-kbt, i )*ra1
502 j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
504 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
509 DO 140 l = kb - k, 1, -1
510 nrt = ( n-j2+ka+l ) / ka1
512 $
CALL zlartv( nrt, ab( l, j2-l+1 ), inca,
513 $ ab( l+1, j2-l+1 ), inca, rwork( j2-ka ),
514 $ work( j2-ka ), ka1 )
516 nr = ( n-j2+ka ) / ka1
517 j1 = j2 + ( nr-1 )*ka1
518 DO 150 j = j1, j2, -ka1
519 work( j ) = work( j-ka )
520 rwork( j ) = rwork( j-ka )
522 DO 160 j = j2, j1, ka1
527 work( j ) = work( j )*ab( 1, j+1 )
528 ab( 1, j+1 ) = rwork( j )*ab( 1, j+1 )
531 IF( i-k.LT.n-ka .AND. k.LE.kbt )
532 $ work( i-k+ka ) = work( i-k )
537 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
538 nr = ( n-j2+ka ) / ka1
539 j1 = j2 + ( nr-1 )*ka1
545 CALL zlargv( nr, ab( 1, j2 ), inca, work( j2 ), ka1,
551 CALL zlartv( nr, ab( ka1-l, j2 ), inca,
552 $ ab( ka-l, j2+1 ), inca, rwork( j2 ),
559 CALL zlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
560 $ ab( ka, j2+1 ), inca, rwork( j2 ),
563 CALL zlacgv( nr, work( j2 ), ka1 )
568 DO 190 l = ka - 1, kb - k + 1, -1
569 nrt = ( n-j2+l ) / ka1
571 $
CALL zlartv( nrt, ab( l, j2+ka1-l ), inca,
572 $ ab( l+1, j2+ka1-l ), inca, rwork( j2 ),
580 DO 200 j = j2, j1, ka1
581 CALL zrot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
582 $ rwork( j ), dconjg( work( j ) ) )
588 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
592 DO 220 l = kb - k, 1, -1
593 nrt = ( n-j2+l ) / ka1
595 $
CALL zlartv( nrt, ab( l, j2+ka1-l ), inca,
596 $ ab( l+1, j2+ka1-l ), inca, rwork( j2-m ),
597 $ work( j2-m ), ka1 )
602 DO 240 j = n - 1, j2 + ka, -1
603 rwork( j-m ) = rwork( j-ka-m )
604 work( j-m ) = work( j-ka-m )
616 bii = dble( bb( 1, i ) )
617 ab( 1, i ) = ( dble( ab( 1, i ) ) / bii ) / bii
619 ab( j-i+1, i ) = ab( j-i+1, i ) / bii
621 DO 260 j = max( 1, i-ka ), i - 1
622 ab( i-j+1, j ) = ab( i-j+1, j ) / bii
624 DO 290 k = i - kbt, i - 1
625 DO 270 j = i - kbt, k
626 ab( k-j+1, j ) = ab( k-j+1, j ) -
627 $ bb( i-j+1, j )*dconjg( ab( i-k+1,
628 $ k ) ) - dconjg( bb( i-k+1, k ) )*
629 $ ab( i-j+1, j ) + dble( ab( 1, i ) )*
630 $ bb( i-j+1, j )*dconjg( bb( i-k+1,
633 DO 280 j = max( 1, i-ka ), i - kbt - 1
634 ab( k-j+1, j ) = ab( k-j+1, j ) -
635 $ dconjg( bb( i-k+1, k ) )*
640 DO 300 k = max( j-ka, i-kbt ), i - 1
641 ab( j-k+1, k ) = ab( j-k+1, k ) -
642 $ bb( i-k+1, k )*ab( j-i+1, i )
650 CALL zdscal( n-m, one / bii, x( m+1, i ), 1 )
652 $
CALL zgeru( n-m, kbt, -cone, x( m+1, i ), 1,
653 $ bb( kbt+1, i-kbt ), ldbb-1,
654 $ x( m+1, i-kbt ), ldx )
659 ra1 = ab( i1-i+1, i )
672 IF( i-k+ka.LT.n .AND. i-k.GT.1 )
THEN
676 CALL zlartg( ab( ka1-k, i ), ra1, rwork( i-k+ka-m ),
677 $ work( i-k+ka-m ), ra )
682 t = -bb( k+1, i-k )*ra1
683 work( i-k ) = rwork( i-k+ka-m )*t -
684 $ dconjg( work( i-k+ka-m ) )*
686 ab( ka1, i-k ) = work( i-k+ka-m )*t +
687 $ rwork( i-k+ka-m )*ab( ka1, i-k )
691 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
692 nr = ( n-j2+ka ) / ka1
693 j1 = j2 + ( nr-1 )*ka1
695 j2t = max( j2, i+2*ka-k+1 )
699 nrt = ( n-j2t+ka ) / ka1
700 DO 320 j = j2t, j1, ka1
705 work( j-m ) = work( j-m )*ab( ka1, j-ka+1 )
706 ab( ka1, j-ka+1 ) = rwork( j-m )*ab( ka1, j-ka+1 )
713 $
CALL zlargv( nrt, ab( ka1, j2t-ka ), inca, work( j2t-m ),
714 $ ka1, rwork( j2t-m ), ka1 )
720 CALL zlartv( nr, ab( l+1, j2-l ), inca,
721 $ ab( l+2, j2-l ), inca, rwork( j2-m ),
722 $ work( j2-m ), ka1 )
728 CALL zlar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2, j2 ),
729 $ inca, rwork( j2-m ), work( j2-m ), ka1 )
731 CALL zlacgv( nr, work( j2-m ), ka1 )
736 DO 340 l = ka - 1, kb - k + 1, -1
737 nrt = ( n-j2+l ) / ka1
739 $
CALL zlartv( nrt, ab( ka1-l+1, j2 ), inca,
740 $ ab( ka1-l, j2+1 ), inca, rwork( j2-m ),
741 $ work( j2-m ), ka1 )
748 DO 350 j = j2, j1, ka1
749 CALL zrot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
750 $ rwork( j-m ), work( j-m ) )
756 IF( i2.LE.n .AND. kbt.GT.0 )
THEN
761 work( i-kbt ) = -bb( kbt+1, i-kbt )*ra1
767 j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
769 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
774 DO 370 l = kb - k, 1, -1
775 nrt = ( n-j2+ka+l ) / ka1
777 $
CALL zlartv( nrt, ab( ka1-l+1, j2-ka ), inca,
778 $ ab( ka1-l, j2-ka+1 ), inca,
779 $ rwork( j2-ka ), work( j2-ka ), ka1 )
781 nr = ( n-j2+ka ) / ka1
782 j1 = j2 + ( nr-1 )*ka1
783 DO 380 j = j1, j2, -ka1
784 work( j ) = work( j-ka )
785 rwork( j ) = rwork( j-ka )
787 DO 390 j = j2, j1, ka1
792 work( j ) = work( j )*ab( ka1, j-ka+1 )
793 ab( ka1, j-ka+1 ) = rwork( j )*ab( ka1, j-ka+1 )
796 IF( i-k.LT.n-ka .AND. k.LE.kbt )
797 $ work( i-k+ka ) = work( i-k )
802 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
803 nr = ( n-j2+ka ) / ka1
804 j1 = j2 + ( nr-1 )*ka1
810 CALL zlargv( nr, ab( ka1, j2-ka ), inca, work( j2 ), ka1,
816 CALL zlartv( nr, ab( l+1, j2-l ), inca,
817 $ ab( l+2, j2-l ), inca, rwork( j2 ),
824 CALL zlar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2, j2 ),
825 $ inca, rwork( j2 ), work( j2 ), ka1 )
827 CALL zlacgv( nr, work( j2 ), ka1 )
832 DO 420 l = ka - 1, kb - k + 1, -1
833 nrt = ( n-j2+l ) / ka1
835 $
CALL zlartv( nrt, ab( ka1-l+1, j2 ), inca,
836 $ ab( ka1-l, j2+1 ), inca, rwork( j2 ),
844 DO 430 j = j2, j1, ka1
845 CALL zrot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
846 $ rwork( j ), work( j ) )
852 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
856 DO 450 l = kb - k, 1, -1
857 nrt = ( n-j2+l ) / ka1
859 $
CALL zlartv( nrt, ab( ka1-l+1, j2 ), inca,
860 $ ab( ka1-l, j2+1 ), inca, rwork( j2-m ),
861 $ work( j2-m ), ka1 )
866 DO 470 j = n - 1, j2 + ka, -1
867 rwork( j-m ) = rwork( j-ka-m )
868 work( j-m ) = work( j-ka-m )
917 IF( i.LT.m-kbt )
THEN
931 bii = dble( bb( kb1, i ) )
932 ab( ka1, i ) = ( dble( ab( ka1, i ) ) / bii ) / bii
934 ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
936 DO 510 j = i + 1, min( n, i+ka )
937 ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
939 DO 540 k = i + 1, i + kbt
940 DO 520 j = k, i + kbt
941 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
943 $ dconjg( ab( i-k+ka1, k ) ) -
944 $ dconjg( bb( i-k+kb1, k ) )*
946 $ dble( ab( ka1, i ) )*
948 $ dconjg( bb( i-k+kb1, k ) )
950 DO 530 j = i + kbt + 1, min( n, i+ka )
951 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
952 $ dconjg( bb( i-k+kb1, k ) )*
957 DO 550 k = i + 1, min( j+ka, i+kbt )
958 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
959 $ bb( i-k+kb1, k )*ab( j-i+ka1, i )
967 CALL zdscal( nx, one / bii, x( 1, i ), 1 )
969 $
CALL zgeru( nx, kbt, -cone, x( 1, i ), 1,
970 $ bb( kb, i+1 ), ldbb-1, x( 1, i+1 ), ldx )
975 ra1 = ab( i1-i+ka1, i )
987 IF( i+k-ka1.GT.0 .AND. i+k.LT.m )
THEN
991 CALL zlartg( ab( k+1, i ), ra1, rwork( i+k-ka ),
992 $ work( i+k-ka ), ra )
997 t = -bb( kb1-k, i+k )*ra1
998 work( m-kb+i+k ) = rwork( i+k-ka )*t -
999 $ dconjg( work( i+k-ka ) )*
1001 ab( 1, i+k ) = work( i+k-ka )*t +
1002 $ rwork( i+k-ka )*ab( 1, i+k )
1006 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1007 nr = ( j2+ka-1 ) / ka1
1008 j1 = j2 - ( nr-1 )*ka1
1010 j2t = min( j2, i-2*ka+k-1 )
1014 nrt = ( j2t+ka-1 ) / ka1
1015 DO 570 j = j1, j2t, ka1
1020 work( j ) = work( j )*ab( 1, j+ka-1 )
1021 ab( 1, j+ka-1 ) = rwork( j )*ab( 1, j+ka-1 )
1028 $
CALL zlargv( nrt, ab( 1, j1+ka ), inca, work( j1 ), ka1,
1029 $ rwork( j1 ), ka1 )
1034 DO 580 l = 1, ka - 1
1035 CALL zlartv( nr, ab( ka1-l, j1+l ), inca,
1036 $ ab( ka-l, j1+l ), inca, rwork( j1 ),
1043 CALL zlar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1044 $ ab( ka, j1 ), inca, rwork( j1 ), work( j1 ),
1047 CALL zlacgv( nr, work( j1 ), ka1 )
1052 DO 590 l = ka - 1, kb - k + 1, -1
1053 nrt = ( j2+l-1 ) / ka1
1054 j1t = j2 - ( nrt-1 )*ka1
1056 $
CALL zlartv( nrt, ab( l, j1t ), inca,
1057 $ ab( l+1, j1t-1 ), inca, rwork( j1t ),
1058 $ work( j1t ), ka1 )
1065 DO 600 j = j1, j2, ka1
1066 CALL zrot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1067 $ rwork( j ), work( j ) )
1073 IF( i2.GT.0 .AND. kbt.GT.0 )
THEN
1078 work( m-kb+i+kbt ) = -bb( kb1-kbt, i+kbt )*ra1
1082 DO 650 k = kb, 1, -1
1084 j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1086 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1091 DO 620 l = kb - k, 1, -1
1092 nrt = ( j2+ka+l-1 ) / ka1
1093 j1t = j2 - ( nrt-1 )*ka1
1095 $
CALL zlartv( nrt, ab( l, j1t+ka ), inca,
1096 $ ab( l+1, j1t+ka-1 ), inca,
1097 $ rwork( m-kb+j1t+ka ),
1098 $ work( m-kb+j1t+ka ), ka1 )
1100 nr = ( j2+ka-1 ) / ka1
1101 j1 = j2 - ( nr-1 )*ka1
1102 DO 630 j = j1, j2, ka1
1103 work( m-kb+j ) = work( m-kb+j+ka )
1104 rwork( m-kb+j ) = rwork( m-kb+j+ka )
1106 DO 640 j = j1, j2, ka1
1111 work( m-kb+j ) = work( m-kb+j )*ab( 1, j+ka-1 )
1112 ab( 1, j+ka-1 ) = rwork( m-kb+j )*ab( 1, j+ka-1 )
1115 IF( i+k.GT.ka1 .AND. k.LE.kbt )
1116 $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1120 DO 690 k = kb, 1, -1
1121 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1122 nr = ( j2+ka-1 ) / ka1
1123 j1 = j2 - ( nr-1 )*ka1
1129 CALL zlargv( nr, ab( 1, j1+ka ), inca, work( m-kb+j1 ),
1130 $ ka1, rwork( m-kb+j1 ), ka1 )
1134 DO 660 l = 1, ka - 1
1135 CALL zlartv( nr, ab( ka1-l, j1+l ), inca,
1136 $ ab( ka-l, j1+l ), inca, rwork( m-kb+j1 ),
1137 $ work( m-kb+j1 ), ka1 )
1143 CALL zlar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1144 $ ab( ka, j1 ), inca, rwork( m-kb+j1 ),
1145 $ work( m-kb+j1 ), ka1 )
1147 CALL zlacgv( nr, work( m-kb+j1 ), ka1 )
1152 DO 670 l = ka - 1, kb - k + 1, -1
1153 nrt = ( j2+l-1 ) / ka1
1154 j1t = j2 - ( nrt-1 )*ka1
1156 $
CALL zlartv( nrt, ab( l, j1t ), inca,
1157 $ ab( l+1, j1t-1 ), inca,
1158 $ rwork( m-kb+j1t ), work( m-kb+j1t ),
1166 DO 680 j = j1, j2, ka1
1167 CALL zrot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1168 $ rwork( m-kb+j ), work( m-kb+j ) )
1173 DO 710 k = 1, kb - 1
1174 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1178 DO 700 l = kb - k, 1, -1
1179 nrt = ( j2+l-1 ) / ka1
1180 j1t = j2 - ( nrt-1 )*ka1
1182 $
CALL zlartv( nrt, ab( l, j1t ), inca,
1183 $ ab( l+1, j1t-1 ), inca, rwork( j1t ),
1184 $ work( j1t ), ka1 )
1189 DO 720 j = 2, i2 - ka
1190 rwork( j ) = rwork( j+ka )
1191 work( j ) = work( j+ka )
1203 bii = dble( bb( 1, i ) )
1204 ab( 1, i ) = ( dble( ab( 1, i ) ) / bii ) / bii
1205 DO 730 j = i1, i - 1
1206 ab( i-j+1, j ) = ab( i-j+1, j ) / bii
1208 DO 740 j = i + 1, min( n, i+ka )
1209 ab( j-i+1, i ) = ab( j-i+1, i ) / bii
1211 DO 770 k = i + 1, i + kbt
1212 DO 750 j = k, i + kbt
1213 ab( j-k+1, k ) = ab( j-k+1, k ) -
1214 $ bb( j-i+1, i )*dconjg( ab( k-i+1,
1215 $ i ) ) - dconjg( bb( k-i+1, i ) )*
1216 $ ab( j-i+1, i ) + dble( ab( 1, i ) )*
1217 $ bb( j-i+1, i )*dconjg( bb( k-i+1,
1220 DO 760 j = i + kbt + 1, min( n, i+ka )
1221 ab( j-k+1, k ) = ab( j-k+1, k ) -
1222 $ dconjg( bb( k-i+1, i ) )*
1227 DO 780 k = i + 1, min( j+ka, i+kbt )
1228 ab( k-j+1, j ) = ab( k-j+1, j ) -
1229 $ bb( k-i+1, i )*ab( i-j+1, j )
1237 CALL zdscal( nx, one / bii, x( 1, i ), 1 )
1239 $
CALL zgerc( nx, kbt, -cone, x( 1, i ), 1, bb( 2, i ),
1240 $ 1, x( 1, i+1 ), ldx )
1245 ra1 = ab( i-i1+1, i1 )
1251 DO 840 k = 1, kb - 1
1257 IF( i+k-ka1.GT.0 .AND. i+k.LT.m )
THEN
1261 CALL zlartg( ab( ka1-k, i+k-ka ), ra1,
1262 $ rwork( i+k-ka ), work( i+k-ka ), ra )
1267 t = -bb( k+1, i )*ra1
1268 work( m-kb+i+k ) = rwork( i+k-ka )*t -
1269 $ dconjg( work( i+k-ka ) )*
1271 ab( ka1, i+k-ka ) = work( i+k-ka )*t +
1272 $ rwork( i+k-ka )*ab( ka1, i+k-ka )
1276 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1277 nr = ( j2+ka-1 ) / ka1
1278 j1 = j2 - ( nr-1 )*ka1
1280 j2t = min( j2, i-2*ka+k-1 )
1284 nrt = ( j2t+ka-1 ) / ka1
1285 DO 800 j = j1, j2t, ka1
1290 work( j ) = work( j )*ab( ka1, j-1 )
1291 ab( ka1, j-1 ) = rwork( j )*ab( ka1, j-1 )
1298 $
CALL zlargv( nrt, ab( ka1, j1 ), inca, work( j1 ), ka1,
1299 $ rwork( j1 ), ka1 )
1304 DO 810 l = 1, ka - 1
1305 CALL zlartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),
1306 $ inca, rwork( j1 ), work( j1 ), ka1 )
1312 CALL zlar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1313 $ ab( 2, j1-1 ), inca, rwork( j1 ),
1316 CALL zlacgv( nr, work( j1 ), ka1 )
1321 DO 820 l = ka - 1, kb - k + 1, -1
1322 nrt = ( j2+l-1 ) / ka1
1323 j1t = j2 - ( nrt-1 )*ka1
1325 $
CALL zlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1326 $ ab( ka1-l, j1t-ka1+l ), inca,
1327 $ rwork( j1t ), work( j1t ), ka1 )
1334 DO 830 j = j1, j2, ka1
1335 CALL zrot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1336 $ rwork( j ), dconjg( work( j ) ) )
1342 IF( i2.GT.0 .AND. kbt.GT.0 )
THEN
1347 work( m-kb+i+kbt ) = -bb( kbt+1, i )*ra1
1351 DO 880 k = kb, 1, -1
1353 j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1355 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1360 DO 850 l = kb - k, 1, -1
1361 nrt = ( j2+ka+l-1 ) / ka1
1362 j1t = j2 - ( nrt-1 )*ka1
1364 $
CALL zlartv( nrt, ab( ka1-l+1, j1t+l-1 ), inca,
1365 $ ab( ka1-l, j1t+l-1 ), inca,
1366 $ rwork( m-kb+j1t+ka ),
1367 $ work( m-kb+j1t+ka ), ka1 )
1369 nr = ( j2+ka-1 ) / ka1
1370 j1 = j2 - ( nr-1 )*ka1
1371 DO 860 j = j1, j2, ka1
1372 work( m-kb+j ) = work( m-kb+j+ka )
1373 rwork( m-kb+j ) = rwork( m-kb+j+ka )
1375 DO 870 j = j1, j2, ka1
1380 work( m-kb+j ) = work( m-kb+j )*ab( ka1, j-1 )
1381 ab( ka1, j-1 ) = rwork( m-kb+j )*ab( ka1, j-1 )
1384 IF( i+k.GT.ka1 .AND. k.LE.kbt )
1385 $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1389 DO 920 k = kb, 1, -1
1390 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1391 nr = ( j2+ka-1 ) / ka1
1392 j1 = j2 - ( nr-1 )*ka1
1398 CALL zlargv( nr, ab( ka1, j1 ), inca, work( m-kb+j1 ),
1399 $ ka1, rwork( m-kb+j1 ), ka1 )
1403 DO 890 l = 1, ka - 1
1404 CALL zlartv( nr, ab( l+1, j1 ), inca, ab( l+2, j1-1 ),
1405 $ inca, rwork( m-kb+j1 ), work( m-kb+j1 ),
1412 CALL zlar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1413 $ ab( 2, j1-1 ), inca, rwork( m-kb+j1 ),
1414 $ work( m-kb+j1 ), ka1 )
1416 CALL zlacgv( nr, work( m-kb+j1 ), ka1 )
1421 DO 900 l = ka - 1, kb - k + 1, -1
1422 nrt = ( j2+l-1 ) / ka1
1423 j1t = j2 - ( nrt-1 )*ka1
1425 $
CALL zlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1426 $ ab( ka1-l, j1t-ka1+l ), inca,
1427 $ rwork( m-kb+j1t ), work( m-kb+j1t ),
1435 DO 910 j = j1, j2, ka1
1436 CALL zrot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1437 $ rwork( m-kb+j ), dconjg( work( m-kb+j ) ) )
1442 DO 940 k = 1, kb - 1
1443 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1447 DO 930 l = kb - k, 1, -1
1448 nrt = ( j2+l-1 ) / ka1
1449 j1t = j2 - ( nrt-1 )*ka1
1451 $
CALL zlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1452 $ ab( ka1-l, j1t-ka1+l ), inca,
1453 $ rwork( j1t ), work( j1t ), ka1 )
1458 DO 950 j = 2, i2 - ka
1459 rwork( j ) = rwork( j+ka )
1460 work( j ) = work( j+ka )
subroutine zlartv(N, X, INCX, Y, INCY, C, S, INCC)
ZLARTV applies a vector of plane rotations with real cosines and complex sines to the elements of a p...
subroutine zgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERC
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
subroutine zlar2v(N, X, Y, Z, INCX, C, S, INCC)
ZLAR2V applies a vector of plane rotations with real cosines and complex sines from both sides to a s...
subroutine xerbla(SRNAME, INFO)
XERBLA
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 zlargv(N, X, INCX, Y, INCY, C, INCC)
ZLARGV generates a vector of plane rotations with real cosines and complex sines. ...
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine zgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERU
subroutine zhbgst(VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X, LDX, WORK, RWORK, INFO)
ZHBGST
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.