LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ cdrgsx()

 subroutine cdrgsx ( integer NSIZE, integer NCMAX, real THRESH, integer NIN, integer NOUT, complex, dimension( lda, * ) A, integer LDA, complex, dimension( lda, * ) B, complex, dimension( lda, * ) AI, complex, dimension( lda, * ) BI, complex, dimension( lda, * ) Z, complex, dimension( lda, * ) Q, complex, dimension( * ) ALPHA, complex, dimension( * ) BETA, complex, dimension( ldc, * ) C, integer LDC, real, dimension( * ) S, complex, dimension( * ) WORK, integer LWORK, real, dimension( * ) RWORK, integer, dimension( * ) IWORK, integer LIWORK, logical, dimension( * ) BWORK, integer INFO )

CDRGSX

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

CGGES factors A and B as Q*S*Z'  and Q*T*Z' , where ' means conjugate
transpose, S and T are  upper triangular (i.e., in generalized Schur
form), and Q and Z are unitary. It also computes the generalized
eigenvalues (alpha(j),beta(j)), j=1,...,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 CDRGSX is called with NSIZE > 0, five (5) types of built-in
matrix pairs are used to test the routine CGGESX.

When CDRGSX is called with NSIZE = 0, it reads in test matrix data
to test CGGESX.
(need more details on what kind of read-in data are needed).

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. triangular form)

(6)   maximum over j of D(j)  where:

|alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
D(j) = ------------------------ + -----------------------
max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)

(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 CGGESX, 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 same 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 CGESVD) 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 CLAKF2``` [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 INFO not equal to 0.)``` [out] A ``` A is COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (LDA, NSIZE) Copy of A, modified by CGGESX.``` [out] BI ``` BI is COMPLEX array, dimension (LDA, NSIZE) Copy of B, modified by CGGESX.``` [out] Z ``` Z is COMPLEX array, dimension (LDA, NSIZE) Z holds the left Schur vectors computed by CGGESX.``` [out] Q ``` Q is COMPLEX array, dimension (LDA, NSIZE) Q holds the right Schur vectors computed by CGGESX.``` [out] ALPHA ` ALPHA is COMPLEX array, dimension (NSIZE)` [out] BETA ``` BETA is COMPLEX array, dimension (NSIZE) On exit, ALPHA/BETA are the eigenvalues.``` [out] C ``` C is COMPLEX array, dimension (LDC, LDC) Store the matrix generated by subroutine CLAKF2, 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 COMPLEX array, dimension (LWORK)` [in] LWORK ``` LWORK is INTEGER The dimension of the array WORK. LWORK >= 3*NSIZE*NSIZE/2``` [out] RWORK ``` RWORK is REAL array, dimension (5*NSIZE*NSIZE/2 - 4)``` [out] IWORK ` IWORK is INTEGER array, dimension (LIWORK)` [in] LIWORK ``` LIWORK is INTEGER The dimension of the array IWORK. LIWORK >= NSIZE + 2.``` [out] BWORK ` BWORK is LOGICAL array, dimension (NSIZE)` [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 346 of file cdrgsx.f.

349*
350* -- LAPACK test routine --
351* -- LAPACK is a software package provided by Univ. of Tennessee, --
352* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
353*
354* .. Scalar Arguments ..
355 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
356 \$ NOUT, NSIZE
357 REAL THRESH
358* ..
359* .. Array Arguments ..
360 LOGICAL BWORK( * )
361 INTEGER IWORK( * )
362 REAL RWORK( * ), S( * )
363 COMPLEX A( LDA, * ), AI( LDA, * ), ALPHA( * ),
364 \$ B( LDA, * ), BETA( * ), BI( LDA, * ),
365 \$ C( LDC, * ), Q( LDA, * ), WORK( * ),
366 \$ Z( LDA, * )
367* ..
368*
369* =====================================================================
370*
371* .. Parameters ..
372 REAL ZERO, ONE, TEN
373 parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1 )
374 COMPLEX CZERO
375 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
376* ..
377* .. Local Scalars ..
379 CHARACTER SENSE
380 INTEGER BDSPAC, I, IFUNC, J, LINFO, MAXWRK, MINWRK, MM,
381 \$ MN2, NERRS, NPTKNT, NTEST, NTESTT, PRTYPE, QBA,
382 \$ QBB
383 REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
384 \$ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
385 COMPLEX X
386* ..
387* .. Local Arrays ..
388 REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 )
389* ..
390* .. External Functions ..
391 LOGICAL CLCTSX
392 INTEGER ILAENV
393 REAL CLANGE, SLAMCH
394 EXTERNAL clctsx, ilaenv, clange, slamch
395* ..
396* .. External Subroutines ..
397 EXTERNAL alasvm, cgesvd, cget51, cggesx, clacpy, clakf2,
399* ..
400* .. Scalars in Common ..
401 LOGICAL FS
402 INTEGER K, M, MPLUSN, N
403* ..
404* .. Common blocks ..
405 COMMON / mn / m, n, mplusn, k, fs
406* ..
407* .. Intrinsic Functions ..
408 INTRINSIC abs, aimag, max, real, sqrt
409* ..
410* .. Statement Functions ..
411 REAL ABS1
412* ..
413* .. Statement Function definitions ..
414 abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
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 = -15
432 ELSE IF( liwork.LT.nsize+2 ) 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 = 3*nsize*nsize / 2
446*
447* workspace for cggesx
448*
449 maxwrk = nsize*( 1+ilaenv( 1, 'CGEQRF', ' ', nsize, 1, nsize,
450 \$ 0 ) )
451 maxwrk = max( maxwrk, nsize*( 1+ilaenv( 1, 'CUNGQR', ' ',
452 \$ nsize, 1, nsize, -1 ) ) )
453*
454* workspace for cgesvd
455*
456 bdspac = 3*nsize*nsize / 2
457 maxwrk = max( maxwrk, nsize*nsize*
458 \$ ( 1+ilaenv( 1, 'CGEBRD', ' ', 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 = -18
469*
470 IF( info.NE.0 ) THEN
471 CALL xerbla( 'CDRGSX', -info )
472 RETURN
473 END IF
474*
475* Important constants
476*
477 ulp = slamch( 'P' )
478 ulpinv = one / ulp
479 smlnum = slamch( 'S' ) / ulp
480 bignum = one / smlnum
481 CALL slabad( 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 CGGESX, 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 claset( 'Full', mplusn, mplusn, czero, czero, ai,
515 \$ lda )
516 CALL claset( 'Full', mplusn, mplusn, czero, czero, bi,
517 \$ lda )
518*
519 CALL clatm5( 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 CLCTSX
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 clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
540 CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
541*
542 CALL cggesx( 'V', 'V', 'S', clctsx, sense, mplusn, ai,
543 \$ lda, bi, lda, mm, alpha, beta, q, lda, z,
544 \$ lda, pl, difest, work, lwork, rwork,
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 )'CGGESX', linfo, mplusn,
550 \$ prtype
551 info = linfo
552 GO TO 30
553 END IF
554*
555* Compute the norm(A, B)
556*
557 CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work,
558 \$ mplusn )
559 CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
560 \$ work( mplusn*mplusn+1 ), mplusn )
561 abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn,
562 \$ rwork )
563*
564* Do tests (1) to (4)
565*
566 result( 2 ) = zero
567 CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
568 \$ lda, work, rwork, result( 1 ) )
569 CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
570 \$ lda, work, rwork, result( 2 ) )
571 CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
572 \$ lda, work, rwork, result( 3 ) )
573 CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
574 \$ lda, work, rwork, 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
586 temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
587 \$ max( smlnum, abs1( alpha( j ) ),
588 \$ abs1( ai( j, j ) ) )+
589 \$ abs1( beta( j )-bi( j, j ) ) /
590 \$ max( smlnum, abs1( beta( j ) ),
591 \$ abs1( 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 temp1 = max( temp1, temp2 )
606 WRITE( nout, fmt = 9997 )j, mplusn, prtype
607 END IF
608 10 CONTINUE
609 result( 6 ) = temp1
610 ntest = ntest + 2
611*
612* Test (7) (if sorting worked)
613*
614 result( 7 ) = zero
615 IF( linfo.EQ.mplusn+3 ) THEN
616 result( 7 ) = ulpinv
617 ELSE IF( mm.NE.n ) THEN
618 result( 7 ) = ulpinv
619 END IF
620 ntest = ntest + 1
621*
622* Test (8): compare the estimated value DIF and its
623* value. first, compute the exact DIF.
624*
625 result( 8 ) = zero
626 mn2 = mm*( mplusn-mm )*2
627 IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
628*
629* Note: for either following two cases, there are
630* almost same number of test cases fail the test.
631*
632 CALL clakf2( mm, mplusn-mm, ai, lda,
633 \$ ai( mm+1, mm+1 ), bi,
634 \$ bi( mm+1, mm+1 ), c, ldc )
635*
636 CALL cgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
637 \$ 1, work( 2 ), 1, work( 3 ), lwork-2,
638 \$ rwork, info )
639 diftru = s( mn2 )
640*
641 IF( difest( 2 ).EQ.zero ) THEN
642 IF( diftru.GT.abnrm*ulp )
643 \$ result( 8 ) = ulpinv
644 ELSE IF( diftru.EQ.zero ) THEN
645 IF( difest( 2 ).GT.abnrm*ulp )
646 \$ result( 8 ) = ulpinv
647 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
648 \$ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
649 result( 8 ) = max( diftru / difest( 2 ),
650 \$ difest( 2 ) / diftru )
651 END IF
652 ntest = ntest + 1
653 END IF
654*
655* Test (9)
656*
657 result( 9 ) = zero
658 IF( linfo.EQ.( mplusn+2 ) ) THEN
659 IF( diftru.GT.abnrm*ulp )
660 \$ result( 9 ) = ulpinv
661 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
662 \$ result( 9 ) = ulpinv
663 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
664 \$ result( 9 ) = ulpinv
665 ntest = ntest + 1
666 END IF
667*
668 ntestt = ntestt + ntest
669*
670* Print out tests which fail.
671*
672 DO 20 j = 1, 9
673 IF( result( j ).GE.thresh ) THEN
674*
675* If this is the first test to fail,
676* print a header to the data file.
677*
678 IF( nerrs.EQ.0 ) THEN
679 WRITE( nout, fmt = 9996 )'CGX'
680*
681* Matrix types
682*
683 WRITE( nout, fmt = 9994 )
684*
685* Tests performed
686*
687 WRITE( nout, fmt = 9993 )'unitary', '''',
688 \$ 'transpose', ( '''', i = 1, 4 )
689*
690 END IF
691 nerrs = nerrs + 1
692 IF( result( j ).LT.10000.0 ) THEN
693 WRITE( nout, fmt = 9992 )mplusn, prtype,
694 \$ weight, m, j, result( j )
695 ELSE
696 WRITE( nout, fmt = 9991 )mplusn, prtype,
697 \$ weight, m, j, result( j )
698 END IF
699 END IF
700 20 CONTINUE
701*
702 30 CONTINUE
703 40 CONTINUE
704 50 CONTINUE
705 60 CONTINUE
706*
707 GO TO 150
708*
709 70 CONTINUE
710*
711* Read in data from file to check accuracy of condition estimation
712* Read input data until N=0
713*
714 nptknt = 0
715*
716 80 CONTINUE
717 READ( nin, fmt = *, END = 140 )mplusn
718 IF( mplusn.EQ.0 )
719 \$ GO TO 140
720 READ( nin, fmt = *, END = 140 )n
721 DO 90 i = 1, mplusn
722 READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
723 90 CONTINUE
724 DO 100 i = 1, mplusn
725 READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
726 100 CONTINUE
727 READ( nin, fmt = * )pltru, diftru
728*
729 nptknt = nptknt + 1
730 fs = .true.
731 k = 0
732 m = mplusn - n
733*
734 CALL clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
735 CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
736*
737* Compute the Schur factorization while swapping the
738* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
739*
740 CALL cggesx( 'V', 'V', 'S', clctsx, 'B', mplusn, ai, lda, bi, lda,
741 \$ mm, alpha, beta, q, lda, z, lda, pl, difest, work,
742 \$ lwork, rwork, iwork, liwork, bwork, linfo )
743*
744 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
745 result( 1 ) = ulpinv
746 WRITE( nout, fmt = 9998 )'CGGESX', linfo, mplusn, nptknt
747 GO TO 130
748 END IF
749*
750* Compute the norm(A, B)
751* (should this be norm of (A,B) or (AI,BI)?)
752*
753 CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
754 CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
755 \$ work( mplusn*mplusn+1 ), mplusn )
756 abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn, rwork )
757*
758* Do tests (1) to (4)
759*
760 CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
761 \$ rwork, result( 1 ) )
762 CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
763 \$ rwork, result( 2 ) )
764 CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
765 \$ rwork, result( 3 ) )
766 CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
767 \$ rwork, result( 4 ) )
768*
769* Do tests (5) and (6): check Schur form of A and compare
770* eigenvalues with diagonals.
771*
772 ntest = 6
773 temp1 = zero
774 result( 5 ) = zero
775 result( 6 ) = zero
776*
777 DO 110 j = 1, mplusn
779 temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
780 \$ max( smlnum, abs1( alpha( j ) ), abs1( ai( j, j ) ) )+
781 \$ abs1( beta( j )-bi( j, j ) ) /
782 \$ max( smlnum, abs1( beta( j ) ), abs1( bi( j, j ) ) ) )
783 \$ / ulp
784 IF( j.LT.mplusn ) THEN
785 IF( ai( j+1, j ).NE.zero ) THEN
787 result( 5 ) = ulpinv
788 END IF
789 END IF
790 IF( j.GT.1 ) THEN
791 IF( ai( j, j-1 ).NE.zero ) THEN
793 result( 5 ) = ulpinv
794 END IF
795 END IF
796 temp1 = max( temp1, temp2 )
798 WRITE( nout, fmt = 9997 )j, mplusn, nptknt
799 END IF
800 110 CONTINUE
801 result( 6 ) = temp1
802*
803* Test (7) (if sorting worked) <--------- need to be checked.
804*
805 ntest = 7
806 result( 7 ) = zero
807 IF( linfo.EQ.mplusn+3 )
808 \$ result( 7 ) = ulpinv
809*
810* Test (8): compare the estimated value of DIF and its true value.
811*
812 ntest = 8
813 result( 8 ) = zero
814 IF( difest( 2 ).EQ.zero ) THEN
815 IF( diftru.GT.abnrm*ulp )
816 \$ result( 8 ) = ulpinv
817 ELSE IF( diftru.EQ.zero ) THEN
818 IF( difest( 2 ).GT.abnrm*ulp )
819 \$ result( 8 ) = ulpinv
820 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
821 \$ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
822 result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
823 END IF
824*
825* Test (9)
826*
827 ntest = 9
828 result( 9 ) = zero
829 IF( linfo.EQ.( mplusn+2 ) ) THEN
830 IF( diftru.GT.abnrm*ulp )
831 \$ result( 9 ) = ulpinv
832 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
833 \$ result( 9 ) = ulpinv
834 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
835 \$ result( 9 ) = ulpinv
836 END IF
837*
838* Test (10): compare the estimated value of PL and it true value.
839*
840 ntest = 10
841 result( 10 ) = zero
842 IF( pl( 1 ).EQ.zero ) THEN
843 IF( pltru.GT.abnrm*ulp )
844 \$ result( 10 ) = ulpinv
845 ELSE IF( pltru.EQ.zero ) THEN
846 IF( pl( 1 ).GT.abnrm*ulp )
847 \$ result( 10 ) = ulpinv
848 ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
849 \$ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
850 result( 10 ) = ulpinv
851 END IF
852*
853 ntestt = ntestt + ntest
854*
855* Print out tests which fail.
856*
857 DO 120 j = 1, ntest
858 IF( result( j ).GE.thresh ) THEN
859*
860* If this is the first test to fail,
861* print a header to the data file.
862*
863 IF( nerrs.EQ.0 ) THEN
864 WRITE( nout, fmt = 9996 )'CGX'
865*
866* Matrix types
867*
868 WRITE( nout, fmt = 9995 )
869*
870* Tests performed
871*
872 WRITE( nout, fmt = 9993 )'unitary', '''', 'transpose',
873 \$ ( '''', i = 1, 4 )
874*
875 END IF
876 nerrs = nerrs + 1
877 IF( result( j ).LT.10000.0 ) THEN
878 WRITE( nout, fmt = 9990 )nptknt, mplusn, j, result( j )
879 ELSE
880 WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
881 END IF
882 END IF
883*
884 120 CONTINUE
885*
886 130 CONTINUE
887 GO TO 80
888 140 CONTINUE
889*
890 150 CONTINUE
891*
892* Summary
893*
894 CALL alasvm( 'CGX', nout, nerrs, ntestt, 0 )
895*
896 work( 1 ) = maxwrk
897*
898 RETURN
899*
900 9999 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
901 \$ i6, ', JTYPE=', i6, ')' )
902*
903 9998 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
904 \$ i6, ', Input Example #', i2, ')' )
905*
906 9997 FORMAT( ' CDRGSX: S not in Schur form at eigenvalue ', i6, '.',
907 \$ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
908*
909 9996 FORMAT( / 1x, a3, ' -- Complex Expert Generalized Schur form',
910 \$ ' problem driver' )
911*
912 9995 FORMAT( 'Input Example' )
913*
914 9994 FORMAT( ' Matrix types: ', /
915 \$ ' 1: A is a block diagonal matrix of Jordan blocks ',
916 \$ 'and B is the identity ', / ' matrix, ',
917 \$ / ' 2: A and B are upper triangular matrices, ',
918 \$ / ' 3: A and B are as type 2, but each second diagonal ',
919 \$ 'block in A_11 and ', /
920 \$ ' each third diaongal block in A_22 are 2x2 blocks,',
921 \$ / ' 4: A and B are block diagonal matrices, ',
922 \$ / ' 5: (A,B) has potentially close or common ',
923 \$ 'eigenvalues.', / )
924*
925 9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
926 \$ 'Q and Z are ', a, ',', / 19x,
927 \$ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
928 \$ / ' 1 = | A - Q S Z', a,
929 \$ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
930 \$ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
931 \$ ' | / ( n ulp ) 4 = | I - ZZ', a,
932 \$ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
933 \$ 'Schur form S', / ' 6 = difference between (alpha,beta)',
934 \$ ' and diagonals of (S,T)', /
935 \$ ' 7 = 1/ULP if SDIM is not the correct number of ',
936 \$ 'selected eigenvalues', /
937 \$ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
938 \$ 'DIFTRU/DIFEST > 10*THRESH',
939 \$ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
940 \$ 'when reordering fails', /
941 \$ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
942 \$ 'PLTRU/PLEST > THRESH', /
943 \$ ' ( Test 10 is only for input examples )', / )
944 9992 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
945 \$ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
946 9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
947 \$ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, e10.3 )
948 9990 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
949 \$ ' result ', i2, ' is', 0p, f8.2 )
950 9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
951 \$ ' result ', i2, ' is', 1p, e10.3 )
952*
953* End of CDRGSX
954*
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 cget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RWORK, RESULT)
CGET51
Definition: cget51.f:155
logical function clctsx(ALPHA, BETA)
CLCTSX
Definition: clctsx.f:57
subroutine clatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
CLATM5
Definition: clatm5.f:268
subroutine clakf2(M, N, A, LDA, B, D, E, Z, LDZ)
CLAKF2
Definition: clakf2.f:105
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:115
subroutine cggesx(JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
CGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: cggesx.f:330
subroutine cgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
CGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: cgesvd.f:214
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
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: