 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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.```

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: