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

SDRGSX

Purpose:
``` SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
problem expert driver SGGESX.

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

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

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 SGGESX, 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 SGESVD) 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 SLAKF2``` [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 IINFO not equal to 0.)``` [out] A ``` A is REAL 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 REAL 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 REAL array, dimension (LDA, NSIZE) Copy of A, modified by SGGESX.``` [out] BI ``` BI is REAL array, dimension (LDA, NSIZE) Copy of B, modified by SGGESX.``` [out] Z ``` Z is REAL array, dimension (LDA, NSIZE) Z holds the left Schur vectors computed by SGGESX.``` [out] Q ``` Q is REAL array, dimension (LDA, NSIZE) Q holds the right Schur vectors computed by SGGESX.``` [out] ALPHAR ` ALPHAR is REAL array, dimension (NSIZE)` [out] ALPHAI ` ALPHAI is REAL array, dimension (NSIZE)` [out] BETA ``` BETA is REAL array, dimension (NSIZE) \verbatim On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.``` [out] C ``` C is REAL array, dimension (LDC, LDC) Store the matrix generated by subroutine SLAKF2, 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 REAL 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.```
Date
June 2016

Definition at line 361 of file sdrgsx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: