LAPACK 3.12.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, 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 ulp = dlamch( 'Precision' )
587 ulpinv = one / ulp
588 rtulp = sqrt( ulp )
589 rtulpi = one / rtulp
590*
591* Loop over sizes, types
592*
593 nerrs = 0
594*
595 DO 140 jsize = 1, nsizes
596 n = nn( jsize )
597 IF( nsizes.NE.1 ) THEN
598 mtypes = min( maxtyp, ntypes )
599 ELSE
600 mtypes = min( maxtyp+1, ntypes )
601 END IF
602*
603 DO 130 jtype = 1, mtypes
604 IF( .NOT.dotype( jtype ) )
605 $ GO TO 130
606*
607* Save ISEED in case of an error.
608*
609 DO 20 j = 1, 4
610 ioldsd( j ) = iseed( j )
611 20 CONTINUE
612*
613* Compute "A"
614*
615* Control parameters:
616*
617* KMAGN KCONDS KMODE KTYPE
618* =1 O(1) 1 clustered 1 zero
619* =2 large large clustered 2 identity
620* =3 small exponential Jordan
621* =4 arithmetic diagonal, (w/ eigenvalues)
622* =5 random log symmetric, w/ eigenvalues
623* =6 random general, w/ eigenvalues
624* =7 random diagonal
625* =8 random symmetric
626* =9 random general
627* =10 random triangular
628*
629 IF( mtypes.GT.maxtyp )
630 $ GO TO 90
631*
632 itype = ktype( jtype )
633 imode = kmode( jtype )
634*
635* Compute norm
636*
637 GO TO ( 30, 40, 50 )kmagn( jtype )
638*
639 30 CONTINUE
640 anorm = one
641 GO TO 60
642*
643 40 CONTINUE
644 anorm = ovfl*ulp
645 GO TO 60
646*
647 50 CONTINUE
648 anorm = unfl*ulpinv
649 GO TO 60
650*
651 60 CONTINUE
652*
653 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
654 iinfo = 0
655 cond = ulpinv
656*
657* Special Matrices -- Identity & Jordan block
658*
659* Zero
660*
661 IF( itype.EQ.1 ) THEN
662 iinfo = 0
663*
664 ELSE IF( itype.EQ.2 ) THEN
665*
666* Identity
667*
668 DO 70 jcol = 1, n
669 a( jcol, jcol ) = anorm
670 70 CONTINUE
671*
672 ELSE IF( itype.EQ.3 ) THEN
673*
674* Jordan Block
675*
676 DO 80 jcol = 1, n
677 a( jcol, jcol ) = anorm
678 IF( jcol.GT.1 )
679 $ a( jcol, jcol-1 ) = one
680 80 CONTINUE
681*
682 ELSE IF( itype.EQ.4 ) THEN
683*
684* Diagonal Matrix, [Eigen]values Specified
685*
686 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
687 $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
688 $ iinfo )
689*
690 ELSE IF( itype.EQ.5 ) THEN
691*
692* Symmetric, eigenvalues specified
693*
694 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
695 $ anorm, n, n, 'N', a, lda, work( n+1 ),
696 $ iinfo )
697*
698 ELSE IF( itype.EQ.6 ) THEN
699*
700* General, eigenvalues specified
701*
702 IF( kconds( jtype ).EQ.1 ) THEN
703 conds = one
704 ELSE IF( kconds( jtype ).EQ.2 ) THEN
705 conds = rtulpi
706 ELSE
707 conds = zero
708 END IF
709*
710 adumma( 1 ) = ' '
711 CALL dlatme( n, 'S', iseed, work, imode, cond, one,
712 $ adumma, 'T', 'T', 'T', work( n+1 ), 4,
713 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
714 $ iinfo )
715*
716 ELSE IF( itype.EQ.7 ) THEN
717*
718* Diagonal, random eigenvalues
719*
720 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
721 $ 'T', 'N', work( n+1 ), 1, one,
722 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
723 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
724*
725 ELSE IF( itype.EQ.8 ) THEN
726*
727* Symmetric, random eigenvalues
728*
729 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
730 $ 'T', 'N', work( n+1 ), 1, one,
731 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
732 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
733*
734 ELSE IF( itype.EQ.9 ) THEN
735*
736* General, random eigenvalues
737*
738 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
739 $ 'T', 'N', work( n+1 ), 1, one,
740 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
741 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
742 IF( n.GE.4 ) THEN
743 CALL dlaset( 'Full', 2, n, zero, zero, a, lda )
744 CALL dlaset( 'Full', n-3, 1, zero, zero, a( 3, 1 ),
745 $ lda )
746 CALL dlaset( 'Full', n-3, 2, zero, zero, a( 3, n-1 ),
747 $ lda )
748 CALL dlaset( 'Full', 1, n, zero, zero, a( n, 1 ),
749 $ lda )
750 END IF
751*
752 ELSE IF( itype.EQ.10 ) THEN
753*
754* Triangular, random eigenvalues
755*
756 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
757 $ 'T', 'N', work( n+1 ), 1, one,
758 $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
759 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
760*
761 ELSE
762*
763 iinfo = 1
764 END IF
765*
766 IF( iinfo.NE.0 ) THEN
767 WRITE( nounit, fmt = 9991 )'Generator', iinfo, n, jtype,
768 $ ioldsd
769 info = abs( iinfo )
770 RETURN
771 END IF
772*
773 90 CONTINUE
774*
775* Test for minimal and generous workspace
776*
777 DO 120 iwk = 1, 2
778 IF( iwk.EQ.1 ) THEN
779 nnwork = 3*n
780 ELSE
781 nnwork = max( 3*n, 2*n*n )
782 END IF
783 nnwork = max( nnwork, 1 )
784*
785 CALL dget24( .false., jtype, thresh, ioldsd, nounit, n,
786 $ a, lda, h, ht, wr, wi, wrt, wit, wrtmp,
787 $ witmp, vs, ldvs, vs1, rcdein, rcdvin, nslct,
788 $ islct, result, work, nnwork, iwork, bwork,
789 $ info )
790*
791* Check for RESULT(j) > THRESH
792*
793 ntest = 0
794 nfail = 0
795 DO 100 j = 1, 15
796 IF( result( j ).GE.zero )
797 $ ntest = ntest + 1
798 IF( result( j ).GE.thresh )
799 $ nfail = nfail + 1
800 100 CONTINUE
801*
802 IF( nfail.GT.0 )
803 $ ntestf = ntestf + 1
804 IF( ntestf.EQ.1 ) THEN
805 WRITE( nounit, fmt = 9999 )path
806 WRITE( nounit, fmt = 9998 )
807 WRITE( nounit, fmt = 9997 )
808 WRITE( nounit, fmt = 9996 )
809 WRITE( nounit, fmt = 9995 )thresh
810 WRITE( nounit, fmt = 9994 )
811 ntestf = 2
812 END IF
813*
814 DO 110 j = 1, 15
815 IF( result( j ).GE.thresh ) THEN
816 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
817 $ j, result( j )
818 END IF
819 110 CONTINUE
820*
821 nerrs = nerrs + nfail
822 ntestt = ntestt + ntest
823*
824 120 CONTINUE
825 130 CONTINUE
826 140 CONTINUE
827*
828 150 CONTINUE
829*
830* Read in data from file to check accuracy of condition estimation
831* Read input data until N=0
832*
833 jtype = 0
834 160 CONTINUE
835 READ( niunit, fmt = *, END = 200 )N, nslct
836 IF( n.EQ.0 )
837 $ GO TO 200
838 jtype = jtype + 1
839 iseed( 1 ) = jtype
840 IF( nslct.GT.0 )
841 $ READ( niunit, fmt = * )( islct( i ), i = 1, nslct )
842 DO 170 i = 1, n
843 READ( niunit, fmt = * )( a( i, j ), j = 1, n )
844 170 CONTINUE
845 READ( niunit, fmt = * )rcdein, rcdvin
846*
847 CALL dget24( .true., 22, thresh, iseed, nounit, n, a, lda, h, ht,
848 $ wr, wi, wrt, wit, wrtmp, witmp, vs, ldvs, vs1,
849 $ rcdein, rcdvin, nslct, islct, result, work, lwork,
850 $ iwork, bwork, info )
851*
852* Check for RESULT(j) > THRESH
853*
854 ntest = 0
855 nfail = 0
856 DO 180 j = 1, 17
857 IF( result( j ).GE.zero )
858 $ ntest = ntest + 1
859 IF( result( j ).GE.thresh )
860 $ nfail = nfail + 1
861 180 CONTINUE
862*
863 IF( nfail.GT.0 )
864 $ ntestf = ntestf + 1
865 IF( ntestf.EQ.1 ) THEN
866 WRITE( nounit, fmt = 9999 )path
867 WRITE( nounit, fmt = 9998 )
868 WRITE( nounit, fmt = 9997 )
869 WRITE( nounit, fmt = 9996 )
870 WRITE( nounit, fmt = 9995 )thresh
871 WRITE( nounit, fmt = 9994 )
872 ntestf = 2
873 END IF
874 DO 190 j = 1, 17
875 IF( result( j ).GE.thresh ) THEN
876 WRITE( nounit, fmt = 9992 )n, jtype, j, result( j )
877 END IF
878 190 CONTINUE
879*
880 nerrs = nerrs + nfail
881 ntestt = ntestt + ntest
882 GO TO 160
883 200 CONTINUE
884*
885* Summary
886*
887 CALL dlasum( path, nounit, nerrs, ntestt )
888*
889 9999 FORMAT( / 1x, a3, ' -- Real Schur Form Decomposition Expert ',
890 $ 'Driver', / ' Matrix types (see DDRVSX for details):' )
891*
892 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
893 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
894 $ / ' 2=Identity matrix. ', ' 6=Diagona',
895 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
896 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
897 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
898 $ 'mall, evenly spaced.' )
899 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
900 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
901 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
902 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
903 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
904 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
905 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
906 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
907 $ ' complx ' )
908 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
909 $ 'with small random entries.', / ' 20=Matrix with large ran',
910 $ 'dom entries. ', / )
911 9995 FORMAT( ' Tests performed with test threshold =', f8.2,
912 $ / ' ( A denotes A on input and T denotes A on output)',
913 $ / / ' 1 = 0 if T in Schur form (no sort), ',
914 $ ' 1/ulp otherwise', /
915 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
916 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
917 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
918 $ ' 1/ulp otherwise', /
919 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
920 $ ' 1/ulp otherwise', /
921 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
922 $ ', 1/ulp otherwise' )
923 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
924 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
925 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
926 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
927 $ ' 1/ulp otherwise', /
928 $ ' 11 = 0 if T same no matter what else computed (sort),',
929 $ ' 1/ulp otherwise', /
930 $ ' 12 = 0 if WR, WI same no matter what else computed ',
931 $ '(sort), 1/ulp otherwise', /
932 $ ' 13 = 0 if sorting successful, 1/ulp otherwise',
933 $ / ' 14 = 0 if RCONDE same no matter what else computed,',
934 $ ' 1/ulp otherwise', /
935 $ ' 15 = 0 if RCONDv same no matter what else computed,',
936 $ ' 1/ulp otherwise', /
937 $ ' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
938 $ / ' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
939 9993 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
940 $ ' type ', i2, ', test(', i2, ')=', g10.3 )
941 9992 FORMAT( ' N=', i5, ', input example =', i3, ', test(', i2, ')=',
942 $ g10.3 )
943 9991 FORMAT( ' DDRVSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
944 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
945*
946 RETURN
947*
948* End of DDRVSX
949*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 dlasum(type, iounit, ie, nrun)
DLASUM
Definition dlasum.f:43
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 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 dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
Definition dlatms.f:321
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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
Here is the call graph for this function:
Here is the caller graph for this function: