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

◆ dchkhs()

subroutine dchkhs ( 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, * )  T1,
double precision, dimension( lda, * )  T2,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldu, * )  Z,
double precision, dimension( ldu, * )  UZ,
double precision, dimension( * )  WR1,
double precision, dimension( * )  WI1,
double precision, dimension( * )  WR2,
double precision, dimension( * )  WI2,
double precision, dimension( * )  WR3,
double precision, dimension( * )  WI3,
double precision, dimension( ldu, * )  EVECTL,
double precision, dimension( ldu, * )  EVECTR,
double precision, dimension( ldu, * )  EVECTY,
double precision, dimension( ldu, * )  EVECTX,
double precision, dimension( ldu, * )  UU,
double precision, dimension( * )  TAU,
double precision, dimension( * )  WORK,
integer  NWORK,
integer, dimension( * )  IWORK,
logical, dimension( * )  SELECT,
double precision, dimension( 14 )  RESULT,
integer  INFO 
)

DCHKHS

Purpose:
    DCHKHS  checks the nonsymmetric eigenvalue problem routines.

            DGEHRD factors A as  U H U' , where ' means transpose,
            H is hessenberg, and U is an orthogonal matrix.

            DORGHR generates the orthogonal matrix U.

            DORMHR multiplies a matrix by the orthogonal matrix U.

            DHSEQR factors H as  Z T Z' , where Z is orthogonal and
            T is "quasi-triangular", and the eigenvalue vector W.

            DTREVC computes the left and right eigenvector matrices
            L and R for T.

            DHSEIN computes the left and right eigenvector matrices
            Y and X for H, using inverse iteration.

    When DCHKHS 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, 14
    tests will be performed:

    (1)     | A - U H U**T | / ( |A| n ulp )

    (2)     | I - UU**T | / ( n ulp )

    (3)     | H - Z T Z**T | / ( |H| n ulp )

    (4)     | I - ZZ**T | / ( n ulp )

    (5)     | A - UZ H (UZ)**T | / ( |A| n ulp )

    (6)     | I - UZ (UZ)**T | / ( n ulp )

    (7)     | T(Z computed) - T(Z not computed) | / ( |T| ulp )

    (8)     | W(Z computed) - W(Z not computed) | / ( |W| ulp )

    (9)     | TR - RW | / ( |T| |R| ulp )

    (10)    | L**H T - W**H L | / ( |T| |L| ulp )

    (11)    | HX - XW | / ( |H| |X| ulp )

    (12)    | Y**H H - W**H Y | / ( |H| |Y| ulp )

    (13)    | AX - XW | / ( |A| |X| ulp )

    (14)    | Y**H A - W**H Y | / ( |A| |Y| ulp )

    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 SQRT( overflow threshold )
    (8)  Same as (4), but multiplied by SQRT( 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 SQRT( overflow threshold )
    (18) Same as (16), but multiplied by SQRT( underflow threshold )

    (19) Nonsymmetric matrix with random entries chosen from (-1,1).
    (20) Same as (19), but multiplied by SQRT( overflow threshold )
    (21) Same as (19), but multiplied by SQRT( underflow threshold )
  NSIZES - INTEGER
           The number of sizes of matrices to use.  If it is zero,
           DCHKHS does nothing.  It must be at least zero.
           Not modified.

  NN     - 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.
           Not modified.

  NTYPES - INTEGER
           The number of elements in DOTYPE.   If it is zero, DCHKHS
           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. .
           Not modified.

  DOTYPE - 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.
           Not modified.

  ISEED  - 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 DCHKHS to continue the same random number
           sequence.
           Modified.

  THRESH - 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.
           Not modified.

  NOUNIT - INTEGER
           The FORTRAN unit number for printing out error messages
           (e.g., if a routine returns IINFO not equal to 0.)
           Not modified.

  A      - 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.
           Modified.

  LDA    - INTEGER
           The leading dimension of A, H, T1 and T2.  It must be at
           least 1 and at least max( NN ).
           Not modified.

  H      - DOUBLE PRECISION array, dimension (LDA,max(NN))
           The upper hessenberg matrix computed by DGEHRD.  On exit,
           H contains the Hessenberg form of the matrix in A.
           Modified.

  T1     - DOUBLE PRECISION array, dimension (LDA,max(NN))
           The Schur (="quasi-triangular") matrix computed by DHSEQR
           if Z is computed.  On exit, T1 contains the Schur form of
           the matrix in A.
           Modified.

  T2     - DOUBLE PRECISION array, dimension (LDA,max(NN))
           The Schur matrix computed by DHSEQR when Z is not computed.
           This should be identical to T1.
           Modified.

  LDU    - INTEGER
           The leading dimension of U, Z, UZ and UU.  It must be at
           least 1 and at least max( NN ).
           Not modified.

  U      - DOUBLE PRECISION array, dimension (LDU,max(NN))
           The orthogonal matrix computed by DGEHRD.
           Modified.

  Z      - DOUBLE PRECISION array, dimension (LDU,max(NN))
           The orthogonal matrix computed by DHSEQR.
           Modified.

  UZ     - DOUBLE PRECISION array, dimension (LDU,max(NN))
           The product of U times Z.
           Modified.

  WR1    - DOUBLE PRECISION array, dimension (max(NN))
  WI1    - DOUBLE PRECISION array, dimension (max(NN))
           The real and imaginary parts of the eigenvalues of A,
           as computed when Z is computed.
           On exit, WR1 + WI1*i are the eigenvalues of the matrix in A.
           Modified.

  WR2    - DOUBLE PRECISION array, dimension (max(NN))
  WI2    - DOUBLE PRECISION array, dimension (max(NN))
           The real and imaginary parts of the eigenvalues of A,
           as computed when T is computed but not Z.
           On exit, WR2 + WI2*i are the eigenvalues of the matrix in A.
           Modified.

  WR3    - DOUBLE PRECISION array, dimension (max(NN))
  WI3    - DOUBLE PRECISION array, dimension (max(NN))
           Like WR1, WI1, these arrays contain the eigenvalues of A,
           but those computed when DHSEQR only computes the
           eigenvalues, i.e., not the Schur vectors and no more of the
           Schur form than is necessary for computing the
           eigenvalues.
           Modified.

  EVECTL - DOUBLE PRECISION array, dimension (LDU,max(NN))
           The (upper triangular) left eigenvector matrix for the
           matrix in T1.  For complex conjugate pairs, the real part
           is stored in one row and the imaginary part in the next.
           Modified.

  EVEZTR - DOUBLE PRECISION array, dimension (LDU,max(NN))
           The (upper triangular) right eigenvector matrix for the
           matrix in T1.  For complex conjugate pairs, the real part
           is stored in one column and the imaginary part in the next.
           Modified.

  EVECTY - DOUBLE PRECISION array, dimension (LDU,max(NN))
           The left eigenvector matrix for the
           matrix in H.  For complex conjugate pairs, the real part
           is stored in one row and the imaginary part in the next.
           Modified.

  EVECTX - DOUBLE PRECISION array, dimension (LDU,max(NN))
           The right eigenvector matrix for the
           matrix in H.  For complex conjugate pairs, the real part
           is stored in one column and the imaginary part in the next.
           Modified.

  UU     - DOUBLE PRECISION array, dimension (LDU,max(NN))
           Details of the orthogonal matrix computed by DGEHRD.
           Modified.

  TAU    - DOUBLE PRECISION array, dimension(max(NN))
           Further details of the orthogonal matrix computed by DGEHRD.
           Modified.

  WORK   - DOUBLE PRECISION array, dimension (NWORK)
           Workspace.
           Modified.

  NWORK  - INTEGER
           The number of entries in WORK.  NWORK >= 4*NN(j)*NN(j) + 2.

  IWORK  - INTEGER array, dimension (max(NN))
           Workspace.
           Modified.

  SELECT - LOGICAL array, dimension (max(NN))
           Workspace.
           Modified.

  RESULT - DOUBLE PRECISION array, dimension (14)
           The values computed by the fourteen tests described above.
           The values are currently limited to 1/ulp, to avoid
           overflow.
           Modified.

  INFO   - 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) ).
           -14: LDU < 1 or LDU < NMAX.
           -28: NWORK too small.
           If  DLATMR, SLATMS, or SLATME returns an error code, the
               absolute value of it is returned.
           If 1, then DHSEQR could not find all the shifts.
           If 2, then the EISPACK code (for small blocks) failed.
           If >2, then 30*N iterations were not enough to find an
               eigenvalue or to decompose the problem.
           Modified.

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

     Some Local Variables and Parameters:
     ---- ----- --------- --- ----------

     ZERO, ONE       Real 0 and 1.
     MAXTYP          The number of types defined.
     MTEST           The number of tests defined: care must be taken
                     that (1) the size of RESULT, (2) the number of
                     tests actually performed, and (3) MTEST agree.
     NTEST           The number of tests performed on this matrix
                     so far.  This should be less than MTEST, and
                     equal to it by the last test.  It will be less
                     if any of the routines being tested indicates
                     that it could not compute the matrices that
                     would be tested.
     NMAX            Largest value in NN.
     NMATS           The number of matrices generated so far.
     NERRS           The number of tests which have exceeded THRESH
                     so far (computed by DLAFTS).
     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.
     RTOVFL, RTUNFL,
     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)       Selects 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 407 of file dchkhs.f.

412*
413* -- LAPACK test routine --
414* -- LAPACK is a software package provided by Univ. of Tennessee, --
415* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
416*
417* .. Scalar Arguments ..
418 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
419 DOUBLE PRECISION THRESH
420* ..
421* .. Array Arguments ..
422 LOGICAL DOTYPE( * ), SELECT( * )
423 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
424 DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ),
425 $ EVECTR( LDU, * ), EVECTX( LDU, * ),
426 $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ),
427 $ T1( LDA, * ), T2( LDA, * ), TAU( * ),
428 $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
429 $ WI1( * ), WI2( * ), WI3( * ), WORK( * ),
430 $ WR1( * ), WR2( * ), WR3( * ), Z( LDU, * )
431* ..
432*
433* =====================================================================
434*
435* .. Parameters ..
436 DOUBLE PRECISION ZERO, ONE
437 parameter( zero = 0.0d0, one = 1.0d0 )
438 INTEGER MAXTYP
439 parameter( maxtyp = 21 )
440* ..
441* .. Local Scalars ..
442 LOGICAL BADNN, MATCH
443 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL,
444 $ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS,
445 $ NMATS, NMAX, NSELC, NSELR, NTEST, NTESTT
446 DOUBLE PRECISION ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP,
447 $ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL
448* ..
449* .. Local Arrays ..
450 CHARACTER ADUMMA( 1 )
451 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
452 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
453 $ KTYPE( MAXTYP )
454 DOUBLE PRECISION DUMMA( 6 )
455* ..
456* .. External Functions ..
457 DOUBLE PRECISION DLAMCH
458 EXTERNAL dlamch
459* ..
460* .. External Subroutines ..
461 EXTERNAL dcopy, dgehrd, dgemm, dget10, dget22, dhsein,
464 $ dtrevc, xerbla
465* ..
466* .. Intrinsic Functions ..
467 INTRINSIC abs, dble, max, min, sqrt
468* ..
469* .. Data statements ..
470 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
471 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
472 $ 3, 1, 2, 3 /
473 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
474 $ 1, 5, 5, 5, 4, 3, 1 /
475 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
476* ..
477* .. Executable Statements ..
478*
479* Check for errors
480*
481 ntestt = 0
482 info = 0
483*
484 badnn = .false.
485 nmax = 0
486 DO 10 j = 1, nsizes
487 nmax = max( nmax, nn( j ) )
488 IF( nn( j ).LT.0 )
489 $ badnn = .true.
490 10 CONTINUE
491*
492* Check for errors
493*
494 IF( nsizes.LT.0 ) THEN
495 info = -1
496 ELSE IF( badnn ) THEN
497 info = -2
498 ELSE IF( ntypes.LT.0 ) THEN
499 info = -3
500 ELSE IF( thresh.LT.zero ) THEN
501 info = -6
502 ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
503 info = -9
504 ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax ) THEN
505 info = -14
506 ELSE IF( 4*nmax*nmax+2.GT.nwork ) THEN
507 info = -28
508 END IF
509*
510 IF( info.NE.0 ) THEN
511 CALL xerbla( 'DCHKHS', -info )
512 RETURN
513 END IF
514*
515* Quick return if possible
516*
517 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
518 $ RETURN
519*
520* More important constants
521*
522 unfl = dlamch( 'Safe minimum' )
523 ovfl = dlamch( 'Overflow' )
524 CALL dlabad( unfl, ovfl )
525 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
526 ulpinv = one / ulp
527 rtunfl = sqrt( unfl )
528 rtovfl = sqrt( ovfl )
529 rtulp = sqrt( ulp )
530 rtulpi = one / rtulp
531*
532* Loop over sizes, types
533*
534 nerrs = 0
535 nmats = 0
536*
537 DO 270 jsize = 1, nsizes
538 n = nn( jsize )
539 IF( n.EQ.0 )
540 $ GO TO 270
541 n1 = max( 1, n )
542 aninv = one / dble( n1 )
543*
544 IF( nsizes.NE.1 ) THEN
545 mtypes = min( maxtyp, ntypes )
546 ELSE
547 mtypes = min( maxtyp+1, ntypes )
548 END IF
549*
550 DO 260 jtype = 1, mtypes
551 IF( .NOT.dotype( jtype ) )
552 $ GO TO 260
553 nmats = nmats + 1
554 ntest = 0
555*
556* Save ISEED in case of an error.
557*
558 DO 20 j = 1, 4
559 ioldsd( j ) = iseed( j )
560 20 CONTINUE
561*
562* Initialize RESULT
563*
564 DO 30 j = 1, 14
565 result( j ) = zero
566 30 CONTINUE
567*
568* Compute "A"
569*
570* Control parameters:
571*
572* KMAGN KCONDS KMODE KTYPE
573* =1 O(1) 1 clustered 1 zero
574* =2 large large clustered 2 identity
575* =3 small exponential Jordan
576* =4 arithmetic diagonal, (w/ eigenvalues)
577* =5 random log symmetric, w/ eigenvalues
578* =6 random general, w/ eigenvalues
579* =7 random diagonal
580* =8 random symmetric
581* =9 random general
582* =10 random triangular
583*
584 IF( mtypes.GT.maxtyp )
585 $ GO TO 100
586*
587 itype = ktype( jtype )
588 imode = kmode( jtype )
589*
590* Compute norm
591*
592 GO TO ( 40, 50, 60 )kmagn( jtype )
593*
594 40 CONTINUE
595 anorm = one
596 GO TO 70
597*
598 50 CONTINUE
599 anorm = ( rtovfl*ulp )*aninv
600 GO TO 70
601*
602 60 CONTINUE
603 anorm = rtunfl*n*ulpinv
604 GO TO 70
605*
606 70 CONTINUE
607*
608 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
609 iinfo = 0
610 cond = ulpinv
611*
612* Special Matrices
613*
614 IF( itype.EQ.1 ) THEN
615*
616* Zero
617*
618 iinfo = 0
619*
620 ELSE IF( itype.EQ.2 ) THEN
621*
622* Identity
623*
624 DO 80 jcol = 1, n
625 a( jcol, jcol ) = anorm
626 80 CONTINUE
627*
628 ELSE IF( itype.EQ.3 ) THEN
629*
630* Jordan Block
631*
632 DO 90 jcol = 1, n
633 a( jcol, jcol ) = anorm
634 IF( jcol.GT.1 )
635 $ a( jcol, jcol-1 ) = one
636 90 CONTINUE
637*
638 ELSE IF( itype.EQ.4 ) THEN
639*
640* Diagonal Matrix, [Eigen]values Specified
641*
642 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
643 $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
644 $ iinfo )
645*
646 ELSE IF( itype.EQ.5 ) THEN
647*
648* Symmetric, eigenvalues specified
649*
650 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
651 $ anorm, n, n, 'N', a, lda, work( n+1 ),
652 $ iinfo )
653*
654 ELSE IF( itype.EQ.6 ) THEN
655*
656* General, eigenvalues specified
657*
658 IF( kconds( jtype ).EQ.1 ) THEN
659 conds = one
660 ELSE IF( kconds( jtype ).EQ.2 ) THEN
661 conds = rtulpi
662 ELSE
663 conds = zero
664 END IF
665*
666 adumma( 1 ) = ' '
667 CALL dlatme( n, 'S', iseed, work, imode, cond, one,
668 $ adumma, 'T', 'T', 'T', work( n+1 ), 4,
669 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
670 $ iinfo )
671*
672 ELSE IF( itype.EQ.7 ) THEN
673*
674* Diagonal, random eigenvalues
675*
676 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
677 $ 'T', 'N', work( n+1 ), 1, one,
678 $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
679 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
680*
681 ELSE IF( itype.EQ.8 ) THEN
682*
683* Symmetric, random eigenvalues
684*
685 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
686 $ 'T', 'N', work( n+1 ), 1, one,
687 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
688 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
689*
690 ELSE IF( itype.EQ.9 ) THEN
691*
692* General, random eigenvalues
693*
694 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
695 $ 'T', 'N', work( n+1 ), 1, one,
696 $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
697 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
698*
699 ELSE IF( itype.EQ.10 ) THEN
700*
701* Triangular, random eigenvalues
702*
703 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
704 $ 'T', 'N', work( n+1 ), 1, one,
705 $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
706 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
707*
708 ELSE
709*
710 iinfo = 1
711 END IF
712*
713 IF( iinfo.NE.0 ) THEN
714 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
715 $ ioldsd
716 info = abs( iinfo )
717 RETURN
718 END IF
719*
720 100 CONTINUE
721*
722* Call DGEHRD to compute H and U, do tests.
723*
724 CALL dlacpy( ' ', n, n, a, lda, h, lda )
725*
726 ntest = 1
727*
728 ilo = 1
729 ihi = n
730*
731 CALL dgehrd( n, ilo, ihi, h, lda, work, work( n+1 ),
732 $ nwork-n, iinfo )
733*
734 IF( iinfo.NE.0 ) THEN
735 result( 1 ) = ulpinv
736 WRITE( nounit, fmt = 9999 )'DGEHRD', iinfo, n, jtype,
737 $ ioldsd
738 info = abs( iinfo )
739 GO TO 250
740 END IF
741*
742 DO 120 j = 1, n - 1
743 uu( j+1, j ) = zero
744 DO 110 i = j + 2, n
745 u( i, j ) = h( i, j )
746 uu( i, j ) = h( i, j )
747 h( i, j ) = zero
748 110 CONTINUE
749 120 CONTINUE
750 CALL dcopy( n-1, work, 1, tau, 1 )
751 CALL dorghr( n, ilo, ihi, u, ldu, work, work( n+1 ),
752 $ nwork-n, iinfo )
753 ntest = 2
754*
755 CALL dhst01( n, ilo, ihi, a, lda, h, lda, u, ldu, work,
756 $ nwork, result( 1 ) )
757*
758* Call DHSEQR to compute T1, T2 and Z, do tests.
759*
760* Eigenvalues only (WR3,WI3)
761*
762 CALL dlacpy( ' ', n, n, h, lda, t2, lda )
763 ntest = 3
764 result( 3 ) = ulpinv
765*
766 CALL dhseqr( 'E', 'N', n, ilo, ihi, t2, lda, wr3, wi3, uz,
767 $ ldu, work, nwork, iinfo )
768 IF( iinfo.NE.0 ) THEN
769 WRITE( nounit, fmt = 9999 )'DHSEQR(E)', iinfo, n, jtype,
770 $ ioldsd
771 IF( iinfo.LE.n+2 ) THEN
772 info = abs( iinfo )
773 GO TO 250
774 END IF
775 END IF
776*
777* Eigenvalues (WR2,WI2) and Full Schur Form (T2)
778*
779 CALL dlacpy( ' ', n, n, h, lda, t2, lda )
780*
781 CALL dhseqr( 'S', 'N', n, ilo, ihi, t2, lda, wr2, wi2, uz,
782 $ ldu, work, nwork, iinfo )
783 IF( iinfo.NE.0 .AND. iinfo.LE.n+2 ) THEN
784 WRITE( nounit, fmt = 9999 )'DHSEQR(S)', iinfo, n, jtype,
785 $ ioldsd
786 info = abs( iinfo )
787 GO TO 250
788 END IF
789*
790* Eigenvalues (WR1,WI1), Schur Form (T1), and Schur vectors
791* (UZ)
792*
793 CALL dlacpy( ' ', n, n, h, lda, t1, lda )
794 CALL dlacpy( ' ', n, n, u, ldu, uz, ldu )
795*
796 CALL dhseqr( 'S', 'V', n, ilo, ihi, t1, lda, wr1, wi1, uz,
797 $ ldu, work, nwork, iinfo )
798 IF( iinfo.NE.0 .AND. iinfo.LE.n+2 ) THEN
799 WRITE( nounit, fmt = 9999 )'DHSEQR(V)', iinfo, n, jtype,
800 $ ioldsd
801 info = abs( iinfo )
802 GO TO 250
803 END IF
804*
805* Compute Z = U' UZ
806*
807 CALL dgemm( 'T', 'N', n, n, n, one, u, ldu, uz, ldu, zero,
808 $ z, ldu )
809 ntest = 8
810*
811* Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
812* and 4: | I - Z Z' | / ( n ulp )
813*
814 CALL dhst01( n, ilo, ihi, h, lda, t1, lda, z, ldu, work,
815 $ nwork, result( 3 ) )
816*
817* Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
818* and 6: | I - UZ (UZ)' | / ( n ulp )
819*
820 CALL dhst01( n, ilo, ihi, a, lda, t1, lda, uz, ldu, work,
821 $ nwork, result( 5 ) )
822*
823* Do Test 7: | T2 - T1 | / ( |T| n ulp )
824*
825 CALL dget10( n, n, t2, lda, t1, lda, work, result( 7 ) )
826*
827* Do Test 8: | W2 - W1 | / ( max(|W1|,|W2|) ulp )
828*
829 temp1 = zero
830 temp2 = zero
831 DO 130 j = 1, n
832 temp1 = max( temp1, abs( wr1( j ) )+abs( wi1( j ) ),
833 $ abs( wr2( j ) )+abs( wi2( j ) ) )
834 temp2 = max( temp2, abs( wr1( j )-wr2( j ) )+
835 & abs( wi1( j )-wi2( j ) ) )
836 130 CONTINUE
837*
838 result( 8 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
839*
840* Compute the Left and Right Eigenvectors of T
841*
842* Compute the Right eigenvector Matrix:
843*
844 ntest = 9
845 result( 9 ) = ulpinv
846*
847* Select last max(N/4,1) real, max(N/4,1) complex eigenvectors
848*
849 nselc = 0
850 nselr = 0
851 j = n
852 140 CONTINUE
853 IF( wi1( j ).EQ.zero ) THEN
854 IF( nselr.LT.max( n / 4, 1 ) ) THEN
855 nselr = nselr + 1
856 SELECT( j ) = .true.
857 ELSE
858 SELECT( j ) = .false.
859 END IF
860 j = j - 1
861 ELSE
862 IF( nselc.LT.max( n / 4, 1 ) ) THEN
863 nselc = nselc + 1
864 SELECT( j ) = .true.
865 SELECT( j-1 ) = .false.
866 ELSE
867 SELECT( j ) = .false.
868 SELECT( j-1 ) = .false.
869 END IF
870 j = j - 2
871 END IF
872 IF( j.GT.0 )
873 $ GO TO 140
874*
875 CALL dtrevc( 'Right', 'All', SELECT, n, t1, lda, dumma, ldu,
876 $ evectr, ldu, n, in, work, iinfo )
877 IF( iinfo.NE.0 ) THEN
878 WRITE( nounit, fmt = 9999 )'DTREVC(R,A)', iinfo, n,
879 $ jtype, ioldsd
880 info = abs( iinfo )
881 GO TO 250
882 END IF
883*
884* Test 9: | TR - RW | / ( |T| |R| ulp )
885*
886 CALL dget22( 'N', 'N', 'N', n, t1, lda, evectr, ldu, wr1,
887 $ wi1, work, dumma( 1 ) )
888 result( 9 ) = dumma( 1 )
889 IF( dumma( 2 ).GT.thresh ) THEN
890 WRITE( nounit, fmt = 9998 )'Right', 'DTREVC',
891 $ dumma( 2 ), n, jtype, ioldsd
892 END IF
893*
894* Compute selected right eigenvectors and confirm that
895* they agree with previous right eigenvectors
896*
897 CALL dtrevc( 'Right', 'Some', SELECT, n, t1, lda, dumma,
898 $ ldu, evectl, ldu, n, in, work, iinfo )
899 IF( iinfo.NE.0 ) THEN
900 WRITE( nounit, fmt = 9999 )'DTREVC(R,S)', iinfo, n,
901 $ jtype, ioldsd
902 info = abs( iinfo )
903 GO TO 250
904 END IF
905*
906 k = 1
907 match = .true.
908 DO 170 j = 1, n
909 IF( SELECT( j ) .AND. wi1( j ).EQ.zero ) THEN
910 DO 150 jj = 1, n
911 IF( evectr( jj, j ).NE.evectl( jj, k ) ) THEN
912 match = .false.
913 GO TO 180
914 END IF
915 150 CONTINUE
916 k = k + 1
917 ELSE IF( SELECT( j ) .AND. wi1( j ).NE.zero ) THEN
918 DO 160 jj = 1, n
919 IF( evectr( jj, j ).NE.evectl( jj, k ) .OR.
920 $ evectr( jj, j+1 ).NE.evectl( jj, k+1 ) ) THEN
921 match = .false.
922 GO TO 180
923 END IF
924 160 CONTINUE
925 k = k + 2
926 END IF
927 170 CONTINUE
928 180 CONTINUE
929 IF( .NOT.match )
930 $ WRITE( nounit, fmt = 9997 )'Right', 'DTREVC', n, jtype,
931 $ ioldsd
932*
933* Compute the Left eigenvector Matrix:
934*
935 ntest = 10
936 result( 10 ) = ulpinv
937 CALL dtrevc( 'Left', 'All', SELECT, n, t1, lda, evectl, ldu,
938 $ dumma, ldu, n, in, work, iinfo )
939 IF( iinfo.NE.0 ) THEN
940 WRITE( nounit, fmt = 9999 )'DTREVC(L,A)', iinfo, n,
941 $ jtype, ioldsd
942 info = abs( iinfo )
943 GO TO 250
944 END IF
945*
946* Test 10: | LT - WL | / ( |T| |L| ulp )
947*
948 CALL dget22( 'Trans', 'N', 'Conj', n, t1, lda, evectl, ldu,
949 $ wr1, wi1, work, dumma( 3 ) )
950 result( 10 ) = dumma( 3 )
951 IF( dumma( 4 ).GT.thresh ) THEN
952 WRITE( nounit, fmt = 9998 )'Left', 'DTREVC', dumma( 4 ),
953 $ n, jtype, ioldsd
954 END IF
955*
956* Compute selected left eigenvectors and confirm that
957* they agree with previous left eigenvectors
958*
959 CALL dtrevc( 'Left', 'Some', SELECT, n, t1, lda, evectr,
960 $ ldu, dumma, ldu, n, in, work, iinfo )
961 IF( iinfo.NE.0 ) THEN
962 WRITE( nounit, fmt = 9999 )'DTREVC(L,S)', iinfo, n,
963 $ jtype, ioldsd
964 info = abs( iinfo )
965 GO TO 250
966 END IF
967*
968 k = 1
969 match = .true.
970 DO 210 j = 1, n
971 IF( SELECT( j ) .AND. wi1( j ).EQ.zero ) THEN
972 DO 190 jj = 1, n
973 IF( evectl( jj, j ).NE.evectr( jj, k ) ) THEN
974 match = .false.
975 GO TO 220
976 END IF
977 190 CONTINUE
978 k = k + 1
979 ELSE IF( SELECT( j ) .AND. wi1( j ).NE.zero ) THEN
980 DO 200 jj = 1, n
981 IF( evectl( jj, j ).NE.evectr( jj, k ) .OR.
982 $ evectl( jj, j+1 ).NE.evectr( jj, k+1 ) ) THEN
983 match = .false.
984 GO TO 220
985 END IF
986 200 CONTINUE
987 k = k + 2
988 END IF
989 210 CONTINUE
990 220 CONTINUE
991 IF( .NOT.match )
992 $ WRITE( nounit, fmt = 9997 )'Left', 'DTREVC', n, jtype,
993 $ ioldsd
994*
995* Call DHSEIN for Right eigenvectors of H, do test 11
996*
997 ntest = 11
998 result( 11 ) = ulpinv
999 DO 230 j = 1, n
1000 SELECT( j ) = .true.
1001 230 CONTINUE
1002*
1003 CALL dhsein( 'Right', 'Qr', 'Ninitv', SELECT, n, h, lda,
1004 $ wr3, wi3, dumma, ldu, evectx, ldu, n1, in,
1005 $ work, iwork, iwork, iinfo )
1006 IF( iinfo.NE.0 ) THEN
1007 WRITE( nounit, fmt = 9999 )'DHSEIN(R)', iinfo, n, jtype,
1008 $ ioldsd
1009 info = abs( iinfo )
1010 IF( iinfo.LT.0 )
1011 $ GO TO 250
1012 ELSE
1013*
1014* Test 11: | HX - XW | / ( |H| |X| ulp )
1015*
1016* (from inverse iteration)
1017*
1018 CALL dget22( 'N', 'N', 'N', n, h, lda, evectx, ldu, wr3,
1019 $ wi3, work, dumma( 1 ) )
1020 IF( dumma( 1 ).LT.ulpinv )
1021 $ result( 11 ) = dumma( 1 )*aninv
1022 IF( dumma( 2 ).GT.thresh ) THEN
1023 WRITE( nounit, fmt = 9998 )'Right', 'DHSEIN',
1024 $ dumma( 2 ), n, jtype, ioldsd
1025 END IF
1026 END IF
1027*
1028* Call DHSEIN for Left eigenvectors of H, do test 12
1029*
1030 ntest = 12
1031 result( 12 ) = ulpinv
1032 DO 240 j = 1, n
1033 SELECT( j ) = .true.
1034 240 CONTINUE
1035*
1036 CALL dhsein( 'Left', 'Qr', 'Ninitv', SELECT, n, h, lda, wr3,
1037 $ wi3, evecty, ldu, dumma, ldu, n1, in, work,
1038 $ iwork, iwork, iinfo )
1039 IF( iinfo.NE.0 ) THEN
1040 WRITE( nounit, fmt = 9999 )'DHSEIN(L)', iinfo, n, jtype,
1041 $ ioldsd
1042 info = abs( iinfo )
1043 IF( iinfo.LT.0 )
1044 $ GO TO 250
1045 ELSE
1046*
1047* Test 12: | YH - WY | / ( |H| |Y| ulp )
1048*
1049* (from inverse iteration)
1050*
1051 CALL dget22( 'C', 'N', 'C', n, h, lda, evecty, ldu, wr3,
1052 $ wi3, work, dumma( 3 ) )
1053 IF( dumma( 3 ).LT.ulpinv )
1054 $ result( 12 ) = dumma( 3 )*aninv
1055 IF( dumma( 4 ).GT.thresh ) THEN
1056 WRITE( nounit, fmt = 9998 )'Left', 'DHSEIN',
1057 $ dumma( 4 ), n, jtype, ioldsd
1058 END IF
1059 END IF
1060*
1061* Call DORMHR for Right eigenvectors of A, do test 13
1062*
1063 ntest = 13
1064 result( 13 ) = ulpinv
1065*
1066 CALL dormhr( 'Left', 'No transpose', n, n, ilo, ihi, uu,
1067 $ ldu, tau, evectx, ldu, work, nwork, iinfo )
1068 IF( iinfo.NE.0 ) THEN
1069 WRITE( nounit, fmt = 9999 )'DORMHR(R)', iinfo, n, jtype,
1070 $ ioldsd
1071 info = abs( iinfo )
1072 IF( iinfo.LT.0 )
1073 $ GO TO 250
1074 ELSE
1075*
1076* Test 13: | AX - XW | / ( |A| |X| ulp )
1077*
1078* (from inverse iteration)
1079*
1080 CALL dget22( 'N', 'N', 'N', n, a, lda, evectx, ldu, wr3,
1081 $ wi3, work, dumma( 1 ) )
1082 IF( dumma( 1 ).LT.ulpinv )
1083 $ result( 13 ) = dumma( 1 )*aninv
1084 END IF
1085*
1086* Call DORMHR for Left eigenvectors of A, do test 14
1087*
1088 ntest = 14
1089 result( 14 ) = ulpinv
1090*
1091 CALL dormhr( 'Left', 'No transpose', n, n, ilo, ihi, uu,
1092 $ ldu, tau, evecty, ldu, work, nwork, iinfo )
1093 IF( iinfo.NE.0 ) THEN
1094 WRITE( nounit, fmt = 9999 )'DORMHR(L)', iinfo, n, jtype,
1095 $ ioldsd
1096 info = abs( iinfo )
1097 IF( iinfo.LT.0 )
1098 $ GO TO 250
1099 ELSE
1100*
1101* Test 14: | YA - WY | / ( |A| |Y| ulp )
1102*
1103* (from inverse iteration)
1104*
1105 CALL dget22( 'C', 'N', 'C', n, a, lda, evecty, ldu, wr3,
1106 $ wi3, work, dumma( 3 ) )
1107 IF( dumma( 3 ).LT.ulpinv )
1108 $ result( 14 ) = dumma( 3 )*aninv
1109 END IF
1110*
1111* End of Loop -- Check for RESULT(j) > THRESH
1112*
1113 250 CONTINUE
1114*
1115 ntestt = ntestt + ntest
1116 CALL dlafts( 'DHS', n, n, jtype, ntest, result, ioldsd,
1117 $ thresh, nounit, nerrs )
1118*
1119 260 CONTINUE
1120 270 CONTINUE
1121*
1122* Summary
1123*
1124 CALL dlasum( 'DHS', nounit, nerrs, ntestt )
1125*
1126 RETURN
1127*
1128 9999 FORMAT( ' DCHKHS: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
1129 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1130 9998 FORMAT( ' DCHKHS: ', a, ' Eigenvectors from ', a, ' incorrectly ',
1131 $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
1132 $ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5,
1133 $ ')' )
1134 9997 FORMAT( ' DCHKHS: Selected ', a, ' Eigenvectors from ', a,
1135 $ ' do not match other eigenvectors ', 9x, 'N=', i6,
1136 $ ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1137*
1138* End of DCHKHS
1139*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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
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 dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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 dget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, WI, WORK, RESULT)
DGET22
Definition: dget22.f:168
subroutine dlafts(TYPE, M, N, IMAT, NTESTS, RESULT, ISEED, THRESH, IOUNIT, IE)
DLAFTS
Definition: dlafts.f:99
subroutine dget10(M, N, A, LDA, B, LDB, WORK, RESULT)
DGET10
Definition: dget10.f:93
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
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:167
subroutine dhsein(SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI, VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL, IFAILR, INFO)
DHSEIN
Definition: dhsein.f:263
subroutine dtrevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTREVC
Definition: dtrevc.f:222
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:316
subroutine dormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMHR
Definition: dormhr.f:178
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: