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

◆ 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( 16 )  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.

            STREVC3 computes left and right eigenvector matrices
            from a Schur matrix T and backtransforms them with Z
            to eigenvector matrices L and R for A. L and R are
            GE matrices.

    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, 16
    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 )

    (15)    | AR - RW | / ( |A| |R| ulp )

    (16)    | LA - WL | / ( |A| |L| 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 (16)
           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 416 of file schkhs.f.

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