284
285
286
287
288
289
290 CHARACTER COMPQ, COMPZ, JOB
291 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
292
293
294 DOUBLE PRECISION RWORK( * )
295 COMPLEX*16 ALPHA( * ), BETA( * ), H( LDH, * ),
296 $ Q( LDQ, * ), T( LDT, * ), WORK( * ),
297 $ Z( LDZ, * )
298
299
300
301
302
303 COMPLEX*16 CZERO, CONE
304 parameter( czero = ( 0.0d+0, 0.0d+0 ),
305 $ cone = ( 1.0d+0, 0.0d+0 ) )
306 DOUBLE PRECISION ZERO, ONE
307 parameter( zero = 0.0d+0, one = 1.0d+0 )
308 DOUBLE PRECISION HALF
309 parameter( half = 0.5d+0 )
310
311
312 LOGICAL ILAZR2, ILAZRO, ILQ, ILSCHR, ILZ, LQUERY
313 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
314 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
315 $ JR, MAXIT
316 DOUBLE PRECISION ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
317 $ C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
318 COMPLEX*16 ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
319 $ CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
320 $ U12, X, ABI12, Y
321
322
323 COMPLEX*16 ZLADIV
324 LOGICAL LSAME
325 DOUBLE PRECISION DLAMCH, ZLANHS
327
328
330
331
332 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min,
333 $ sqrt
334
335
336 DOUBLE PRECISION ABS1
337
338
339 abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
340
341
342
343
344
345 IF(
lsame( job,
'E' ) )
THEN
346 ilschr = .false.
347 ischur = 1
348 ELSE IF(
lsame( job,
'S' ) )
THEN
349 ilschr = .true.
350 ischur = 2
351 ELSE
352 ilschr = .true.
353 ischur = 0
354 END IF
355
356 IF(
lsame( compq,
'N' ) )
THEN
357 ilq = .false.
358 icompq = 1
359 ELSE IF(
lsame( compq,
'V' ) )
THEN
360 ilq = .true.
361 icompq = 2
362 ELSE IF(
lsame( compq,
'I' ) )
THEN
363 ilq = .true.
364 icompq = 3
365 ELSE
366 ilq = .true.
367 icompq = 0
368 END IF
369
370 IF(
lsame( compz,
'N' ) )
THEN
371 ilz = .false.
372 icompz = 1
373 ELSE IF(
lsame( compz,
'V' ) )
THEN
374 ilz = .true.
375 icompz = 2
376 ELSE IF(
lsame( compz,
'I' ) )
THEN
377 ilz = .true.
378 icompz = 3
379 ELSE
380 ilz = .true.
381 icompz = 0
382 END IF
383
384
385
386 info = 0
387 work( 1 ) = max( 1, n )
388 lquery = ( lwork.EQ.-1 )
389 IF( ischur.EQ.0 ) THEN
390 info = -1
391 ELSE IF( icompq.EQ.0 ) THEN
392 info = -2
393 ELSE IF( icompz.EQ.0 ) THEN
394 info = -3
395 ELSE IF( n.LT.0 ) THEN
396 info = -4
397 ELSE IF( ilo.LT.1 ) THEN
398 info = -5
399 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
400 info = -6
401 ELSE IF( ldh.LT.n ) THEN
402 info = -8
403 ELSE IF( ldt.LT.n ) THEN
404 info = -10
405 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
406 info = -14
407 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
408 info = -16
409 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
410 info = -18
411 END IF
412 IF( info.NE.0 ) THEN
413 CALL xerbla(
'ZHGEQZ', -info )
414 RETURN
415 ELSE IF( lquery ) THEN
416 RETURN
417 END IF
418
419
420
421
422 IF( n.LE.0 ) THEN
423 work( 1 ) = dcmplx( 1 )
424 RETURN
425 END IF
426
427
428
429 IF( icompq.EQ.3 )
430 $
CALL zlaset(
'Full', n, n, czero, cone, q, ldq )
431 IF( icompz.EQ.3 )
432 $
CALL zlaset(
'Full', n, n, czero, cone, z, ldz )
433
434
435
436 in = ihi + 1 - ilo
439 anorm =
zlanhs(
'F', in, h( ilo, ilo ), ldh, rwork )
440 bnorm =
zlanhs(
'F', in, t( ilo, ilo ), ldt, rwork )
441 atol = max( safmin, ulp*anorm )
442 btol = max( safmin, ulp*bnorm )
443 ascale = one / max( safmin, anorm )
444 bscale = one / max( safmin, bnorm )
445
446
447
448
449 DO 10 j = ihi + 1, n
450 absb = abs( t( j, j ) )
451 IF( absb.GT.safmin ) THEN
452 signbc = dconjg( t( j, j ) / absb )
453 t( j, j ) = absb
454 IF( ilschr ) THEN
455 CALL zscal( j-1, signbc, t( 1, j ), 1 )
456 CALL zscal( j, signbc, h( 1, j ), 1 )
457 ELSE
458 CALL zscal( 1, signbc, h( j, j ), 1 )
459 END IF
460 IF( ilz )
461 $
CALL zscal( n, signbc, z( 1, j ), 1 )
462 ELSE
463 t( j, j ) = czero
464 END IF
465 alpha( j ) = h( j, j )
466 beta( j ) = t( j, j )
467 10 CONTINUE
468
469
470
471 IF( ihi.LT.ilo )
472 $ GO TO 190
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489 ilast = ihi
490 IF( ilschr ) THEN
491 ifrstm = 1
492 ilastm = n
493 ELSE
494 ifrstm = ilo
495 ilastm = ihi
496 END IF
497 iiter = 0
498 eshift = czero
499 maxit = 30*( ihi-ilo+1 )
500
501 DO 170 jiter = 1, maxit
502
503
504
505 IF( jiter.GT.maxit )
506 $ GO TO 180
507
508
509
510
511
512
513
514
515
516 IF( ilast.EQ.ilo ) THEN
517 GO TO 60
518 ELSE
519 IF( abs1( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*(
520 $ abs1( h( ilast, ilast ) ) + abs1( h( ilast-1, ilast-1 )
521 $ ) ) ) ) THEN
522 h( ilast, ilast-1 ) = czero
523 GO TO 60
524 END IF
525 END IF
526
527 IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
528 t( ilast, ilast ) = czero
529 GO TO 50
530 END IF
531
532
533
534 DO 40 j = ilast - 1, ilo, -1
535
536
537
538 IF( j.EQ.ilo ) THEN
539 ilazro = .true.
540 ELSE
541 IF( abs1( h( j, j-1 ) ).LE.max( safmin, ulp*(
542 $ abs1( h( j, j ) ) + abs1( h( j-1, j-1 ) )
543 $ ) ) ) THEN
544 h( j, j-1 ) = czero
545 ilazro = .true.
546 ELSE
547 ilazro = .false.
548 END IF
549 END IF
550
551
552
553 IF( abs( t( j, j ) ).LT.btol ) THEN
554 t( j, j ) = czero
555
556
557
558 ilazr2 = .false.
559 IF( .NOT.ilazro ) THEN
560 IF( abs1( h( j, j-1 ) )*( ascale*abs1( h( j+1,
561 $ j ) ) ).LE.abs1( h( j, j ) )*( ascale*atol ) )
562 $ ilazr2 = .true.
563 END IF
564
565
566
567
568
569
570
571 IF( ilazro .OR. ilazr2 ) THEN
572 DO 20 jch = j, ilast - 1
573 ctemp = h( jch, jch )
574 CALL zlartg( ctemp, h( jch+1, jch ), c, s,
575 $ h( jch, jch ) )
576 h( jch+1, jch ) = czero
577 CALL zrot( ilastm-jch, h( jch, jch+1 ), ldh,
578 $ h( jch+1, jch+1 ), ldh, c, s )
579 CALL zrot( ilastm-jch, t( jch, jch+1 ), ldt,
580 $ t( jch+1, jch+1 ), ldt, c, s )
581 IF( ilq )
582 $
CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
583 $ c, dconjg( s ) )
584 IF( ilazr2 )
585 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
586 ilazr2 = .false.
587 IF( abs1( t( jch+1, jch+1 ) ).GE.btol ) THEN
588 IF( jch+1.GE.ilast ) THEN
589 GO TO 60
590 ELSE
591 ifirst = jch + 1
592 GO TO 70
593 END IF
594 END IF
595 t( jch+1, jch+1 ) = czero
596 20 CONTINUE
597 GO TO 50
598 ELSE
599
600
601
602
603 DO 30 jch = j, ilast - 1
604 ctemp = t( jch, jch+1 )
605 CALL zlartg( ctemp, t( jch+1, jch+1 ), c, s,
606 $ t( jch, jch+1 ) )
607 t( jch+1, jch+1 ) = czero
608 IF( jch.LT.ilastm-1 )
609 $
CALL zrot( ilastm-jch-1, t( jch, jch+2 ), ldt,
610 $ t( jch+1, jch+2 ), ldt, c, s )
611 CALL zrot( ilastm-jch+2, h( jch, jch-1 ), ldh,
612 $ h( jch+1, jch-1 ), ldh, c, s )
613 IF( ilq )
614 $
CALL zrot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
615 $ c, dconjg( s ) )
616 ctemp = h( jch+1, jch )
617 CALL zlartg( ctemp, h( jch+1, jch-1 ), c, s,
618 $ h( jch+1, jch ) )
619 h( jch+1, jch-1 ) = czero
620 CALL zrot( jch+1-ifrstm, h( ifrstm, jch ), 1,
621 $ h( ifrstm, jch-1 ), 1, c, s )
622 CALL zrot( jch-ifrstm, t( ifrstm, jch ), 1,
623 $ t( ifrstm, jch-1 ), 1, c, s )
624 IF( ilz )
625 $
CALL zrot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
626 $ c, s )
627 30 CONTINUE
628 GO TO 50
629 END IF
630 ELSE IF( ilazro ) THEN
631
632
633
634 ifirst = j
635 GO TO 70
636 END IF
637
638
639
640 40 CONTINUE
641
642
643
644 info = 2*n + 1
645 GO TO 210
646
647
648
649
650 50 CONTINUE
651 ctemp = h( ilast, ilast )
652 CALL zlartg( ctemp, h( ilast, ilast-1 ), c, s,
653 $ h( ilast, ilast ) )
654 h( ilast, ilast-1 ) = czero
655 CALL zrot( ilast-ifrstm, h( ifrstm, ilast ), 1,
656 $ h( ifrstm, ilast-1 ), 1, c, s )
657 CALL zrot( ilast-ifrstm, t( ifrstm, ilast ), 1,
658 $ t( ifrstm, ilast-1 ), 1, c, s )
659 IF( ilz )
660 $
CALL zrot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
661
662
663
664 60 CONTINUE
665 absb = abs( t( ilast, ilast ) )
666 IF( absb.GT.safmin ) THEN
667 signbc = dconjg( t( ilast, ilast ) / absb )
668 t( ilast, ilast ) = absb
669 IF( ilschr ) THEN
670 CALL zscal( ilast-ifrstm, signbc, t( ifrstm, ilast ), 1 )
671 CALL zscal( ilast+1-ifrstm, signbc, h( ifrstm, ilast ),
672 $ 1 )
673 ELSE
674 CALL zscal( 1, signbc, h( ilast, ilast ), 1 )
675 END IF
676 IF( ilz )
677 $
CALL zscal( n, signbc, z( 1, ilast ), 1 )
678 ELSE
679 t( ilast, ilast ) = czero
680 END IF
681 alpha( ilast ) = h( ilast, ilast )
682 beta( ilast ) = t( ilast, ilast )
683
684
685
686 ilast = ilast - 1
687 IF( ilast.LT.ilo )
688 $ GO TO 190
689
690
691
692 iiter = 0
693 eshift = czero
694 IF( .NOT.ilschr ) THEN
695 ilastm = ilast
696 IF( ifrstm.GT.ilast )
697 $ ifrstm = ilo
698 END IF
699 GO TO 160
700
701
702
703
704
705
706 70 CONTINUE
707 iiter = iiter + 1
708 IF( .NOT.ilschr ) THEN
709 ifrstm = ifirst
710 END IF
711
712
713
714
715
716
717
718 IF( ( iiter / 10 )*10.NE.iiter ) THEN
719
720
721
722
723
724
725
726
727 u12 = ( bscale*t( ilast-1, ilast ) ) /
728 $ ( bscale*t( ilast, ilast ) )
729 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
730 $ ( bscale*t( ilast-1, ilast-1 ) )
731 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
732 $ ( bscale*t( ilast-1, ilast-1 ) )
733 ad12 = ( ascale*h( ilast-1, ilast ) ) /
734 $ ( bscale*t( ilast, ilast ) )
735 ad22 = ( ascale*h( ilast, ilast ) ) /
736 $ ( bscale*t( ilast, ilast ) )
737 abi22 = ad22 - u12*ad21
738 abi12 = ad12 - u12*ad11
739
740 shift = abi22
741 ctemp = sqrt( abi12 )*sqrt( ad21 )
742 temp = abs1( ctemp )
743 IF( ctemp.NE.zero ) THEN
744 x = half*( ad11-shift )
745 temp2 = abs1( x )
746 temp = max( temp, abs1( x ) )
747 y = temp*sqrt( ( x / temp )**2+( ctemp / temp )**2 )
748 IF( temp2.GT.zero ) THEN
749 IF( dble( x / temp2 )*dble( y )+
750 $ dimag( x / temp2 )*dimag( y ).LT.zero )y = -y
751 END IF
752 shift = shift - ctemp*
zladiv( ctemp, ( x+y ) )
753 END IF
754 ELSE
755
756
757
758 IF( ( iiter / 20 )*20.EQ.iiter .AND.
759 $ bscale*abs1(t( ilast, ilast )).GT.safmin ) THEN
760 eshift = eshift + ( ascale*h( ilast,
761 $ ilast ) )/( bscale*t( ilast, ilast ) )
762 ELSE
763 eshift = eshift + ( ascale*h( ilast,
764 $ ilast-1 ) )/( bscale*t( ilast-1, ilast-1 ) )
765 END IF
766 shift = eshift
767 END IF
768
769
770
771 DO 80 j = ilast - 1, ifirst + 1, -1
772 istart = j
773 ctemp = ascale*h( j, j ) - shift*( bscale*t( j, j ) )
774 temp = abs1( ctemp )
775 temp2 = ascale*abs1( h( j+1, j ) )
776 tempr = max( temp, temp2 )
777 IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
778 temp = temp / tempr
779 temp2 = temp2 / tempr
780 END IF
781 IF( abs1( h( j, j-1 ) )*temp2.LE.temp*atol )
782 $ GO TO 90
783 80 CONTINUE
784
785 istart = ifirst
786 ctemp = ascale*h( ifirst, ifirst ) -
787 $ shift*( bscale*t( ifirst, ifirst ) )
788 90 CONTINUE
789
790
791
792
793
794 ctemp2 = ascale*h( istart+1, istart )
795 CALL zlartg( ctemp, ctemp2, c, s, ctemp3 )
796
797
798
799 DO 150 j = istart, ilast - 1
800 IF( j.GT.istart ) THEN
801 ctemp = h( j, j-1 )
802 CALL zlartg( ctemp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
803 h( j+1, j-1 ) = czero
804 END IF
805
806 DO 100 jc = j, ilastm
807 ctemp = c*h( j, jc ) + s*h( j+1, jc )
808 h( j+1, jc ) = -dconjg( s )*h( j, jc ) + c*h( j+1, jc )
809 h( j, jc ) = ctemp
810 ctemp2 = c*t( j, jc ) + s*t( j+1, jc )
811 t( j+1, jc ) = -dconjg( s )*t( j, jc ) + c*t( j+1, jc )
812 t( j, jc ) = ctemp2
813 100 CONTINUE
814 IF( ilq ) THEN
815 DO 110 jr = 1, n
816 ctemp = c*q( jr, j ) + dconjg( s )*q( jr, j+1 )
817 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
818 q( jr, j ) = ctemp
819 110 CONTINUE
820 END IF
821
822 ctemp = t( j+1, j+1 )
823 CALL zlartg( ctemp, t( j+1, j ), c, s, t( j+1, j+1 ) )
824 t( j+1, j ) = czero
825
826 DO 120 jr = ifrstm, min( j+2, ilast )
827 ctemp = c*h( jr, j+1 ) + s*h( jr, j )
828 h( jr, j ) = -dconjg( s )*h( jr, j+1 ) + c*h( jr, j )
829 h( jr, j+1 ) = ctemp
830 120 CONTINUE
831 DO 130 jr = ifrstm, j
832 ctemp = c*t( jr, j+1 ) + s*t( jr, j )
833 t( jr, j ) = -dconjg( s )*t( jr, j+1 ) + c*t( jr, j )
834 t( jr, j+1 ) = ctemp
835 130 CONTINUE
836 IF( ilz ) THEN
837 DO 140 jr = 1, n
838 ctemp = c*z( jr, j+1 ) + s*z( jr, j )
839 z( jr, j ) = -dconjg( s )*z( jr, j+1 ) + c*z( jr, j )
840 z( jr, j+1 ) = ctemp
841 140 CONTINUE
842 END IF
843 150 CONTINUE
844
845 160 CONTINUE
846
847 170 CONTINUE
848
849
850
851 180 CONTINUE
852 info = ilast
853 GO TO 210
854
855
856
857 190 CONTINUE
858
859
860
861 DO 200 j = 1, ilo - 1
862 absb = abs( t( j, j ) )
863 IF( absb.GT.safmin ) THEN
864 signbc = dconjg( t( j, j ) / absb )
865 t( j, j ) = absb
866 IF( ilschr ) THEN
867 CALL zscal( j-1, signbc, t( 1, j ), 1 )
868 CALL zscal( j, signbc, h( 1, j ), 1 )
869 ELSE
870 CALL zscal( 1, signbc, h( j, j ), 1 )
871 END IF
872 IF( ilz )
873 $
CALL zscal( n, signbc, z( 1, j ), 1 )
874 ELSE
875 t( j, j ) = czero
876 END IF
877 alpha( j ) = h( j, j )
878 beta( j ) = t( j, j )
879 200 CONTINUE
880
881
882
883 info = 0
884
885
886
887 210 CONTINUE
888 work( 1 ) = dcmplx( n )
889 RETURN
890
891
892
subroutine xerbla(srname, info)
complex *16 function zladiv(x, y)
ZLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
double precision function dlamch(cmach)
DLAMCH
double precision function zlanhs(norm, n, a, lda, work)
ZLANHS returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
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 zscal(n, za, zx, incx)
ZSCAL