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