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

◆ ddrves()

subroutine ddrves ( integer  nsizes,
integer, dimension( * )  nn,
integer  ntypes,
logical, dimension( * )  dotype,
integer, dimension( 4 )  iseed,
double precision  thresh,
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( ldvs, * )  vs,
integer  ldvs,
double precision, dimension( 13 )  result,
double precision, dimension( * )  work,
integer  nwork,
integer, dimension( * )  iwork,
logical, dimension( * )  bwork,
integer  info 
)

DDRVES

Purpose:
    DDRVES checks the nonsymmetric eigenvalue (Schur form) problem
    driver DGEES.

    When DDRVES 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, 13
    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
            (with sorting of eigenvalues)

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

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

    (13)    if sorting worked and SDIM is the number of
            eigenvalues which were SELECTed

    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
Parameters
[in]NSIZES
          NSIZES is INTEGER
          The number of sizes of matrices to use.  If it is zero,
          DDRVES does nothing.  It must be at least zero.
[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.   If it is zero, DDRVES
          does nothing.  It must be at least zero.  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 DDRVES 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]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 DGEES.
[out]HT
          HT is DOUBLE PRECISION array, dimension (LDA, max(NN))
          Yet another copy of the test matrix A, modified by DGEES.
[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 DGEES only computes a partial
          eigendecomposition, i.e. not Schur vectors
[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]RESULT
          RESULT is DOUBLE PRECISION array, dimension (13)
          The values computed by the 13 tests described above.
          The values are currently limited to 1/ulp, to avoid overflow.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (NWORK)
[in]NWORK
          NWORK is INTEGER
          The number of entries in WORK.  This must be at least
          5*NN(j)+2*NN(j)**2 for all j.
[out]IWORK
          IWORK is INTEGER array, dimension (max(NN))
[out]BWORK
          BWORK is LOGICAL array, dimension (max(NN))
[out]INFO
          INFO is INTEGER
          If 0, then everything ran OK.
           -1: NSIZES < 0
           -2: Some NN(j) < 0
           -3: NTYPES < 0
           -6: THRESH < 0
           -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
          -17: LDVS < 1 or LDVS < NMAX, where NMAX is max( NN(j) ).
          -20: NWORK too small.
          If  DLATMR, SLATMS, SLATME or DGEES returns an error code,
              the absolute value of it is returned.

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

     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 385 of file ddrves.f.

388*
389* -- LAPACK test routine --
390* -- LAPACK is a software package provided by Univ. of Tennessee, --
391* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
392*
393* .. Scalar Arguments ..
394 INTEGER INFO, LDA, LDVS, NOUNIT, NSIZES, NTYPES, NWORK
395 DOUBLE PRECISION THRESH
396* ..
397* .. Array Arguments ..
398 LOGICAL BWORK( * ), DOTYPE( * )
399 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
400 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ),
401 $ RESULT( 13 ), VS( LDVS, * ), WI( * ), WIT( * ),
402 $ WORK( * ), WR( * ), WRT( * )
403* ..
404*
405* =====================================================================
406*
407* .. Parameters ..
408 DOUBLE PRECISION ZERO, ONE
409 parameter( zero = 0.0d0, one = 1.0d0 )
410 INTEGER MAXTYP
411 parameter( maxtyp = 21 )
412* ..
413* .. Local Scalars ..
414 LOGICAL BADNN
415 CHARACTER SORT
416 CHARACTER*3 PATH
417 INTEGER I, IINFO, IMODE, ISORT, ITYPE, IWK, J, JCOL,
418 $ JSIZE, JTYPE, KNTEIG, LWORK, MTYPES, N, NERRS,
419 $ NFAIL, NMAX, NNWORK, NTEST, NTESTF, NTESTT,
420 $ RSUB, SDIM
421 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RTULP, RTULPI, TMP,
422 $ ULP, ULPINV, UNFL
423* ..
424* .. Local Arrays ..
425 CHARACTER ADUMMA( 1 )
426 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
427 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
428 $ KTYPE( MAXTYP )
429 DOUBLE PRECISION RES( 2 )
430* ..
431* .. Arrays in Common ..
432 LOGICAL SELVAL( 20 )
433 DOUBLE PRECISION SELWI( 20 ), SELWR( 20 )
434* ..
435* .. Scalars in Common ..
436 INTEGER SELDIM, SELOPT
437* ..
438* .. Common blocks ..
439 COMMON / sslct / selopt, seldim, selval, selwr, selwi
440* ..
441* .. External Functions ..
442 LOGICAL DSLECT
443 DOUBLE PRECISION DLAMCH
444 EXTERNAL dslect, dlamch
445* ..
446* .. External Subroutines ..
447 EXTERNAL dgees, dhst01, dlacpy, dlaset, dlasum, dlatme,
449* ..
450* .. Intrinsic Functions ..
451 INTRINSIC abs, max, sign, sqrt
452* ..
453* .. Data statements ..
454 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
455 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
456 $ 3, 1, 2, 3 /
457 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
458 $ 1, 5, 5, 5, 4, 3, 1 /
459 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
460* ..
461* .. Executable Statements ..
462*
463 path( 1: 1 ) = 'Double precision'
464 path( 2: 3 ) = 'ES'
465*
466* Check for errors
467*
468 ntestt = 0
469 ntestf = 0
470 info = 0
471 selopt = 0
472*
473* Important constants
474*
475 badnn = .false.
476 nmax = 0
477 DO 10 j = 1, nsizes
478 nmax = max( nmax, nn( j ) )
479 IF( nn( j ).LT.0 )
480 $ badnn = .true.
481 10 CONTINUE
482*
483* Check for errors
484*
485 IF( nsizes.LT.0 ) THEN
486 info = -1
487 ELSE IF( badnn ) THEN
488 info = -2
489 ELSE IF( ntypes.LT.0 ) THEN
490 info = -3
491 ELSE IF( thresh.LT.zero ) THEN
492 info = -6
493 ELSE IF( nounit.LE.0 ) THEN
494 info = -7
495 ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
496 info = -9
497 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax ) THEN
498 info = -17
499 ELSE IF( 5*nmax+2*nmax**2.GT.nwork ) THEN
500 info = -20
501 END IF
502*
503 IF( info.NE.0 ) THEN
504 CALL xerbla( 'DDRVES', -info )
505 RETURN
506 END IF
507*
508* Quick return if nothing to do
509*
510 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
511 $ RETURN
512*
513* More Important constants
514*
515 unfl = dlamch( 'Safe minimum' )
516 ovfl = one / unfl
517 ulp = dlamch( 'Precision' )
518 ulpinv = one / ulp
519 rtulp = sqrt( ulp )
520 rtulpi = one / rtulp
521*
522* Loop over sizes, types
523*
524 nerrs = 0
525*
526 DO 270 jsize = 1, nsizes
527 n = nn( jsize )
528 mtypes = maxtyp
529 IF( nsizes.EQ.1 .AND. ntypes.EQ.maxtyp+1 )
530 $ mtypes = mtypes + 1
531*
532 DO 260 jtype = 1, mtypes
533 IF( .NOT.dotype( jtype ) )
534 $ GO TO 260
535*
536* Save ISEED in case of an error.
537*
538 DO 20 j = 1, 4
539 ioldsd( j ) = iseed( j )
540 20 CONTINUE
541*
542* Compute "A"
543*
544* Control parameters:
545*
546* KMAGN KCONDS KMODE KTYPE
547* =1 O(1) 1 clustered 1 zero
548* =2 large large clustered 2 identity
549* =3 small exponential Jordan
550* =4 arithmetic diagonal, (w/ eigenvalues)
551* =5 random log symmetric, w/ eigenvalues
552* =6 random general, w/ eigenvalues
553* =7 random diagonal
554* =8 random symmetric
555* =9 random general
556* =10 random triangular
557*
558 IF( mtypes.GT.maxtyp )
559 $ GO TO 90
560*
561 itype = ktype( jtype )
562 imode = kmode( jtype )
563*
564* Compute norm
565*
566 GO TO ( 30, 40, 50 )kmagn( jtype )
567*
568 30 CONTINUE
569 anorm = one
570 GO TO 60
571*
572 40 CONTINUE
573 anorm = ovfl*ulp
574 GO TO 60
575*
576 50 CONTINUE
577 anorm = unfl*ulpinv
578 GO TO 60
579*
580 60 CONTINUE
581*
582 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
583 iinfo = 0
584 cond = ulpinv
585*
586* Special Matrices -- Identity & Jordan block
587*
588* Zero
589*
590 IF( itype.EQ.1 ) THEN
591 iinfo = 0
592*
593 ELSE IF( itype.EQ.2 ) THEN
594*
595* Identity
596*
597 DO 70 jcol = 1, n
598 a( jcol, jcol ) = anorm
599 70 CONTINUE
600*
601 ELSE IF( itype.EQ.3 ) THEN
602*
603* Jordan Block
604*
605 DO 80 jcol = 1, n
606 a( jcol, jcol ) = anorm
607 IF( jcol.GT.1 )
608 $ a( jcol, jcol-1 ) = one
609 80 CONTINUE
610*
611 ELSE IF( itype.EQ.4 ) THEN
612*
613* Diagonal Matrix, [Eigen]values Specified
614*
615 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
616 $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
617 $ iinfo )
618*
619 ELSE IF( itype.EQ.5 ) THEN
620*
621* Symmetric, eigenvalues specified
622*
623 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
624 $ anorm, n, n, 'N', a, lda, work( n+1 ),
625 $ iinfo )
626*
627 ELSE IF( itype.EQ.6 ) THEN
628*
629* General, eigenvalues specified
630*
631 IF( kconds( jtype ).EQ.1 ) THEN
632 conds = one
633 ELSE IF( kconds( jtype ).EQ.2 ) THEN
634 conds = rtulpi
635 ELSE
636 conds = zero
637 END IF
638*
639 adumma( 1 ) = ' '
640 CALL dlatme( n, 'S', iseed, work, imode, cond, one,
641 $ adumma, 'T', 'T', 'T', work( n+1 ), 4,
642 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
643 $ iinfo )
644*
645 ELSE IF( itype.EQ.7 ) THEN
646*
647* Diagonal, random eigenvalues
648*
649 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
650 $ 'T', 'N', work( n+1 ), 1, one,
651 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
652 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
653*
654 ELSE IF( itype.EQ.8 ) THEN
655*
656* Symmetric, random eigenvalues
657*
658 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
659 $ 'T', 'N', work( n+1 ), 1, one,
660 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
661 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
662*
663 ELSE IF( itype.EQ.9 ) THEN
664*
665* General, random eigenvalues
666*
667 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
668 $ 'T', 'N', work( n+1 ), 1, one,
669 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
670 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
671 IF( n.GE.4 ) THEN
672 CALL dlaset( 'Full', 2, n, zero, zero, a, lda )
673 CALL dlaset( 'Full', n-3, 1, zero, zero, a( 3, 1 ),
674 $ lda )
675 CALL dlaset( 'Full', n-3, 2, zero, zero, a( 3, n-1 ),
676 $ lda )
677 CALL dlaset( 'Full', 1, n, zero, zero, a( n, 1 ),
678 $ lda )
679 END IF
680*
681 ELSE IF( itype.EQ.10 ) THEN
682*
683* Triangular, random eigenvalues
684*
685 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
686 $ 'T', 'N', work( n+1 ), 1, one,
687 $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
688 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
689*
690 ELSE
691*
692 iinfo = 1
693 END IF
694*
695 IF( iinfo.NE.0 ) THEN
696 WRITE( nounit, fmt = 9992 )'Generator', iinfo, n, jtype,
697 $ ioldsd
698 info = abs( iinfo )
699 RETURN
700 END IF
701*
702 90 CONTINUE
703*
704* Test for minimal and generous workspace
705*
706 DO 250 iwk = 1, 2
707 IF( iwk.EQ.1 ) THEN
708 nnwork = 3*n
709 ELSE
710 nnwork = 5*n + 2*n**2
711 END IF
712 nnwork = max( nnwork, 1 )
713*
714* Initialize RESULT
715*
716 DO 100 j = 1, 13
717 result( j ) = -one
718 100 CONTINUE
719*
720* Test with and without sorting of eigenvalues
721*
722 DO 210 isort = 0, 1
723 IF( isort.EQ.0 ) THEN
724 sort = 'N'
725 rsub = 0
726 ELSE
727 sort = 'S'
728 rsub = 6
729 END IF
730*
731* Compute Schur form and Schur vectors, and test them
732*
733 CALL dlacpy( 'F', n, n, a, lda, h, lda )
734 CALL dgees( 'V', sort, dslect, n, h, lda, sdim, wr,
735 $ wi, vs, ldvs, work, nnwork, bwork, iinfo )
736 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
737 result( 1+rsub ) = ulpinv
738 WRITE( nounit, fmt = 9992 )'DGEES1', iinfo, n,
739 $ jtype, ioldsd
740 info = abs( iinfo )
741 GO TO 220
742 END IF
743*
744* Do Test (1) or Test (7)
745*
746 result( 1+rsub ) = zero
747 DO 120 j = 1, n - 2
748 DO 110 i = j + 2, n
749 IF( h( i, j ).NE.zero )
750 $ result( 1+rsub ) = ulpinv
751 110 CONTINUE
752 120 CONTINUE
753 DO 130 i = 1, n - 2
754 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.
755 $ zero )result( 1+rsub ) = ulpinv
756 130 CONTINUE
757 DO 140 i = 1, n - 1
758 IF( h( i+1, i ).NE.zero ) THEN
759 IF( h( i, i ).NE.h( i+1, i+1 ) .OR.
760 $ h( i, i+1 ).EQ.zero .OR.
761 $ sign( one, h( i+1, i ) ).EQ.
762 $ sign( one, h( i, i+1 ) ) )result( 1+rsub )
763 $ = ulpinv
764 END IF
765 140 CONTINUE
766*
767* Do Tests (2) and (3) or Tests (8) and (9)
768*
769 lwork = max( 1, 2*n*n )
770 CALL dhst01( n, 1, n, a, lda, h, lda, vs, ldvs, work,
771 $ lwork, res )
772 result( 2+rsub ) = res( 1 )
773 result( 3+rsub ) = res( 2 )
774*
775* Do Test (4) or Test (10)
776*
777 result( 4+rsub ) = zero
778 DO 150 i = 1, n
779 IF( h( i, i ).NE.wr( i ) )
780 $ result( 4+rsub ) = ulpinv
781 150 CONTINUE
782 IF( n.GT.1 ) THEN
783 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
784 $ result( 4+rsub ) = ulpinv
785 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
786 $ result( 4+rsub ) = ulpinv
787 END IF
788 DO 160 i = 1, n - 1
789 IF( h( i+1, i ).NE.zero ) THEN
790 tmp = sqrt( abs( h( i+1, i ) ) )*
791 $ sqrt( abs( h( i, i+1 ) ) )
792 result( 4+rsub ) = max( result( 4+rsub ),
793 $ abs( wi( i )-tmp ) /
794 $ max( ulp*tmp, unfl ) )
795 result( 4+rsub ) = max( result( 4+rsub ),
796 $ abs( wi( i+1 )+tmp ) /
797 $ max( ulp*tmp, unfl ) )
798 ELSE IF( i.GT.1 ) THEN
799 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.
800 $ zero .AND. wi( i ).NE.zero )result( 4+rsub )
801 $ = ulpinv
802 END IF
803 160 CONTINUE
804*
805* Do Test (5) or Test (11)
806*
807 CALL dlacpy( 'F', n, n, a, lda, ht, lda )
808 CALL dgees( 'N', sort, dslect, n, ht, lda, sdim, wrt,
809 $ wit, vs, ldvs, work, nnwork, bwork,
810 $ iinfo )
811 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
812 result( 5+rsub ) = ulpinv
813 WRITE( nounit, fmt = 9992 )'DGEES2', iinfo, n,
814 $ jtype, ioldsd
815 info = abs( iinfo )
816 GO TO 220
817 END IF
818*
819 result( 5+rsub ) = zero
820 DO 180 j = 1, n
821 DO 170 i = 1, n
822 IF( h( i, j ).NE.ht( i, j ) )
823 $ result( 5+rsub ) = ulpinv
824 170 CONTINUE
825 180 CONTINUE
826*
827* Do Test (6) or Test (12)
828*
829 result( 6+rsub ) = zero
830 DO 190 i = 1, n
831 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
832 $ result( 6+rsub ) = ulpinv
833 190 CONTINUE
834*
835* Do Test (13)
836*
837 IF( isort.EQ.1 ) THEN
838 result( 13 ) = zero
839 knteig = 0
840 DO 200 i = 1, n
841 IF( dslect( wr( i ), wi( i ) ) .OR.
842 $ dslect( wr( i ), -wi( i ) ) )
843 $ knteig = knteig + 1
844 IF( i.LT.n ) THEN
845 IF( ( dslect( wr( i+1 ),
846 $ wi( i+1 ) ) .OR. dslect( wr( i+1 ),
847 $ -wi( i+1 ) ) ) .AND.
848 $ ( .NOT.( dslect( wr( i ),
849 $ wi( i ) ) .OR. dslect( wr( i ),
850 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )
851 $ result( 13 ) = ulpinv
852 END IF
853 200 CONTINUE
854 IF( sdim.NE.knteig ) THEN
855 result( 13 ) = ulpinv
856 END IF
857 END IF
858*
859 210 CONTINUE
860*
861* End of Loop -- Check for RESULT(j) > THRESH
862*
863 220 CONTINUE
864*
865 ntest = 0
866 nfail = 0
867 DO 230 j = 1, 13
868 IF( result( j ).GE.zero )
869 $ ntest = ntest + 1
870 IF( result( j ).GE.thresh )
871 $ nfail = nfail + 1
872 230 CONTINUE
873*
874 IF( nfail.GT.0 )
875 $ ntestf = ntestf + 1
876 IF( ntestf.EQ.1 ) THEN
877 WRITE( nounit, fmt = 9999 )path
878 WRITE( nounit, fmt = 9998 )
879 WRITE( nounit, fmt = 9997 )
880 WRITE( nounit, fmt = 9996 )
881 WRITE( nounit, fmt = 9995 )thresh
882 WRITE( nounit, fmt = 9994 )
883 ntestf = 2
884 END IF
885*
886 DO 240 j = 1, 13
887 IF( result( j ).GE.thresh ) THEN
888 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
889 $ j, result( j )
890 END IF
891 240 CONTINUE
892*
893 nerrs = nerrs + nfail
894 ntestt = ntestt + ntest
895*
896 250 CONTINUE
897 260 CONTINUE
898 270 CONTINUE
899*
900* Summary
901*
902 CALL dlasum( path, nounit, nerrs, ntestt )
903*
904 9999 FORMAT( / 1x, a3, ' -- Real Schur Form Decomposition Driver',
905 $ / ' Matrix types (see DDRVES for details): ' )
906*
907 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
908 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
909 $ / ' 2=Identity matrix. ', ' 6=Diagona',
910 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
911 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
912 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
913 $ 'mall, evenly spaced.' )
914 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
915 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
916 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
917 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
918 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
919 $ 'lex ', / ' 12=Well-cond., random complex ', 6x, ' ',
920 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
921 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
922 $ ' complx ' )
923 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
924 $ 'with small random entries.', / ' 20=Matrix with large ran',
925 $ 'dom entries. ', / )
926 9995 FORMAT( ' Tests performed with test threshold =', f8.2,
927 $ / ' ( A denotes A on input and T denotes A on output)',
928 $ / / ' 1 = 0 if T in Schur form (no sort), ',
929 $ ' 1/ulp otherwise', /
930 $ ' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
931 $ / ' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
932 $ ' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
933 $ ' 1/ulp otherwise', /
934 $ ' 5 = 0 if T same no matter if VS computed (no sort),',
935 $ ' 1/ulp otherwise', /
936 $ ' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
937 $ ', 1/ulp otherwise' )
938 9994 FORMAT( ' 7 = 0 if T in Schur form (sort), ', ' 1/ulp otherwise',
939 $ / ' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
940 $ / ' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
941 $ / ' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
942 $ ' 1/ulp otherwise', /
943 $ ' 11 = 0 if T same no matter if VS computed (sort),',
944 $ ' 1/ulp otherwise', /
945 $ ' 12 = 0 if WR, WI same no matter if VS computed (sort),',
946 $ ' 1/ulp otherwise', /
947 $ ' 13 = 0 if sorting successful, 1/ulp otherwise', / )
948 9993 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
949 $ ' type ', i2, ', test(', i2, ')=', g10.3 )
950 9992 FORMAT( ' DDRVES: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
951 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
952*
953 RETURN
954*
955* End of DDRVES
956*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dhst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
DHST01
Definition dhst01.f:134
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
logical function dslect(zr, zi)
DSLECT
Definition dslect.f:62
subroutine dgees(jobvs, sort, select, n, a, lda, sdim, wr, wi, vs, ldvs, work, lwork, bwork, info)
DGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE m...
Definition dgees.f:216
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
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: