LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dget24 ( logical  COMP,
integer  JTYPE,
double precision  THRESH,
integer, dimension( 4 )  ISEED,
integer  NOUNIT,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( lda, * )  H,
double precision, dimension( lda, * )  HT,
double precision, dimension( * )  WR,
double precision, dimension( * )  WI,
double precision, dimension( * )  WRT,
double precision, dimension( * )  WIT,
double precision, dimension( * )  WRTMP,
double precision, dimension( * )  WITMP,
double precision, dimension( ldvs, * )  VS,
integer  LDVS,
double precision, dimension( ldvs, * )  VS1,
double precision  RCDEIN,
double precision  RCDVIN,
integer  NSLCT,
integer, dimension( * )  ISLCT,
double precision, dimension( 17 )  RESULT,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

DGET24

Purpose:
    DGET24 checks the nonsymmetric eigenvalue (Schur form) problem
    expert driver DGEESX.

    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 DGEESX 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 DGEESX 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 DOUBLE PRECISION
          A test will count as "failed" if the "error", computed as
          described above, exceeds THRESH.  Note that the error
          is scaled to be O(1), so THRESH should be a reasonably
          small multiple of 1, e.g., 10 or 100.  In particular,
          it should not depend on the precision (single vs. double)
          or the size of the matrix.  It must be at least zero.
[in]NOUNIT
          NOUNIT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns INFO not equal to 0.)
[in]N
          N is INTEGER
          The dimension of A. N must be at least 0.
[in,out]A
          A is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDA, N)
          Another copy of the test matrix A, modified by DGEESX.
[out]HT
          HT is DOUBLE PRECISION array, dimension (LDA, N)
          Yet another copy of the test matrix A, modified by DGEESX.
[out]WR
          WR is DOUBLE PRECISION array, dimension (N)
[out]WI
          WI is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
[out]WIT
          WIT is DOUBLE PRECISION array, dimension (N)

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

          Like WR, WI, these arrays contain the eigenvalues of A,
          but sorted by increasing real part.
[out]VS
          VS is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDVS, N)
          VS1 holds another copy of the computed Schur vectors.
[in]RCDEIN
          RCDEIN is DOUBLE PRECISION
          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
          condition number for the average of selected eigenvalues.
[in]RCDVIN
          RCDVIN is DOUBLE PRECISION
          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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK to be passed to DGEESX. 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, DGEESX 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.
Date
November 2011

Definition at line 345 of file dget24.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: