164
165
166
167
168
169
170 CHARACTER UPLO, VECT
171 INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
172
173
174 DOUBLE PRECISION RWORK( * )
175 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
176 $ X( LDX, * )
177
178
179
180
181
182 COMPLEX*16 CZERO, CONE
183 DOUBLE PRECISION ONE
184 parameter( czero = ( 0.0d+0, 0.0d+0 ),
185 $ cone = ( 1.0d+0, 0.0d+0 ), one = 1.0d+0 )
186
187
188 LOGICAL UPDATE, UPPER, WANTX
189 INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
190 $ KA1, KB1, KBT, L, M, NR, NRT, NX
191 DOUBLE PRECISION BII
192 COMPLEX*16 RA, RA1, T
193
194
195 LOGICAL LSAME
197
198
202
203
204 INTRINSIC dble, dconjg, max, min
205
206
207
208
209
210 wantx =
lsame( vect,
'V' )
211 upper =
lsame( uplo,
'U' )
212 ka1 = ka + 1
213 kb1 = kb + 1
214 info = 0
215 IF( .NOT.wantx .AND. .NOT.
lsame( vect,
'N' ) )
THEN
216 info = -1
217 ELSE IF( .NOT.upper .AND. .NOT.
lsame( uplo,
'L' ) )
THEN
218 info = -2
219 ELSE IF( n.LT.0 ) THEN
220 info = -3
221 ELSE IF( ka.LT.0 ) THEN
222 info = -4
223 ELSE IF( kb.LT.0 .OR. kb.GT.ka ) THEN
224 info = -5
225 ELSE IF( ldab.LT.ka+1 ) THEN
226 info = -7
227 ELSE IF( ldbb.LT.kb+1 ) THEN
228 info = -9
229 ELSE IF( ldx.LT.1 .OR. wantx .AND. ldx.LT.max( 1, n ) ) THEN
230 info = -11
231 END IF
232 IF( info.NE.0 ) THEN
233 CALL xerbla(
'ZHBGST', -info )
234 RETURN
235 END IF
236
237
238
239 IF( n.EQ.0 )
240 $ RETURN
241
242 inca = ldab*ka1
243
244
245
246 IF( wantx )
247 $
CALL zlaset(
'Full', n, n, czero, cone, x, ldx )
248
249
250
251
252
253 m = ( n+kb ) / 2
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314 update = .true.
315 i = n + 1
316 10 CONTINUE
317 IF( update ) THEN
318 i = i - 1
319 kbt = min( kb, i-1 )
320 i0 = i - 1
321 i1 = min( n, i+ka )
322 i2 = i - kbt + ka1
323 IF( i.LT.m+1 ) THEN
324 update = .false.
325 i = i + 1
326 i0 = m
327 IF( ka.EQ.0 )
328 $ GO TO 480
329 GO TO 10
330 END IF
331 ELSE
332 i = i + ka
333 IF( i.GT.n-1 )
334 $ GO TO 480
335 END IF
336
337 IF( upper ) THEN
338
339
340
341 IF( update ) THEN
342
343
344
345 bii = dble( bb( kb1, i ) )
346 ab( ka1, i ) = ( dble( ab( ka1, i ) ) / bii ) / bii
347 DO 20 j = i + 1, i1
348 ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
349 20 CONTINUE
350 DO 30 j = max( 1, i-ka ), i - 1
351 ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
352 30 CONTINUE
353 DO 60 k = i - kbt, i - 1
354 DO 40 j = i - kbt, k
355 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
356 $ bb( j-i+kb1, i )*
357 $ dconjg( ab( k-i+ka1, i ) ) -
358 $ dconjg( bb( k-i+kb1, i ) )*
359 $ ab( j-i+ka1, i ) +
360 $ dble( ab( ka1, i ) )*
361 $ bb( j-i+kb1, i )*
362 $ dconjg( bb( k-i+kb1, i ) )
363 40 CONTINUE
364 DO 50 j = max( 1, i-ka ), i - kbt - 1
365 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
366 $ dconjg( bb( k-i+kb1, i ) )*
367 $ ab( j-i+ka1, i )
368 50 CONTINUE
369 60 CONTINUE
370 DO 80 j = i, i1
371 DO 70 k = max( j-ka, i-kbt ), i - 1
372 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
373 $ bb( k-i+kb1, i )*ab( i-j+ka1, j )
374 70 CONTINUE
375 80 CONTINUE
376
377 IF( wantx ) THEN
378
379
380
381 CALL zdscal( n-m, one / bii, x( m+1, i ), 1 )
382 IF( kbt.GT.0 )
383 $
CALL zgerc( n-m, kbt, -cone, x( m+1, i ), 1,
384 $ bb( kb1-kbt, i ), 1, x( m+1, i-kbt ),
385 $ ldx )
386 END IF
387
388
389
390 ra1 = ab( i-i1+ka1, i1 )
391 END IF
392
393
394
395
396
397 DO 130 k = 1, kb - 1
398 IF( update ) THEN
399
400
401
402
403 IF( i-k+ka.LT.n .AND. i-k.GT.1 ) THEN
404
405
406
407 CALL zlartg( ab( k+1, i-k+ka ), ra1,
408 $ rwork( i-k+ka-m ), work( i-k+ka-m ), ra )
409
410
411
412
413 t = -bb( kb1-k, i )*ra1
414 work( i-k ) = rwork( i-k+ka-m )*t -
415 $ dconjg( work( i-k+ka-m ) )*
416 $ ab( 1, i-k+ka )
417 ab( 1, i-k+ka ) = work( i-k+ka-m )*t +
418 $ rwork( i-k+ka-m )*ab( 1, i-k+ka )
419 ra1 = ra
420 END IF
421 END IF
422 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
423 nr = ( n-j2+ka ) / ka1
424 j1 = j2 + ( nr-1 )*ka1
425 IF( update ) THEN
426 j2t = max( j2, i+2*ka-k+1 )
427 ELSE
428 j2t = j2
429 END IF
430 nrt = ( n-j2t+ka ) / ka1
431 DO 90 j = j2t, j1, ka1
432
433
434
435
436 work( j-m ) = work( j-m )*ab( 1, j+1 )
437 ab( 1, j+1 ) = rwork( j-m )*ab( 1, j+1 )
438 90 CONTINUE
439
440
441
442
443 IF( nrt.GT.0 )
444 $
CALL zlargv( nrt, ab( 1, j2t ), inca, work( j2t-m ),
445 $ ka1,
446 $ rwork( j2t-m ), ka1 )
447 IF( nr.GT.0 ) THEN
448
449
450
451 DO 100 l = 1, ka - 1
452 CALL zlartv( nr, ab( ka1-l, j2 ), inca,
453 $ ab( ka-l, j2+1 ), inca, rwork( j2-m ),
454 $ work( j2-m ), ka1 )
455 100 CONTINUE
456
457
458
459
460 CALL zlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
461 $ ab( ka, j2+1 ), inca, rwork( j2-m ),
462 $ work( j2-m ), ka1 )
463
464 CALL zlacgv( nr, work( j2-m ), ka1 )
465 END IF
466
467
468
469 DO 110 l = ka - 1, kb - k + 1, -1
470 nrt = ( n-j2+l ) / ka1
471 IF( nrt.GT.0 )
472 $
CALL zlartv( nrt, ab( l, j2+ka1-l ), inca,
473 $ ab( l+1, j2+ka1-l ), inca, rwork( j2-m ),
474 $ work( j2-m ), ka1 )
475 110 CONTINUE
476
477 IF( wantx ) THEN
478
479
480
481 DO 120 j = j2, j1, ka1
482 CALL zrot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
483 $ rwork( j-m ), dconjg( work( j-m ) ) )
484 120 CONTINUE
485 END IF
486 130 CONTINUE
487
488 IF( update ) THEN
489 IF( i2.LE.n .AND. kbt.GT.0 ) THEN
490
491
492
493
494 work( i-kbt ) = -bb( kb1-kbt, i )*ra1
495 END IF
496 END IF
497
498 DO 170 k = kb, 1, -1
499 IF( update ) THEN
500 j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
501 ELSE
502 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
503 END IF
504
505
506
507 DO 140 l = kb - k, 1, -1
508 nrt = ( n-j2+ka+l ) / ka1
509 IF( nrt.GT.0 )
510 $
CALL zlartv( nrt, ab( l, j2-l+1 ), inca,
511 $ ab( l+1, j2-l+1 ), inca, rwork( j2-ka ),
512 $ work( j2-ka ), ka1 )
513 140 CONTINUE
514 nr = ( n-j2+ka ) / ka1
515 j1 = j2 + ( nr-1 )*ka1
516 DO 150 j = j1, j2, -ka1
517 work( j ) = work( j-ka )
518 rwork( j ) = rwork( j-ka )
519 150 CONTINUE
520 DO 160 j = j2, j1, ka1
521
522
523
524
525 work( j ) = work( j )*ab( 1, j+1 )
526 ab( 1, j+1 ) = rwork( j )*ab( 1, j+1 )
527 160 CONTINUE
528 IF( update ) THEN
529 IF( i-k.LT.n-ka .AND. k.LE.kbt )
530 $ work( i-k+ka ) = work( i-k )
531 END IF
532 170 CONTINUE
533
534 DO 210 k = kb, 1, -1
535 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
536 nr = ( n-j2+ka ) / ka1
537 j1 = j2 + ( nr-1 )*ka1
538 IF( nr.GT.0 ) THEN
539
540
541
542
543 CALL zlargv( nr, ab( 1, j2 ), inca, work( j2 ), ka1,
544 $ rwork( j2 ), ka1 )
545
546
547
548 DO 180 l = 1, ka - 1
549 CALL zlartv( nr, ab( ka1-l, j2 ), inca,
550 $ ab( ka-l, j2+1 ), inca, rwork( j2 ),
551 $ work( j2 ), ka1 )
552 180 CONTINUE
553
554
555
556
557 CALL zlar2v( nr, ab( ka1, j2 ), ab( ka1, j2+1 ),
558 $ ab( ka, j2+1 ), inca, rwork( j2 ),
559 $ work( j2 ), ka1 )
560
561 CALL zlacgv( nr, work( j2 ), ka1 )
562 END IF
563
564
565
566 DO 190 l = ka - 1, kb - k + 1, -1
567 nrt = ( n-j2+l ) / ka1
568 IF( nrt.GT.0 )
569 $
CALL zlartv( nrt, ab( l, j2+ka1-l ), inca,
570 $ ab( l+1, j2+ka1-l ), inca, rwork( j2 ),
571 $ work( j2 ), ka1 )
572 190 CONTINUE
573
574 IF( wantx ) THEN
575
576
577
578 DO 200 j = j2, j1, ka1
579 CALL zrot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
580 $ rwork( j ), dconjg( work( j ) ) )
581 200 CONTINUE
582 END IF
583 210 CONTINUE
584
585 DO 230 k = 1, kb - 1
586 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
587
588
589
590 DO 220 l = kb - k, 1, -1
591 nrt = ( n-j2+l ) / ka1
592 IF( nrt.GT.0 )
593 $
CALL zlartv( nrt, ab( l, j2+ka1-l ), inca,
594 $ ab( l+1, j2+ka1-l ), inca, rwork( j2-m ),
595 $ work( j2-m ), ka1 )
596 220 CONTINUE
597 230 CONTINUE
598
599 IF( kb.GT.1 ) THEN
600 DO 240 j = n - 1, j2 + ka, -1
601 rwork( j-m ) = rwork( j-ka-m )
602 work( j-m ) = work( j-ka-m )
603 240 CONTINUE
604 END IF
605
606 ELSE
607
608
609
610 IF( update ) THEN
611
612
613
614 bii = dble( bb( 1, i ) )
615 ab( 1, i ) = ( dble( ab( 1, i ) ) / bii ) / bii
616 DO 250 j = i + 1, i1
617 ab( j-i+1, i ) = ab( j-i+1, i ) / bii
618 250 CONTINUE
619 DO 260 j = max( 1, i-ka ), i - 1
620 ab( i-j+1, j ) = ab( i-j+1, j ) / bii
621 260 CONTINUE
622 DO 290 k = i - kbt, i - 1
623 DO 270 j = i - kbt, k
624 ab( k-j+1, j ) = ab( k-j+1, j ) -
625 $ bb( i-j+1, j )*dconjg( ab( i-k+1,
626 $ k ) ) - dconjg( bb( i-k+1, k ) )*
627 $ ab( i-j+1, j ) + dble( ab( 1, i ) )*
628 $ bb( i-j+1, j )*dconjg( bb( i-k+1,
629 $ k ) )
630 270 CONTINUE
631 DO 280 j = max( 1, i-ka ), i - kbt - 1
632 ab( k-j+1, j ) = ab( k-j+1, j ) -
633 $ dconjg( bb( i-k+1, k ) )*
634 $ ab( i-j+1, j )
635 280 CONTINUE
636 290 CONTINUE
637 DO 310 j = i, i1
638 DO 300 k = max( j-ka, i-kbt ), i - 1
639 ab( j-k+1, k ) = ab( j-k+1, k ) -
640 $ bb( i-k+1, k )*ab( j-i+1, i )
641 300 CONTINUE
642 310 CONTINUE
643
644 IF( wantx ) THEN
645
646
647
648 CALL zdscal( n-m, one / bii, x( m+1, i ), 1 )
649 IF( kbt.GT.0 )
650 $
CALL zgeru( n-m, kbt, -cone, x( m+1, i ), 1,
651 $ bb( kbt+1, i-kbt ), ldbb-1,
652 $ x( m+1, i-kbt ), ldx )
653 END IF
654
655
656
657 ra1 = ab( i1-i+1, i )
658 END IF
659
660
661
662
663
664 DO 360 k = 1, kb - 1
665 IF( update ) THEN
666
667
668
669
670 IF( i-k+ka.LT.n .AND. i-k.GT.1 ) THEN
671
672
673
674 CALL zlartg( ab( ka1-k, i ), ra1,
675 $ rwork( i-k+ka-m ),
676 $ work( i-k+ka-m ), ra )
677
678
679
680
681 t = -bb( k+1, i-k )*ra1
682 work( i-k ) = rwork( i-k+ka-m )*t -
683 $ dconjg( work( i-k+ka-m ) )*
684 $ ab( ka1, i-k )
685 ab( ka1, i-k ) = work( i-k+ka-m )*t +
686 $ rwork( i-k+ka-m )*ab( ka1, i-k )
687 ra1 = ra
688 END IF
689 END IF
690 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
691 nr = ( n-j2+ka ) / ka1
692 j1 = j2 + ( nr-1 )*ka1
693 IF( update ) THEN
694 j2t = max( j2, i+2*ka-k+1 )
695 ELSE
696 j2t = j2
697 END IF
698 nrt = ( n-j2t+ka ) / ka1
699 DO 320 j = j2t, j1, ka1
700
701
702
703
704 work( j-m ) = work( j-m )*ab( ka1, j-ka+1 )
705 ab( ka1, j-ka+1 ) = rwork( j-m )*ab( ka1, j-ka+1 )
706 320 CONTINUE
707
708
709
710
711 IF( nrt.GT.0 )
712 $
CALL zlargv( nrt, ab( ka1, j2t-ka ), inca,
713 $ work( j2t-m ),
714 $ ka1, rwork( j2t-m ), ka1 )
715 IF( nr.GT.0 ) THEN
716
717
718
719 DO 330 l = 1, ka - 1
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 )
723 330 CONTINUE
724
725
726
727
728 CALL zlar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2,
729 $ j2 ),
730 $ inca, rwork( j2-m ), work( j2-m ), ka1 )
731
732 CALL zlacgv( nr, work( j2-m ), ka1 )
733 END IF
734
735
736
737 DO 340 l = ka - 1, kb - k + 1, -1
738 nrt = ( n-j2+l ) / ka1
739 IF( nrt.GT.0 )
740 $
CALL zlartv( nrt, ab( ka1-l+1, j2 ), inca,
741 $ ab( ka1-l, j2+1 ), inca, rwork( j2-m ),
742 $ work( j2-m ), ka1 )
743 340 CONTINUE
744
745 IF( wantx ) THEN
746
747
748
749 DO 350 j = j2, j1, ka1
750 CALL zrot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
751 $ rwork( j-m ), work( j-m ) )
752 350 CONTINUE
753 END IF
754 360 CONTINUE
755
756 IF( update ) THEN
757 IF( i2.LE.n .AND. kbt.GT.0 ) THEN
758
759
760
761
762 work( i-kbt ) = -bb( kbt+1, i-kbt )*ra1
763 END IF
764 END IF
765
766 DO 400 k = kb, 1, -1
767 IF( update ) THEN
768 j2 = i - k - 1 + max( 2, k-i0+1 )*ka1
769 ELSE
770 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
771 END IF
772
773
774
775 DO 370 l = kb - k, 1, -1
776 nrt = ( n-j2+ka+l ) / ka1
777 IF( nrt.GT.0 )
778 $
CALL zlartv( nrt, ab( ka1-l+1, j2-ka ), inca,
779 $ ab( ka1-l, j2-ka+1 ), inca,
780 $ rwork( j2-ka ), work( j2-ka ), ka1 )
781 370 CONTINUE
782 nr = ( n-j2+ka ) / ka1
783 j1 = j2 + ( nr-1 )*ka1
784 DO 380 j = j1, j2, -ka1
785 work( j ) = work( j-ka )
786 rwork( j ) = rwork( j-ka )
787 380 CONTINUE
788 DO 390 j = j2, j1, ka1
789
790
791
792
793 work( j ) = work( j )*ab( ka1, j-ka+1 )
794 ab( ka1, j-ka+1 ) = rwork( j )*ab( ka1, j-ka+1 )
795 390 CONTINUE
796 IF( update ) THEN
797 IF( i-k.LT.n-ka .AND. k.LE.kbt )
798 $ work( i-k+ka ) = work( i-k )
799 END IF
800 400 CONTINUE
801
802 DO 440 k = kb, 1, -1
803 j2 = i - k - 1 + max( 1, k-i0+1 )*ka1
804 nr = ( n-j2+ka ) / ka1
805 j1 = j2 + ( nr-1 )*ka1
806 IF( nr.GT.0 ) THEN
807
808
809
810
811 CALL zlargv( nr, ab( ka1, j2-ka ), inca, work( j2 ),
812 $ ka1,
813 $ rwork( j2 ), ka1 )
814
815
816
817 DO 410 l = 1, ka - 1
818 CALL zlartv( nr, ab( l+1, j2-l ), inca,
819 $ ab( l+2, j2-l ), inca, rwork( j2 ),
820 $ work( j2 ), ka1 )
821 410 CONTINUE
822
823
824
825
826 CALL zlar2v( nr, ab( 1, j2 ), ab( 1, j2+1 ), ab( 2,
827 $ j2 ),
828 $ inca, rwork( j2 ), work( j2 ), ka1 )
829
830 CALL zlacgv( nr, work( j2 ), ka1 )
831 END IF
832
833
834
835 DO 420 l = ka - 1, kb - k + 1, -1
836 nrt = ( n-j2+l ) / ka1
837 IF( nrt.GT.0 )
838 $
CALL zlartv( nrt, ab( ka1-l+1, j2 ), inca,
839 $ ab( ka1-l, j2+1 ), inca, rwork( j2 ),
840 $ work( j2 ), ka1 )
841 420 CONTINUE
842
843 IF( wantx ) THEN
844
845
846
847 DO 430 j = j2, j1, ka1
848 CALL zrot( n-m, x( m+1, j ), 1, x( m+1, j+1 ), 1,
849 $ rwork( j ), work( j ) )
850 430 CONTINUE
851 END IF
852 440 CONTINUE
853
854 DO 460 k = 1, kb - 1
855 j2 = i - k - 1 + max( 1, k-i0+2 )*ka1
856
857
858
859 DO 450 l = kb - k, 1, -1
860 nrt = ( n-j2+l ) / ka1
861 IF( nrt.GT.0 )
862 $
CALL zlartv( nrt, ab( ka1-l+1, j2 ), inca,
863 $ ab( ka1-l, j2+1 ), inca, rwork( j2-m ),
864 $ work( j2-m ), ka1 )
865 450 CONTINUE
866 460 CONTINUE
867
868 IF( kb.GT.1 ) THEN
869 DO 470 j = n - 1, j2 + ka, -1
870 rwork( j-m ) = rwork( j-ka-m )
871 work( j-m ) = work( j-ka-m )
872 470 CONTINUE
873 END IF
874
875 END IF
876
877 GO TO 10
878
879 480 CONTINUE
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897 update = .true.
898 i = 0
899 490 CONTINUE
900 IF( update ) THEN
901 i = i + 1
902 kbt = min( kb, m-i )
903 i0 = i + 1
904 i1 = max( 1, i-ka )
905 i2 = i + kbt - ka1
906 IF( i.GT.m ) THEN
907 update = .false.
908 i = i - 1
909 i0 = m + 1
910 IF( ka.EQ.0 )
911 $ RETURN
912 GO TO 490
913 END IF
914 ELSE
915 i = i - ka
916 IF( i.LT.2 )
917 $ RETURN
918 END IF
919
920 IF( i.LT.m-kbt ) THEN
921 nx = m
922 ELSE
923 nx = n
924 END IF
925
926 IF( upper ) THEN
927
928
929
930 IF( update ) THEN
931
932
933
934 bii = dble( bb( kb1, i ) )
935 ab( ka1, i ) = ( dble( ab( ka1, i ) ) / bii ) / bii
936 DO 500 j = i1, i - 1
937 ab( j-i+ka1, i ) = ab( j-i+ka1, i ) / bii
938 500 CONTINUE
939 DO 510 j = i + 1, min( n, i+ka )
940 ab( i-j+ka1, j ) = ab( i-j+ka1, j ) / bii
941 510 CONTINUE
942 DO 540 k = i + 1, i + kbt
943 DO 520 j = k, i + kbt
944 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
945 $ bb( i-j+kb1, j )*
946 $ dconjg( ab( i-k+ka1, k ) ) -
947 $ dconjg( bb( i-k+kb1, k ) )*
948 $ ab( i-j+ka1, j ) +
949 $ dble( ab( ka1, i ) )*
950 $ bb( i-j+kb1, j )*
951 $ dconjg( bb( i-k+kb1, k ) )
952 520 CONTINUE
953 DO 530 j = i + kbt + 1, min( n, i+ka )
954 ab( k-j+ka1, j ) = ab( k-j+ka1, j ) -
955 $ dconjg( bb( i-k+kb1, k ) )*
956 $ ab( i-j+ka1, j )
957 530 CONTINUE
958 540 CONTINUE
959 DO 560 j = i1, i
960 DO 550 k = i + 1, min( j+ka, i+kbt )
961 ab( j-k+ka1, k ) = ab( j-k+ka1, k ) -
962 $ bb( i-k+kb1, k )*ab( j-i+ka1, i )
963 550 CONTINUE
964 560 CONTINUE
965
966 IF( wantx ) THEN
967
968
969
970 CALL zdscal( nx, one / bii, x( 1, i ), 1 )
971 IF( kbt.GT.0 )
972 $
CALL zgeru( nx, kbt, -cone, x( 1, i ), 1,
973 $ bb( kb, i+1 ), ldbb-1, x( 1, i+1 ), ldx )
974 END IF
975
976
977
978 ra1 = ab( i1-i+ka1, i )
979 END IF
980
981
982
983
984 DO 610 k = 1, kb - 1
985 IF( update ) THEN
986
987
988
989
990 IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
991
992
993
994 CALL zlartg( ab( k+1, i ), ra1, rwork( i+k-ka ),
995 $ work( i+k-ka ), ra )
996
997
998
999
1000 t = -bb( kb1-k, i+k )*ra1
1001 work( m-kb+i+k ) = rwork( i+k-ka )*t -
1002 $ dconjg( work( i+k-ka ) )*
1003 $ ab( 1, i+k )
1004 ab( 1, i+k ) = work( i+k-ka )*t +
1005 $ rwork( i+k-ka )*ab( 1, i+k )
1006 ra1 = ra
1007 END IF
1008 END IF
1009 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1010 nr = ( j2+ka-1 ) / ka1
1011 j1 = j2 - ( nr-1 )*ka1
1012 IF( update ) THEN
1013 j2t = min( j2, i-2*ka+k-1 )
1014 ELSE
1015 j2t = j2
1016 END IF
1017 nrt = ( j2t+ka-1 ) / ka1
1018 DO 570 j = j1, j2t, ka1
1019
1020
1021
1022
1023 work( j ) = work( j )*ab( 1, j+ka-1 )
1024 ab( 1, j+ka-1 ) = rwork( j )*ab( 1, j+ka-1 )
1025 570 CONTINUE
1026
1027
1028
1029
1030 IF( nrt.GT.0 )
1031 $
CALL zlargv( nrt, ab( 1, j1+ka ), inca, work( j1 ),
1032 $ ka1,
1033 $ rwork( j1 ), ka1 )
1034 IF( nr.GT.0 ) THEN
1035
1036
1037
1038 DO 580 l = 1, ka - 1
1039 CALL zlartv( nr, ab( ka1-l, j1+l ), inca,
1040 $ ab( ka-l, j1+l ), inca, rwork( j1 ),
1041 $ work( j1 ), ka1 )
1042 580 CONTINUE
1043
1044
1045
1046
1047 CALL zlar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1048 $ ab( ka, j1 ), inca, rwork( j1 ), work( j1 ),
1049 $ ka1 )
1050
1051 CALL zlacgv( nr, work( j1 ), ka1 )
1052 END IF
1053
1054
1055
1056 DO 590 l = ka - 1, kb - k + 1, -1
1057 nrt = ( j2+l-1 ) / ka1
1058 j1t = j2 - ( nrt-1 )*ka1
1059 IF( nrt.GT.0 )
1060 $
CALL zlartv( nrt, ab( l, j1t ), inca,
1061 $ ab( l+1, j1t-1 ), inca, rwork( j1t ),
1062 $ work( j1t ), ka1 )
1063 590 CONTINUE
1064
1065 IF( wantx ) THEN
1066
1067
1068
1069 DO 600 j = j1, j2, ka1
1070 CALL zrot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1071 $ rwork( j ), work( j ) )
1072 600 CONTINUE
1073 END IF
1074 610 CONTINUE
1075
1076 IF( update ) THEN
1077 IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1078
1079
1080
1081
1082 work( m-kb+i+kbt ) = -bb( kb1-kbt, i+kbt )*ra1
1083 END IF
1084 END IF
1085
1086 DO 650 k = kb, 1, -1
1087 IF( update ) THEN
1088 j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1089 ELSE
1090 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1091 END IF
1092
1093
1094
1095 DO 620 l = kb - k, 1, -1
1096 nrt = ( j2+ka+l-1 ) / ka1
1097 j1t = j2 - ( nrt-1 )*ka1
1098 IF( nrt.GT.0 )
1099 $
CALL zlartv( nrt, ab( l, j1t+ka ), inca,
1100 $ ab( l+1, j1t+ka-1 ), inca,
1101 $ rwork( m-kb+j1t+ka ),
1102 $ work( m-kb+j1t+ka ), ka1 )
1103 620 CONTINUE
1104 nr = ( j2+ka-1 ) / ka1
1105 j1 = j2 - ( nr-1 )*ka1
1106 DO 630 j = j1, j2, ka1
1107 work( m-kb+j ) = work( m-kb+j+ka )
1108 rwork( m-kb+j ) = rwork( m-kb+j+ka )
1109 630 CONTINUE
1110 DO 640 j = j1, j2, ka1
1111
1112
1113
1114
1115 work( m-kb+j ) = work( m-kb+j )*ab( 1, j+ka-1 )
1116 ab( 1, j+ka-1 ) = rwork( m-kb+j )*ab( 1, j+ka-1 )
1117 640 CONTINUE
1118 IF( update ) THEN
1119 IF( i+k.GT.ka1 .AND. k.LE.kbt )
1120 $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1121 END IF
1122 650 CONTINUE
1123
1124 DO 690 k = kb, 1, -1
1125 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1126 nr = ( j2+ka-1 ) / ka1
1127 j1 = j2 - ( nr-1 )*ka1
1128 IF( nr.GT.0 ) THEN
1129
1130
1131
1132
1133 CALL zlargv( nr, ab( 1, j1+ka ), inca,
1134 $ work( m-kb+j1 ),
1135 $ ka1, rwork( m-kb+j1 ), ka1 )
1136
1137
1138
1139 DO 660 l = 1, ka - 1
1140 CALL zlartv( nr, ab( ka1-l, j1+l ), inca,
1141 $ ab( ka-l, j1+l ), inca, rwork( m-kb+j1 ),
1142 $ work( m-kb+j1 ), ka1 )
1143 660 CONTINUE
1144
1145
1146
1147
1148 CALL zlar2v( nr, ab( ka1, j1 ), ab( ka1, j1-1 ),
1149 $ ab( ka, j1 ), inca, rwork( m-kb+j1 ),
1150 $ work( m-kb+j1 ), ka1 )
1151
1152 CALL zlacgv( nr, work( m-kb+j1 ), ka1 )
1153 END IF
1154
1155
1156
1157 DO 670 l = ka - 1, kb - k + 1, -1
1158 nrt = ( j2+l-1 ) / ka1
1159 j1t = j2 - ( nrt-1 )*ka1
1160 IF( nrt.GT.0 )
1161 $
CALL zlartv( nrt, ab( l, j1t ), inca,
1162 $ ab( l+1, j1t-1 ), inca,
1163 $ rwork( m-kb+j1t ), work( m-kb+j1t ),
1164 $ ka1 )
1165 670 CONTINUE
1166
1167 IF( wantx ) THEN
1168
1169
1170
1171 DO 680 j = j1, j2, ka1
1172 CALL zrot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1173 $ rwork( m-kb+j ), work( m-kb+j ) )
1174 680 CONTINUE
1175 END IF
1176 690 CONTINUE
1177
1178 DO 710 k = 1, kb - 1
1179 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1180
1181
1182
1183 DO 700 l = kb - k, 1, -1
1184 nrt = ( j2+l-1 ) / ka1
1185 j1t = j2 - ( nrt-1 )*ka1
1186 IF( nrt.GT.0 )
1187 $
CALL zlartv( nrt, ab( l, j1t ), inca,
1188 $ ab( l+1, j1t-1 ), inca, rwork( j1t ),
1189 $ work( j1t ), ka1 )
1190 700 CONTINUE
1191 710 CONTINUE
1192
1193 IF( kb.GT.1 ) THEN
1194 DO 720 j = 2, i2 - ka
1195 rwork( j ) = rwork( j+ka )
1196 work( j ) = work( j+ka )
1197 720 CONTINUE
1198 END IF
1199
1200 ELSE
1201
1202
1203
1204 IF( update ) THEN
1205
1206
1207
1208 bii = dble( bb( 1, i ) )
1209 ab( 1, i ) = ( dble( ab( 1, i ) ) / bii ) / bii
1210 DO 730 j = i1, i - 1
1211 ab( i-j+1, j ) = ab( i-j+1, j ) / bii
1212 730 CONTINUE
1213 DO 740 j = i + 1, min( n, i+ka )
1214 ab( j-i+1, i ) = ab( j-i+1, i ) / bii
1215 740 CONTINUE
1216 DO 770 k = i + 1, i + kbt
1217 DO 750 j = k, i + kbt
1218 ab( j-k+1, k ) = ab( j-k+1, k ) -
1219 $ bb( j-i+1, i )*dconjg( ab( k-i+1,
1220 $ i ) ) - dconjg( bb( k-i+1, i ) )*
1221 $ ab( j-i+1, i ) + dble( ab( 1, i ) )*
1222 $ bb( j-i+1, i )*dconjg( bb( k-i+1,
1223 $ i ) )
1224 750 CONTINUE
1225 DO 760 j = i + kbt + 1, min( n, i+ka )
1226 ab( j-k+1, k ) = ab( j-k+1, k ) -
1227 $ dconjg( bb( k-i+1, i ) )*
1228 $ ab( j-i+1, i )
1229 760 CONTINUE
1230 770 CONTINUE
1231 DO 790 j = i1, i
1232 DO 780 k = i + 1, min( j+ka, i+kbt )
1233 ab( k-j+1, j ) = ab( k-j+1, j ) -
1234 $ bb( k-i+1, i )*ab( i-j+1, j )
1235 780 CONTINUE
1236 790 CONTINUE
1237
1238 IF( wantx ) THEN
1239
1240
1241
1242 CALL zdscal( nx, one / bii, x( 1, i ), 1 )
1243 IF( kbt.GT.0 )
1244 $
CALL zgerc( nx, kbt, -cone, x( 1, i ), 1, bb( 2,
1245 $ i ),
1246 $ 1, x( 1, i+1 ), ldx )
1247 END IF
1248
1249
1250
1251 ra1 = ab( i-i1+1, i1 )
1252 END IF
1253
1254
1255
1256
1257 DO 840 k = 1, kb - 1
1258 IF( update ) THEN
1259
1260
1261
1262
1263 IF( i+k-ka1.GT.0 .AND. i+k.LT.m ) THEN
1264
1265
1266
1267 CALL zlartg( ab( ka1-k, i+k-ka ), ra1,
1268 $ rwork( i+k-ka ), work( i+k-ka ), ra )
1269
1270
1271
1272
1273 t = -bb( k+1, i )*ra1
1274 work( m-kb+i+k ) = rwork( i+k-ka )*t -
1275 $ dconjg( work( i+k-ka ) )*
1276 $ ab( ka1, i+k-ka )
1277 ab( ka1, i+k-ka ) = work( i+k-ka )*t +
1278 $ rwork( i+k-ka )*ab( ka1, i+k-ka )
1279 ra1 = ra
1280 END IF
1281 END IF
1282 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1283 nr = ( j2+ka-1 ) / ka1
1284 j1 = j2 - ( nr-1 )*ka1
1285 IF( update ) THEN
1286 j2t = min( j2, i-2*ka+k-1 )
1287 ELSE
1288 j2t = j2
1289 END IF
1290 nrt = ( j2t+ka-1 ) / ka1
1291 DO 800 j = j1, j2t, ka1
1292
1293
1294
1295
1296 work( j ) = work( j )*ab( ka1, j-1 )
1297 ab( ka1, j-1 ) = rwork( j )*ab( ka1, j-1 )
1298 800 CONTINUE
1299
1300
1301
1302
1303 IF( nrt.GT.0 )
1304 $
CALL zlargv( nrt, ab( ka1, j1 ), inca, work( j1 ),
1305 $ ka1,
1306 $ rwork( j1 ), ka1 )
1307 IF( nr.GT.0 ) THEN
1308
1309
1310
1311 DO 810 l = 1, ka - 1
1312 CALL zlartv( nr, ab( l+1, j1 ), inca, ab( l+2,
1313 $ j1-1 ),
1314 $ inca, rwork( j1 ), work( j1 ), ka1 )
1315 810 CONTINUE
1316
1317
1318
1319
1320 CALL zlar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1321 $ ab( 2, j1-1 ), inca, rwork( j1 ),
1322 $ work( j1 ), ka1 )
1323
1324 CALL zlacgv( nr, work( j1 ), ka1 )
1325 END IF
1326
1327
1328
1329 DO 820 l = ka - 1, kb - k + 1, -1
1330 nrt = ( j2+l-1 ) / ka1
1331 j1t = j2 - ( nrt-1 )*ka1
1332 IF( nrt.GT.0 )
1333 $
CALL zlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1334 $ ab( ka1-l, j1t-ka1+l ), inca,
1335 $ rwork( j1t ), work( j1t ), ka1 )
1336 820 CONTINUE
1337
1338 IF( wantx ) THEN
1339
1340
1341
1342 DO 830 j = j1, j2, ka1
1343 CALL zrot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1344 $ rwork( j ), dconjg( work( j ) ) )
1345 830 CONTINUE
1346 END IF
1347 840 CONTINUE
1348
1349 IF( update ) THEN
1350 IF( i2.GT.0 .AND. kbt.GT.0 ) THEN
1351
1352
1353
1354
1355 work( m-kb+i+kbt ) = -bb( kbt+1, i )*ra1
1356 END IF
1357 END IF
1358
1359 DO 880 k = kb, 1, -1
1360 IF( update ) THEN
1361 j2 = i + k + 1 - max( 2, k+i0-m )*ka1
1362 ELSE
1363 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1364 END IF
1365
1366
1367
1368 DO 850 l = kb - k, 1, -1
1369 nrt = ( j2+ka+l-1 ) / ka1
1370 j1t = j2 - ( nrt-1 )*ka1
1371 IF( nrt.GT.0 )
1372 $
CALL zlartv( nrt, ab( ka1-l+1, j1t+l-1 ), inca,
1373 $ ab( ka1-l, j1t+l-1 ), inca,
1374 $ rwork( m-kb+j1t+ka ),
1375 $ work( m-kb+j1t+ka ), ka1 )
1376 850 CONTINUE
1377 nr = ( j2+ka-1 ) / ka1
1378 j1 = j2 - ( nr-1 )*ka1
1379 DO 860 j = j1, j2, ka1
1380 work( m-kb+j ) = work( m-kb+j+ka )
1381 rwork( m-kb+j ) = rwork( m-kb+j+ka )
1382 860 CONTINUE
1383 DO 870 j = j1, j2, ka1
1384
1385
1386
1387
1388 work( m-kb+j ) = work( m-kb+j )*ab( ka1, j-1 )
1389 ab( ka1, j-1 ) = rwork( m-kb+j )*ab( ka1, j-1 )
1390 870 CONTINUE
1391 IF( update ) THEN
1392 IF( i+k.GT.ka1 .AND. k.LE.kbt )
1393 $ work( m-kb+i+k-ka ) = work( m-kb+i+k )
1394 END IF
1395 880 CONTINUE
1396
1397 DO 920 k = kb, 1, -1
1398 j2 = i + k + 1 - max( 1, k+i0-m )*ka1
1399 nr = ( j2+ka-1 ) / ka1
1400 j1 = j2 - ( nr-1 )*ka1
1401 IF( nr.GT.0 ) THEN
1402
1403
1404
1405
1406 CALL zlargv( nr, ab( ka1, j1 ), inca, work( m-kb+j1 ),
1407 $ ka1, rwork( m-kb+j1 ), ka1 )
1408
1409
1410
1411 DO 890 l = 1, ka - 1
1412 CALL zlartv( nr, ab( l+1, j1 ), inca, ab( l+2,
1413 $ j1-1 ),
1414 $ inca, rwork( m-kb+j1 ), work( m-kb+j1 ),
1415 $ ka1 )
1416 890 CONTINUE
1417
1418
1419
1420
1421 CALL zlar2v( nr, ab( 1, j1 ), ab( 1, j1-1 ),
1422 $ ab( 2, j1-1 ), inca, rwork( m-kb+j1 ),
1423 $ work( m-kb+j1 ), ka1 )
1424
1425 CALL zlacgv( nr, work( m-kb+j1 ), ka1 )
1426 END IF
1427
1428
1429
1430 DO 900 l = ka - 1, kb - k + 1, -1
1431 nrt = ( j2+l-1 ) / ka1
1432 j1t = j2 - ( nrt-1 )*ka1
1433 IF( nrt.GT.0 )
1434 $
CALL zlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1435 $ ab( ka1-l, j1t-ka1+l ), inca,
1436 $ rwork( m-kb+j1t ), work( m-kb+j1t ),
1437 $ ka1 )
1438 900 CONTINUE
1439
1440 IF( wantx ) THEN
1441
1442
1443
1444 DO 910 j = j1, j2, ka1
1445 CALL zrot( nx, x( 1, j ), 1, x( 1, j-1 ), 1,
1446 $ rwork( m-kb+j ), dconjg( work( m-kb+j ) ) )
1447 910 CONTINUE
1448 END IF
1449 920 CONTINUE
1450
1451 DO 940 k = 1, kb - 1
1452 j2 = i + k + 1 - max( 1, k+i0-m+1 )*ka1
1453
1454
1455
1456 DO 930 l = kb - k, 1, -1
1457 nrt = ( j2+l-1 ) / ka1
1458 j1t = j2 - ( nrt-1 )*ka1
1459 IF( nrt.GT.0 )
1460 $
CALL zlartv( nrt, ab( ka1-l+1, j1t-ka1+l ), inca,
1461 $ ab( ka1-l, j1t-ka1+l ), inca,
1462 $ rwork( j1t ), work( j1t ), ka1 )
1463 930 CONTINUE
1464 940 CONTINUE
1465
1466 IF( kb.GT.1 ) THEN
1467 DO 950 j = 2, i2 - ka
1468 rwork( j ) = rwork( j+ka )
1469 work( j ) = work( j+ka )
1470 950 CONTINUE
1471 END IF
1472
1473 END IF
1474
1475 GO TO 490
1476
1477
1478
subroutine xerbla(srname, info)
subroutine zgeru(m, n, alpha, x, incx, y, incy, a, lda)
ZGERU
subroutine zgerc(m, n, alpha, x, incx, y, incy, a, lda)
ZGERC
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
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 zlargv(n, x, incx, y, incy, c, incc)
ZLARGV generates a vector of plane rotations with real cosines and complex sines.
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
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 zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME
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 zdscal(n, da, zx, incx)
ZDSCAL