LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zdrgsx()

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

ZDRGSX

Purpose:
 ZDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
 problem expert driver ZGGESX.

 ZGGES 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 ZDRGSX is called with NSIZE > 0, five (5) types of built-in
 matrix pairs are used to test the routine ZGGESX.

 When ZDRGSX is called with NSIZE = 0, it reads in test matrix data
 to test ZGGESX.
 (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 ZGGESX, 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 ZGESVD) 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 ZLAKF2
[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 INFO not equal to 0.)
[out]A
          A is COMPLEX*16 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*16 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*16 array, dimension (LDA, NSIZE)
          Copy of A, modified by ZGGESX.
[out]BI
          BI is COMPLEX*16 array, dimension (LDA, NSIZE)
          Copy of B, modified by ZGGESX.
[out]Z
          Z is COMPLEX*16 array, dimension (LDA, NSIZE)
          Z holds the left Schur vectors computed by ZGGESX.
[out]Q
          Q is COMPLEX*16 array, dimension (LDA, NSIZE)
          Q holds the right Schur vectors computed by ZGGESX.
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (NSIZE)
[out]BETA
          BETA is COMPLEX*16 array, dimension (NSIZE)

          On exit, ALPHA/BETA are the eigenvalues.
[out]C
          C is COMPLEX*16 array, dimension (LDC, LDC)
          Store the matrix generated by subroutine ZLAKF2, 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 COMPLEX*16 array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= 3*NSIZE*NSIZE/2
[out]RWORK
          RWORK is DOUBLE PRECISION 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.
Date
June 2016

Definition at line 351 of file zdrgsx.f.

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