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,
'.', / 9
x,
'N=',
952 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
953 9998
FORMAT(
' ZGET24: ', a,
' returned INFO=', i6,
'.', / 9
x,
'N=',
954 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zunt01(ROWCOL, M, N, U, LDU, WORK, LWORK, RWORK, RESID)
ZUNT01
DOUBLE PRECISION function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine xerbla(SRNAME, INFO)
XERBLA
LOGICAL function zslect(Z)
ZSLECT
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zget24(COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK, LWORK, RWORK, BWORK, INFO)
ZGET24
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine zgeesx(JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, W, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, RWORK, BWORK, INFO)
ZGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...