LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ ddrgsx()

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.

Definition at line 356 of file ddrgsx.f.

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