LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ddrgsx ( integer  NSIZE,
integer  NCMAX,
double precision  THRESH,
integer  NIN,
integer  NOUT,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( lda, * )  B,
double precision, dimension( lda, * )  AI,
double precision, dimension( lda, * )  BI,
double precision, dimension( lda, * )  Z,
double precision, dimension( lda, * )  Q,
double precision, dimension( * )  ALPHAR,
double precision, dimension( * )  ALPHAI,
double precision, dimension( * )  BETA,
double precision, dimension( ldc, * )  C,
integer  LDC,
double precision, dimension( * )  S,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

DDRGSX

Purpose:
 DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
 problem expert driver DGGESX.

 DGGESX factors A and B as Q S Z' and Q T Z', where ' means
 transpose, T is upper triangular, S is in generalized Schur form
 (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
 the 2x2 blocks corresponding to complex conjugate pairs of
 generalized eigenvalues), and Q and Z are orthogonal.  It also
 computes the generalized eigenvalues (alpha(1),beta(1)), ...,
 (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the
 characteristic equation

     det( A - w(j) B ) = 0

 Optionally it also reorders the eigenvalues so that a selected
 cluster of eigenvalues appears in the leading diagonal block of the
 Schur forms; computes a reciprocal condition number for the average
 of the selected eigenvalues; and computes a reciprocal condition
 number for the right and left deflating subspaces corresponding to
 the selected eigenvalues.

 When DDRGSX is called with NSIZE > 0, five (5) types of built-in
 matrix pairs are used to test the routine DGGESX.

 When DDRGSX is called with NSIZE = 0, it reads in test matrix data
 to test DGGESX.

 For each matrix pair, the following tests will be performed and
 compared with the threshold THRESH except for the tests (7) and (9):

 (1)   | A - Q S Z' | / ( |A| n ulp )

 (2)   | B - Q T Z' | / ( |B| n ulp )

 (3)   | I - QQ' | / ( n ulp )

 (4)   | I - ZZ' | / ( n ulp )

 (5)   if A is in Schur form (i.e. quasi-triangular form)

 (6)   maximum over j of D(j)  where:

       if alpha(j) is real:
                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
           D(j) = ------------------------ + -----------------------
                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)

       if alpha(j) is complex:
                                 | det( s S - w T ) |
           D(j) = ---------------------------------------------------
                  ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )

           and S and T are here the 2 x 2 diagonal blocks of S and T
           corresponding to the j-th and j+1-th eigenvalues.

 (7)   if sorting worked and SDIM is the number of eigenvalues
       which were selected.

 (8)   the estimated value DIF does not differ from the true values of
       Difu and Difl more than a factor 10*THRESH. If the estimate DIF
       equals zero the corresponding true values of Difu and Difl
       should be less than EPS*norm(A, B). If the true value of Difu
       and Difl equal zero, the estimate DIF should be less than
       EPS*norm(A, B).

 (9)   If INFO = N+3 is returned by DGGESX, the reordering "failed"
       and we check that DIF = PL = PR = 0 and that the true value of
       Difu and Difl is < EPS*norm(A, B). We count the events when
       INFO=N+3.

 For read-in test matrices, the above tests are run except that the
 exact value for DIF (and PL) is input data.  Additionally, there is
 one more test run for read-in test matrices:

 (10)  the estimated value PL does not differ from the true value of
       PLTRU more than a factor THRESH. If the estimate PL equals
       zero the corresponding true value of PLTRU should be less than
       EPS*norm(A, B). If the true value of PLTRU equal zero, the
       estimate PL should be less than EPS*norm(A, B).

 Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
 matrix pairs are generated and tested. NSIZE should be kept small.

 SVD (routine DGESVD) is used for computing the true value of DIF_u
 and DIF_l when testing the built-in test problems.

 Built-in Test Matrices
 ======================

 All built-in test matrices are the 2 by 2 block of triangular
 matrices

          A = [ A11 A12 ]    and      B = [ B11 B12 ]
              [     A22 ]                 [     B22 ]

 where for different type of A11 and A22 are given as the following.
 A12 and B12 are chosen so that the generalized Sylvester equation

          A11*R - L*A22 = -A12
          B11*R - L*B22 = -B12

 have prescribed solution R and L.

 Type 1:  A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
          B11 = I_m, B22 = I_k
          where J_k(a,b) is the k-by-k Jordan block with ``a'' on
          diagonal and ``b'' on superdiagonal.

 Type 2:  A11 = (a_ij) = ( 2(.5-sin(i)) ) and
          B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
          A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
          B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k

 Type 3:  A11, A22 and B11, B22 are chosen as for Type 2, but each
          second diagonal block in A_11 and each third diagonal block
          in A_22 are made as 2 by 2 blocks.

 Type 4:  A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
             for i=1,...,m,  j=1,...,m and
          A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
             for i=m+1,...,k,  j=m+1,...,k

 Type 5:  (A,B) and have potentially close or common eigenvalues and
          very large departure from block diagonality A_11 is chosen
          as the m x m leading submatrix of A_1:
                  |  1  b                            |
                  | -b  1                            |
                  |        1+d  b                    |
                  |         -b 1+d                   |
           A_1 =  |                  d  1            |
                  |                 -1  d            |
                  |                        -d  1     |
                  |                        -1 -d     |
                  |                               1  |
          and A_22 is chosen as the k x k leading submatrix of A_2:
                  | -1  b                            |
                  | -b -1                            |
                  |       1-d  b                     |
                  |       -b  1-d                    |
           A_2 =  |                 d 1+b            |
                  |               -1-b d             |
                  |                       -d  1+b    |
                  |                      -1+b  -d    |
                  |                              1-d |
          and matrix B are chosen as identity matrices (see DLATM5).
Parameters
[in]NSIZE
          NSIZE is INTEGER
          The maximum size of the matrices to use. NSIZE >= 0.
          If NSIZE = 0, no built-in tests matrices are used, but
          read-in test matrices are used to test DGGESX.
[in]NCMAX
          NCMAX is INTEGER
          Maximum allowable NMAX for generating Kroneker matrix
          in call to DLAKF2
[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.  THRESH >= 0.
[in]NIN
          NIN is INTEGER
          The FORTRAN unit number for reading in the data file of
          problems to solve.
[in]NOUT
          NOUT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns IINFO not equal to 0.)
[out]A
          A is DOUBLE PRECISION array, dimension (LDA, NSIZE)
          Used to store the matrix whose eigenvalues are to be
          computed.  On exit, A contains the last matrix actually used.
[in]LDA
          LDA is INTEGER
          The leading dimension of A, B, AI, BI, Z and Q,
          LDA >= max( 1, NSIZE ). For the read-in test,
          LDA >= max( 1, N ), N is the size of the test matrices.
[out]B
          B is DOUBLE PRECISION array, dimension (LDA, NSIZE)
          Used to store the matrix whose eigenvalues are to be
          computed.  On exit, B contains the last matrix actually used.
[out]AI
          AI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
          Copy of A, modified by DGGESX.
[out]BI
          BI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
          Copy of B, modified by DGGESX.
[out]Z
          Z is DOUBLE PRECISION array, dimension (LDA, NSIZE)
          Z holds the left Schur vectors computed by DGGESX.
[out]Q
          Q is DOUBLE PRECISION array, dimension (LDA, NSIZE)
          Q holds the right Schur vectors computed by DGGESX.
[out]ALPHAR
          ALPHAR is DOUBLE PRECISION array, dimension (NSIZE)
[out]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (NSIZE)
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (NSIZE)

          On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
[out]C
          C is DOUBLE PRECISION array, dimension (LDC, LDC)
          Store the matrix generated by subroutine DLAKF2, this is the
          matrix formed by Kronecker products used for estimating
          DIF.
[in]LDC
          LDC is INTEGER
          The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
[out]S
          S is DOUBLE PRECISION array, dimension (LDC)
          Singular values of C
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
[out]IWORK
          IWORK is INTEGER array, dimension (LIWORK)
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK. LIWORK >= NSIZE + 6.
[out]BWORK
          BWORK is LOGICAL array, dimension (LDA)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  A routine returned an error code.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
June 2016

Definition at line 361 of file ddrgsx.f.

361 *
362 * -- LAPACK test routine (version 3.6.1) --
363 * -- LAPACK is a software package provided by Univ. of Tennessee, --
364 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
365 * June 2016
366 *
367 * .. Scalar Arguments ..
368  INTEGER info, lda, ldc, liwork, lwork, ncmax, nin,
369  $ nout, nsize
370  DOUBLE PRECISION thresh
371 * ..
372 * .. Array Arguments ..
373  LOGICAL bwork( * )
374  INTEGER iwork( * )
375  DOUBLE PRECISION a( lda, * ), ai( lda, * ), alphai( * ),
376  $ alphar( * ), b( lda, * ), beta( * ),
377  $ bi( lda, * ), c( ldc, * ), q( lda, * ), s( * ),
378  $ work( * ), z( lda, * )
379 * ..
380 *
381 * =====================================================================
382 *
383 * .. Parameters ..
384  DOUBLE PRECISION zero, one, ten
385  parameter ( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1 )
386 * ..
387 * .. Local Scalars ..
388  LOGICAL ilabad
389  CHARACTER sense
390  INTEGER bdspac, i, i1, ifunc, iinfo, j, linfo, maxwrk,
391  $ minwrk, mm, mn2, nerrs, nptknt, ntest, ntestt,
392  $ prtype, qba, qbb
393  DOUBLE PRECISION abnrm, bignum, diftru, pltru, smlnum, temp1,
394  $ temp2, thrsh2, ulp, ulpinv, weight
395 * ..
396 * .. Local Arrays ..
397  DOUBLE PRECISION difest( 2 ), pl( 2 ), result( 10 )
398 * ..
399 * .. External Functions ..
400  LOGICAL dlctsx
401  INTEGER ilaenv
402  DOUBLE PRECISION dlamch, dlange
403  EXTERNAL dlctsx, ilaenv, dlamch, dlange
404 * ..
405 * .. External Subroutines ..
406  EXTERNAL alasvm, dgesvd, dget51, dget53, dggesx, dlabad,
408 * ..
409 * .. Intrinsic Functions ..
410  INTRINSIC abs, max, sqrt
411 * ..
412 * .. Scalars in Common ..
413  LOGICAL fs
414  INTEGER k, m, mplusn, n
415 * ..
416 * .. Common blocks ..
417  COMMON / mn / m, n, mplusn, k, fs
418 * ..
419 * .. Executable Statements ..
420 *
421 * Check for errors
422 *
423  IF( nsize.LT.0 ) THEN
424  info = -1
425  ELSE IF( thresh.LT.zero ) THEN
426  info = -2
427  ELSE IF( nin.LE.0 ) THEN
428  info = -3
429  ELSE IF( nout.LE.0 ) THEN
430  info = -4
431  ELSE IF( lda.LT.1 .OR. lda.LT.nsize ) THEN
432  info = -6
433  ELSE IF( ldc.LT.1 .OR. ldc.LT.nsize*nsize / 2 ) THEN
434  info = -17
435  ELSE IF( liwork.LT.nsize+6 ) THEN
436  info = -21
437  END IF
438 *
439 * Compute workspace
440 * (Note: Comments in the code beginning "Workspace:" describe the
441 * minimal amount of workspace needed at that point in the code,
442 * as well as the preferred amount for good performance.
443 * NB refers to the optimal block size for the immediately
444 * following subroutine, as returned by ILAENV.)
445 *
446  minwrk = 1
447  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
448  minwrk = max( 10*( nsize+1 ), 5*nsize*nsize / 2 )
449 *
450 * workspace for sggesx
451 *
452  maxwrk = 9*( nsize+1 ) + nsize*
453  $ ilaenv( 1, 'DGEQRF', ' ', nsize, 1, nsize, 0 )
454  maxwrk = max( maxwrk, 9*( nsize+1 )+nsize*
455  $ ilaenv( 1, 'DORGQR', ' ', nsize, 1, nsize, -1 ) )
456 *
457 * workspace for dgesvd
458 *
459  bdspac = 5*nsize*nsize / 2
460  maxwrk = max( maxwrk, 3*nsize*nsize / 2+nsize*nsize*
461  $ ilaenv( 1, 'DGEBRD', ' ', nsize*nsize / 2,
462  $ nsize*nsize / 2, -1, -1 ) )
463  maxwrk = max( maxwrk, bdspac )
464 *
465  maxwrk = max( maxwrk, minwrk )
466 *
467  work( 1 ) = maxwrk
468  END IF
469 *
470  IF( lwork.LT.minwrk )
471  $ info = -19
472 *
473  IF( info.NE.0 ) THEN
474  CALL xerbla( 'DDRGSX', -info )
475  RETURN
476  END IF
477 *
478 * Important constants
479 *
480  ulp = dlamch( 'P' )
481  ulpinv = one / ulp
482  smlnum = dlamch( 'S' ) / ulp
483  bignum = one / smlnum
484  CALL dlabad( smlnum, bignum )
485  thrsh2 = ten*thresh
486  ntestt = 0
487  nerrs = 0
488 *
489 * Go to the tests for read-in matrix pairs
490 *
491  ifunc = 0
492  IF( nsize.EQ.0 )
493  $ GO TO 70
494 *
495 * Test the built-in matrix pairs.
496 * Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
497 * of test matrices, different size (M+N)
498 *
499  prtype = 0
500  qba = 3
501  qbb = 4
502  weight = sqrt( ulp )
503 *
504  DO 60 ifunc = 0, 3
505  DO 50 prtype = 1, 5
506  DO 40 m = 1, nsize - 1
507  DO 30 n = 1, nsize - m
508 *
509  weight = one / weight
510  mplusn = m + n
511 *
512 * Generate test matrices
513 *
514  fs = .true.
515  k = 0
516 *
517  CALL dlaset( 'Full', mplusn, mplusn, zero, zero, ai,
518  $ lda )
519  CALL dlaset( 'Full', mplusn, mplusn, zero, zero, bi,
520  $ lda )
521 *
522  CALL dlatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
523  $ lda, ai( 1, m+1 ), lda, bi, lda,
524  $ bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
525  $ q, lda, z, lda, weight, qba, qbb )
526 *
527 * Compute the Schur factorization and swapping the
528 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
529 * Swapping is accomplished via the function DLCTSX
530 * which is supplied below.
531 *
532  IF( ifunc.EQ.0 ) THEN
533  sense = 'N'
534  ELSE IF( ifunc.EQ.1 ) THEN
535  sense = 'E'
536  ELSE IF( ifunc.EQ.2 ) THEN
537  sense = 'V'
538  ELSE IF( ifunc.EQ.3 ) THEN
539  sense = 'B'
540  END IF
541 *
542  CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
543  CALL dlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
544 *
545  CALL dggesx( 'V', 'V', 'S', dlctsx, sense, mplusn, ai,
546  $ lda, bi, lda, mm, alphar, alphai, beta,
547  $ q, lda, z, lda, pl, difest, work, lwork,
548  $ iwork, liwork, bwork, linfo )
549 *
550  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
551  result( 1 ) = ulpinv
552  WRITE( nout, fmt = 9999 )'DGGESX', linfo, mplusn,
553  $ prtype
554  info = linfo
555  GO TO 30
556  END IF
557 *
558 * Compute the norm(A, B)
559 *
560  CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, work,
561  $ mplusn )
562  CALL dlacpy( 'Full', mplusn, mplusn, bi, lda,
563  $ work( mplusn*mplusn+1 ), mplusn )
564  abnrm = dlange( 'Fro', mplusn, 2*mplusn, work, mplusn,
565  $ work )
566 *
567 * Do tests (1) to (4)
568 *
569  CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
570  $ lda, work, result( 1 ) )
571  CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
572  $ lda, work, result( 2 ) )
573  CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
574  $ lda, work, result( 3 ) )
575  CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
576  $ lda, work, result( 4 ) )
577  ntest = 4
578 *
579 * Do tests (5) and (6): check Schur form of A and
580 * compare eigenvalues with diagonals.
581 *
582  temp1 = zero
583  result( 5 ) = zero
584  result( 6 ) = zero
585 *
586  DO 10 j = 1, mplusn
587  ilabad = .false.
588  IF( alphai( j ).EQ.zero ) THEN
589  temp2 = ( abs( alphar( j )-ai( j, j ) ) /
590  $ max( smlnum, abs( alphar( j ) ),
591  $ abs( ai( j, j ) ) )+
592  $ abs( beta( j )-bi( j, j ) ) /
593  $ max( smlnum, abs( beta( j ) ),
594  $ abs( bi( j, j ) ) ) ) / ulp
595  IF( j.LT.mplusn ) THEN
596  IF( ai( j+1, j ).NE.zero ) THEN
597  ilabad = .true.
598  result( 5 ) = ulpinv
599  END IF
600  END IF
601  IF( j.GT.1 ) THEN
602  IF( ai( j, j-1 ).NE.zero ) THEN
603  ilabad = .true.
604  result( 5 ) = ulpinv
605  END IF
606  END IF
607  ELSE
608  IF( alphai( j ).GT.zero ) THEN
609  i1 = j
610  ELSE
611  i1 = j - 1
612  END IF
613  IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
614  ilabad = .true.
615  ELSE IF( i1.LT.mplusn-1 ) THEN
616  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
617  ilabad = .true.
618  result( 5 ) = ulpinv
619  END IF
620  ELSE IF( i1.GT.1 ) THEN
621  IF( ai( i1, i1-1 ).NE.zero ) THEN
622  ilabad = .true.
623  result( 5 ) = ulpinv
624  END IF
625  END IF
626  IF( .NOT.ilabad ) THEN
627  CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ),
628  $ lda, beta( j ), alphar( j ),
629  $ alphai( j ), temp2, iinfo )
630  IF( iinfo.GE.3 ) THEN
631  WRITE( nout, fmt = 9997 )iinfo, j,
632  $ mplusn, prtype
633  info = abs( iinfo )
634  END IF
635  ELSE
636  temp2 = ulpinv
637  END IF
638  END IF
639  temp1 = max( temp1, temp2 )
640  IF( ilabad ) THEN
641  WRITE( nout, fmt = 9996 )j, mplusn, prtype
642  END IF
643  10 CONTINUE
644  result( 6 ) = temp1
645  ntest = ntest + 2
646 *
647 * Test (7) (if sorting worked)
648 *
649  result( 7 ) = zero
650  IF( linfo.EQ.mplusn+3 ) THEN
651  result( 7 ) = ulpinv
652  ELSE IF( mm.NE.n ) THEN
653  result( 7 ) = ulpinv
654  END IF
655  ntest = ntest + 1
656 *
657 * Test (8): compare the estimated value DIF and its
658 * value. first, compute the exact DIF.
659 *
660  result( 8 ) = zero
661  mn2 = mm*( mplusn-mm )*2
662  IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
663 *
664 * Note: for either following two causes, there are
665 * almost same number of test cases fail the test.
666 *
667  CALL dlakf2( mm, mplusn-mm, ai, lda,
668  $ ai( mm+1, mm+1 ), bi,
669  $ bi( mm+1, mm+1 ), c, ldc )
670 *
671  CALL dgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
672  $ 1, work( 2 ), 1, work( 3 ), lwork-2,
673  $ info )
674  diftru = s( mn2 )
675 *
676  IF( difest( 2 ).EQ.zero ) THEN
677  IF( diftru.GT.abnrm*ulp )
678  $ result( 8 ) = ulpinv
679  ELSE IF( diftru.EQ.zero ) THEN
680  IF( difest( 2 ).GT.abnrm*ulp )
681  $ result( 8 ) = ulpinv
682  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
683  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
684  result( 8 ) = max( diftru / difest( 2 ),
685  $ difest( 2 ) / diftru )
686  END IF
687  ntest = ntest + 1
688  END IF
689 *
690 * Test (9)
691 *
692  result( 9 ) = zero
693  IF( linfo.EQ.( mplusn+2 ) ) THEN
694  IF( diftru.GT.abnrm*ulp )
695  $ result( 9 ) = ulpinv
696  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
697  $ result( 9 ) = ulpinv
698  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
699  $ result( 9 ) = ulpinv
700  ntest = ntest + 1
701  END IF
702 *
703  ntestt = ntestt + ntest
704 *
705 * Print out tests which fail.
706 *
707  DO 20 j = 1, 9
708  IF( result( j ).GE.thresh ) THEN
709 *
710 * If this is the first test to fail,
711 * print a header to the data file.
712 *
713  IF( nerrs.EQ.0 ) THEN
714  WRITE( nout, fmt = 9995 )'DGX'
715 *
716 * Matrix types
717 *
718  WRITE( nout, fmt = 9993 )
719 *
720 * Tests performed
721 *
722  WRITE( nout, fmt = 9992 )'orthogonal', '''',
723  $ 'transpose', ( '''', i = 1, 4 )
724 *
725  END IF
726  nerrs = nerrs + 1
727  IF( result( j ).LT.10000.0d0 ) THEN
728  WRITE( nout, fmt = 9991 )mplusn, prtype,
729  $ weight, m, j, result( j )
730  ELSE
731  WRITE( nout, fmt = 9990 )mplusn, prtype,
732  $ weight, m, j, result( j )
733  END IF
734  END IF
735  20 CONTINUE
736 *
737  30 CONTINUE
738  40 CONTINUE
739  50 CONTINUE
740  60 CONTINUE
741 *
742  GO TO 150
743 *
744  70 CONTINUE
745 *
746 * Read in data from file to check accuracy of condition estimation
747 * Read input data until N=0
748 *
749  nptknt = 0
750 *
751  80 CONTINUE
752  READ( nin, fmt = *, end = 140 )mplusn
753  IF( mplusn.EQ.0 )
754  $ GO TO 140
755  READ( nin, fmt = *, end = 140 )n
756  DO 90 i = 1, mplusn
757  READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
758  90 CONTINUE
759  DO 100 i = 1, mplusn
760  READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
761  100 CONTINUE
762  READ( nin, fmt = * )pltru, diftru
763 *
764  nptknt = nptknt + 1
765  fs = .true.
766  k = 0
767  m = mplusn - n
768 *
769  CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
770  CALL dlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
771 *
772 * Compute the Schur factorization while swaping the
773 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
774 *
775  CALL dggesx( 'V', 'V', 'S', dlctsx, 'B', mplusn, ai, lda, bi, lda,
776  $ mm, alphar, alphai, beta, q, lda, z, lda, pl, difest,
777  $ work, lwork, iwork, liwork, bwork, linfo )
778 *
779  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
780  result( 1 ) = ulpinv
781  WRITE( nout, fmt = 9998 )'DGGESX', linfo, mplusn, nptknt
782  GO TO 130
783  END IF
784 *
785 * Compute the norm(A, B)
786 * (should this be norm of (A,B) or (AI,BI)?)
787 *
788  CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
789  CALL dlacpy( 'Full', mplusn, mplusn, bi, lda,
790  $ work( mplusn*mplusn+1 ), mplusn )
791  abnrm = dlange( 'Fro', mplusn, 2*mplusn, work, mplusn, work )
792 *
793 * Do tests (1) to (4)
794 *
795  CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
796  $ result( 1 ) )
797  CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
798  $ result( 2 ) )
799  CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
800  $ result( 3 ) )
801  CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
802  $ result( 4 ) )
803 *
804 * Do tests (5) and (6): check Schur form of A and compare
805 * eigenvalues with diagonals.
806 *
807  ntest = 6
808  temp1 = zero
809  result( 5 ) = zero
810  result( 6 ) = zero
811 *
812  DO 110 j = 1, mplusn
813  ilabad = .false.
814  IF( alphai( j ).EQ.zero ) THEN
815  temp2 = ( abs( alphar( j )-ai( j, j ) ) /
816  $ max( smlnum, abs( alphar( j ) ), abs( ai( j,
817  $ j ) ) )+abs( beta( j )-bi( j, j ) ) /
818  $ max( smlnum, abs( beta( j ) ), abs( bi( j, j ) ) ) )
819  $ / ulp
820  IF( j.LT.mplusn ) THEN
821  IF( ai( j+1, j ).NE.zero ) THEN
822  ilabad = .true.
823  result( 5 ) = ulpinv
824  END IF
825  END IF
826  IF( j.GT.1 ) THEN
827  IF( ai( j, j-1 ).NE.zero ) THEN
828  ilabad = .true.
829  result( 5 ) = ulpinv
830  END IF
831  END IF
832  ELSE
833  IF( alphai( j ).GT.zero ) THEN
834  i1 = j
835  ELSE
836  i1 = j - 1
837  END IF
838  IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
839  ilabad = .true.
840  ELSE IF( i1.LT.mplusn-1 ) THEN
841  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
842  ilabad = .true.
843  result( 5 ) = ulpinv
844  END IF
845  ELSE IF( i1.GT.1 ) THEN
846  IF( ai( i1, i1-1 ).NE.zero ) THEN
847  ilabad = .true.
848  result( 5 ) = ulpinv
849  END IF
850  END IF
851  IF( .NOT.ilabad ) THEN
852  CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ), lda,
853  $ beta( j ), alphar( j ), alphai( j ), temp2,
854  $ iinfo )
855  IF( iinfo.GE.3 ) THEN
856  WRITE( nout, fmt = 9997 )iinfo, j, mplusn, nptknt
857  info = abs( iinfo )
858  END IF
859  ELSE
860  temp2 = ulpinv
861  END IF
862  END IF
863  temp1 = max( temp1, temp2 )
864  IF( ilabad ) THEN
865  WRITE( nout, fmt = 9996 )j, mplusn, nptknt
866  END IF
867  110 CONTINUE
868  result( 6 ) = temp1
869 *
870 * Test (7) (if sorting worked) <--------- need to be checked.
871 *
872  ntest = 7
873  result( 7 ) = zero
874  IF( linfo.EQ.mplusn+3 )
875  $ result( 7 ) = ulpinv
876 *
877 * Test (8): compare the estimated value of DIF and its true value.
878 *
879  ntest = 8
880  result( 8 ) = zero
881  IF( difest( 2 ).EQ.zero ) THEN
882  IF( diftru.GT.abnrm*ulp )
883  $ result( 8 ) = ulpinv
884  ELSE IF( diftru.EQ.zero ) THEN
885  IF( difest( 2 ).GT.abnrm*ulp )
886  $ result( 8 ) = ulpinv
887  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
888  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
889  result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
890  END IF
891 *
892 * Test (9)
893 *
894  ntest = 9
895  result( 9 ) = zero
896  IF( linfo.EQ.( mplusn+2 ) ) THEN
897  IF( diftru.GT.abnrm*ulp )
898  $ result( 9 ) = ulpinv
899  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
900  $ result( 9 ) = ulpinv
901  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
902  $ result( 9 ) = ulpinv
903  END IF
904 *
905 * Test (10): compare the estimated value of PL and it true value.
906 *
907  ntest = 10
908  result( 10 ) = zero
909  IF( pl( 1 ).EQ.zero ) THEN
910  IF( pltru.GT.abnrm*ulp )
911  $ result( 10 ) = ulpinv
912  ELSE IF( pltru.EQ.zero ) THEN
913  IF( pl( 1 ).GT.abnrm*ulp )
914  $ result( 10 ) = ulpinv
915  ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
916  $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
917  result( 10 ) = ulpinv
918  END IF
919 *
920  ntestt = ntestt + ntest
921 *
922 * Print out tests which fail.
923 *
924  DO 120 j = 1, ntest
925  IF( result( j ).GE.thresh ) THEN
926 *
927 * If this is the first test to fail,
928 * print a header to the data file.
929 *
930  IF( nerrs.EQ.0 ) THEN
931  WRITE( nout, fmt = 9995 )'DGX'
932 *
933 * Matrix types
934 *
935  WRITE( nout, fmt = 9994 )
936 *
937 * Tests performed
938 *
939  WRITE( nout, fmt = 9992 )'orthogonal', '''',
940  $ 'transpose', ( '''', i = 1, 4 )
941 *
942  END IF
943  nerrs = nerrs + 1
944  IF( result( j ).LT.10000.0d0 ) THEN
945  WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
946  ELSE
947  WRITE( nout, fmt = 9988 )nptknt, mplusn, j, result( j )
948  END IF
949  END IF
950 *
951  120 CONTINUE
952 *
953  130 CONTINUE
954  GO TO 80
955  140 CONTINUE
956 *
957  150 CONTINUE
958 *
959 * Summary
960 *
961  CALL alasvm( 'DGX', nout, nerrs, ntestt, 0 )
962 *
963  work( 1 ) = maxwrk
964 *
965  RETURN
966 *
967  9999 FORMAT( ' DDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
968  $ i6, ', JTYPE=', i6, ')' )
969 *
970  9998 FORMAT( ' DDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
971  $ i6, ', Input Example #', i2, ')' )
972 *
973  9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', i1, ' for eigenvalue ',
974  $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
975 *
976  9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', i6, '.',
977  $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
978 *
979  9995 FORMAT( / 1x, a3, ' -- Real Expert Generalized Schur form',
980  $ ' problem driver' )
981 *
982  9994 FORMAT( 'Input Example' )
983 *
984  9993 FORMAT( ' Matrix types: ', /
985  $ ' 1: A is a block diagonal matrix of Jordan blocks ',
986  $ 'and B is the identity ', / ' matrix, ',
987  $ / ' 2: A and B are upper triangular matrices, ',
988  $ / ' 3: A and B are as type 2, but each second diagonal ',
989  $ 'block in A_11 and ', /
990  $ ' each third diaongal block in A_22 are 2x2 blocks,',
991  $ / ' 4: A and B are block diagonal matrices, ',
992  $ / ' 5: (A,B) has potentially close or common ',
993  $ 'eigenvalues.', / )
994 *
995  9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
996  $ 'Q and Z are ', a, ',', / 19x,
997  $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
998  $ / ' 1 = | A - Q S Z', a,
999  $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
1000  $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
1001  $ ' | / ( n ulp ) 4 = | I - ZZ', a,
1002  $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
1003  $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
1004  $ ' and diagonals of (S,T)', /
1005  $ ' 7 = 1/ULP if SDIM is not the correct number of ',
1006  $ 'selected eigenvalues', /
1007  $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
1008  $ 'DIFTRU/DIFEST > 10*THRESH',
1009  $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
1010  $ 'when reordering fails', /
1011  $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
1012  $ 'PLTRU/PLEST > THRESH', /
1013  $ ' ( Test 10 is only for input examples )', / )
1014  9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
1015  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
1016  9990 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
1017  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, d10.3 )
1018  9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1019  $ ' result ', i2, ' is', 0p, f8.2 )
1020  9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1021  $ ' result ', i2, ' is', 1p, d10.3 )
1022 *
1023 * End of DDRGSX
1024 *
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
DLAKF2
Definition: dlakf2.f:107
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 dlatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
DLATM5
Definition: dlatm5.f:270
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
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 dlctsx(AR, AI, BETA)
DLCTSX
Definition: dlctsx.f:67
subroutine dget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
DGET51
Definition: dget51.f:151
subroutine dget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
DGET53
Definition: dget53.f:128
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine dggesx(JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
DGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
Definition: dggesx.f:367
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
DGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: dgesvd.f:213

Here is the call graph for this function:

Here is the caller graph for this function: