 LAPACK  3.10.0 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.

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 ..
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
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
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
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
612  ELSE IF( i1.LT.mplusn-1 ) THEN
613  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
615  result( 5 ) = ulpinv
616  END IF
617  ELSE IF( i1.GT.1 ) THEN
618  IF( ai( i1, i1-1 ).NE.zero ) THEN
620  result( 5 ) = ulpinv
621  END IF
622  END IF
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 )
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
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
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
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
837  ELSE IF( i1.LT.mplusn-1 ) THEN
838  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
840  result( 5 ) = ulpinv
841  END IF
842  ELSE IF( i1.GT.1 ) THEN
843  IF( ai( i1, i1-1 ).NE.zero ) THEN
845  result( 5 ) = ulpinv
846  END IF
847  END IF
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 )
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 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: