LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cdrgsx()

subroutine cdrgsx ( integer  NSIZE,
integer  NCMAX,
real  THRESH,
integer  NIN,
integer  NOUT,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( lda, * )  B,
complex, dimension( lda, * )  AI,
complex, dimension( lda, * )  BI,
complex, dimension( lda, * )  Z,
complex, dimension( lda, * )  Q,
complex, dimension( * )  ALPHA,
complex, dimension( * )  BETA,
complex, dimension( ldc, * )  C,
integer  LDC,
real, dimension( * )  S,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

CDRGSX

Purpose:
 CDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
 problem expert driver CGGESX.

 CGGES factors A and B as Q*S*Z'  and Q*T*Z' , where ' means conjugate
 transpose, S and T are  upper triangular (i.e., in generalized Schur
 form), and Q and Z are unitary. It also computes the generalized
 eigenvalues (alpha(j),beta(j)), j=1,...,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 CDRGSX is called with NSIZE > 0, five (5) types of built-in
 matrix pairs are used to test the routine CGGESX.

 When CDRGSX is called with NSIZE = 0, it reads in test matrix data
 to test CGGESX.
 (need more details on what kind of read-in data are needed).

 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. triangular form)

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

                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
           D(j) = ------------------------ + -----------------------
                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)

 (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 CGGESX, 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 same 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 CGESVD) 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 SLATM5).
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 SGGESX.
[in]NCMAX
          NCMAX is INTEGER
          Maximum allowable NMAX for generating Kroneker matrix
          in call to CLAKF2
[in]THRESH
          THRESH is REAL
          A test will count as "failed" if the "error", computed as
          described above, exceeds THRESH.  Note that the error
          is scaled to be O(1), so THRESH should be a reasonably
          small multiple of 1, e.g., 10 or 100.  In particular,
          it should not depend on the precision (single vs. double)
          or the size of the matrix.  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 INFO not equal to 0.)
[out]A
          A is COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (LDA, NSIZE)
          Copy of A, modified by CGGESX.
[out]BI
          BI is COMPLEX array, dimension (LDA, NSIZE)
          Copy of B, modified by CGGESX.
[out]Z
          Z is COMPLEX array, dimension (LDA, NSIZE)
          Z holds the left Schur vectors computed by CGGESX.
[out]Q
          Q is COMPLEX array, dimension (LDA, NSIZE)
          Q holds the right Schur vectors computed by CGGESX.
[out]ALPHA
          ALPHA is COMPLEX array, dimension (NSIZE)
[out]BETA
          BETA is COMPLEX array, dimension (NSIZE)

          On exit, ALPHA/BETA are the eigenvalues.
[out]C
          C is COMPLEX array, dimension (LDC, LDC)
          Store the matrix generated by subroutine CLAKF2, 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 REAL array, dimension (LDC)
          Singular values of C
[out]WORK
          WORK is COMPLEX array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= 3*NSIZE*NSIZE/2
[out]RWORK
          RWORK is REAL array,
                                 dimension (5*NSIZE*NSIZE/2 - 4)
[out]IWORK
          IWORK is INTEGER array, dimension (LIWORK)
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK. LIWORK >= NSIZE + 2.
[out]BWORK
          BWORK is LOGICAL array, dimension (NSIZE)
[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 346 of file cdrgsx.f.

349 *
350 * -- LAPACK test routine --
351 * -- LAPACK is a software package provided by Univ. of Tennessee, --
352 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
353 *
354 * .. Scalar Arguments ..
355  INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
356  $ NOUT, NSIZE
357  REAL THRESH
358 * ..
359 * .. Array Arguments ..
360  LOGICAL BWORK( * )
361  INTEGER IWORK( * )
362  REAL RWORK( * ), S( * )
363  COMPLEX A( LDA, * ), AI( LDA, * ), ALPHA( * ),
364  $ B( LDA, * ), BETA( * ), BI( LDA, * ),
365  $ C( LDC, * ), Q( LDA, * ), WORK( * ),
366  $ Z( LDA, * )
367 * ..
368 *
369 * =====================================================================
370 *
371 * .. Parameters ..
372  REAL ZERO, ONE, TEN
373  parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1 )
374  COMPLEX CZERO
375  parameter( czero = ( 0.0e+0, 0.0e+0 ) )
376 * ..
377 * .. Local Scalars ..
378  LOGICAL ILABAD
379  CHARACTER SENSE
380  INTEGER BDSPAC, I, IFUNC, J, LINFO, MAXWRK, MINWRK, MM,
381  $ MN2, NERRS, NPTKNT, NTEST, NTESTT, PRTYPE, QBA,
382  $ QBB
383  REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
384  $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
385  COMPLEX X
386 * ..
387 * .. Local Arrays ..
388  REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 )
389 * ..
390 * .. External Functions ..
391  LOGICAL CLCTSX
392  INTEGER ILAENV
393  REAL CLANGE, SLAMCH
394  EXTERNAL clctsx, ilaenv, clange, slamch
395 * ..
396 * .. External Subroutines ..
397  EXTERNAL alasvm, cgesvd, cget51, cggesx, clacpy, clakf2,
399 * ..
400 * .. Scalars in Common ..
401  LOGICAL FS
402  INTEGER K, M, MPLUSN, N
403 * ..
404 * .. Common blocks ..
405  COMMON / mn / m, n, mplusn, k, fs
406 * ..
407 * .. Intrinsic Functions ..
408  INTRINSIC abs, aimag, max, real, sqrt
409 * ..
410 * .. Statement Functions ..
411  REAL ABS1
412 * ..
413 * .. Statement Function definitions ..
414  abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
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 = -15
432  ELSE IF( liwork.LT.nsize+2 ) 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 = 3*nsize*nsize / 2
446 *
447 * workspace for cggesx
448 *
449  maxwrk = nsize*( 1+ilaenv( 1, 'CGEQRF', ' ', nsize, 1, nsize,
450  $ 0 ) )
451  maxwrk = max( maxwrk, nsize*( 1+ilaenv( 1, 'CUNGQR', ' ',
452  $ nsize, 1, nsize, -1 ) ) )
453 *
454 * workspace for cgesvd
455 *
456  bdspac = 3*nsize*nsize / 2
457  maxwrk = max( maxwrk, nsize*nsize*
458  $ ( 1+ilaenv( 1, 'CGEBRD', ' ', 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 = -18
469 *
470  IF( info.NE.0 ) THEN
471  CALL xerbla( 'CDRGSX', -info )
472  RETURN
473  END IF
474 *
475 * Important constants
476 *
477  ulp = slamch( 'P' )
478  ulpinv = one / ulp
479  smlnum = slamch( 'S' ) / ulp
480  bignum = one / smlnum
481  CALL slabad( 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 CGGESX, 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 claset( 'Full', mplusn, mplusn, czero, czero, ai,
515  $ lda )
516  CALL claset( 'Full', mplusn, mplusn, czero, czero, bi,
517  $ lda )
518 *
519  CALL clatm5( 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 CLCTSX
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 clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
540  CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
541 *
542  CALL cggesx( 'V', 'V', 'S', clctsx, sense, mplusn, ai,
543  $ lda, bi, lda, mm, alpha, beta, q, lda, z,
544  $ lda, pl, difest, work, lwork, rwork,
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 )'CGGESX', linfo, mplusn,
550  $ prtype
551  info = linfo
552  GO TO 30
553  END IF
554 *
555 * Compute the norm(A, B)
556 *
557  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work,
558  $ mplusn )
559  CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
560  $ work( mplusn*mplusn+1 ), mplusn )
561  abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn,
562  $ rwork )
563 *
564 * Do tests (1) to (4)
565 *
566  result( 2 ) = zero
567  CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
568  $ lda, work, rwork, result( 1 ) )
569  CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
570  $ lda, work, rwork, result( 2 ) )
571  CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
572  $ lda, work, rwork, result( 3 ) )
573  CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
574  $ lda, work, rwork, result( 4 ) )
575  ntest = 4
576 *
577 * Do tests (5) and (6): check Schur form of A and
578 * compare eigenvalues with diagonals.
579 *
580  temp1 = zero
581  result( 5 ) = zero
582  result( 6 ) = zero
583 *
584  DO 10 j = 1, mplusn
585  ilabad = .false.
586  temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
587  $ max( smlnum, abs1( alpha( j ) ),
588  $ abs1( ai( j, j ) ) )+
589  $ abs1( beta( j )-bi( j, j ) ) /
590  $ max( smlnum, abs1( beta( j ) ),
591  $ abs1( 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  temp1 = max( temp1, temp2 )
605  IF( ilabad ) THEN
606  WRITE( nout, fmt = 9997 )j, mplusn, prtype
607  END IF
608  10 CONTINUE
609  result( 6 ) = temp1
610  ntest = ntest + 2
611 *
612 * Test (7) (if sorting worked)
613 *
614  result( 7 ) = zero
615  IF( linfo.EQ.mplusn+3 ) THEN
616  result( 7 ) = ulpinv
617  ELSE IF( mm.NE.n ) THEN
618  result( 7 ) = ulpinv
619  END IF
620  ntest = ntest + 1
621 *
622 * Test (8): compare the estimated value DIF and its
623 * value. first, compute the exact DIF.
624 *
625  result( 8 ) = zero
626  mn2 = mm*( mplusn-mm )*2
627  IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
628 *
629 * Note: for either following two cases, there are
630 * almost same number of test cases fail the test.
631 *
632  CALL clakf2( mm, mplusn-mm, ai, lda,
633  $ ai( mm+1, mm+1 ), bi,
634  $ bi( mm+1, mm+1 ), c, ldc )
635 *
636  CALL cgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
637  $ 1, work( 2 ), 1, work( 3 ), lwork-2,
638  $ rwork, info )
639  diftru = s( mn2 )
640 *
641  IF( difest( 2 ).EQ.zero ) THEN
642  IF( diftru.GT.abnrm*ulp )
643  $ result( 8 ) = ulpinv
644  ELSE IF( diftru.EQ.zero ) THEN
645  IF( difest( 2 ).GT.abnrm*ulp )
646  $ result( 8 ) = ulpinv
647  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
648  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
649  result( 8 ) = max( diftru / difest( 2 ),
650  $ difest( 2 ) / diftru )
651  END IF
652  ntest = ntest + 1
653  END IF
654 *
655 * Test (9)
656 *
657  result( 9 ) = zero
658  IF( linfo.EQ.( mplusn+2 ) ) THEN
659  IF( diftru.GT.abnrm*ulp )
660  $ result( 9 ) = ulpinv
661  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
662  $ result( 9 ) = ulpinv
663  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
664  $ result( 9 ) = ulpinv
665  ntest = ntest + 1
666  END IF
667 *
668  ntestt = ntestt + ntest
669 *
670 * Print out tests which fail.
671 *
672  DO 20 j = 1, 9
673  IF( result( j ).GE.thresh ) THEN
674 *
675 * If this is the first test to fail,
676 * print a header to the data file.
677 *
678  IF( nerrs.EQ.0 ) THEN
679  WRITE( nout, fmt = 9996 )'CGX'
680 *
681 * Matrix types
682 *
683  WRITE( nout, fmt = 9994 )
684 *
685 * Tests performed
686 *
687  WRITE( nout, fmt = 9993 )'unitary', '''',
688  $ 'transpose', ( '''', i = 1, 4 )
689 *
690  END IF
691  nerrs = nerrs + 1
692  IF( result( j ).LT.10000.0 ) THEN
693  WRITE( nout, fmt = 9992 )mplusn, prtype,
694  $ weight, m, j, result( j )
695  ELSE
696  WRITE( nout, fmt = 9991 )mplusn, prtype,
697  $ weight, m, j, result( j )
698  END IF
699  END IF
700  20 CONTINUE
701 *
702  30 CONTINUE
703  40 CONTINUE
704  50 CONTINUE
705  60 CONTINUE
706 *
707  GO TO 150
708 *
709  70 CONTINUE
710 *
711 * Read in data from file to check accuracy of condition estimation
712 * Read input data until N=0
713 *
714  nptknt = 0
715 *
716  80 CONTINUE
717  READ( nin, fmt = *, END = 140 )mplusn
718  IF( mplusn.EQ.0 )
719  $ GO TO 140
720  READ( nin, fmt = *, END = 140 )n
721  DO 90 i = 1, mplusn
722  READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
723  90 CONTINUE
724  DO 100 i = 1, mplusn
725  READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
726  100 CONTINUE
727  READ( nin, fmt = * )pltru, diftru
728 *
729  nptknt = nptknt + 1
730  fs = .true.
731  k = 0
732  m = mplusn - n
733 *
734  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
735  CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
736 *
737 * Compute the Schur factorization while swapping the
738 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
739 *
740  CALL cggesx( 'V', 'V', 'S', clctsx, 'B', mplusn, ai, lda, bi, lda,
741  $ mm, alpha, beta, q, lda, z, lda, pl, difest, work,
742  $ lwork, rwork, iwork, liwork, bwork, linfo )
743 *
744  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
745  result( 1 ) = ulpinv
746  WRITE( nout, fmt = 9998 )'CGGESX', linfo, mplusn, nptknt
747  GO TO 130
748  END IF
749 *
750 * Compute the norm(A, B)
751 * (should this be norm of (A,B) or (AI,BI)?)
752 *
753  CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
754  CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
755  $ work( mplusn*mplusn+1 ), mplusn )
756  abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn, rwork )
757 *
758 * Do tests (1) to (4)
759 *
760  CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
761  $ rwork, result( 1 ) )
762  CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
763  $ rwork, result( 2 ) )
764  CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
765  $ rwork, result( 3 ) )
766  CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
767  $ rwork, result( 4 ) )
768 *
769 * Do tests (5) and (6): check Schur form of A and compare
770 * eigenvalues with diagonals.
771 *
772  ntest = 6
773  temp1 = zero
774  result( 5 ) = zero
775  result( 6 ) = zero
776 *
777  DO 110 j = 1, mplusn
778  ilabad = .false.
779  temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
780  $ max( smlnum, abs1( alpha( j ) ), abs1( ai( j, j ) ) )+
781  $ abs1( beta( j )-bi( j, j ) ) /
782  $ max( smlnum, abs1( beta( j ) ), abs1( bi( j, j ) ) ) )
783  $ / ulp
784  IF( j.LT.mplusn ) THEN
785  IF( ai( j+1, j ).NE.zero ) THEN
786  ilabad = .true.
787  result( 5 ) = ulpinv
788  END IF
789  END IF
790  IF( j.GT.1 ) THEN
791  IF( ai( j, j-1 ).NE.zero ) THEN
792  ilabad = .true.
793  result( 5 ) = ulpinv
794  END IF
795  END IF
796  temp1 = max( temp1, temp2 )
797  IF( ilabad ) THEN
798  WRITE( nout, fmt = 9997 )j, mplusn, nptknt
799  END IF
800  110 CONTINUE
801  result( 6 ) = temp1
802 *
803 * Test (7) (if sorting worked) <--------- need to be checked.
804 *
805  ntest = 7
806  result( 7 ) = zero
807  IF( linfo.EQ.mplusn+3 )
808  $ result( 7 ) = ulpinv
809 *
810 * Test (8): compare the estimated value of DIF and its true value.
811 *
812  ntest = 8
813  result( 8 ) = zero
814  IF( difest( 2 ).EQ.zero ) THEN
815  IF( diftru.GT.abnrm*ulp )
816  $ result( 8 ) = ulpinv
817  ELSE IF( diftru.EQ.zero ) THEN
818  IF( difest( 2 ).GT.abnrm*ulp )
819  $ result( 8 ) = ulpinv
820  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
821  $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
822  result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
823  END IF
824 *
825 * Test (9)
826 *
827  ntest = 9
828  result( 9 ) = zero
829  IF( linfo.EQ.( mplusn+2 ) ) THEN
830  IF( diftru.GT.abnrm*ulp )
831  $ result( 9 ) = ulpinv
832  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
833  $ result( 9 ) = ulpinv
834  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
835  $ result( 9 ) = ulpinv
836  END IF
837 *
838 * Test (10): compare the estimated value of PL and it true value.
839 *
840  ntest = 10
841  result( 10 ) = zero
842  IF( pl( 1 ).EQ.zero ) THEN
843  IF( pltru.GT.abnrm*ulp )
844  $ result( 10 ) = ulpinv
845  ELSE IF( pltru.EQ.zero ) THEN
846  IF( pl( 1 ).GT.abnrm*ulp )
847  $ result( 10 ) = ulpinv
848  ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
849  $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
850  result( 10 ) = ulpinv
851  END IF
852 *
853  ntestt = ntestt + ntest
854 *
855 * Print out tests which fail.
856 *
857  DO 120 j = 1, ntest
858  IF( result( j ).GE.thresh ) THEN
859 *
860 * If this is the first test to fail,
861 * print a header to the data file.
862 *
863  IF( nerrs.EQ.0 ) THEN
864  WRITE( nout, fmt = 9996 )'CGX'
865 *
866 * Matrix types
867 *
868  WRITE( nout, fmt = 9995 )
869 *
870 * Tests performed
871 *
872  WRITE( nout, fmt = 9993 )'unitary', '''', 'transpose',
873  $ ( '''', i = 1, 4 )
874 *
875  END IF
876  nerrs = nerrs + 1
877  IF( result( j ).LT.10000.0 ) THEN
878  WRITE( nout, fmt = 9990 )nptknt, mplusn, j, result( j )
879  ELSE
880  WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
881  END IF
882  END IF
883 *
884  120 CONTINUE
885 *
886  130 CONTINUE
887  GO TO 80
888  140 CONTINUE
889 *
890  150 CONTINUE
891 *
892 * Summary
893 *
894  CALL alasvm( 'CGX', nout, nerrs, ntestt, 0 )
895 *
896  work( 1 ) = maxwrk
897 *
898  RETURN
899 *
900  9999 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
901  $ i6, ', JTYPE=', i6, ')' )
902 *
903  9998 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
904  $ i6, ', Input Example #', i2, ')' )
905 *
906  9997 FORMAT( ' CDRGSX: S not in Schur form at eigenvalue ', i6, '.',
907  $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
908 *
909  9996 FORMAT( / 1x, a3, ' -- Complex Expert Generalized Schur form',
910  $ ' problem driver' )
911 *
912  9995 FORMAT( 'Input Example' )
913 *
914  9994 FORMAT( ' Matrix types: ', /
915  $ ' 1: A is a block diagonal matrix of Jordan blocks ',
916  $ 'and B is the identity ', / ' matrix, ',
917  $ / ' 2: A and B are upper triangular matrices, ',
918  $ / ' 3: A and B are as type 2, but each second diagonal ',
919  $ 'block in A_11 and ', /
920  $ ' each third diaongal block in A_22 are 2x2 blocks,',
921  $ / ' 4: A and B are block diagonal matrices, ',
922  $ / ' 5: (A,B) has potentially close or common ',
923  $ 'eigenvalues.', / )
924 *
925  9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
926  $ 'Q and Z are ', a, ',', / 19x,
927  $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
928  $ / ' 1 = | A - Q S Z', a,
929  $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
930  $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
931  $ ' | / ( n ulp ) 4 = | I - ZZ', a,
932  $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
933  $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
934  $ ' and diagonals of (S,T)', /
935  $ ' 7 = 1/ULP if SDIM is not the correct number of ',
936  $ 'selected eigenvalues', /
937  $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
938  $ 'DIFTRU/DIFEST > 10*THRESH',
939  $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
940  $ 'when reordering fails', /
941  $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
942  $ 'PLTRU/PLEST > THRESH', /
943  $ ' ( Test 10 is only for input examples )', / )
944  9992 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
945  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
946  9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
947  $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, e10.3 )
948  9990 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
949  $ ' result ', i2, ' is', 0p, f8.2 )
950  9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
951  $ ' result ', i2, ' is', 1p, e10.3 )
952 *
953 * End of CDRGSX
954 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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
subroutine cget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RWORK, RESULT)
CGET51
Definition: cget51.f:155
logical function clctsx(ALPHA, BETA)
CLCTSX
Definition: clctsx.f:57
subroutine clatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
CLATM5
Definition: clatm5.f:268
subroutine clakf2(M, N, A, LDA, B, D, E, Z, LDZ)
CLAKF2
Definition: clakf2.f:105
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:115
subroutine cggesx(JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
CGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: cggesx.f:330
subroutine cgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
CGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: cgesvd.f:214
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: