 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ dget24()

 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.

Definition at line 339 of file dget24.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  DOUBLE PRECISION RCDEIN, RCDVIN, THRESH
352 * ..
353 * .. Array Arguments ..
354  LOGICAL BWORK( * )
355  INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
356  DOUBLE PRECISION 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  DOUBLE PRECISION ZERO, ONE
366  parameter( zero = 0.0d0, one = 1.0d0 )
367  DOUBLE PRECISION EPSIN
368  parameter( epsin = 5.9605d-8 )
369 * ..
370 * .. Local Scalars ..
371  CHARACTER SORT
372  INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
373  \$ RSUB, SDIM, SDIM1
374  DOUBLE PRECISION 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  DOUBLE PRECISION 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 DSLECT
393  DOUBLE PRECISION DLAMCH, DLANGE
394  EXTERNAL dslect, dlamch, dlange
395 * ..
396 * .. External Subroutines ..
397  EXTERNAL dcopy, dgeesx, dgemm, dlacpy, dort01, xerbla
398 * ..
399 * .. Intrinsic Functions ..
400  INTRINSIC abs, dble, max, min, 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( 'DGET24', -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 = dlamch( 'Safe minimum' )
438  ulp = dlamch( '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 dlacpy( 'F', n, n, a, lda, h, lda )
457  CALL dgeesx( 'V', sort, dslect, '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 )'DGEESX1', iinfo, n, jtype,
464  \$ iseed
465  ELSE
466  WRITE( nounit, fmt = 9999 )'DGEESX1', 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 dcopy( n, wr, 1, wrtmp, 1 )
474  CALL dcopy( 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 dlacpy( ' ', n, n, a, lda, vs1, ldvs )
503 *
504 * Compute Q*H and store in HT.
505 *
506  CALL dgemm( '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 dgemm( 'No transpose', 'Transpose', n, n, n, -one, ht,
512  \$ lda, vs, ldvs, one, vs1, ldvs )
513 *
514  anorm = max( dlange( '1', n, n, a, lda, work ), smlnum )
515  wnorm = dlange( '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, dble( 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 dort01( '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 dlacpy( 'F', n, n, a, lda, ht, lda )
566  CALL dgeesx( 'N', sort, dslect, '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 )'DGEESX2', iinfo, n, jtype,
573  \$ iseed
574  ELSE
575  WRITE( nounit, fmt = 9999 )'DGEESX2', 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( dslect( wr( i ), wi( i ) ) .OR.
605  \$ dslect( wr( i ), -wi( i ) ) )knteig = knteig + 1
606  IF( i.LT.n ) THEN
607  IF( ( dslect( wr( i+1 ), wi( i+1 ) ) .OR.
608  \$ dslect( wr( i+1 ), -wi( i+1 ) ) ) .AND.
609  \$ ( .NOT.( dslect( wr( i ),
610  \$ wi( i ) ) .OR. dslect( 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 dlacpy( 'F', n, n, a, lda, ht, lda )
632  CALL dgeesx( 'V', sort, dslect, '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 )'DGEESX3', iinfo, n, jtype,
640  \$ iseed
641  ELSE
642  WRITE( nounit, fmt = 9999 )'DGEESX3', 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 dlacpy( 'F', n, n, a, lda, ht, lda )
667  CALL dgeesx( 'N', sort, dslect, '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 )'DGEESX4', iinfo, n, jtype,
675  \$ iseed
676  ELSE
677  WRITE( nounit, fmt = 9999 )'DGEESX4', 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 dlacpy( 'F', n, n, a, lda, ht, lda )
709  CALL dgeesx( 'V', sort, dslect, '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 )'DGEESX5', iinfo, n, jtype,
716  \$ iseed
717  ELSE
718  WRITE( nounit, fmt = 9999 )'DGEESX5', 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 dlacpy( 'F', n, n, a, lda, ht, lda )
748  CALL dgeesx( 'N', sort, dslect, '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 )'DGEESX6', iinfo, n, jtype,
755  \$ iseed
756  ELSE
757  WRITE( nounit, fmt = 9999 )'DGEESX6', 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 dlacpy( 'F', n, n, a, lda, ht, lda )
787  CALL dgeesx( 'V', sort, dslect, '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 )'DGEESX7', iinfo, n, jtype,
794  \$ iseed
795  ELSE
796  WRITE( nounit, fmt = 9999 )'DGEESX7', 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 dlacpy( 'F', n, n, a, lda, ht, lda )
826  CALL dgeesx( 'N', sort, dslect, '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 )'DGEESX8', iinfo, n, jtype,
833  \$ iseed
834  ELSE
835  WRITE( nounit, fmt = 9999 )'DGEESX8', 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 DSLECT 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 dlacpy( 'F', n, n, a, lda, ht, lda )
910  CALL dgeesx( 'N', 'S', dslect, '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 )'DGEESX9', 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 = dlange( '1', n, n, a, lda, work )
925  v = max( dble( 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( ' DGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
984  \$ i6, ', INPUT EXAMPLE NUMBER = ', i4 )
985  9998 FORMAT( ' DGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
986  \$ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
987 *
988  RETURN
989 *
990 * End of DGET24
991 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
DORT01
Definition: dort01.f:116
logical function dslect(ZR, ZI)
DSLECT
Definition: dslect.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:114
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:281
Here is the call graph for this function:
Here is the caller graph for this function: