333 SUBROUTINE zget24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
334 $ h, ht, w, wt, wtmp, vs, ldvs, vs1, rcdein,
335 $ rcdvin, nslct, islct, isrt, result, work,
336 $ lwork, rwork, bwork, info )
345 INTEGER info, isrt, jtype, lda, ldvs, lwork, n, nounit,
347 DOUBLE PRECISION rcdein, rcdvin, thresh
351 INTEGER iseed( 4 ), islct( * )
352 DOUBLE PRECISION result( 17 ), rwork( * )
353 COMPLEX*16 a( lda, * ), h( lda, * ), ht( lda, * ),
354 $ vs( ldvs, * ), vs1( ldvs, * ), w( * ),
355 $ work( * ), wt( * ), wtmp( * )
361 COMPLEX*16 czero, cone
362 parameter( czero = ( 0.0d+0, 0.0d+0 ),
363 $ cone = ( 1.0d+0, 0.0d+0 ) )
364 DOUBLE PRECISION zero, one
365 parameter( zero = 0.0d+0, one = 1.0d+0 )
366 DOUBLE PRECISION epsin
367 parameter( epsin = 5.9605d-8 )
371 INTEGER i, iinfo, isort, itmp, j, kmin, knteig, rsub,
373 DOUBLE PRECISION anorm, eps, rcnde1, rcndv1, rconde, rcondv,
374 $ smlnum, tol, tolin, ulp, ulpinv, v, vricmp,
390 INTRINSIC abs, dble, dimag, max, min
394 DOUBLE PRECISION selwi( 20 ), selwr( 20 )
397 INTEGER seldim, selopt
400 common / sslct / selopt, seldim, selval, selwr, selwi
407 IF( thresh.LT.zero )
THEN
409 ELSE IF( nounit.LE.0 )
THEN
411 ELSE IF( n.LT.0 )
THEN
413 ELSE IF( lda.LT.1 .OR. lda.LT.n )
THEN
415 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n )
THEN
417 ELSE IF( lwork.LT.2*n )
THEN
422 CALL
xerbla(
'ZGET24', -info )
437 smlnum =
dlamch(
'Safe minimum' )
438 ulp =
dlamch(
'Precision' )
445 IF( isort.EQ.0 )
THEN
455 CALL
zlacpy(
'F', n, n, a, lda, h, lda )
456 CALL
zgeesx(
'V', sort,
zslect,
'N', n, h, lda, sdim, w, vs,
457 $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
459 IF( iinfo.NE.0 )
THEN
460 result( 1+rsub ) = ulpinv
461 IF( jtype.NE.22 )
THEN
462 WRITE( nounit, fmt = 9998 )
'ZGEESX1', iinfo, n, jtype,
465 WRITE( nounit, fmt = 9999 )
'ZGEESX1', iinfo, n,
471 IF( isort.EQ.0 )
THEN
472 CALL
zcopy( n, w, 1, wtmp, 1 )
477 result( 1+rsub ) = zero
480 IF( h( i, j ).NE.czero )
481 $ result( 1+rsub ) = ulpinv
489 CALL
zlacpy(
' ', n, n, a, lda, vs1, ldvs )
493 CALL
zgemm(
'No transpose',
'No transpose', n, n, n, cone, vs,
494 $ ldvs, h, lda, czero, ht, lda )
498 CALL
zgemm(
'No transpose',
'Conjugate transpose', n, n, n,
499 $ -cone, ht, lda, vs, ldvs, cone, vs1, ldvs )
501 anorm = max(
zlange(
'1', n, n, a, lda, rwork ), smlnum )
502 wnorm =
zlange(
'1', n, n, vs1, ldvs, rwork )
504 IF( anorm.GT.wnorm )
THEN
505 result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
507 IF( anorm.LT.one )
THEN
508 result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
511 result( 2+rsub ) = min( wnorm / anorm, dble( n ) ) /
518 CALL
zunt01(
'Columns', n, n, vs, ldvs, work, lwork, rwork,
523 result( 4+rsub ) = zero
525 IF( h( i, i ).NE.w( i ) )
526 $ result( 4+rsub ) = ulpinv
531 CALL
zlacpy(
'F', n, n, a, lda, ht, lda )
532 CALL
zgeesx(
'N', sort,
zslect,
'N', n, ht, lda, sdim, wt, vs,
533 $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
535 IF( iinfo.NE.0 )
THEN
536 result( 5+rsub ) = ulpinv
537 IF( jtype.NE.22 )
THEN
538 WRITE( nounit, fmt = 9998 )
'ZGEESX2', iinfo, n, jtype,
541 WRITE( nounit, fmt = 9999 )
'ZGEESX2', iinfo, n,
548 result( 5+rsub ) = zero
551 IF( h( i, j ).NE.ht( i, j ) )
552 $ result( 5+rsub ) = ulpinv
558 result( 6+rsub ) = zero
560 IF( w( i ).NE.wt( i ) )
561 $ result( 6+rsub ) = ulpinv
566 IF( isort.EQ.1 )
THEN
571 $ knteig = knteig + 1
573 IF(
zslect( w( i+1 ) ) .AND.
574 $ ( .NOT.
zslect( w( i ) ) ) )result( 13 ) = ulpinv
578 $ result( 13 ) = ulpinv
586 IF( lwork.GE.( n*( n+1 ) ) / 2 )
THEN
593 CALL
zlacpy(
'F', n, n, a, lda, ht, lda )
594 CALL
zgeesx(
'V', sort,
zslect,
'B', n, ht, lda, sdim1, wt,
595 $ vs1, ldvs, rconde, rcondv, work, lwork, rwork,
597 IF( iinfo.NE.0 )
THEN
598 result( 14 ) = ulpinv
599 result( 15 ) = ulpinv
600 IF( jtype.NE.22 )
THEN
601 WRITE( nounit, fmt = 9998 )
'ZGEESX3', iinfo, n, jtype,
604 WRITE( nounit, fmt = 9999 )
'ZGEESX3', iinfo, n,
614 IF( w( i ).NE.wt( i ) )
615 $ result( 10 ) = ulpinv
617 IF( h( i, j ).NE.ht( i, j ) )
618 $ result( 11 ) = ulpinv
619 IF( vs( i, j ).NE.vs1( i, j ) )
620 $ result( 12 ) = ulpinv
624 $ result( 13 ) = ulpinv
628 CALL
zlacpy(
'F', n, n, a, lda, ht, lda )
629 CALL
zgeesx(
'N', sort,
zslect,
'B', n, ht, lda, sdim1, wt,
630 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
632 IF( iinfo.NE.0 )
THEN
633 result( 14 ) = ulpinv
634 result( 15 ) = ulpinv
635 IF( jtype.NE.22 )
THEN
636 WRITE( nounit, fmt = 9998 )
'ZGEESX4', iinfo, n, jtype,
639 WRITE( nounit, fmt = 9999 )
'ZGEESX4', iinfo, n,
648 IF( rcnde1.NE.rconde )
649 $ result( 14 ) = ulpinv
650 IF( rcndv1.NE.rcondv )
651 $ result( 15 ) = ulpinv
656 IF( w( i ).NE.wt( i ) )
657 $ result( 10 ) = ulpinv
659 IF( h( i, j ).NE.ht( i, j ) )
660 $ result( 11 ) = ulpinv
661 IF( vs( i, j ).NE.vs1( i, j ) )
662 $ result( 12 ) = ulpinv
666 $ result( 13 ) = ulpinv
670 CALL
zlacpy(
'F', n, n, a, lda, ht, lda )
671 CALL
zgeesx(
'V', sort,
zslect,
'E', n, ht, lda, sdim1, wt,
672 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
674 IF( iinfo.NE.0 )
THEN
675 result( 14 ) = ulpinv
676 IF( jtype.NE.22 )
THEN
677 WRITE( nounit, fmt = 9998 )
'ZGEESX5', iinfo, n, jtype,
680 WRITE( nounit, fmt = 9999 )
'ZGEESX5', iinfo, n,
689 IF( rcnde1.NE.rconde )
690 $ result( 14 ) = ulpinv
695 IF( w( i ).NE.wt( i ) )
696 $ result( 10 ) = ulpinv
698 IF( h( i, j ).NE.ht( i, j ) )
699 $ result( 11 ) = ulpinv
700 IF( vs( i, j ).NE.vs1( i, j ) )
701 $ result( 12 ) = ulpinv
705 $ result( 13 ) = ulpinv
709 CALL
zlacpy(
'F', n, n, a, lda, ht, lda )
710 CALL
zgeesx(
'N', sort,
zslect,
'E', n, ht, lda, sdim1, wt,
711 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
713 IF( iinfo.NE.0 )
THEN
714 result( 14 ) = ulpinv
715 IF( jtype.NE.22 )
THEN
716 WRITE( nounit, fmt = 9998 )
'ZGEESX6', iinfo, n, jtype,
719 WRITE( nounit, fmt = 9999 )
'ZGEESX6', iinfo, n,
728 IF( rcnde1.NE.rconde )
729 $ result( 14 ) = ulpinv
734 IF( w( i ).NE.wt( i ) )
735 $ result( 10 ) = ulpinv
737 IF( h( i, j ).NE.ht( i, j ) )
738 $ result( 11 ) = ulpinv
739 IF( vs( i, j ).NE.vs1( i, j ) )
740 $ result( 12 ) = ulpinv
744 $ result( 13 ) = ulpinv
748 CALL
zlacpy(
'F', n, n, a, lda, ht, lda )
749 CALL
zgeesx(
'V', sort,
zslect,
'V', n, ht, lda, sdim1, wt,
750 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
752 IF( iinfo.NE.0 )
THEN
753 result( 15 ) = ulpinv
754 IF( jtype.NE.22 )
THEN
755 WRITE( nounit, fmt = 9998 )
'ZGEESX7', iinfo, n, jtype,
758 WRITE( nounit, fmt = 9999 )
'ZGEESX7', iinfo, n,
767 IF( rcndv1.NE.rcondv )
768 $ result( 15 ) = ulpinv
773 IF( w( i ).NE.wt( i ) )
774 $ result( 10 ) = ulpinv
776 IF( h( i, j ).NE.ht( i, j ) )
777 $ result( 11 ) = ulpinv
778 IF( vs( i, j ).NE.vs1( i, j ) )
779 $ result( 12 ) = ulpinv
783 $ result( 13 ) = ulpinv
787 CALL
zlacpy(
'F', n, n, a, lda, ht, lda )
788 CALL
zgeesx(
'N', sort,
zslect,
'V', n, ht, lda, sdim1, wt,
789 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
791 IF( iinfo.NE.0 )
THEN
792 result( 15 ) = ulpinv
793 IF( jtype.NE.22 )
THEN
794 WRITE( nounit, fmt = 9998 )
'ZGEESX8', iinfo, n, jtype,
797 WRITE( nounit, fmt = 9999 )
'ZGEESX8', iinfo, n,
806 IF( rcndv1.NE.rcondv )
807 $ result( 15 ) = ulpinv
812 IF( w( i ).NE.wt( i ) )
813 $ result( 10 ) = ulpinv
815 IF( h( i, j ).NE.ht( i, j ) )
816 $ result( 11 ) = ulpinv
817 IF( vs( i, j ).NE.vs1( i, j ) )
818 $ result( 12 ) = ulpinv
822 $ result( 13 ) = ulpinv
839 eps = max( ulp, epsin )
842 selval( i ) = .false.
843 selwr( i ) = dble( wtmp( i ) )
844 selwi( i ) = dimag( wtmp( i ) )
849 vrimin = dble( wtmp( i ) )
851 vrimin = dimag( wtmp( i ) )
855 vricmp = dble( wtmp( j ) )
857 vricmp = dimag( wtmp( j ) )
859 IF( vricmp.LT.vrimin )
THEN
865 wtmp( kmin ) = wtmp( i )
868 ipnt( i ) = ipnt( kmin )
872 selval( ipnt( islct( i ) ) ) = .true.
877 CALL
zlacpy(
'F', n, n, a, lda, ht, lda )
878 CALL
zgeesx(
'N',
'S',
zslect,
'B', n, ht, lda, sdim1, wt, vs1,
879 $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
881 IF( iinfo.NE.0 )
THEN
882 result( 16 ) = ulpinv
883 result( 17 ) = ulpinv
884 WRITE( nounit, fmt = 9999 )
'ZGEESX9', iinfo, n, iseed( 1 )
892 anorm =
zlange(
'1', n, n, a, lda, rwork )
893 v = max( dble( n )*eps*anorm, smlnum )
896 IF( v.GT.rcondv )
THEN
901 IF( v.GT.rcdvin )
THEN
906 tol = max( tol, smlnum / eps )
907 tolin = max( tolin, smlnum / eps )
908 IF( eps*( rcdein-tolin ).GT.rconde+tol )
THEN
909 result( 16 ) = ulpinv
910 ELSE IF( rcdein-tolin.GT.rconde+tol )
THEN
911 result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
912 ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) )
THEN
913 result( 16 ) = ulpinv
914 ELSE IF( rcdein+tolin.LT.rconde-tol )
THEN
915 result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
923 IF( v.GT.rcondv*rconde )
THEN
928 IF( v.GT.rcdvin*rcdein )
THEN
933 tol = max( tol, smlnum / eps )
934 tolin = max( tolin, smlnum / eps )
935 IF( eps*( rcdvin-tolin ).GT.rcondv+tol )
THEN
936 result( 17 ) = ulpinv
937 ELSE IF( rcdvin-tolin.GT.rcondv+tol )
THEN
938 result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
939 ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) )
THEN
940 result( 17 ) = ulpinv
941 ELSE IF( rcdvin+tolin.LT.rcondv-tol )
THEN
942 result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
951 9999 format(
' ZGET24: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
952 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
953 9998 format(
' ZGET24: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
954 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )