LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ sdrgsx()

 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.```

Definition at line 356 of file sdrgsx.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  REAL THRESH
368 * ..
369 * .. Array Arguments ..
370  LOGICAL BWORK( * )
371  INTEGER IWORK( * )
372  REAL 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  REAL ZERO, ONE, TEN
382  parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+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  REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
391  \$ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
392 * ..
393 * .. Local Arrays ..
394  REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 )
395 * ..
396 * .. External Functions ..
397  LOGICAL SLCTSX
398  INTEGER ILAENV
399  REAL SLAMCH, SLANGE
400  EXTERNAL slctsx, ilaenv, slamch, slange
401 * ..
402 * .. External Subroutines ..
403  EXTERNAL alasvm, sgesvd, sget51, sget53, sggesx, slabad,
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 c MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 )
446  minwrk = max( 10*( nsize+1 ), 5*nsize*nsize / 2 )
447 *
448 * workspace for sggesx
449 *
450  maxwrk = 9*( nsize+1 ) + nsize*
451  \$ ilaenv( 1, 'SGEQRF', ' ', nsize, 1, nsize, 0 )
452  maxwrk = max( maxwrk, 9*( nsize+1 )+nsize*
453  \$ ilaenv( 1, 'SORGQR', ' ', nsize, 1, nsize, -1 ) )
454 *
455 * workspace for sgesvd
456 *
457  bdspac = 5*nsize*nsize / 2
458  maxwrk = max( maxwrk, 3*nsize*nsize / 2+nsize*nsize*
459  \$ ilaenv( 1, 'SGEBRD', ' ', nsize*nsize / 2,
460  \$ nsize*nsize / 2, -1, -1 ) )
461  maxwrk = max( maxwrk, bdspac )
462 *
463  maxwrk = max( maxwrk, minwrk )
464 *
465  work( 1 ) = maxwrk
466  END IF
467 *
468  IF( lwork.LT.minwrk )
469  \$ info = -19
470 *
471  IF( info.NE.0 ) THEN
472  CALL xerbla( 'SDRGSX', -info )
473  RETURN
474  END IF
475 *
476 * Important constants
477 *
478  ulp = slamch( 'P' )
479  ulpinv = one / ulp
480  smlnum = slamch( 'S' ) / ulp
481  bignum = one / smlnum
482  CALL slabad( smlnum, bignum )
483  thrsh2 = ten*thresh
484  ntestt = 0
485  nerrs = 0
486 *
487 * Go to the tests for read-in matrix pairs
488 *
489  ifunc = 0
490  IF( nsize.EQ.0 )
491  \$ GO TO 70
492 *
493 * Test the built-in matrix pairs.
494 * Loop over different functions (IFUNC) of SGGESX, types (PRTYPE)
495 * of test matrices, different size (M+N)
496 *
497  prtype = 0
498  qba = 3
499  qbb = 4
500  weight = sqrt( ulp )
501 *
502  DO 60 ifunc = 0, 3
503  DO 50 prtype = 1, 5
504  DO 40 m = 1, nsize - 1
505  DO 30 n = 1, nsize - m
506 *
507  weight = one / weight
508  mplusn = m + n
509 *
510 * Generate test matrices
511 *
512  fs = .true.
513  k = 0
514 *
515  CALL slaset( 'Full', mplusn, mplusn, zero, zero, ai,
516  \$ lda )
517  CALL slaset( 'Full', mplusn, mplusn, zero, zero, bi,
518  \$ lda )
519 *
520  CALL slatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
521  \$ lda, ai( 1, m+1 ), lda, bi, lda,
522  \$ bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
523  \$ q, lda, z, lda, weight, qba, qbb )
524 *
525 * Compute the Schur factorization and swapping the
526 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
527 * Swapping is accomplished via the function SLCTSX
528 * which is supplied below.
529 *
530  IF( ifunc.EQ.0 ) THEN
531  sense = 'N'
532  ELSE IF( ifunc.EQ.1 ) THEN
533  sense = 'E'
534  ELSE IF( ifunc.EQ.2 ) THEN
535  sense = 'V'
536  ELSE IF( ifunc.EQ.3 ) THEN
537  sense = 'B'
538  END IF
539 *
540  CALL slacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
541  CALL slacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
542 *
543  CALL sggesx( 'V', 'V', 'S', slctsx, sense, mplusn, ai,
544  \$ lda, bi, lda, mm, alphar, alphai, beta,
545  \$ q, lda, z, lda, pl, difest, work, lwork,
546  \$ iwork, liwork, bwork, linfo )
547 *
548  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
549  result( 1 ) = ulpinv
550  WRITE( nout, fmt = 9999 )'SGGESX', linfo, mplusn,
551  \$ prtype
552  info = linfo
553  GO TO 30
554  END IF
555 *
556 * Compute the norm(A, B)
557 *
558  CALL slacpy( 'Full', mplusn, mplusn, ai, lda, work,
559  \$ mplusn )
560  CALL slacpy( 'Full', mplusn, mplusn, bi, lda,
561  \$ work( mplusn*mplusn+1 ), mplusn )
562  abnrm = slange( 'Fro', mplusn, 2*mplusn, work, mplusn,
563  \$ work )
564 *
565 * Do tests (1) to (4)
566 *
567  CALL sget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
568  \$ lda, work, result( 1 ) )
569  CALL sget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
570  \$ lda, work, result( 2 ) )
571  CALL sget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
572  \$ lda, work, result( 3 ) )
573  CALL sget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
574  \$ lda, work, 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  IF( alphai( j ).EQ.zero ) THEN
587  temp2 = ( abs( alphar( j )-ai( j, j ) ) /
588  \$ max( smlnum, abs( alphar( j ) ),
589  \$ abs( ai( j, j ) ) )+
590  \$ abs( beta( j )-bi( j, j ) ) /
591  \$ max( smlnum, abs( beta( j ) ),
592  \$ abs( bi( j, j ) ) ) ) / ulp
593  IF( j.LT.mplusn ) THEN
594  IF( ai( j+1, j ).NE.zero ) THEN
595  ilabad = .true.
596  result( 5 ) = ulpinv
597  END IF
598  END IF
599  IF( j.GT.1 ) THEN
600  IF( ai( j, j-1 ).NE.zero ) THEN
601  ilabad = .true.
602  result( 5 ) = ulpinv
603  END IF
604  END IF
605  ELSE
606  IF( alphai( j ).GT.zero ) THEN
607  i1 = j
608  ELSE
609  i1 = j - 1
610  END IF
611  IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
612  ilabad = .true.
613  ELSE IF( i1.LT.mplusn-1 ) THEN
614  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
615  ilabad = .true.
616  result( 5 ) = ulpinv
617  END IF
618  ELSE IF( i1.GT.1 ) THEN
619  IF( ai( i1, i1-1 ).NE.zero ) THEN
620  ilabad = .true.
621  result( 5 ) = ulpinv
622  END IF
623  END IF
624  IF( .NOT.ilabad ) THEN
625  CALL sget53( ai( i1, i1 ), lda, bi( i1, i1 ),
626  \$ lda, beta( j ), alphar( j ),
627  \$ alphai( j ), temp2, iinfo )
628  IF( iinfo.GE.3 ) THEN
629  WRITE( nout, fmt = 9997 )iinfo, j,
630  \$ mplusn, prtype
631  info = abs( iinfo )
632  END IF
633  ELSE
634  temp2 = ulpinv
635  END IF
636  END IF
637  temp1 = max( temp1, temp2 )
638  IF( ilabad ) THEN
639  WRITE( nout, fmt = 9996 )j, mplusn, prtype
640  END IF
641  10 CONTINUE
642  result( 6 ) = temp1
643  ntest = ntest + 2
644 *
645 * Test (7) (if sorting worked)
646 *
647  result( 7 ) = zero
648  IF( linfo.EQ.mplusn+3 ) THEN
649  result( 7 ) = ulpinv
650  ELSE IF( mm.NE.n ) THEN
651  result( 7 ) = ulpinv
652  END IF
653  ntest = ntest + 1
654 *
655 * Test (8): compare the estimated value DIF and its
656 * value. first, compute the exact DIF.
657 *
658  result( 8 ) = zero
659  mn2 = mm*( mplusn-mm )*2
660  IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
661 *
662 * Note: for either following two causes, there are
663 * almost same number of test cases fail the test.
664 *
665  CALL slakf2( mm, mplusn-mm, ai, lda,
666  \$ ai( mm+1, mm+1 ), bi,
667  \$ bi( mm+1, mm+1 ), c, ldc )
668 *
669  CALL sgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
670  \$ 1, work( 2 ), 1, work( 3 ), lwork-2,
671  \$ info )
672  diftru = s( mn2 )
673 *
674  IF( difest( 2 ).EQ.zero ) THEN
675  IF( diftru.GT.abnrm*ulp )
676  \$ result( 8 ) = ulpinv
677  ELSE IF( diftru.EQ.zero ) THEN
678  IF( difest( 2 ).GT.abnrm*ulp )
679  \$ result( 8 ) = ulpinv
680  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
681  \$ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
682  result( 8 ) = max( diftru / difest( 2 ),
683  \$ difest( 2 ) / diftru )
684  END IF
685  ntest = ntest + 1
686  END IF
687 *
688 * Test (9)
689 *
690  result( 9 ) = zero
691  IF( linfo.EQ.( mplusn+2 ) ) THEN
692  IF( diftru.GT.abnrm*ulp )
693  \$ result( 9 ) = ulpinv
694  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
695  \$ result( 9 ) = ulpinv
696  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
697  \$ result( 9 ) = ulpinv
698  ntest = ntest + 1
699  END IF
700 *
701  ntestt = ntestt + ntest
702 *
703 * Print out tests which fail.
704 *
705  DO 20 j = 1, 9
706  IF( result( j ).GE.thresh ) THEN
707 *
708 * If this is the first test to fail,
709 * print a header to the data file.
710 *
711  IF( nerrs.EQ.0 ) THEN
712  WRITE( nout, fmt = 9995 )'SGX'
713 *
714 * Matrix types
715 *
716  WRITE( nout, fmt = 9993 )
717 *
718 * Tests performed
719 *
720  WRITE( nout, fmt = 9992 )'orthogonal', '''',
721  \$ 'transpose', ( '''', i = 1, 4 )
722 *
723  END IF
724  nerrs = nerrs + 1
725  IF( result( j ).LT.10000.0 ) THEN
726  WRITE( nout, fmt = 9991 )mplusn, prtype,
727  \$ weight, m, j, result( j )
728  ELSE
729  WRITE( nout, fmt = 9990 )mplusn, prtype,
730  \$ weight, m, j, result( j )
731  END IF
732  END IF
733  20 CONTINUE
734 *
735  30 CONTINUE
736  40 CONTINUE
737  50 CONTINUE
738  60 CONTINUE
739 *
740  GO TO 150
741 *
742  70 CONTINUE
743 *
744 * Read in data from file to check accuracy of condition estimation
745 * Read input data until N=0
746 *
747  nptknt = 0
748 *
749  80 CONTINUE
750  READ( nin, fmt = *, END = 140 )mplusn
751  IF( mplusn.EQ.0 )
752  \$ GO TO 140
753  READ( nin, fmt = *, END = 140 )n
754  DO 90 i = 1, mplusn
755  READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
756  90 CONTINUE
757  DO 100 i = 1, mplusn
758  READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
759  100 CONTINUE
760  READ( nin, fmt = * )pltru, diftru
761 *
762  nptknt = nptknt + 1
763  fs = .true.
764  k = 0
765  m = mplusn - n
766 *
767  CALL slacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
768  CALL slacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
769 *
770 * Compute the Schur factorization while swapping the
771 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
772 *
773  CALL sggesx( 'V', 'V', 'S', slctsx, 'B', mplusn, ai, lda, bi, lda,
774  \$ mm, alphar, alphai, beta, q, lda, z, lda, pl, difest,
775  \$ work, lwork, iwork, liwork, bwork, linfo )
776 *
777  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
778  result( 1 ) = ulpinv
779  WRITE( nout, fmt = 9998 )'SGGESX', linfo, mplusn, nptknt
780  GO TO 130
781  END IF
782 *
783 * Compute the norm(A, B)
784 * (should this be norm of (A,B) or (AI,BI)?)
785 *
786  CALL slacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
787  CALL slacpy( 'Full', mplusn, mplusn, bi, lda,
788  \$ work( mplusn*mplusn+1 ), mplusn )
789  abnrm = slange( 'Fro', mplusn, 2*mplusn, work, mplusn, work )
790 *
791 * Do tests (1) to (4)
792 *
793  CALL sget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
794  \$ result( 1 ) )
795  CALL sget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
796  \$ result( 2 ) )
797  CALL sget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
798  \$ result( 3 ) )
799  CALL sget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
800  \$ result( 4 ) )
801 *
802 * Do tests (5) and (6): check Schur form of A and compare
803 * eigenvalues with diagonals.
804 *
805  ntest = 6
806  temp1 = zero
807  result( 5 ) = zero
808  result( 6 ) = zero
809 *
810  DO 110 j = 1, mplusn
811  ilabad = .false.
812  IF( alphai( j ).EQ.zero ) THEN
813  temp2 = ( abs( alphar( j )-ai( j, j ) ) /
814  \$ max( smlnum, abs( alphar( j ) ), abs( ai( j,
815  \$ j ) ) )+abs( beta( j )-bi( j, j ) ) /
816  \$ max( smlnum, abs( beta( j ) ), abs( bi( j, j ) ) ) )
817  \$ / ulp
818  IF( j.LT.mplusn ) THEN
819  IF( ai( j+1, j ).NE.zero ) THEN
820  ilabad = .true.
821  result( 5 ) = ulpinv
822  END IF
823  END IF
824  IF( j.GT.1 ) THEN
825  IF( ai( j, j-1 ).NE.zero ) THEN
826  ilabad = .true.
827  result( 5 ) = ulpinv
828  END IF
829  END IF
830  ELSE
831  IF( alphai( j ).GT.zero ) THEN
832  i1 = j
833  ELSE
834  i1 = j - 1
835  END IF
836  IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
837  ilabad = .true.
838  ELSE IF( i1.LT.mplusn-1 ) THEN
839  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
840  ilabad = .true.
841  result( 5 ) = ulpinv
842  END IF
843  ELSE IF( i1.GT.1 ) THEN
844  IF( ai( i1, i1-1 ).NE.zero ) THEN
845  ilabad = .true.
846  result( 5 ) = ulpinv
847  END IF
848  END IF
849  IF( .NOT.ilabad ) THEN
850  CALL sget53( ai( i1, i1 ), lda, bi( i1, i1 ), lda,
851  \$ beta( j ), alphar( j ), alphai( j ), temp2,
852  \$ iinfo )
853  IF( iinfo.GE.3 ) THEN
854  WRITE( nout, fmt = 9997 )iinfo, j, mplusn, nptknt
855  info = abs( iinfo )
856  END IF
857  ELSE
858  temp2 = ulpinv
859  END IF
860  END IF
861  temp1 = max( temp1, temp2 )
862  IF( ilabad ) THEN
863  WRITE( nout, fmt = 9996 )j, mplusn, nptknt
864  END IF
865  110 CONTINUE
866  result( 6 ) = temp1
867 *
868 * Test (7) (if sorting worked) <--------- need to be checked.
869 *
870  ntest = 7
871  result( 7 ) = zero
872  IF( linfo.EQ.mplusn+3 )
873  \$ result( 7 ) = ulpinv
874 *
875 * Test (8): compare the estimated value of DIF and its true value.
876 *
877  ntest = 8
878  result( 8 ) = zero
879  IF( difest( 2 ).EQ.zero ) THEN
880  IF( diftru.GT.abnrm*ulp )
881  \$ result( 8 ) = ulpinv
882  ELSE IF( diftru.EQ.zero ) THEN
883  IF( difest( 2 ).GT.abnrm*ulp )
884  \$ result( 8 ) = ulpinv
885  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
886  \$ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
887  result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
888  END IF
889 *
890 * Test (9)
891 *
892  ntest = 9
893  result( 9 ) = zero
894  IF( linfo.EQ.( mplusn+2 ) ) THEN
895  IF( diftru.GT.abnrm*ulp )
896  \$ result( 9 ) = ulpinv
897  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
898  \$ result( 9 ) = ulpinv
899  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
900  \$ result( 9 ) = ulpinv
901  END IF
902 *
903 * Test (10): compare the estimated value of PL and it true value.
904 *
905  ntest = 10
906  result( 10 ) = zero
907  IF( pl( 1 ).EQ.zero ) THEN
908  IF( pltru.GT.abnrm*ulp )
909  \$ result( 10 ) = ulpinv
910  ELSE IF( pltru.EQ.zero ) THEN
911  IF( pl( 1 ).GT.abnrm*ulp )
912  \$ result( 10 ) = ulpinv
913  ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
914  \$ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
915  result( 10 ) = ulpinv
916  END IF
917 *
918  ntestt = ntestt + ntest
919 *
920 * Print out tests which fail.
921 *
922  DO 120 j = 1, ntest
923  IF( result( j ).GE.thresh ) THEN
924 *
925 * If this is the first test to fail,
926 * print a header to the data file.
927 *
928  IF( nerrs.EQ.0 ) THEN
929  WRITE( nout, fmt = 9995 )'SGX'
930 *
931 * Matrix types
932 *
933  WRITE( nout, fmt = 9994 )
934 *
935 * Tests performed
936 *
937  WRITE( nout, fmt = 9992 )'orthogonal', '''',
938  \$ 'transpose', ( '''', i = 1, 4 )
939 *
940  END IF
941  nerrs = nerrs + 1
942  IF( result( j ).LT.10000.0 ) THEN
943  WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
944  ELSE
945  WRITE( nout, fmt = 9988 )nptknt, mplusn, j, result( j )
946  END IF
947  END IF
948 *
949  120 CONTINUE
950 *
951  130 CONTINUE
952  GO TO 80
953  140 CONTINUE
954 *
955  150 CONTINUE
956 *
957 * Summary
958 *
959  CALL alasvm( 'SGX', nout, nerrs, ntestt, 0 )
960 *
961  work( 1 ) = maxwrk
962 *
963  RETURN
964 *
965  9999 FORMAT( ' SDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
966  \$ i6, ', JTYPE=', i6, ')' )
967 *
968  9998 FORMAT( ' SDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
969  \$ i6, ', Input Example #', i2, ')' )
970 *
971  9997 FORMAT( ' SDRGSX: SGET53 returned INFO=', i1, ' for eigenvalue ',
972  \$ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
973 *
974  9996 FORMAT( ' SDRGSX: S not in Schur form at eigenvalue ', i6, '.',
975  \$ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
976 *
977  9995 FORMAT( / 1x, a3, ' -- Real Expert Generalized Schur form',
978  \$ ' problem driver' )
979 *
980  9994 FORMAT( 'Input Example' )
981 *
982  9993 FORMAT( ' Matrix types: ', /
983  \$ ' 1: A is a block diagonal matrix of Jordan blocks ',
984  \$ 'and B is the identity ', / ' matrix, ',
985  \$ / ' 2: A and B are upper triangular matrices, ',
986  \$ / ' 3: A and B are as type 2, but each second diagonal ',
987  \$ 'block in A_11 and ', /
988  \$ ' each third diaongal block in A_22 are 2x2 blocks,',
989  \$ / ' 4: A and B are block diagonal matrices, ',
990  \$ / ' 5: (A,B) has potentially close or common ',
991  \$ 'eigenvalues.', / )
992 *
993  9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
994  \$ 'Q and Z are ', a, ',', / 19x,
995  \$ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
996  \$ / ' 1 = | A - Q S Z', a,
997  \$ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
998  \$ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
999  \$ ' | / ( n ulp ) 4 = | I - ZZ', a,
1000  \$ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
1001  \$ 'Schur form S', / ' 6 = difference between (alpha,beta)',
1002  \$ ' and diagonals of (S,T)', /
1003  \$ ' 7 = 1/ULP if SDIM is not the correct number of ',
1004  \$ 'selected eigenvalues', /
1005  \$ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
1006  \$ 'DIFTRU/DIFEST > 10*THRESH',
1007  \$ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
1008  \$ 'when reordering fails', /
1009  \$ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
1010  \$ 'PLTRU/PLEST > THRESH', /
1011  \$ ' ( Test 10 is only for input examples )', / )
1012  9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
1013  \$ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
1014  9990 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
1015  \$ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, e10.3 )
1016  9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1017  \$ ' result ', i2, ' is', 0p, f8.2 )
1018  9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1019  \$ ' result ', i2, ' is', 1p, e10.3 )
1020 *
1021 * End of SDRGSX
1022 *
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:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
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 slakf2(M, N, A, LDA, B, D, E, Z, LDZ)
SLAKF2
Definition: slakf2.f:105
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:268
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:114
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:365
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:211
subroutine sget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
SGET53
Definition: sget53.f:126
logical function slctsx(AR, AI, BETA)
SLCTSX
Definition: slctsx.f:65
subroutine sget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
SGET51
Definition: sget51.f:149
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: