LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ddrvsx()

subroutine ddrvsx ( integer  NSIZES,
integer, dimension( * )  NN,
integer  NTYPES,
logical, dimension( * )  DOTYPE,
integer, dimension( 4 )  ISEED,
double precision  THRESH,
integer  NIUNIT,
integer  NOUNIT,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( lda, * )  H,
double precision, dimension( lda, * )  HT,
double precision, dimension( * )  WR,
double precision, dimension( * )  WI,
double precision, dimension( * )  WRT,
double precision, dimension( * )  WIT,
double precision, dimension( * )  WRTMP,
double precision, dimension( * )  WITMP,
double precision, dimension( ldvs, * )  VS,
integer  LDVS,
double precision, dimension( ldvs, * )  VS1,
double precision, dimension( 17 )  RESULT,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

DDRVSX

Purpose:
    DDRVSX checks the nonsymmetric eigenvalue (Schur form) problem
    expert driver DGEESX.

    DDRVSX uses both test matrices generated randomly depending on
    data supplied in the calling sequence, as well as on data
    read from an input file and including precomputed condition
    numbers to which it compares the ones it computes.

    When DDRVSX is called, a number of matrix "sizes" ("n's") and a
    number of matrix "types" are specified.  For each size ("n")
    and each type of matrix, one matrix will be generated and used
    to test the nonsymmetric eigenroutines.  For each matrix, 15
    tests will be performed:

    (1)     0 if T is in Schur form, 1/ulp otherwise
           (no sorting of eigenvalues)

    (2)     | A - VS T VS' | / ( n |A| ulp )

      Here VS is the matrix of Schur eigenvectors, and T is in Schur
      form  (no sorting of eigenvalues).

    (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).

    (4)     0     if WR+sqrt(-1)*WI are eigenvalues of T
            1/ulp otherwise
            (no sorting of eigenvalues)

    (5)     0     if T(with VS) = T(without VS),
            1/ulp otherwise
            (no sorting of eigenvalues)

    (6)     0     if eigenvalues(with VS) = eigenvalues(without VS),
            1/ulp otherwise
            (no sorting of eigenvalues)

    (7)     0 if T is in Schur form, 1/ulp otherwise
            (with sorting of eigenvalues)

    (8)     | A - VS T VS' | / ( n |A| ulp )

      Here VS is the matrix of Schur eigenvectors, and T is in Schur
      form  (with sorting of eigenvalues).

    (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).

    (10)    0     if WR+sqrt(-1)*WI are eigenvalues of T
            1/ulp otherwise
            If workspace sufficient, also compare WR, WI with and
            without reciprocal condition numbers
            (with sorting of eigenvalues)

    (11)    0     if T(with VS) = T(without VS),
            1/ulp otherwise
            If workspace sufficient, also compare T with and without
            reciprocal condition numbers
            (with sorting of eigenvalues)

    (12)    0     if eigenvalues(with VS) = eigenvalues(without VS),
            1/ulp otherwise
            If workspace sufficient, also compare VS with and without
            reciprocal condition numbers
            (with sorting of eigenvalues)

    (13)    if sorting worked and SDIM is the number of
            eigenvalues which were SELECTed
            If workspace sufficient, also compare SDIM with and
            without reciprocal condition numbers

    (14)    if RCONDE the same no matter if VS and/or RCONDV computed

    (15)    if RCONDV the same no matter if VS and/or RCONDE computed

    The "sizes" are specified by an array NN(1:NSIZES); the value of
    each element NN(j) specifies one size.
    The "types" are specified by a logical array DOTYPE( 1:NTYPES );
    if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
    Currently, the list of possible types is:

    (1)  The zero matrix.
    (2)  The identity matrix.
    (3)  A (transposed) Jordan block, with 1's on the diagonal.

    (4)  A diagonal matrix with evenly spaced entries
         1, ..., ULP  and random signs.
         (ULP = (first number larger than 1) - 1 )
    (5)  A diagonal matrix with geometrically spaced entries
         1, ..., ULP  and random signs.
    (6)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
         and random signs.

    (7)  Same as (4), but multiplied by a constant near
         the overflow threshold
    (8)  Same as (4), but multiplied by a constant near
         the underflow threshold

    (9)  A matrix of the form  U' T U, where U is orthogonal and
         T has evenly spaced entries 1, ..., ULP with random signs
         on the diagonal and random O(1) entries in the upper
         triangle.

    (10) A matrix of the form  U' T U, where U is orthogonal and
         T has geometrically spaced entries 1, ..., ULP with random
         signs on the diagonal and random O(1) entries in the upper
         triangle.

    (11) A matrix of the form  U' T U, where U is orthogonal and
         T has "clustered" entries 1, ULP,..., ULP with random
         signs on the diagonal and random O(1) entries in the upper
         triangle.

    (12) A matrix of the form  U' T U, where U is orthogonal and
         T has real or complex conjugate paired eigenvalues randomly
         chosen from ( ULP, 1 ) and random O(1) entries in the upper
         triangle.

    (13) A matrix of the form  X' T X, where X has condition
         SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
         with random signs on the diagonal and random O(1) entries
         in the upper triangle.

    (14) A matrix of the form  X' T X, where X has condition
         SQRT( ULP ) and T has geometrically spaced entries
         1, ..., ULP with random signs on the diagonal and random
         O(1) entries in the upper triangle.

    (15) A matrix of the form  X' T X, where X has condition
         SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
         with random signs on the diagonal and random O(1) entries
         in the upper triangle.

    (16) A matrix of the form  X' T X, where X has condition
         SQRT( ULP ) and T has real or complex conjugate paired
         eigenvalues randomly chosen from ( ULP, 1 ) and random
         O(1) entries in the upper triangle.

    (17) Same as (16), but multiplied by a constant
         near the overflow threshold
    (18) Same as (16), but multiplied by a constant
         near the underflow threshold

    (19) Nonsymmetric matrix with random entries chosen from (-1,1).
         If N is at least 4, all entries in first two rows and last
         row, and first column and last two columns are zero.
    (20) Same as (19), but multiplied by a constant
         near the overflow threshold
    (21) Same as (19), but multiplied by a constant
         near the underflow threshold

    In addition, an input file will be read from logical unit number
    NIUNIT. The file contains matrices along with precomputed
    eigenvalues and reciprocal condition numbers for the eigenvalue
    average and right invariant subspace. For these matrices, in
    addition to tests (1) to (15) we will compute the following two
    tests:

   (16)  |RCONDE - RCDEIN| / cond(RCONDE)

      RCONDE is the reciprocal average eigenvalue condition number
      computed by DGEESX and RCDEIN (the precomputed true value)
      is supplied as input.  cond(RCONDE) is the condition number
      of RCONDE, and takes errors in computing RCONDE into account,
      so that the resulting quantity should be O(ULP). cond(RCONDE)
      is essentially given by norm(A)/RCONDV.

   (17)  |RCONDV - RCDVIN| / cond(RCONDV)

      RCONDV is the reciprocal right invariant subspace condition
      number computed by DGEESX and RCDVIN (the precomputed true
      value) is supplied as input. cond(RCONDV) is the condition
      number of RCONDV, and takes errors in computing RCONDV into
      account, so that the resulting quantity should be O(ULP).
      cond(RCONDV) is essentially given by norm(A)/RCONDE.
Parameters
[in]NSIZES
          NSIZES is INTEGER
          The number of sizes of matrices to use.  NSIZES must be at
          least zero. If it is zero, no randomly generated matrices
          are tested, but any test matrices read from NIUNIT will be
          tested.
[in]NN
          NN is INTEGER array, dimension (NSIZES)
          An array containing the sizes to be used for the matrices.
          Zero values will be skipped.  The values must be at least
          zero.
[in]NTYPES
          NTYPES is INTEGER
          The number of elements in DOTYPE. NTYPES must be at least
          zero. If it is zero, no randomly generated test matrices
          are tested, but and test matrices read from NIUNIT will be
          tested. If it is MAXTYP+1 and NSIZES is 1, then an
          additional type, MAXTYP+1 is defined, which is to use
          whatever matrix is in A.  This is only useful if
          DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .
[in]DOTYPE
          DOTYPE is LOGICAL array, dimension (NTYPES)
          If DOTYPE(j) is .TRUE., then for each size in NN a
          matrix of that size and of type j will be generated.
          If NTYPES is smaller than the maximum number of types
          defined (PARAMETER MAXTYP), then types NTYPES+1 through
          MAXTYP will not be generated.  If NTYPES is larger
          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
          will be ignored.
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry ISEED specifies the seed of the random number
          generator. The array elements should be between 0 and 4095;
          if not they will be reduced mod 4096.  Also, ISEED(4) must
          be odd.  The random number generator uses a linear
          congruential sequence limited to small integers, and so
          should produce machine independent random numbers. The
          values of ISEED are changed on exit, and can be used in the
          next call to DDRVSX to continue the same random number
          sequence.
[in]THRESH
          THRESH is DOUBLE PRECISION
          A test will count as "failed" if the "error", computed as
          described above, exceeds THRESH.  Note that the error
          is scaled to be O(1), so THRESH should be a reasonably
          small multiple of 1, e.g., 10 or 100.  In particular,
          it should not depend on the precision (single vs. double)
          or the size of the matrix.  It must be at least zero.
[in]NIUNIT
          NIUNIT is INTEGER
          The FORTRAN unit number for reading in the data file of
          problems to solve.
[in]NOUNIT
          NOUNIT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns INFO not equal to 0.)
[out]A
          A is DOUBLE PRECISION array, dimension (LDA, max(NN))
          Used to hold the matrix whose eigenvalues are to be
          computed.  On exit, A contains the last matrix actually used.
[in]LDA
          LDA is INTEGER
          The leading dimension of A, and H. LDA must be at
          least 1 and at least max( NN ).
[out]H
          H is DOUBLE PRECISION array, dimension (LDA, max(NN))
          Another copy of the test matrix A, modified by DGEESX.
[out]HT
          HT is DOUBLE PRECISION array, dimension (LDA, max(NN))
          Yet another copy of the test matrix A, modified by DGEESX.
[out]WR
          WR is DOUBLE PRECISION array, dimension (max(NN))
[out]WI
          WI is DOUBLE PRECISION array, dimension (max(NN))

          The real and imaginary parts of the eigenvalues of A.
          On exit, WR + WI*i are the eigenvalues of the matrix in A.
[out]WRT
          WRT is DOUBLE PRECISION array, dimension (max(NN))
[out]WIT
          WIT is DOUBLE PRECISION array, dimension (max(NN))

          Like WR, WI, these arrays contain the eigenvalues of A,
          but those computed when DGEESX only computes a partial
          eigendecomposition, i.e. not Schur vectors
[out]WRTMP
          WRTMP is DOUBLE PRECISION array, dimension (max(NN))
[out]WITMP
          WITMP is DOUBLE PRECISION array, dimension (max(NN))

          More temporary storage for eigenvalues.
[out]VS
          VS is DOUBLE PRECISION array, dimension (LDVS, max(NN))
          VS holds the computed Schur vectors.
[in]LDVS
          LDVS is INTEGER
          Leading dimension of VS. Must be at least max(1,max(NN)).
[out]VS1
          VS1 is DOUBLE PRECISION array, dimension (LDVS, max(NN))
          VS1 holds another copy of the computed Schur vectors.
[out]RESULT
          RESULT is DOUBLE PRECISION array, dimension (17)
          The values computed by the 17 tests described above.
          The values are currently limited to 1/ulp, to avoid overflow.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK.  This must be at least
          max(3*NN(j),2*NN(j)**2) for all j.
[out]IWORK
          IWORK is INTEGER array, dimension (max(NN)*max(NN))
[out]BWORK
          BWORK is LOGICAL array, dimension (max(NN))
[out]INFO
          INFO is INTEGER
          If 0,  successful exit.
            <0,  input parameter -INFO is incorrect
            >0,  DLATMR, SLATMS, SLATME or DGET24 returned an error
                 code and INFO is its absolute value

-----------------------------------------------------------------------

     Some Local Variables and Parameters:
     ---- ----- --------- --- ----------
     ZERO, ONE       Real 0 and 1.
     MAXTYP          The number of types defined.
     NMAX            Largest value in NN.
     NERRS           The number of tests which have exceeded THRESH
     COND, CONDS,
     IMODE           Values to be passed to the matrix generators.
     ANORM           Norm of A; passed to matrix generators.

     OVFL, UNFL      Overflow and underflow thresholds.
     ULP, ULPINV     Finest relative precision and its inverse.
     RTULP, RTULPI   Square roots of the previous 4 values.
             The following four arrays decode JTYPE:
     KTYPE(j)        The general type (1-10) for type "j".
     KMODE(j)        The MODE value to be passed to the matrix
                     generator for type "j".
     KMAGN(j)        The order of magnitude ( O(1),
                     O(overflow^(1/2) ), O(underflow^(1/2) )
     KCONDS(j)       Selectw whether CONDS is to be 1 or
                     1/sqrt(ulp).  (0 means irrelevant.)
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 450 of file ddrvsx.f.

454*
455* -- LAPACK test routine --
456* -- LAPACK is a software package provided by Univ. of Tennessee, --
457* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
458*
459* .. Scalar Arguments ..
460 INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
461 $ NTYPES
462 DOUBLE PRECISION THRESH
463* ..
464* .. Array Arguments ..
465 LOGICAL BWORK( * ), DOTYPE( * )
466 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
467 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ),
468 $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
469 $ WI( * ), WIT( * ), WITMP( * ), WORK( * ),
470 $ WR( * ), WRT( * ), WRTMP( * )
471* ..
472*
473* =====================================================================
474*
475* .. Parameters ..
476 DOUBLE PRECISION ZERO, ONE
477 parameter( zero = 0.0d0, one = 1.0d0 )
478 INTEGER MAXTYP
479 parameter( maxtyp = 21 )
480* ..
481* .. Local Scalars ..
482 LOGICAL BADNN
483 CHARACTER*3 PATH
484 INTEGER I, IINFO, IMODE, ITYPE, IWK, J, JCOL, JSIZE,
485 $ JTYPE, MTYPES, N, NERRS, NFAIL, NMAX, NNWORK,
486 $ NSLCT, NTEST, NTESTF, NTESTT
487 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
488 $ RTULP, RTULPI, ULP, ULPINV, UNFL
489* ..
490* .. Local Arrays ..
491 CHARACTER ADUMMA( 1 )
492 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
493 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ),
494 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
495* ..
496* .. Arrays in Common ..
497 LOGICAL SELVAL( 20 )
498 DOUBLE PRECISION SELWI( 20 ), SELWR( 20 )
499* ..
500* .. Scalars in Common ..
501 INTEGER SELDIM, SELOPT
502* ..
503* .. Common blocks ..
504 COMMON / sslct / selopt, seldim, selval, selwr, selwi
505* ..
506* .. External Functions ..
507 DOUBLE PRECISION DLAMCH
508 EXTERNAL dlamch
509* ..
510* .. External Subroutines ..
511 EXTERNAL dget24, dlabad, dlaset, dlasum, dlatme, dlatmr,
512 $ dlatms, xerbla
513* ..
514* .. Intrinsic Functions ..
515 INTRINSIC abs, max, min, sqrt
516* ..
517* .. Data statements ..
518 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
519 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
520 $ 3, 1, 2, 3 /
521 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
522 $ 1, 5, 5, 5, 4, 3, 1 /
523 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
524* ..
525* .. Executable Statements ..
526*
527 path( 1: 1 ) = 'Double precision'
528 path( 2: 3 ) = 'SX'
529*
530* Check for errors
531*
532 ntestt = 0
533 ntestf = 0
534 info = 0
535*
536* Important constants
537*
538 badnn = .false.
539*
540* 12 is the largest dimension in the input file of precomputed
541* problems
542*
543 nmax = 12
544 DO 10 j = 1, nsizes
545 nmax = max( nmax, nn( j ) )
546 IF( nn( j ).LT.0 )
547 $ badnn = .true.
548 10 CONTINUE
549*
550* Check for errors
551*
552 IF( nsizes.LT.0 ) THEN
553 info = -1
554 ELSE IF( badnn ) THEN
555 info = -2
556 ELSE IF( ntypes.LT.0 ) THEN
557 info = -3
558 ELSE IF( thresh.LT.zero ) THEN
559 info = -6
560 ELSE IF( niunit.LE.0 ) THEN
561 info = -7
562 ELSE IF( nounit.LE.0 ) THEN
563 info = -8
564 ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
565 info = -10
566 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax ) THEN
567 info = -20
568 ELSE IF( max( 3*nmax, 2*nmax**2 ).GT.lwork ) THEN
569 info = -24
570 END IF
571*
572 IF( info.NE.0 ) THEN
573 CALL xerbla( 'DDRVSX', -info )
574 RETURN
575 END IF
576*
577* If nothing to do check on NIUNIT
578*
579 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
580 $ GO TO 150
581*
582* More Important constants
583*
584 unfl = dlamch( 'Safe minimum' )
585 ovfl = one / unfl
586 CALL dlabad( unfl, ovfl )
587 ulp = dlamch( 'Precision' )
588 ulpinv = one / ulp
589 rtulp = sqrt( ulp )
590 rtulpi = one / rtulp
591*
592* Loop over sizes, types
593*
594 nerrs = 0
595*
596 DO 140 jsize = 1, nsizes
597 n = nn( jsize )
598 IF( nsizes.NE.1 ) THEN
599 mtypes = min( maxtyp, ntypes )
600 ELSE
601 mtypes = min( maxtyp+1, ntypes )
602 END IF
603*
604 DO 130 jtype = 1, mtypes
605 IF( .NOT.dotype( jtype ) )
606 $ GO TO 130
607*
608* Save ISEED in case of an error.
609*
610 DO 20 j = 1, 4
611 ioldsd( j ) = iseed( j )
612 20 CONTINUE
613*
614* Compute "A"
615*
616* Control parameters:
617*
618* KMAGN KCONDS KMODE KTYPE
619* =1 O(1) 1 clustered 1 zero
620* =2 large large clustered 2 identity
621* =3 small exponential Jordan
622* =4 arithmetic diagonal, (w/ eigenvalues)
623* =5 random log symmetric, w/ eigenvalues
624* =6 random general, w/ eigenvalues
625* =7 random diagonal
626* =8 random symmetric
627* =9 random general
628* =10 random triangular
629*
630 IF( mtypes.GT.maxtyp )
631 $ GO TO 90
632*
633 itype = ktype( jtype )
634 imode = kmode( jtype )
635*
636* Compute norm
637*
638 GO TO ( 30, 40, 50 )kmagn( jtype )
639*
640 30 CONTINUE
641 anorm = one
642 GO TO 60
643*
644 40 CONTINUE
645 anorm = ovfl*ulp
646 GO TO 60
647*
648 50 CONTINUE
649 anorm = unfl*ulpinv
650 GO TO 60
651*
652 60 CONTINUE
653*
654 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
655 iinfo = 0
656 cond = ulpinv
657*
658* Special Matrices -- Identity & Jordan block
659*
660* Zero
661*
662 IF( itype.EQ.1 ) THEN
663 iinfo = 0
664*
665 ELSE IF( itype.EQ.2 ) THEN
666*
667* Identity
668*
669 DO 70 jcol = 1, n
670 a( jcol, jcol ) = anorm
671 70 CONTINUE
672*
673 ELSE IF( itype.EQ.3 ) THEN
674*
675* Jordan Block
676*
677 DO 80 jcol = 1, n
678 a( jcol, jcol ) = anorm
679 IF( jcol.GT.1 )
680 $ a( jcol, jcol-1 ) = one
681 80 CONTINUE
682*
683 ELSE IF( itype.EQ.4 ) THEN
684*
685* Diagonal Matrix, [Eigen]values Specified
686*
687 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
688 $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
689 $ iinfo )
690*
691 ELSE IF( itype.EQ.5 ) THEN
692*
693* Symmetric, eigenvalues specified
694*
695 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
696 $ anorm, n, n, 'N', a, lda, work( n+1 ),
697 $ iinfo )
698*
699 ELSE IF( itype.EQ.6 ) THEN
700*
701* General, eigenvalues specified
702*
703 IF( kconds( jtype ).EQ.1 ) THEN
704 conds = one
705 ELSE IF( kconds( jtype ).EQ.2 ) THEN
706 conds = rtulpi
707 ELSE
708 conds = zero
709 END IF
710*
711 adumma( 1 ) = ' '
712 CALL dlatme( n, 'S', iseed, work, imode, cond, one,
713 $ adumma, 'T', 'T', 'T', work( n+1 ), 4,
714 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
715 $ iinfo )
716*
717 ELSE IF( itype.EQ.7 ) THEN
718*
719* Diagonal, random eigenvalues
720*
721 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
722 $ 'T', 'N', work( n+1 ), 1, one,
723 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
724 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
725*
726 ELSE IF( itype.EQ.8 ) THEN
727*
728* Symmetric, random eigenvalues
729*
730 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
731 $ 'T', 'N', work( n+1 ), 1, one,
732 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
733 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
734*
735 ELSE IF( itype.EQ.9 ) THEN
736*
737* General, random eigenvalues
738*
739 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
740 $ 'T', 'N', work( n+1 ), 1, one,
741 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
742 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
743 IF( n.GE.4 ) THEN
744 CALL dlaset( 'Full', 2, n, zero, zero, a, lda )
745 CALL dlaset( 'Full', n-3, 1, zero, zero, a( 3, 1 ),
746 $ lda )
747 CALL dlaset( 'Full', n-3, 2, zero, zero, a( 3, n-1 ),
748 $ lda )
749 CALL dlaset( 'Full', 1, n, zero, zero, a( n, 1 ),
750 $ lda )
751 END IF
752*
753 ELSE IF( itype.EQ.10 ) THEN
754*
755* Triangular, random eigenvalues
756*
757 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
758 $ 'T', 'N', work( n+1 ), 1, one,
759 $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
760 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
761*
762 ELSE
763*
764 iinfo = 1
765 END IF
766*
767 IF( iinfo.NE.0 ) THEN
768 WRITE( nounit, fmt = 9991 )'Generator', iinfo, n, jtype,
769 $ ioldsd
770 info = abs( iinfo )
771 RETURN
772 END IF
773*
774 90 CONTINUE
775*
776* Test for minimal and generous workspace
777*
778 DO 120 iwk = 1, 2
779 IF( iwk.EQ.1 ) THEN
780 nnwork = 3*n
781 ELSE
782 nnwork = max( 3*n, 2*n*n )
783 END IF
784 nnwork = max( nnwork, 1 )
785*
786 CALL dget24( .false., jtype, thresh, ioldsd, nounit, n,
787 $ a, lda, h, ht, wr, wi, wrt, wit, wrtmp,
788 $ witmp, vs, ldvs, vs1, rcdein, rcdvin, nslct,
789 $ islct, result, work, nnwork, iwork, bwork,
790 $ info )
791*
792* Check for RESULT(j) > THRESH
793*
794 ntest = 0
795 nfail = 0
796 DO 100 j = 1, 15
797 IF( result( j ).GE.zero )
798 $ ntest = ntest + 1
799 IF( result( j ).GE.thresh )
800 $ nfail = nfail + 1
801 100 CONTINUE
802*
803 IF( nfail.GT.0 )
804 $ ntestf = ntestf + 1
805 IF( ntestf.EQ.1 ) THEN
806 WRITE( nounit, fmt = 9999 )path
807 WRITE( nounit, fmt = 9998 )
808 WRITE( nounit, fmt = 9997 )
809 WRITE( nounit, fmt = 9996 )
810 WRITE( nounit, fmt = 9995 )thresh
811 WRITE( nounit, fmt = 9994 )
812 ntestf = 2
813 END IF
814*
815 DO 110 j = 1, 15
816 IF( result( j ).GE.thresh ) THEN
817 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
818 $ j, result( j )
819 END IF
820 110 CONTINUE
821*
822 nerrs = nerrs + nfail
823 ntestt = ntestt + ntest
824*
825 120 CONTINUE
826 130 CONTINUE
827 140 CONTINUE
828*
829 150 CONTINUE
830*
831* Read in data from file to check accuracy of condition estimation
832* Read input data until N=0
833*
834 jtype = 0
835 160 CONTINUE
836 READ( niunit, fmt = *, END = 200 )N, nslct
837 IF( n.EQ.0 )
838 $ GO TO 200
839 jtype = jtype + 1
840 iseed( 1 ) = jtype
841 IF( nslct.GT.0 )
842 $ READ( niunit, fmt = * )( islct( i ), i = 1, nslct )
843 DO 170 i = 1, n
844 READ( niunit, fmt = * )( a( i, j ), j = 1, n )
845 170 CONTINUE
846 READ( niunit, fmt = * )rcdein, rcdvin
847*
848 CALL dget24( .true., 22, thresh, iseed, nounit, n, a, lda, h, ht,
849 $ wr, wi, wrt, wit, wrtmp, witmp, vs, ldvs, vs1,
850 $ rcdein, rcdvin, nslct, islct, result, work, lwork,
851 $ iwork, bwork, info )
852*
853* Check for RESULT(j) > THRESH
854*
855 ntest = 0
856 nfail = 0
857 DO 180 j = 1, 17
858 IF( result( j ).GE.zero )
859 $ ntest = ntest + 1
860 IF( result( j ).GE.thresh )
861 $ nfail = nfail + 1
862 180 CONTINUE
863*
864 IF( nfail.GT.0 )
865 $ ntestf = ntestf + 1
866 IF( ntestf.EQ.1 ) THEN
867 WRITE( nounit, fmt = 9999 )path
868 WRITE( nounit, fmt = 9998 )
869 WRITE( nounit, fmt = 9997 )
870 WRITE( nounit, fmt = 9996 )
871 WRITE( nounit, fmt = 9995 )thresh
872 WRITE( nounit, fmt = 9994 )
873 ntestf = 2
874 END IF
875 DO 190 j = 1, 17
876 IF( result( j ).GE.thresh ) THEN
877 WRITE( nounit, fmt = 9992 )n, jtype, j, result( j )
878 END IF
879 190 CONTINUE
880*
881 nerrs = nerrs + nfail
882 ntestt = ntestt + ntest
883 GO TO 160
884 200 CONTINUE
885*
886* Summary
887*
888 CALL dlasum( path, nounit, nerrs, ntestt )
889*
890 9999 FORMAT( / 1x, a3, ' -- Real Schur Form Decomposition Expert ',
891 $ 'Driver', / ' Matrix types (see DDRVSX for details):' )
892*
893 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
894 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
895 $ / ' 2=Identity matrix. ', ' 6=Diagona',
896 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
897 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
898 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
899 $ 'mall, evenly spaced.' )
900 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
901 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
902 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
903 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
904 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
905 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
906 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
907 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
908 $ ' complx ' )
909 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
910 $ 'with small random entries.', / ' 20=Matrix with large ran',
911 $ 'dom entries. ', / )
912 9995 FORMAT( ' Tests performed with test threshold =', f8.2,
913 $ / ' ( A denotes A on input and T denotes A on output)',
914 $ / / ' 1 = 0 if T in Schur form (no sort), ',
915 $ ' 1/ulp otherwise', /
916 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
917 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
918 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
919 $ ' 1/ulp otherwise', /
920 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
921 $ ' 1/ulp otherwise', /
922 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
923 $ ', 1/ulp otherwise' )
924 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
925 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
926 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
927 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
928 $ ' 1/ulp otherwise', /
929 $ ' 11 = 0 if T same no matter what else computed (sort),',
930 $ ' 1/ulp otherwise', /
931 $ ' 12 = 0 if WR, WI same no matter what else computed ',
932 $ '(sort), 1/ulp otherwise', /
933 $ ' 13 = 0 if sorting successful, 1/ulp otherwise',
934 $ / ' 14 = 0 if RCONDE same no matter what else computed,',
935 $ ' 1/ulp otherwise', /
936 $ ' 15 = 0 if RCONDv same no matter what else computed,',
937 $ ' 1/ulp otherwise', /
938 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
939 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
940 9993 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
941 $ ' type ', i2, ', test(', i2, ')=', g10.3 )
942 9992 FORMAT( ' N=', i5, ', input example =', i3, ', test(', i2, ')=',
943 $ g10.3 )
944 9991 FORMAT( ' DDRVSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
945 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
946*
947 RETURN
948*
949* End of DDRVSX
950*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
Definition: dlasum.f:43
subroutine dget24(COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK, IWORK, BWORK, INFO)
DGET24
Definition: dget24.f:343
subroutine dlatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
DLATMR
Definition: dlatmr.f:471
subroutine dlatme(N, DIST, ISEED, D, MODE, COND, DMAX, EI, RSIGN, UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, A, LDA, WORK, INFO)
DLATME
Definition: dlatme.f:332
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:321
Here is the call graph for this function:
Here is the caller graph for this function: