LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sget24()

subroutine sget24 ( logical  COMP,
integer  JTYPE,
real  THRESH,
integer, dimension( 4 )  ISEED,
integer  NOUNIT,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( lda, * )  H,
real, dimension( lda, * )  HT,
real, dimension( * )  WR,
real, dimension( * )  WI,
real, dimension( * )  WRT,
real, dimension( * )  WIT,
real, dimension( * )  WRTMP,
real, dimension( * )  WITMP,
real, dimension( ldvs, * )  VS,
integer  LDVS,
real, dimension( ldvs, * )  VS1,
real  RCDEIN,
real  RCDVIN,
integer  NSLCT,
integer, dimension( * )  ISLCT,
real, dimension( 17 )  RESULT,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

SGET24

Purpose:
    SGET24 checks the nonsymmetric eigenvalue (Schur form) problem
    expert driver SGEESX.

    If COMP = .FALSE., the first 13 of the following tests will be
    be performed on the input matrix A, and also tests 14 and 15
    if LWORK is sufficiently large.
    If COMP = .TRUE., all 17 test will be performed.

    (1)     0 if T is in Schur form, 1/ulp otherwise
           (no sorting of eigenvalues)

    (2)     | A - VS T VS' | / ( n |A| ulp )

      Here VS is the matrix of Schur eigenvectors, and T is in Schur
      form  (no sorting of eigenvalues).

    (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).

    (4)     0     if WR+sqrt(-1)*WI are eigenvalues of T
            1/ulp otherwise
            (no sorting of eigenvalues)

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

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

    (7)     0 if T is in Schur form, 1/ulp otherwise
            (with sorting of eigenvalues)

    (8)     | A - VS T VS' | / ( n |A| ulp )

      Here VS is the matrix of Schur eigenvectors, and T is in Schur
      form  (with sorting of eigenvalues).

    (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).

    (10)    0     if WR+sqrt(-1)*WI are eigenvalues of T
            1/ulp otherwise
            If workspace sufficient, also compare WR, WI with and
            without reciprocal condition numbers
            (with sorting of eigenvalues)

    (11)    0     if T(with VS) = T(without VS),
            1/ulp otherwise
            If workspace sufficient, also compare T with and without
            reciprocal condition numbers
            (with sorting of eigenvalues)

    (12)    0     if eigenvalues(with VS) = eigenvalues(without VS),
            1/ulp otherwise
            If workspace sufficient, also compare VS with and without
            reciprocal condition numbers
            (with sorting of eigenvalues)

    (13)    if sorting worked and SDIM is the number of
            eigenvalues which were SELECTed
            If workspace sufficient, also compare SDIM with and
            without reciprocal condition numbers

    (14)    if RCONDE the same no matter if VS and/or RCONDV computed

    (15)    if RCONDV the same no matter if VS and/or RCONDE computed

    (16)  |RCONDE - RCDEIN| / cond(RCONDE)

       RCONDE is the reciprocal average eigenvalue condition number
       computed by SGEESX and RCDEIN (the precomputed true value)
       is supplied as input.  cond(RCONDE) is the condition number
       of RCONDE, and takes errors in computing RCONDE into account,
       so that the resulting quantity should be O(ULP). cond(RCONDE)
       is essentially given by norm(A)/RCONDV.

    (17)  |RCONDV - RCDVIN| / cond(RCONDV)

       RCONDV is the reciprocal right invariant subspace condition
       number computed by SGEESX and RCDVIN (the precomputed true
       value) is supplied as input. cond(RCONDV) is the condition
       number of RCONDV, and takes errors in computing RCONDV into
       account, so that the resulting quantity should be O(ULP).
       cond(RCONDV) is essentially given by norm(A)/RCONDE.
Parameters
[in]COMP
          COMP is LOGICAL
          COMP describes which input tests to perform:
            = .FALSE. if the computed condition numbers are not to
                      be tested against RCDVIN and RCDEIN
            = .TRUE.  if they are to be compared
[in]JTYPE
          JTYPE is INTEGER
          Type of input matrix. Used to label output if error occurs.
[in]ISEED
          ISEED is INTEGER array, dimension (4)
          If COMP = .FALSE., the random number generator seed
          used to produce matrix.
          If COMP = .TRUE., ISEED(1) = the number of the example.
          Used to label output if error occurs.
[in]THRESH
          THRESH is 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.
[in]NOUNIT
          NOUNIT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns INFO not equal to 0.)
[in]N
          N is INTEGER
          The dimension of A. N must be at least 0.
[in,out]A
          A is REAL array, dimension (LDA, N)
          Used to hold the matrix whose eigenvalues are to be
          computed.
[in]LDA
          LDA is INTEGER
          The leading dimension of A, and H. LDA must be at
          least 1 and at least N.
[out]H
          H is REAL array, dimension (LDA, N)
          Another copy of the test matrix A, modified by SGEESX.
[out]HT
          HT is REAL array, dimension (LDA, N)
          Yet another copy of the test matrix A, modified by SGEESX.
[out]WR
          WR is REAL array, dimension (N)
[out]WI
          WI is REAL array, dimension (N)

          The real and imaginary parts of the eigenvalues of A.
          On exit, WR + WI*i are the eigenvalues of the matrix in A.
[out]WRT
          WRT is REAL array, dimension (N)
[out]WIT
          WIT is REAL array, dimension (N)

          Like WR, WI, these arrays contain the eigenvalues of A,
          but those computed when SGEESX only computes a partial
          eigendecomposition, i.e. not Schur vectors
[out]WRTMP
          WRTMP is REAL array, dimension (N)
[out]WITMP
          WITMP is REAL array, dimension (N)

          Like WR, WI, these arrays contain the eigenvalues of A,
          but sorted by increasing real part.
[out]VS
          VS is REAL array, dimension (LDVS, N)
          VS holds the computed Schur vectors.
[in]LDVS
          LDVS is INTEGER
          Leading dimension of VS. Must be at least max(1, N).
[out]VS1
          VS1 is REAL array, dimension (LDVS, N)
          VS1 holds another copy of the computed Schur vectors.
[in]RCDEIN
          RCDEIN is REAL
          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
          condition number for the average of selected eigenvalues.
[in]RCDVIN
          RCDVIN is REAL
          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
          condition number for the selected right invariant subspace.
[in]NSLCT
          NSLCT is INTEGER
          When COMP = .TRUE. the number of selected eigenvalues
          corresponding to the precomputed values RCDEIN and RCDVIN.
[in]ISLCT
          ISLCT is INTEGER array, dimension (NSLCT)
          When COMP = .TRUE. ISLCT selects the eigenvalues of the
          input matrix corresponding to the precomputed values RCDEIN
          and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
          eigenvalue with the J-th largest real part is selected.
          Not referenced if COMP = .FALSE.
[out]RESULT
          RESULT is REAL array, dimension (17)
          The values computed by the 17 tests described above.
          The values are currently limited to 1/ulp, to avoid
          overflow.
[out]WORK
          WORK is REAL array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK to be passed to SGEESX. This
          must be at least 3*N, and N+N**2 if tests 14--16 are to
          be performed.
[out]IWORK
          IWORK is INTEGER array, dimension (N*N)
[out]BWORK
          BWORK is LOGICAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          If 0,  successful exit.
          If <0, input parameter -INFO had an incorrect value.
          If >0, SGEESX returned an error code, the absolute
                 value of which is returned.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 339 of file sget24.f.

343 *
344 * -- LAPACK test routine --
345 * -- LAPACK is a software package provided by Univ. of Tennessee, --
346 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
347 *
348 * .. Scalar Arguments ..
349  LOGICAL COMP
350  INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
351  REAL RCDEIN, RCDVIN, THRESH
352 * ..
353 * .. Array Arguments ..
354  LOGICAL BWORK( * )
355  INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
356  REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ),
357  $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
358  $ WI( * ), WIT( * ), WITMP( * ), WORK( * ),
359  $ WR( * ), WRT( * ), WRTMP( * )
360 * ..
361 *
362 * =====================================================================
363 *
364 * .. Parameters ..
365  REAL ZERO, ONE
366  parameter( zero = 0.0e0, one = 1.0e0 )
367  REAL EPSIN
368  parameter( epsin = 5.9605e-8 )
369 * ..
370 * .. Local Scalars ..
371  CHARACTER SORT
372  INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
373  $ RSUB, SDIM, SDIM1
374  REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
375  $ SMLNUM, TMP, TOL, TOLIN, ULP, ULPINV, V, VIMIN,
376  $ VRMIN, WNORM
377 * ..
378 * .. Local Arrays ..
379  INTEGER IPNT( 20 )
380 * ..
381 * .. Arrays in Common ..
382  LOGICAL SELVAL( 20 )
383  REAL SELWI( 20 ), SELWR( 20 )
384 * ..
385 * .. Scalars in Common ..
386  INTEGER SELDIM, SELOPT
387 * ..
388 * .. Common blocks ..
389  COMMON / sslct / selopt, seldim, selval, selwr, selwi
390 * ..
391 * .. External Functions ..
392  LOGICAL SSLECT
393  REAL SLAMCH, SLANGE
394  EXTERNAL sslect, slamch, slange
395 * ..
396 * .. External Subroutines ..
397  EXTERNAL scopy, sgeesx, sgemm, slacpy, sort01, xerbla
398 * ..
399 * .. Intrinsic Functions ..
400  INTRINSIC abs, max, min, real, sign, sqrt
401 * ..
402 * .. Executable Statements ..
403 *
404 * Check for errors
405 *
406  info = 0
407  IF( thresh.LT.zero ) THEN
408  info = -3
409  ELSE IF( nounit.LE.0 ) THEN
410  info = -5
411  ELSE IF( n.LT.0 ) THEN
412  info = -6
413  ELSE IF( lda.LT.1 .OR. lda.LT.n ) THEN
414  info = -8
415  ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n ) THEN
416  info = -18
417  ELSE IF( lwork.LT.3*n ) THEN
418  info = -26
419  END IF
420 *
421  IF( info.NE.0 ) THEN
422  CALL xerbla( 'SGET24', -info )
423  RETURN
424  END IF
425 *
426 * Quick return if nothing to do
427 *
428  DO 10 i = 1, 17
429  result( i ) = -one
430  10 CONTINUE
431 *
432  IF( n.EQ.0 )
433  $ RETURN
434 *
435 * Important constants
436 *
437  smlnum = slamch( 'Safe minimum' )
438  ulp = slamch( 'Precision' )
439  ulpinv = one / ulp
440 *
441 * Perform tests (1)-(13)
442 *
443  selopt = 0
444  liwork = n*n
445  DO 120 isort = 0, 1
446  IF( isort.EQ.0 ) THEN
447  sort = 'N'
448  rsub = 0
449  ELSE
450  sort = 'S'
451  rsub = 6
452  END IF
453 *
454 * Compute Schur form and Schur vectors, and test them
455 *
456  CALL slacpy( 'F', n, n, a, lda, h, lda )
457  CALL sgeesx( 'V', sort, sslect, 'N', n, h, lda, sdim, wr, wi,
458  $ vs, ldvs, rconde, rcondv, work, lwork, iwork,
459  $ liwork, bwork, iinfo )
460  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
461  result( 1+rsub ) = ulpinv
462  IF( jtype.NE.22 ) THEN
463  WRITE( nounit, fmt = 9998 )'SGEESX1', iinfo, n, jtype,
464  $ iseed
465  ELSE
466  WRITE( nounit, fmt = 9999 )'SGEESX1', iinfo, n,
467  $ iseed( 1 )
468  END IF
469  info = abs( iinfo )
470  RETURN
471  END IF
472  IF( isort.EQ.0 ) THEN
473  CALL scopy( n, wr, 1, wrtmp, 1 )
474  CALL scopy( n, wi, 1, witmp, 1 )
475  END IF
476 *
477 * Do Test (1) or Test (7)
478 *
479  result( 1+rsub ) = zero
480  DO 30 j = 1, n - 2
481  DO 20 i = j + 2, n
482  IF( h( i, j ).NE.zero )
483  $ result( 1+rsub ) = ulpinv
484  20 CONTINUE
485  30 CONTINUE
486  DO 40 i = 1, n - 2
487  IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.zero )
488  $ result( 1+rsub ) = ulpinv
489  40 CONTINUE
490  DO 50 i = 1, n - 1
491  IF( h( i+1, i ).NE.zero ) THEN
492  IF( h( i, i ).NE.h( i+1, i+1 ) .OR. h( i, i+1 ).EQ.
493  $ zero .OR. sign( one, h( i+1, i ) ).EQ.
494  $ sign( one, h( i, i+1 ) ) )result( 1+rsub ) = ulpinv
495  END IF
496  50 CONTINUE
497 *
498 * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
499 *
500 * Copy A to VS1, used as workspace
501 *
502  CALL slacpy( ' ', n, n, a, lda, vs1, ldvs )
503 *
504 * Compute Q*H and store in HT.
505 *
506  CALL sgemm( 'No transpose', 'No transpose', n, n, n, one, vs,
507  $ ldvs, h, lda, zero, ht, lda )
508 *
509 * Compute A - Q*H*Q'
510 *
511  CALL sgemm( 'No transpose', 'Transpose', n, n, n, -one, ht,
512  $ lda, vs, ldvs, one, vs1, ldvs )
513 *
514  anorm = max( slange( '1', n, n, a, lda, work ), smlnum )
515  wnorm = slange( '1', n, n, vs1, ldvs, work )
516 *
517  IF( anorm.GT.wnorm ) THEN
518  result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
519  ELSE
520  IF( anorm.LT.one ) THEN
521  result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
522  $ ( n*ulp )
523  ELSE
524  result( 2+rsub ) = min( wnorm / anorm, real( n ) ) /
525  $ ( n*ulp )
526  END IF
527  END IF
528 *
529 * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
530 *
531  CALL sort01( 'Columns', n, n, vs, ldvs, work, lwork,
532  $ result( 3+rsub ) )
533 *
534 * Do Test (4) or Test (10)
535 *
536  result( 4+rsub ) = zero
537  DO 60 i = 1, n
538  IF( h( i, i ).NE.wr( i ) )
539  $ result( 4+rsub ) = ulpinv
540  60 CONTINUE
541  IF( n.GT.1 ) THEN
542  IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
543  $ result( 4+rsub ) = ulpinv
544  IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
545  $ result( 4+rsub ) = ulpinv
546  END IF
547  DO 70 i = 1, n - 1
548  IF( h( i+1, i ).NE.zero ) THEN
549  tmp = sqrt( abs( h( i+1, i ) ) )*
550  $ sqrt( abs( h( i, i+1 ) ) )
551  result( 4+rsub ) = max( result( 4+rsub ),
552  $ abs( wi( i )-tmp ) /
553  $ max( ulp*tmp, smlnum ) )
554  result( 4+rsub ) = max( result( 4+rsub ),
555  $ abs( wi( i+1 )+tmp ) /
556  $ max( ulp*tmp, smlnum ) )
557  ELSE IF( i.GT.1 ) THEN
558  IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.zero .AND.
559  $ wi( i ).NE.zero )result( 4+rsub ) = ulpinv
560  END IF
561  70 CONTINUE
562 *
563 * Do Test (5) or Test (11)
564 *
565  CALL slacpy( 'F', n, n, a, lda, ht, lda )
566  CALL sgeesx( 'N', sort, sslect, 'N', n, ht, lda, sdim, wrt,
567  $ wit, vs, ldvs, rconde, rcondv, work, lwork, iwork,
568  $ liwork, bwork, iinfo )
569  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
570  result( 5+rsub ) = ulpinv
571  IF( jtype.NE.22 ) THEN
572  WRITE( nounit, fmt = 9998 )'SGEESX2', iinfo, n, jtype,
573  $ iseed
574  ELSE
575  WRITE( nounit, fmt = 9999 )'SGEESX2', iinfo, n,
576  $ iseed( 1 )
577  END IF
578  info = abs( iinfo )
579  GO TO 250
580  END IF
581 *
582  result( 5+rsub ) = zero
583  DO 90 j = 1, n
584  DO 80 i = 1, n
585  IF( h( i, j ).NE.ht( i, j ) )
586  $ result( 5+rsub ) = ulpinv
587  80 CONTINUE
588  90 CONTINUE
589 *
590 * Do Test (6) or Test (12)
591 *
592  result( 6+rsub ) = zero
593  DO 100 i = 1, n
594  IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
595  $ result( 6+rsub ) = ulpinv
596  100 CONTINUE
597 *
598 * Do Test (13)
599 *
600  IF( isort.EQ.1 ) THEN
601  result( 13 ) = zero
602  knteig = 0
603  DO 110 i = 1, n
604  IF( sslect( wr( i ), wi( i ) ) .OR.
605  $ sslect( wr( i ), -wi( i ) ) )knteig = knteig + 1
606  IF( i.LT.n ) THEN
607  IF( ( sslect( wr( i+1 ), wi( i+1 ) ) .OR.
608  $ sslect( wr( i+1 ), -wi( i+1 ) ) ) .AND.
609  $ ( .NOT.( sslect( wr( i ),
610  $ wi( i ) ) .OR. sslect( wr( i ),
611  $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )result( 13 )
612  $ = ulpinv
613  END IF
614  110 CONTINUE
615  IF( sdim.NE.knteig )
616  $ result( 13 ) = ulpinv
617  END IF
618 *
619  120 CONTINUE
620 *
621 * If there is enough workspace, perform tests (14) and (15)
622 * as well as (10) through (13)
623 *
624  IF( lwork.GE.n+( n*n ) / 2 ) THEN
625 *
626 * Compute both RCONDE and RCONDV with VS
627 *
628  sort = 'S'
629  result( 14 ) = zero
630  result( 15 ) = zero
631  CALL slacpy( 'F', n, n, a, lda, ht, lda )
632  CALL sgeesx( 'V', sort, sslect, 'B', n, ht, lda, sdim1, wrt,
633  $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
634  $ iwork, liwork, bwork, iinfo )
635  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
636  result( 14 ) = ulpinv
637  result( 15 ) = ulpinv
638  IF( jtype.NE.22 ) THEN
639  WRITE( nounit, fmt = 9998 )'SGEESX3', iinfo, n, jtype,
640  $ iseed
641  ELSE
642  WRITE( nounit, fmt = 9999 )'SGEESX3', iinfo, n,
643  $ iseed( 1 )
644  END IF
645  info = abs( iinfo )
646  GO TO 250
647  END IF
648 *
649 * Perform tests (10), (11), (12), and (13)
650 *
651  DO 140 i = 1, n
652  IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
653  $ result( 10 ) = ulpinv
654  DO 130 j = 1, n
655  IF( h( i, j ).NE.ht( i, j ) )
656  $ result( 11 ) = ulpinv
657  IF( vs( i, j ).NE.vs1( i, j ) )
658  $ result( 12 ) = ulpinv
659  130 CONTINUE
660  140 CONTINUE
661  IF( sdim.NE.sdim1 )
662  $ result( 13 ) = ulpinv
663 *
664 * Compute both RCONDE and RCONDV without VS, and compare
665 *
666  CALL slacpy( 'F', n, n, a, lda, ht, lda )
667  CALL sgeesx( 'N', sort, sslect, 'B', n, ht, lda, sdim1, wrt,
668  $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
669  $ iwork, liwork, bwork, iinfo )
670  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
671  result( 14 ) = ulpinv
672  result( 15 ) = ulpinv
673  IF( jtype.NE.22 ) THEN
674  WRITE( nounit, fmt = 9998 )'SGEESX4', iinfo, n, jtype,
675  $ iseed
676  ELSE
677  WRITE( nounit, fmt = 9999 )'SGEESX4', iinfo, n,
678  $ iseed( 1 )
679  END IF
680  info = abs( iinfo )
681  GO TO 250
682  END IF
683 *
684 * Perform tests (14) and (15)
685 *
686  IF( rcnde1.NE.rconde )
687  $ result( 14 ) = ulpinv
688  IF( rcndv1.NE.rcondv )
689  $ result( 15 ) = ulpinv
690 *
691 * Perform tests (10), (11), (12), and (13)
692 *
693  DO 160 i = 1, n
694  IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
695  $ result( 10 ) = ulpinv
696  DO 150 j = 1, n
697  IF( h( i, j ).NE.ht( i, j ) )
698  $ result( 11 ) = ulpinv
699  IF( vs( i, j ).NE.vs1( i, j ) )
700  $ result( 12 ) = ulpinv
701  150 CONTINUE
702  160 CONTINUE
703  IF( sdim.NE.sdim1 )
704  $ result( 13 ) = ulpinv
705 *
706 * Compute RCONDE with VS, and compare
707 *
708  CALL slacpy( 'F', n, n, a, lda, ht, lda )
709  CALL sgeesx( 'V', sort, sslect, 'E', n, ht, lda, sdim1, wrt,
710  $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
711  $ iwork, liwork, bwork, iinfo )
712  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
713  result( 14 ) = ulpinv
714  IF( jtype.NE.22 ) THEN
715  WRITE( nounit, fmt = 9998 )'SGEESX5', iinfo, n, jtype,
716  $ iseed
717  ELSE
718  WRITE( nounit, fmt = 9999 )'SGEESX5', iinfo, n,
719  $ iseed( 1 )
720  END IF
721  info = abs( iinfo )
722  GO TO 250
723  END IF
724 *
725 * Perform test (14)
726 *
727  IF( rcnde1.NE.rconde )
728  $ result( 14 ) = ulpinv
729 *
730 * Perform tests (10), (11), (12), and (13)
731 *
732  DO 180 i = 1, n
733  IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
734  $ result( 10 ) = ulpinv
735  DO 170 j = 1, n
736  IF( h( i, j ).NE.ht( i, j ) )
737  $ result( 11 ) = ulpinv
738  IF( vs( i, j ).NE.vs1( i, j ) )
739  $ result( 12 ) = ulpinv
740  170 CONTINUE
741  180 CONTINUE
742  IF( sdim.NE.sdim1 )
743  $ result( 13 ) = ulpinv
744 *
745 * Compute RCONDE without VS, and compare
746 *
747  CALL slacpy( 'F', n, n, a, lda, ht, lda )
748  CALL sgeesx( 'N', sort, sslect, 'E', n, ht, lda, sdim1, wrt,
749  $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
750  $ iwork, liwork, bwork, iinfo )
751  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
752  result( 14 ) = ulpinv
753  IF( jtype.NE.22 ) THEN
754  WRITE( nounit, fmt = 9998 )'SGEESX6', iinfo, n, jtype,
755  $ iseed
756  ELSE
757  WRITE( nounit, fmt = 9999 )'SGEESX6', iinfo, n,
758  $ iseed( 1 )
759  END IF
760  info = abs( iinfo )
761  GO TO 250
762  END IF
763 *
764 * Perform test (14)
765 *
766  IF( rcnde1.NE.rconde )
767  $ result( 14 ) = ulpinv
768 *
769 * Perform tests (10), (11), (12), and (13)
770 *
771  DO 200 i = 1, n
772  IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
773  $ result( 10 ) = ulpinv
774  DO 190 j = 1, n
775  IF( h( i, j ).NE.ht( i, j ) )
776  $ result( 11 ) = ulpinv
777  IF( vs( i, j ).NE.vs1( i, j ) )
778  $ result( 12 ) = ulpinv
779  190 CONTINUE
780  200 CONTINUE
781  IF( sdim.NE.sdim1 )
782  $ result( 13 ) = ulpinv
783 *
784 * Compute RCONDV with VS, and compare
785 *
786  CALL slacpy( 'F', n, n, a, lda, ht, lda )
787  CALL sgeesx( 'V', sort, sslect, 'V', n, ht, lda, sdim1, wrt,
788  $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
789  $ iwork, liwork, bwork, iinfo )
790  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
791  result( 15 ) = ulpinv
792  IF( jtype.NE.22 ) THEN
793  WRITE( nounit, fmt = 9998 )'SGEESX7', iinfo, n, jtype,
794  $ iseed
795  ELSE
796  WRITE( nounit, fmt = 9999 )'SGEESX7', iinfo, n,
797  $ iseed( 1 )
798  END IF
799  info = abs( iinfo )
800  GO TO 250
801  END IF
802 *
803 * Perform test (15)
804 *
805  IF( rcndv1.NE.rcondv )
806  $ result( 15 ) = ulpinv
807 *
808 * Perform tests (10), (11), (12), and (13)
809 *
810  DO 220 i = 1, n
811  IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
812  $ result( 10 ) = ulpinv
813  DO 210 j = 1, n
814  IF( h( i, j ).NE.ht( i, j ) )
815  $ result( 11 ) = ulpinv
816  IF( vs( i, j ).NE.vs1( i, j ) )
817  $ result( 12 ) = ulpinv
818  210 CONTINUE
819  220 CONTINUE
820  IF( sdim.NE.sdim1 )
821  $ result( 13 ) = ulpinv
822 *
823 * Compute RCONDV without VS, and compare
824 *
825  CALL slacpy( 'F', n, n, a, lda, ht, lda )
826  CALL sgeesx( 'N', sort, sslect, 'V', n, ht, lda, sdim1, wrt,
827  $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
828  $ iwork, liwork, bwork, iinfo )
829  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
830  result( 15 ) = ulpinv
831  IF( jtype.NE.22 ) THEN
832  WRITE( nounit, fmt = 9998 )'SGEESX8', iinfo, n, jtype,
833  $ iseed
834  ELSE
835  WRITE( nounit, fmt = 9999 )'SGEESX8', iinfo, n,
836  $ iseed( 1 )
837  END IF
838  info = abs( iinfo )
839  GO TO 250
840  END IF
841 *
842 * Perform test (15)
843 *
844  IF( rcndv1.NE.rcondv )
845  $ result( 15 ) = ulpinv
846 *
847 * Perform tests (10), (11), (12), and (13)
848 *
849  DO 240 i = 1, n
850  IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
851  $ result( 10 ) = ulpinv
852  DO 230 j = 1, n
853  IF( h( i, j ).NE.ht( i, j ) )
854  $ result( 11 ) = ulpinv
855  IF( vs( i, j ).NE.vs1( i, j ) )
856  $ result( 12 ) = ulpinv
857  230 CONTINUE
858  240 CONTINUE
859  IF( sdim.NE.sdim1 )
860  $ result( 13 ) = ulpinv
861 *
862  END IF
863 *
864  250 CONTINUE
865 *
866 * If there are precomputed reciprocal condition numbers, compare
867 * computed values with them.
868 *
869  IF( comp ) THEN
870 *
871 * First set up SELOPT, SELDIM, SELVAL, SELWR, and SELWI so that
872 * the logical function SSLECT selects the eigenvalues specified
873 * by NSLCT and ISLCT.
874 *
875  seldim = n
876  selopt = 1
877  eps = max( ulp, epsin )
878  DO 260 i = 1, n
879  ipnt( i ) = i
880  selval( i ) = .false.
881  selwr( i ) = wrtmp( i )
882  selwi( i ) = witmp( i )
883  260 CONTINUE
884  DO 280 i = 1, n - 1
885  kmin = i
886  vrmin = wrtmp( i )
887  vimin = witmp( i )
888  DO 270 j = i + 1, n
889  IF( wrtmp( j ).LT.vrmin ) THEN
890  kmin = j
891  vrmin = wrtmp( j )
892  vimin = witmp( j )
893  END IF
894  270 CONTINUE
895  wrtmp( kmin ) = wrtmp( i )
896  witmp( kmin ) = witmp( i )
897  wrtmp( i ) = vrmin
898  witmp( i ) = vimin
899  itmp = ipnt( i )
900  ipnt( i ) = ipnt( kmin )
901  ipnt( kmin ) = itmp
902  280 CONTINUE
903  DO 290 i = 1, nslct
904  selval( ipnt( islct( i ) ) ) = .true.
905  290 CONTINUE
906 *
907 * Compute condition numbers
908 *
909  CALL slacpy( 'F', n, n, a, lda, ht, lda )
910  CALL sgeesx( 'N', 'S', sslect, 'B', n, ht, lda, sdim1, wrt,
911  $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
912  $ iwork, liwork, bwork, iinfo )
913  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
914  result( 16 ) = ulpinv
915  result( 17 ) = ulpinv
916  WRITE( nounit, fmt = 9999 )'SGEESX9', iinfo, n, iseed( 1 )
917  info = abs( iinfo )
918  GO TO 300
919  END IF
920 *
921 * Compare condition number for average of selected eigenvalues
922 * taking its condition number into account
923 *
924  anorm = slange( '1', n, n, a, lda, work )
925  v = max( real( n )*eps*anorm, smlnum )
926  IF( anorm.EQ.zero )
927  $ v = one
928  IF( v.GT.rcondv ) THEN
929  tol = one
930  ELSE
931  tol = v / rcondv
932  END IF
933  IF( v.GT.rcdvin ) THEN
934  tolin = one
935  ELSE
936  tolin = v / rcdvin
937  END IF
938  tol = max( tol, smlnum / eps )
939  tolin = max( tolin, smlnum / eps )
940  IF( eps*( rcdein-tolin ).GT.rconde+tol ) THEN
941  result( 16 ) = ulpinv
942  ELSE IF( rcdein-tolin.GT.rconde+tol ) THEN
943  result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
944  ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) ) THEN
945  result( 16 ) = ulpinv
946  ELSE IF( rcdein+tolin.LT.rconde-tol ) THEN
947  result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
948  ELSE
949  result( 16 ) = one
950  END IF
951 *
952 * Compare condition numbers for right invariant subspace
953 * taking its condition number into account
954 *
955  IF( v.GT.rcondv*rconde ) THEN
956  tol = rcondv
957  ELSE
958  tol = v / rconde
959  END IF
960  IF( v.GT.rcdvin*rcdein ) THEN
961  tolin = rcdvin
962  ELSE
963  tolin = v / rcdein
964  END IF
965  tol = max( tol, smlnum / eps )
966  tolin = max( tolin, smlnum / eps )
967  IF( eps*( rcdvin-tolin ).GT.rcondv+tol ) THEN
968  result( 17 ) = ulpinv
969  ELSE IF( rcdvin-tolin.GT.rcondv+tol ) THEN
970  result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
971  ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) ) THEN
972  result( 17 ) = ulpinv
973  ELSE IF( rcdvin+tolin.LT.rcondv-tol ) THEN
974  result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
975  ELSE
976  result( 17 ) = one
977  END IF
978 *
979  300 CONTINUE
980 *
981  END IF
982 *
983  9999 FORMAT( ' SGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
984  $ i6, ', INPUT EXAMPLE NUMBER = ', i4 )
985  9998 FORMAT( ' SGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
986  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
987 *
988  RETURN
989 *
990 * End of SGET24
991 *
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
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
subroutine sgeesx(JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
SGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: sgeesx.f:281
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 sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
Definition: sort01.f:116
logical function sslect(ZR, ZI)
SSLECT
Definition: sslect.f:62
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: