331 SUBROUTINE cget24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
332 $ H, HT, W, WT, WTMP, VS, LDVS, VS1, RCDEIN,
333 $ RCDVIN, NSLCT, ISLCT, ISRT, RESULT, WORK,
334 $ LWORK, RWORK, BWORK, INFO )
342 INTEGER INFO, ISRT, JTYPE, LDA, LDVS, LWORK, N, NOUNIT,
344 REAL RCDEIN, RCDVIN, THRESH
348 INTEGER ISEED( 4 ), ISLCT( * )
349 REAL RESULT( 17 ), RWORK( * )
350 COMPLEX A( LDA, * ), H( LDA, * ), HT( LDA, * ),
351 $ vs( ldvs, * ), vs1( ldvs, * ), w( * ),
352 $ work( * ), wt( * ), wtmp( * )
359 PARAMETER ( CZERO = ( 0.0e+0, 0.0e+0 ),
360 $ cone = ( 1.0e+0, 0.0e+0 ) )
362 parameter( zero = 0.0e+0, one = 1.0e+0 )
364 parameter( epsin = 5.9605e-8 )
368 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, RSUB,
370 REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
371 $ smlnum, tol, tolin, ulp, ulpinv, v, vricmp,
381 EXTERNAL CSLECT, CLANGE, SLAMCH
387 INTRINSIC abs, aimag, max, min, real
391 REAL SELWI( 20 ), SELWR( 20 )
394 INTEGER SELDIM, SELOPT
397 COMMON / sslct / selopt, seldim, selval, selwr, selwi
404 IF( thresh.LT.zero )
THEN
406 ELSE IF( nounit.LE.0 )
THEN
408 ELSE IF( n.LT.0 )
THEN
410 ELSE IF( lda.LT.1 .OR. lda.LT.n )
THEN
412 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n )
THEN
414 ELSE IF( lwork.LT.2*n )
THEN
419 CALL xerbla(
'CGET24', -info )
434 smlnum = slamch(
'Safe minimum' )
435 ulp = slamch(
'Precision' )
442 IF( isort.EQ.0 )
THEN
452 CALL clacpy(
'F', n, n, a, lda, h, lda )
453 CALL cgeesx(
'V', sort, cslect,
'N', n, h, lda, sdim, w, vs,
454 $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
456 IF( iinfo.NE.0 )
THEN
457 result( 1+rsub ) = ulpinv
458 IF( jtype.NE.22 )
THEN
459 WRITE( nounit, fmt = 9998 )
'CGEESX1', iinfo, n, jtype,
462 WRITE( nounit, fmt = 9999 )
'CGEESX1', iinfo, n,
468 IF( isort.EQ.0 )
THEN
469 CALL ccopy( n, w, 1, wtmp, 1 )
474 result( 1+rsub ) = zero
477 IF( h( i, j ).NE.czero )
478 $ result( 1+rsub ) = ulpinv
486 CALL clacpy(
' ', n, n, a, lda, vs1, ldvs )
490 CALL cgemm(
'No transpose',
'No transpose', n, n, n, cone, vs,
491 $ ldvs, h, lda, czero, ht, lda )
495 CALL cgemm(
'No transpose',
'Conjugate transpose', n, n, n,
496 $ -cone, ht, lda, vs, ldvs, cone, vs1, ldvs )
498 anorm = max( clange(
'1', n, n, a, lda, rwork ), smlnum )
499 wnorm = clange(
'1', n, n, vs1, ldvs, rwork )
501 IF( anorm.GT.wnorm )
THEN
502 result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
504 IF( anorm.LT.one )
THEN
505 result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
508 result( 2+rsub ) = min( wnorm / anorm, real( n ) ) /
515 CALL cunt01(
'Columns', n, n, vs, ldvs, work, lwork, rwork,
520 result( 4+rsub ) = zero
522 IF( h( i, i ).NE.w( i ) )
523 $ result( 4+rsub ) = ulpinv
528 CALL clacpy(
'F', n, n, a, lda, ht, lda )
529 CALL cgeesx(
'N', sort, cslect,
'N', n, ht, lda, sdim, wt, vs,
530 $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
532 IF( iinfo.NE.0 )
THEN
533 result( 5+rsub ) = ulpinv
534 IF( jtype.NE.22 )
THEN
535 WRITE( nounit, fmt = 9998 )
'CGEESX2', iinfo, n, jtype,
538 WRITE( nounit, fmt = 9999 )
'CGEESX2', iinfo, n,
545 result( 5+rsub ) = zero
548 IF( h( i, j ).NE.ht( i, j ) )
549 $ result( 5+rsub ) = ulpinv
555 result( 6+rsub ) = zero
557 IF( w( i ).NE.wt( i ) )
558 $ result( 6+rsub ) = ulpinv
563 IF( isort.EQ.1 )
THEN
567 IF( cslect( w( i ) ) )
568 $ knteig = knteig + 1
570 IF( cslect( w( i+1 ) ) .AND.
571 $ ( .NOT.cslect( w( i ) ) ) )result( 13 ) = ulpinv
575 $ result( 13 ) = ulpinv
583 IF( lwork.GE.( n*( n+1 ) ) / 2 )
THEN
590 CALL clacpy(
'F', n, n, a, lda, ht, lda )
591 CALL cgeesx(
'V', sort, cslect,
'B', n, ht, lda, sdim1, wt,
592 $ vs1, ldvs, rconde, rcondv, work, lwork, rwork,
594 IF( iinfo.NE.0 )
THEN
595 result( 14 ) = ulpinv
596 result( 15 ) = ulpinv
597 IF( jtype.NE.22 )
THEN
598 WRITE( nounit, fmt = 9998 )
'CGEESX3', iinfo, n, jtype,
601 WRITE( nounit, fmt = 9999 )
'CGEESX3', iinfo, n,
611 IF( w( i ).NE.wt( i ) )
612 $ result( 10 ) = ulpinv
614 IF( h( i, j ).NE.ht( i, j ) )
615 $ result( 11 ) = ulpinv
616 IF( vs( i, j ).NE.vs1( i, j ) )
617 $ result( 12 ) = ulpinv
621 $ result( 13 ) = ulpinv
625 CALL clacpy(
'F', n, n, a, lda, ht, lda )
626 CALL cgeesx(
'N', sort, cslect,
'B', n, ht, lda, sdim1, wt,
627 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
629 IF( iinfo.NE.0 )
THEN
630 result( 14 ) = ulpinv
631 result( 15 ) = ulpinv
632 IF( jtype.NE.22 )
THEN
633 WRITE( nounit, fmt = 9998 )
'CGEESX4', iinfo, n, jtype,
636 WRITE( nounit, fmt = 9999 )
'CGEESX4', iinfo, n,
645 IF( rcnde1.NE.rconde )
646 $ result( 14 ) = ulpinv
647 IF( rcndv1.NE.rcondv )
648 $ result( 15 ) = ulpinv
653 IF( w( i ).NE.wt( i ) )
654 $ result( 10 ) = ulpinv
656 IF( h( i, j ).NE.ht( i, j ) )
657 $ result( 11 ) = ulpinv
658 IF( vs( i, j ).NE.vs1( i, j ) )
659 $ result( 12 ) = ulpinv
663 $ result( 13 ) = ulpinv
667 CALL clacpy(
'F', n, n, a, lda, ht, lda )
668 CALL cgeesx(
'V', sort, cslect,
'E', n, ht, lda, sdim1, wt,
669 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
671 IF( iinfo.NE.0 )
THEN
672 result( 14 ) = ulpinv
673 IF( jtype.NE.22 )
THEN
674 WRITE( nounit, fmt = 9998 )
'CGEESX5', iinfo, n, jtype,
677 WRITE( nounit, fmt = 9999 )
'CGEESX5', iinfo, n,
686 IF( rcnde1.NE.rconde )
687 $ result( 14 ) = ulpinv
692 IF( w( i ).NE.wt( i ) )
693 $ result( 10 ) = ulpinv
695 IF( h( i, j ).NE.ht( i, j ) )
696 $ result( 11 ) = ulpinv
697 IF( vs( i, j ).NE.vs1( i, j ) )
698 $ result( 12 ) = ulpinv
702 $ result( 13 ) = ulpinv
706 CALL clacpy(
'F', n, n, a, lda, ht, lda )
707 CALL cgeesx(
'N', sort, cslect,
'E', n, ht, lda, sdim1, wt,
708 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
710 IF( iinfo.NE.0 )
THEN
711 result( 14 ) = ulpinv
712 IF( jtype.NE.22 )
THEN
713 WRITE( nounit, fmt = 9998 )
'CGEESX6', iinfo, n, jtype,
716 WRITE( nounit, fmt = 9999 )
'CGEESX6', iinfo, n,
725 IF( rcnde1.NE.rconde )
726 $ result( 14 ) = ulpinv
731 IF( w( i ).NE.wt( i ) )
732 $ result( 10 ) = ulpinv
734 IF( h( i, j ).NE.ht( i, j ) )
735 $ result( 11 ) = ulpinv
736 IF( vs( i, j ).NE.vs1( i, j ) )
737 $ result( 12 ) = ulpinv
741 $ result( 13 ) = ulpinv
745 CALL clacpy(
'F', n, n, a, lda, ht, lda )
746 CALL cgeesx(
'V', sort, cslect,
'V', n, ht, lda, sdim1, wt,
747 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
749 IF( iinfo.NE.0 )
THEN
750 result( 15 ) = ulpinv
751 IF( jtype.NE.22 )
THEN
752 WRITE( nounit, fmt = 9998 )
'CGEESX7', iinfo, n, jtype,
755 WRITE( nounit, fmt = 9999 )
'CGEESX7', iinfo, n,
764 IF( rcndv1.NE.rcondv )
765 $ result( 15 ) = ulpinv
770 IF( w( i ).NE.wt( i ) )
771 $ result( 10 ) = ulpinv
773 IF( h( i, j ).NE.ht( i, j ) )
774 $ result( 11 ) = ulpinv
775 IF( vs( i, j ).NE.vs1( i, j ) )
776 $ result( 12 ) = ulpinv
780 $ result( 13 ) = ulpinv
784 CALL clacpy(
'F', n, n, a, lda, ht, lda )
785 CALL cgeesx(
'N', sort, cslect,
'V', n, ht, lda, sdim1, wt,
786 $ vs1, ldvs, rcnde1, rcndv1, work, lwork, rwork,
788 IF( iinfo.NE.0 )
THEN
789 result( 15 ) = ulpinv
790 IF( jtype.NE.22 )
THEN
791 WRITE( nounit, fmt = 9998 )
'CGEESX8', iinfo, n, jtype,
794 WRITE( nounit, fmt = 9999 )
'CGEESX8', iinfo, n,
803 IF( rcndv1.NE.rcondv )
804 $ result( 15 ) = ulpinv
809 IF( w( i ).NE.wt( i ) )
810 $ result( 10 ) = ulpinv
812 IF( h( i, j ).NE.ht( i, j ) )
813 $ result( 11 ) = ulpinv
814 IF( vs( i, j ).NE.vs1( i, j ) )
815 $ result( 12 ) = ulpinv
819 $ result( 13 ) = ulpinv
836 eps = max( ulp, epsin )
839 selval( i ) = .false.
840 selwr( i ) = real( wtmp( i ) )
841 selwi( i ) = aimag( wtmp( i ) )
846 vrimin = real( wtmp( i ) )
848 vrimin = aimag( wtmp( i ) )
852 vricmp = real( wtmp( j ) )
854 vricmp = aimag( wtmp( j ) )
856 IF( vricmp.LT.vrimin )
THEN
862 wtmp( kmin ) = wtmp( i )
865 ipnt( i ) = ipnt( kmin )
869 selval( ipnt( islct( i ) ) ) = .true.
874 CALL clacpy(
'F', n, n, a, lda, ht, lda )
875 CALL cgeesx(
'N',
'S', cslect,
'B', n, ht, lda, sdim1, wt, vs1,
876 $ ldvs, rconde, rcondv, work, lwork, rwork, bwork,
878 IF( iinfo.NE.0 )
THEN
879 result( 16 ) = ulpinv
880 result( 17 ) = ulpinv
881 WRITE( nounit, fmt = 9999 )
'CGEESX9', iinfo, n, iseed( 1 )
889 anorm = clange(
'1', n, n, a, lda, rwork )
890 v = max( real( n )*eps*anorm, smlnum )
893 IF( v.GT.rcondv )
THEN
898 IF( v.GT.rcdvin )
THEN
903 tol = max( tol, smlnum / eps )
904 tolin = max( tolin, smlnum / eps )
905 IF( eps*( rcdein-tolin ).GT.rconde+tol )
THEN
906 result( 16 ) = ulpinv
907 ELSE IF( rcdein-tolin.GT.rconde+tol )
THEN
908 result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
909 ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) )
THEN
910 result( 16 ) = ulpinv
911 ELSE IF( rcdein+tolin.LT.rconde-tol )
THEN
912 result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
920 IF( v.GT.rcondv*rconde )
THEN
925 IF( v.GT.rcdvin*rcdein )
THEN
930 tol = max( tol, smlnum / eps )
931 tolin = max( tolin, smlnum / eps )
932 IF( eps*( rcdvin-tolin ).GT.rcondv+tol )
THEN
933 result( 17 ) = ulpinv
934 ELSE IF( rcdvin-tolin.GT.rcondv+tol )
THEN
935 result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
936 ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) )
THEN
937 result( 17 ) = ulpinv
938 ELSE IF( rcdvin+tolin.LT.rcondv-tol )
THEN
939 result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
948 9999
FORMAT(
' CGET24: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
949 $ i6,
', INPUT EXAMPLE NUMBER = ', i4 )
950 9998
FORMAT(
' CGET24: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
951 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
subroutine xerbla(srname, info)
subroutine cget24(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)
CGET24
subroutine cunt01(rowcol, m, n, u, ldu, work, lwork, rwork, resid)
CUNT01
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
subroutine cgeesx(jobvs, sort, select, sense, n, a, lda, sdim, w, vs, ldvs, rconde, rcondv, work, lwork, rwork, bwork, info)
CGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.