LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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 ..
385 LOGICAL ILABAD
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,
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 thrsh2 = ten*thresh
482 ntestt = 0
483 nerrs = 0
484*
485* Go to the tests for read-in matrix pairs
486*
487 ifunc = 0
488 IF( nsize.EQ.0 )
489 $ GO TO 70
490*
491* Test the built-in matrix pairs.
492* Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
493* of test matrices, different size (M+N)
494*
495 prtype = 0
496 qba = 3
497 qbb = 4
498 weight = sqrt( ulp )
499*
500 DO 60 ifunc = 0, 3
501 DO 50 prtype = 1, 5
502 DO 40 m = 1, nsize - 1
503 DO 30 n = 1, nsize - m
504*
505 weight = one / weight
506 mplusn = m + n
507*
508* Generate test matrices
509*
510 fs = .true.
511 k = 0
512*
513 CALL dlaset( 'Full', mplusn, mplusn, zero, zero, ai,
514 $ lda )
515 CALL dlaset( 'Full', mplusn, mplusn, zero, zero, bi,
516 $ lda )
517*
518 CALL dlatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
519 $ lda, ai( 1, m+1 ), lda, bi, lda,
520 $ bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
521 $ q, lda, z, lda, weight, qba, qbb )
522*
523* Compute the Schur factorization and swapping the
524* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
525* Swapping is accomplished via the function DLCTSX
526* which is supplied below.
527*
528 IF( ifunc.EQ.0 ) THEN
529 sense = 'N'
530 ELSE IF( ifunc.EQ.1 ) THEN
531 sense = 'E'
532 ELSE IF( ifunc.EQ.2 ) THEN
533 sense = 'V'
534 ELSE IF( ifunc.EQ.3 ) THEN
535 sense = 'B'
536 END IF
537*
538 CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
539 CALL dlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
540*
541 CALL dggesx( 'V', 'V', 'S', dlctsx, sense, mplusn, ai,
542 $ lda, bi, lda, mm, alphar, alphai, beta,
543 $ q, lda, z, lda, pl, difest, work, lwork,
544 $ iwork, liwork, bwork, linfo )
545*
546 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
547 result( 1 ) = ulpinv
548 WRITE( nout, fmt = 9999 )'DGGESX', linfo, mplusn,
549 $ prtype
550 info = linfo
551 GO TO 30
552 END IF
553*
554* Compute the norm(A, B)
555*
556 CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, work,
557 $ mplusn )
558 CALL dlacpy( 'Full', mplusn, mplusn, bi, lda,
559 $ work( mplusn*mplusn+1 ), mplusn )
560 abnrm = dlange( 'Fro', mplusn, 2*mplusn, work, mplusn,
561 $ work )
562*
563* Do tests (1) to (4)
564*
565 CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
566 $ lda, work, result( 1 ) )
567 CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
568 $ lda, work, result( 2 ) )
569 CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
570 $ lda, work, result( 3 ) )
571 CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
572 $ lda, work, result( 4 ) )
573 ntest = 4
574*
575* Do tests (5) and (6): check Schur form of A and
576* compare eigenvalues with diagonals.
577*
578 temp1 = zero
579 result( 5 ) = zero
580 result( 6 ) = zero
581*
582 DO 10 j = 1, mplusn
583 ilabad = .false.
584 IF( alphai( j ).EQ.zero ) THEN
585 temp2 = ( abs( alphar( j )-ai( j, j ) ) /
586 $ max( smlnum, abs( alphar( j ) ),
587 $ abs( ai( j, j ) ) )+
588 $ abs( beta( j )-bi( j, j ) ) /
589 $ max( smlnum, abs( beta( j ) ),
590 $ abs( bi( j, j ) ) ) ) / ulp
591 IF( j.LT.mplusn ) THEN
592 IF( ai( j+1, j ).NE.zero ) THEN
593 ilabad = .true.
594 result( 5 ) = ulpinv
595 END IF
596 END IF
597 IF( j.GT.1 ) THEN
598 IF( ai( j, j-1 ).NE.zero ) THEN
599 ilabad = .true.
600 result( 5 ) = ulpinv
601 END IF
602 END IF
603 ELSE
604 IF( alphai( j ).GT.zero ) THEN
605 i1 = j
606 ELSE
607 i1 = j - 1
608 END IF
609 IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
610 ilabad = .true.
611 ELSE IF( i1.LT.mplusn-1 ) THEN
612 IF( ai( i1+2, i1+1 ).NE.zero ) THEN
613 ilabad = .true.
614 result( 5 ) = ulpinv
615 END IF
616 ELSE IF( i1.GT.1 ) THEN
617 IF( ai( i1, i1-1 ).NE.zero ) THEN
618 ilabad = .true.
619 result( 5 ) = ulpinv
620 END IF
621 END IF
622 IF( .NOT.ilabad ) THEN
623 CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ),
624 $ lda, beta( j ), alphar( j ),
625 $ alphai( j ), temp2, iinfo )
626 IF( iinfo.GE.3 ) THEN
627 WRITE( nout, fmt = 9997 )iinfo, j,
628 $ mplusn, prtype
629 info = abs( iinfo )
630 END IF
631 ELSE
632 temp2 = ulpinv
633 END IF
634 END IF
635 temp1 = max( temp1, temp2 )
636 IF( ilabad ) THEN
637 WRITE( nout, fmt = 9996 )j, mplusn, prtype
638 END IF
639 10 CONTINUE
640 result( 6 ) = temp1
641 ntest = ntest + 2
642*
643* Test (7) (if sorting worked)
644*
645 result( 7 ) = zero
646 IF( linfo.EQ.mplusn+3 ) THEN
647 result( 7 ) = ulpinv
648 ELSE IF( mm.NE.n ) THEN
649 result( 7 ) = ulpinv
650 END IF
651 ntest = ntest + 1
652*
653* Test (8): compare the estimated value DIF and its
654* value. first, compute the exact DIF.
655*
656 result( 8 ) = zero
657 mn2 = mm*( mplusn-mm )*2
658 IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
659*
660* Note: for either following two causes, there are
661* almost same number of test cases fail the test.
662*
663 CALL dlakf2( mm, mplusn-mm, ai, lda,
664 $ ai( mm+1, mm+1 ), bi,
665 $ bi( mm+1, mm+1 ), c, ldc )
666*
667 CALL dgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
668 $ 1, work( 2 ), 1, work( 3 ), lwork-2,
669 $ info )
670 diftru = s( mn2 )
671*
672 IF( difest( 2 ).EQ.zero ) THEN
673 IF( diftru.GT.abnrm*ulp )
674 $ result( 8 ) = ulpinv
675 ELSE IF( diftru.EQ.zero ) THEN
676 IF( difest( 2 ).GT.abnrm*ulp )
677 $ result( 8 ) = ulpinv
678 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
679 $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
680 result( 8 ) = max( diftru / difest( 2 ),
681 $ difest( 2 ) / diftru )
682 END IF
683 ntest = ntest + 1
684 END IF
685*
686* Test (9)
687*
688 result( 9 ) = zero
689 IF( linfo.EQ.( mplusn+2 ) ) THEN
690 IF( diftru.GT.abnrm*ulp )
691 $ result( 9 ) = ulpinv
692 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
693 $ result( 9 ) = ulpinv
694 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
695 $ result( 9 ) = ulpinv
696 ntest = ntest + 1
697 END IF
698*
699 ntestt = ntestt + ntest
700*
701* Print out tests which fail.
702*
703 DO 20 j = 1, 9
704 IF( result( j ).GE.thresh ) THEN
705*
706* If this is the first test to fail,
707* print a header to the data file.
708*
709 IF( nerrs.EQ.0 ) THEN
710 WRITE( nout, fmt = 9995 )'DGX'
711*
712* Matrix types
713*
714 WRITE( nout, fmt = 9993 )
715*
716* Tests performed
717*
718 WRITE( nout, fmt = 9992 )'orthogonal', '''',
719 $ 'transpose', ( '''', i = 1, 4 )
720*
721 END IF
722 nerrs = nerrs + 1
723 IF( result( j ).LT.10000.0d0 ) THEN
724 WRITE( nout, fmt = 9991 )mplusn, prtype,
725 $ weight, m, j, result( j )
726 ELSE
727 WRITE( nout, fmt = 9990 )mplusn, prtype,
728 $ weight, m, j, result( j )
729 END IF
730 END IF
731 20 CONTINUE
732*
733 30 CONTINUE
734 40 CONTINUE
735 50 CONTINUE
736 60 CONTINUE
737*
738 GO TO 150
739*
740 70 CONTINUE
741*
742* Read in data from file to check accuracy of condition estimation
743* Read input data until N=0
744*
745 nptknt = 0
746*
747 80 CONTINUE
748 READ( nin, fmt = *, END = 140 )mplusn
749 IF( mplusn.EQ.0 )
750 $ GO TO 140
751 READ( nin, fmt = *, END = 140 )n
752 DO 90 i = 1, mplusn
753 READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
754 90 CONTINUE
755 DO 100 i = 1, mplusn
756 READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
757 100 CONTINUE
758 READ( nin, fmt = * )pltru, diftru
759*
760 nptknt = nptknt + 1
761 fs = .true.
762 k = 0
763 m = mplusn - n
764*
765 CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
766 CALL dlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
767*
768* Compute the Schur factorization while swapping the
769* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
770*
771 CALL dggesx( 'V', 'V', 'S', dlctsx, 'B', mplusn, ai, lda, bi, lda,
772 $ mm, alphar, alphai, beta, q, lda, z, lda, pl, difest,
773 $ work, lwork, iwork, liwork, bwork, linfo )
774*
775 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
776 result( 1 ) = ulpinv
777 WRITE( nout, fmt = 9998 )'DGGESX', linfo, mplusn, nptknt
778 GO TO 130
779 END IF
780*
781* Compute the norm(A, B)
782* (should this be norm of (A,B) or (AI,BI)?)
783*
784 CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
785 CALL dlacpy( 'Full', mplusn, mplusn, bi, lda,
786 $ work( mplusn*mplusn+1 ), mplusn )
787 abnrm = dlange( 'Fro', mplusn, 2*mplusn, work, mplusn, work )
788*
789* Do tests (1) to (4)
790*
791 CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
792 $ result( 1 ) )
793 CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
794 $ result( 2 ) )
795 CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
796 $ result( 3 ) )
797 CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
798 $ result( 4 ) )
799*
800* Do tests (5) and (6): check Schur form of A and compare
801* eigenvalues with diagonals.
802*
803 ntest = 6
804 temp1 = zero
805 result( 5 ) = zero
806 result( 6 ) = zero
807*
808 DO 110 j = 1, mplusn
809 ilabad = .false.
810 IF( alphai( j ).EQ.zero ) THEN
811 temp2 = ( abs( alphar( j )-ai( j, j ) ) /
812 $ max( smlnum, abs( alphar( j ) ), abs( ai( j,
813 $ j ) ) )+abs( beta( j )-bi( j, j ) ) /
814 $ max( smlnum, abs( beta( j ) ), abs( bi( j, j ) ) ) )
815 $ / ulp
816 IF( j.LT.mplusn ) THEN
817 IF( ai( j+1, j ).NE.zero ) THEN
818 ilabad = .true.
819 result( 5 ) = ulpinv
820 END IF
821 END IF
822 IF( j.GT.1 ) THEN
823 IF( ai( j, j-1 ).NE.zero ) THEN
824 ilabad = .true.
825 result( 5 ) = ulpinv
826 END IF
827 END IF
828 ELSE
829 IF( alphai( j ).GT.zero ) THEN
830 i1 = j
831 ELSE
832 i1 = j - 1
833 END IF
834 IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
835 ilabad = .true.
836 ELSE IF( i1.LT.mplusn-1 ) THEN
837 IF( ai( i1+2, i1+1 ).NE.zero ) THEN
838 ilabad = .true.
839 result( 5 ) = ulpinv
840 END IF
841 ELSE IF( i1.GT.1 ) THEN
842 IF( ai( i1, i1-1 ).NE.zero ) THEN
843 ilabad = .true.
844 result( 5 ) = ulpinv
845 END IF
846 END IF
847 IF( .NOT.ilabad ) THEN
848 CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ), lda,
849 $ beta( j ), alphar( j ), alphai( j ), temp2,
850 $ iinfo )
851 IF( iinfo.GE.3 ) THEN
852 WRITE( nout, fmt = 9997 )iinfo, j, mplusn, nptknt
853 info = abs( iinfo )
854 END IF
855 ELSE
856 temp2 = ulpinv
857 END IF
858 END IF
859 temp1 = max( temp1, temp2 )
860 IF( ilabad ) THEN
861 WRITE( nout, fmt = 9996 )j, mplusn, nptknt
862 END IF
863 110 CONTINUE
864 result( 6 ) = temp1
865*
866* Test (7) (if sorting worked) <--------- need to be checked.
867*
868 ntest = 7
869 result( 7 ) = zero
870 IF( linfo.EQ.mplusn+3 )
871 $ result( 7 ) = ulpinv
872*
873* Test (8): compare the estimated value of DIF and its true value.
874*
875 ntest = 8
876 result( 8 ) = zero
877 IF( difest( 2 ).EQ.zero ) THEN
878 IF( diftru.GT.abnrm*ulp )
879 $ result( 8 ) = ulpinv
880 ELSE IF( diftru.EQ.zero ) THEN
881 IF( difest( 2 ).GT.abnrm*ulp )
882 $ result( 8 ) = ulpinv
883 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
884 $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
885 result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
886 END IF
887*
888* Test (9)
889*
890 ntest = 9
891 result( 9 ) = zero
892 IF( linfo.EQ.( mplusn+2 ) ) THEN
893 IF( diftru.GT.abnrm*ulp )
894 $ result( 9 ) = ulpinv
895 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
896 $ result( 9 ) = ulpinv
897 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
898 $ result( 9 ) = ulpinv
899 END IF
900*
901* Test (10): compare the estimated value of PL and it true value.
902*
903 ntest = 10
904 result( 10 ) = zero
905 IF( pl( 1 ).EQ.zero ) THEN
906 IF( pltru.GT.abnrm*ulp )
907 $ result( 10 ) = ulpinv
908 ELSE IF( pltru.EQ.zero ) THEN
909 IF( pl( 1 ).GT.abnrm*ulp )
910 $ result( 10 ) = ulpinv
911 ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
912 $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
913 result( 10 ) = ulpinv
914 END IF
915*
916 ntestt = ntestt + ntest
917*
918* Print out tests which fail.
919*
920 DO 120 j = 1, ntest
921 IF( result( j ).GE.thresh ) THEN
922*
923* If this is the first test to fail,
924* print a header to the data file.
925*
926 IF( nerrs.EQ.0 ) THEN
927 WRITE( nout, fmt = 9995 )'DGX'
928*
929* Matrix types
930*
931 WRITE( nout, fmt = 9994 )
932*
933* Tests performed
934*
935 WRITE( nout, fmt = 9992 )'orthogonal', '''',
936 $ 'transpose', ( '''', i = 1, 4 )
937*
938 END IF
939 nerrs = nerrs + 1
940 IF( result( j ).LT.10000.0d0 ) THEN
941 WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
942 ELSE
943 WRITE( nout, fmt = 9988 )nptknt, mplusn, j, result( j )
944 END IF
945 END IF
946*
947 120 CONTINUE
948*
949 130 CONTINUE
950 GO TO 80
951 140 CONTINUE
952*
953 150 CONTINUE
954*
955* Summary
956*
957 CALL alasvm( 'DGX', nout, nerrs, ntestt, 0 )
958*
959 work( 1 ) = maxwrk
960*
961 RETURN
962*
963 9999 FORMAT( ' DDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
964 $ i6, ', JTYPE=', i6, ')' )
965*
966 9998 FORMAT( ' DDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
967 $ i6, ', Input Example #', i2, ')' )
968*
969 9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', i1, ' for eigenvalue ',
970 $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
971*
972 9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', i6, '.',
973 $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
974*
975 9995 FORMAT( / 1x, a3, ' -- Real Expert Generalized Schur form',
976 $ ' problem driver' )
977*
978 9994 FORMAT( 'Input Example' )
979*
980 9993 FORMAT( ' Matrix types: ', /
981 $ ' 1: A is a block diagonal matrix of Jordan blocks ',
982 $ 'and B is the identity ', / ' matrix, ',
983 $ / ' 2: A and B are upper triangular matrices, ',
984 $ / ' 3: A and B are as type 2, but each second diagonal ',
985 $ 'block in A_11 and ', /
986 $ ' each third diagonal block in A_22 are 2x2 blocks,',
987 $ / ' 4: A and B are block diagonal matrices, ',
988 $ / ' 5: (A,B) has potentially close or common ',
989 $ 'eigenvalues.', / )
990*
991 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
992 $ 'Q and Z are ', a, ',', / 19x,
993 $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
994 $ / ' 1 = | A - Q S Z', a,
995 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
996 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
997 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
998 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
999 $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
1000 $ ' and diagonals of (S,T)', /
1001 $ ' 7 = 1/ULP if SDIM is not the correct number of ',
1002 $ 'selected eigenvalues', /
1003 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
1004 $ 'DIFTRU/DIFEST > 10*THRESH',
1005 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
1006 $ 'when reordering fails', /
1007 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
1008 $ 'PLTRU/PLEST > THRESH', /
1009 $ ' ( Test 10 is only for input examples )', / )
1010 9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
1011 $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
1012 9990 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
1013 $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, d10.3 )
1014 9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1015 $ ' result ', i2, ' is', 0p, f8.2 )
1016 9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1017 $ ' result ', i2, ' is', 1p, d10.3 )
1018*
1019* End of DDRGSX
1020*
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dget51(itype, n, a, lda, b, ldb, u, ldu, v, ldv, work, result)
DGET51
Definition dget51.f:149
subroutine dget53(a, lda, b, ldb, scale, wr, wi, result, info)
DGET53
Definition dget53.f:126
subroutine dlakf2(m, n, a, lda, b, d, e, z, ldz)
DLAKF2
Definition dlakf2.f:105
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
logical function dlctsx(ar, ai, beta)
DLCTSX
Definition dlctsx.f:65
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
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
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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 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
Here is the call graph for this function:
Here is the caller graph for this function: