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 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 zgerc(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERC
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 xerbla(SRNAME, INFO)
XERBLA
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
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 zlargv(N, X, INCX, Y, INCY, C, INCC)
ZLARGV generates a vector of plane rotations with real cosines and complex sines. ...
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 zgeru(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
ZGERU
logical function lsame(CA, CB)
LSAME
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
subroutine zlartg(F, G, CS, SN, R)
ZLARTG generates a plane rotation with real cosine and complex sine.