LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ schkhs()

subroutine schkhs ( integer  NSIZES,
integer, dimension( * )  NN,
integer  NTYPES,
logical, dimension( * )  DOTYPE,
integer, dimension( 4 )  ISEED,
real  THRESH,
integer  NOUNIT,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( lda, * )  H,
real, dimension( lda, * )  T1,
real, dimension( lda, * )  T2,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( ldu, * )  Z,
real, dimension( ldu, * )  UZ,
real, dimension( * )  WR1,
real, dimension( * )  WI1,
real, dimension( * )  WR2,
real, dimension( * )  WI2,
real, dimension( * )  WR3,
real, dimension( * )  WI3,
real, dimension( ldu, * )  EVECTL,
real, dimension( ldu, * )  EVECTR,
real, dimension( ldu, * )  EVECTY,
real, dimension( ldu, * )  EVECTX,
real, dimension( ldu, * )  UU,
real, dimension( * )  TAU,
real, dimension( * )  WORK,
integer  NWORK,
integer, dimension( * )  IWORK,
logical, dimension( * )  SELECT,
real, dimension( 14 )  RESULT,
integer  INFO 
)

SCHKHS

Purpose:
    SCHKHS  checks the nonsymmetric eigenvalue problem routines.

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

            SORGHR generates the orthogonal matrix U.

            SORMHR multiplies a matrix by the orthogonal matrix U.

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

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

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

    When SCHKHS 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,
           SCHKHS 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, SCHKHS
           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 SCHKHS to continue the same random number
           sequence.
           Modified.

  THRESH - REAL
           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      - REAL 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      - REAL array, dimension (LDA,max(NN))
           The upper hessenberg matrix computed by SGEHRD.  On exit,
           H contains the Hessenberg form of the matrix in A.
           Modified.

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

  T2     - REAL array, dimension (LDA,max(NN))
           The Schur matrix computed by SHSEQR 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      - REAL array, dimension (LDU,max(NN))
           The orthogonal matrix computed by SGEHRD.
           Modified.

  Z      - REAL array, dimension (LDU,max(NN))
           The orthogonal matrix computed by SHSEQR.
           Modified.

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

  WR1    - REAL array, dimension (max(NN))
  WI1    - REAL 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    - REAL array, dimension (max(NN))
  WI2    - REAL 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    - REAL array, dimension (max(NN))
  WI3    - REAL array, dimension (max(NN))
           Like WR1, WI1, these arrays contain the eigenvalues of A,
           but those computed when SHSEQR 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 - REAL 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.

  EVECTR - REAL 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 - REAL 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 - REAL 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     - REAL array, dimension (LDU,max(NN))
           Details of the orthogonal matrix computed by SGEHRD.
           Modified.

  TAU    - REAL array, dimension(max(NN))
           Further details of the orthogonal matrix computed by SGEHRD.
           Modified.

  WORK   - REAL 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 - REAL 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  SLATMR, SLATMS, or SLATME returns an error code, the
               absolute value of it is returned.
           If 1, then SHSEQR 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 SLAFTS).
     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 schkhs.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  REAL THRESH
420 * ..
421 * .. Array Arguments ..
422  LOGICAL DOTYPE( * ), SELECT( * )
423  INTEGER ISEED( 4 ), IWORK( * ), NN( * )
424  REAL 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  REAL ZERO, ONE
437  parameter( zero = 0.0, one = 1.0 )
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  REAL 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  REAL DUMMA( 6 )
455 * ..
456 * .. External Functions ..
457  REAL SLAMCH
458  EXTERNAL slamch
459 * ..
460 * .. External Subroutines ..
461  EXTERNAL scopy, sgehrd, sgemm, sget10, sget22, shsein,
464  $ strevc, xerbla
465 * ..
466 * .. Intrinsic Functions ..
467  INTRINSIC abs, max, min, real, 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( 'SCHKHS', -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 = slamch( 'Safe minimum' )
523  ovfl = slamch( 'Overflow' )
524  CALL slabad( unfl, ovfl )
525  ulp = slamch( 'Epsilon' )*slamch( '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 / real( 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 slaset( '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 slatms( 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 slatms( 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 slatme( 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 slatmr( 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 slatmr( 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 slatmr( 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 slatmr( 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 SGEHRD to compute H and U, do tests.
723 *
724  CALL slacpy( ' ', n, n, a, lda, h, lda )
725 *
726  ntest = 1
727 *
728  ilo = 1
729  ihi = n
730 *
731  CALL sgehrd( 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 )'SGEHRD', 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 scopy( n-1, work, 1, tau, 1 )
751  CALL sorghr( n, ilo, ihi, u, ldu, work, work( n+1 ),
752  $ nwork-n, iinfo )
753  ntest = 2
754 *
755  CALL shst01( n, ilo, ihi, a, lda, h, lda, u, ldu, work,
756  $ nwork, result( 1 ) )
757 *
758 * Call SHSEQR to compute T1, T2 and Z, do tests.
759 *
760 * Eigenvalues only (WR3,WI3)
761 *
762  CALL slacpy( ' ', n, n, h, lda, t2, lda )
763  ntest = 3
764  result( 3 ) = ulpinv
765 *
766  CALL shseqr( '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 )'SHSEQR(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 slacpy( ' ', n, n, h, lda, t2, lda )
780 *
781  CALL shseqr( '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 )'SHSEQR(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 slacpy( ' ', n, n, h, lda, t1, lda )
794  CALL slacpy( ' ', n, n, u, ldu, uz, ldu )
795 *
796  CALL shseqr( '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 )'SHSEQR(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 sgemm( '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 shst01( 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 shst01( 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 sget10( 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 strevc( '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 )'STREVC(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 sget22( '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', 'STREVC',
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 strevc( '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 )'STREVC(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', 'STREVC', n, jtype,
931  $ ioldsd
932 *
933 * Compute the Left eigenvector Matrix:
934 *
935  ntest = 10
936  result( 10 ) = ulpinv
937  CALL strevc( '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 )'STREVC(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 sget22( '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', 'STREVC', 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 strevc( '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 )'STREVC(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', 'STREVC', n, jtype,
993  $ ioldsd
994 *
995 * Call SHSEIN 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 shsein( '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 )'SHSEIN(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 sget22( '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', 'SHSEIN',
1024  $ dumma( 2 ), n, jtype, ioldsd
1025  END IF
1026  END IF
1027 *
1028 * Call SHSEIN 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 shsein( '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 )'SHSEIN(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 sget22( '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', 'SHSEIN',
1057  $ dumma( 4 ), n, jtype, ioldsd
1058  END IF
1059  END IF
1060 *
1061 * Call SORMHR for Right eigenvectors of A, do test 13
1062 *
1063  ntest = 13
1064  result( 13 ) = ulpinv
1065 *
1066  CALL sormhr( '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 )'SORMHR(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 sget22( '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 SORMHR for Left eigenvectors of A, do test 14
1087 *
1088  ntest = 14
1089  result( 14 ) = ulpinv
1090 *
1091  CALL sormhr( '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 )'SORMHR(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 sget22( '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 slafts( 'SHS', n, n, jtype, ntest, result, ioldsd,
1117  $ thresh, nounit, nerrs )
1118 *
1119  260 CONTINUE
1120  270 CONTINUE
1121 *
1122 * Summary
1123 *
1124  CALL slasum( 'SHS', nounit, nerrs, ntestt )
1125 *
1126  RETURN
1127 *
1128  9999 FORMAT( ' SCHKHS: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
1129  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1130  9998 FORMAT( ' SCHKHS: ', 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( ' SCHKHS: 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 SCHKHS
1139 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:321
subroutine slatme(N, DIST, ISEED, D, MODE, COND, DMAX, EI, RSIGN, UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, A, LDA, WORK, INFO)
SLATME
Definition: slatme.f:332
subroutine slatmr(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)
SLATMR
Definition: slatmr.f:471
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:167
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:126
subroutine sormhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMHR
Definition: sormhr.f:179
subroutine strevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STREVC
Definition: strevc.f:222
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:316
subroutine shsein(SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI, VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL, IFAILR, INFO)
SHSEIN
Definition: shsein.f:263
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
subroutine sget10(M, N, A, LDA, B, LDB, WORK, RESULT)
SGET10
Definition: sget10.f:93
subroutine sget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, WI, WORK, RESULT)
SGET22
Definition: sget22.f:168
subroutine slafts(TYPE, M, N, IMAT, NTESTS, RESULT, ISEED, THRESH, IOUNIT, IE)
SLAFTS
Definition: slafts.f:99
subroutine shst01(N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK, LWORK, RESULT)
SHST01
Definition: shst01.f:134
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
subroutine slasum(TYPE, IOUNIT, IE, NRUN)
SLASUM
Definition: slasum.f:41
Here is the call graph for this function:
Here is the caller graph for this function: