LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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 ..
378 LOGICAL ILABAD
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 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 CGGESX, 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 claset( 'Full', mplusn, mplusn, czero, czero, ai,
514 $ lda )
515 CALL claset( 'Full', mplusn, mplusn, czero, czero, bi,
516 $ lda )
517*
518 CALL clatm5( 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 CLCTSX
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 clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
539 CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
540*
541 CALL cggesx( 'V', 'V', 'S', clctsx, sense, mplusn, ai,
542 $ lda, bi, lda, mm, alpha, beta, q, lda, z,
543 $ lda, pl, difest, work, lwork, rwork,
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 )'CGGESX', linfo, mplusn,
549 $ prtype
550 info = linfo
551 GO TO 30
552 END IF
553*
554* Compute the norm(A, B)
555*
556 CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work,
557 $ mplusn )
558 CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
559 $ work( mplusn*mplusn+1 ), mplusn )
560 abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn,
561 $ rwork )
562*
563* Do tests (1) to (4)
564*
565 result( 2 ) = zero
566 CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
567 $ lda, work, rwork, result( 1 ) )
568 CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
569 $ lda, work, rwork, result( 2 ) )
570 CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
571 $ lda, work, rwork, result( 3 ) )
572 CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
573 $ lda, work, rwork, 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
584 ilabad = .false.
585 temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
586 $ max( smlnum, abs1( alpha( j ) ),
587 $ abs1( ai( j, j ) ) )+
588 $ abs1( beta( j )-bi( j, j ) ) /
589 $ max( smlnum, abs1( beta( j ) ),
590 $ abs1( 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 temp1 = max( temp1, temp2 )
604 IF( ilabad ) THEN
605 WRITE( nout, fmt = 9997 )j, mplusn, prtype
606 END IF
607 10 CONTINUE
608 result( 6 ) = temp1
609 ntest = ntest + 2
610*
611* Test (7) (if sorting worked)
612*
613 result( 7 ) = zero
614 IF( linfo.EQ.mplusn+3 ) THEN
615 result( 7 ) = ulpinv
616 ELSE IF( mm.NE.n ) THEN
617 result( 7 ) = ulpinv
618 END IF
619 ntest = ntest + 1
620*
621* Test (8): compare the estimated value DIF and its
622* value. first, compute the exact DIF.
623*
624 result( 8 ) = zero
625 mn2 = mm*( mplusn-mm )*2
626 IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
627*
628* Note: for either following two cases, there are
629* almost same number of test cases fail the test.
630*
631 CALL clakf2( mm, mplusn-mm, ai, lda,
632 $ ai( mm+1, mm+1 ), bi,
633 $ bi( mm+1, mm+1 ), c, ldc )
634*
635 CALL cgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
636 $ 1, work( 2 ), 1, work( 3 ), lwork-2,
637 $ rwork, info )
638 diftru = s( mn2 )
639*
640 IF( difest( 2 ).EQ.zero ) THEN
641 IF( diftru.GT.abnrm*ulp )
642 $ result( 8 ) = ulpinv
643 ELSE IF( diftru.EQ.zero ) THEN
644 IF( difest( 2 ).GT.abnrm*ulp )
645 $ result( 8 ) = ulpinv
646 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
647 $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
648 result( 8 ) = max( diftru / difest( 2 ),
649 $ difest( 2 ) / diftru )
650 END IF
651 ntest = ntest + 1
652 END IF
653*
654* Test (9)
655*
656 result( 9 ) = zero
657 IF( linfo.EQ.( mplusn+2 ) ) THEN
658 IF( diftru.GT.abnrm*ulp )
659 $ result( 9 ) = ulpinv
660 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
661 $ result( 9 ) = ulpinv
662 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
663 $ result( 9 ) = ulpinv
664 ntest = ntest + 1
665 END IF
666*
667 ntestt = ntestt + ntest
668*
669* Print out tests which fail.
670*
671 DO 20 j = 1, 9
672 IF( result( j ).GE.thresh ) THEN
673*
674* If this is the first test to fail,
675* print a header to the data file.
676*
677 IF( nerrs.EQ.0 ) THEN
678 WRITE( nout, fmt = 9996 )'CGX'
679*
680* Matrix types
681*
682 WRITE( nout, fmt = 9994 )
683*
684* Tests performed
685*
686 WRITE( nout, fmt = 9993 )'unitary', '''',
687 $ 'transpose', ( '''', i = 1, 4 )
688*
689 END IF
690 nerrs = nerrs + 1
691 IF( result( j ).LT.10000.0 ) THEN
692 WRITE( nout, fmt = 9992 )mplusn, prtype,
693 $ weight, m, j, result( j )
694 ELSE
695 WRITE( nout, fmt = 9991 )mplusn, prtype,
696 $ weight, m, j, result( j )
697 END IF
698 END IF
699 20 CONTINUE
700*
701 30 CONTINUE
702 40 CONTINUE
703 50 CONTINUE
704 60 CONTINUE
705*
706 GO TO 150
707*
708 70 CONTINUE
709*
710* Read in data from file to check accuracy of condition estimation
711* Read input data until N=0
712*
713 nptknt = 0
714*
715 80 CONTINUE
716 READ( nin, fmt = *, END = 140 )mplusn
717 IF( mplusn.EQ.0 )
718 $ GO TO 140
719 READ( nin, fmt = *, END = 140 )n
720 DO 90 i = 1, mplusn
721 READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
722 90 CONTINUE
723 DO 100 i = 1, mplusn
724 READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
725 100 CONTINUE
726 READ( nin, fmt = * )pltru, diftru
727*
728 nptknt = nptknt + 1
729 fs = .true.
730 k = 0
731 m = mplusn - n
732*
733 CALL clacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
734 CALL clacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
735*
736* Compute the Schur factorization while swapping the
737* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
738*
739 CALL cggesx( 'V', 'V', 'S', clctsx, 'B', mplusn, ai, lda, bi, lda,
740 $ mm, alpha, beta, q, lda, z, lda, pl, difest, work,
741 $ lwork, rwork, iwork, liwork, bwork, linfo )
742*
743 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
744 result( 1 ) = ulpinv
745 WRITE( nout, fmt = 9998 )'CGGESX', linfo, mplusn, nptknt
746 GO TO 130
747 END IF
748*
749* Compute the norm(A, B)
750* (should this be norm of (A,B) or (AI,BI)?)
751*
752 CALL clacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
753 CALL clacpy( 'Full', mplusn, mplusn, bi, lda,
754 $ work( mplusn*mplusn+1 ), mplusn )
755 abnrm = clange( 'Fro', mplusn, 2*mplusn, work, mplusn, rwork )
756*
757* Do tests (1) to (4)
758*
759 CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
760 $ rwork, result( 1 ) )
761 CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
762 $ rwork, result( 2 ) )
763 CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
764 $ rwork, result( 3 ) )
765 CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
766 $ rwork, result( 4 ) )
767*
768* Do tests (5) and (6): check Schur form of A and compare
769* eigenvalues with diagonals.
770*
771 ntest = 6
772 temp1 = zero
773 result( 5 ) = zero
774 result( 6 ) = zero
775*
776 DO 110 j = 1, mplusn
777 ilabad = .false.
778 temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
779 $ max( smlnum, abs1( alpha( j ) ), abs1( ai( j, j ) ) )+
780 $ abs1( beta( j )-bi( j, j ) ) /
781 $ max( smlnum, abs1( beta( j ) ), abs1( bi( j, j ) ) ) )
782 $ / ulp
783 IF( j.LT.mplusn ) THEN
784 IF( ai( j+1, j ).NE.zero ) THEN
785 ilabad = .true.
786 result( 5 ) = ulpinv
787 END IF
788 END IF
789 IF( j.GT.1 ) THEN
790 IF( ai( j, j-1 ).NE.zero ) THEN
791 ilabad = .true.
792 result( 5 ) = ulpinv
793 END IF
794 END IF
795 temp1 = max( temp1, temp2 )
796 IF( ilabad ) THEN
797 WRITE( nout, fmt = 9997 )j, mplusn, nptknt
798 END IF
799 110 CONTINUE
800 result( 6 ) = temp1
801*
802* Test (7) (if sorting worked) <--------- need to be checked.
803*
804 ntest = 7
805 result( 7 ) = zero
806 IF( linfo.EQ.mplusn+3 )
807 $ result( 7 ) = ulpinv
808*
809* Test (8): compare the estimated value of DIF and its true value.
810*
811 ntest = 8
812 result( 8 ) = zero
813 IF( difest( 2 ).EQ.zero ) THEN
814 IF( diftru.GT.abnrm*ulp )
815 $ result( 8 ) = ulpinv
816 ELSE IF( diftru.EQ.zero ) THEN
817 IF( difest( 2 ).GT.abnrm*ulp )
818 $ result( 8 ) = ulpinv
819 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
820 $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
821 result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
822 END IF
823*
824* Test (9)
825*
826 ntest = 9
827 result( 9 ) = zero
828 IF( linfo.EQ.( mplusn+2 ) ) THEN
829 IF( diftru.GT.abnrm*ulp )
830 $ result( 9 ) = ulpinv
831 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
832 $ result( 9 ) = ulpinv
833 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
834 $ result( 9 ) = ulpinv
835 END IF
836*
837* Test (10): compare the estimated value of PL and it true value.
838*
839 ntest = 10
840 result( 10 ) = zero
841 IF( pl( 1 ).EQ.zero ) THEN
842 IF( pltru.GT.abnrm*ulp )
843 $ result( 10 ) = ulpinv
844 ELSE IF( pltru.EQ.zero ) THEN
845 IF( pl( 1 ).GT.abnrm*ulp )
846 $ result( 10 ) = ulpinv
847 ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
848 $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
849 result( 10 ) = ulpinv
850 END IF
851*
852 ntestt = ntestt + ntest
853*
854* Print out tests which fail.
855*
856 DO 120 j = 1, ntest
857 IF( result( j ).GE.thresh ) THEN
858*
859* If this is the first test to fail,
860* print a header to the data file.
861*
862 IF( nerrs.EQ.0 ) THEN
863 WRITE( nout, fmt = 9996 )'CGX'
864*
865* Matrix types
866*
867 WRITE( nout, fmt = 9995 )
868*
869* Tests performed
870*
871 WRITE( nout, fmt = 9993 )'unitary', '''', 'transpose',
872 $ ( '''', i = 1, 4 )
873*
874 END IF
875 nerrs = nerrs + 1
876 IF( result( j ).LT.10000.0 ) THEN
877 WRITE( nout, fmt = 9990 )nptknt, mplusn, j, result( j )
878 ELSE
879 WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
880 END IF
881 END IF
882*
883 120 CONTINUE
884*
885 130 CONTINUE
886 GO TO 80
887 140 CONTINUE
888*
889 150 CONTINUE
890*
891* Summary
892*
893 CALL alasvm( 'CGX', nout, nerrs, ntestt, 0 )
894*
895 work( 1 ) = maxwrk
896*
897 RETURN
898*
899 9999 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
900 $ i6, ', JTYPE=', i6, ')' )
901*
902 9998 FORMAT( ' CDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
903 $ i6, ', Input Example #', i2, ')' )
904*
905 9997 FORMAT( ' CDRGSX: S not in Schur form at eigenvalue ', i6, '.',
906 $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
907*
908 9996 FORMAT( / 1x, a3, ' -- Complex Expert Generalized Schur form',
909 $ ' problem driver' )
910*
911 9995 FORMAT( 'Input Example' )
912*
913 9994 FORMAT( ' Matrix types: ', /
914 $ ' 1: A is a block diagonal matrix of Jordan blocks ',
915 $ 'and B is the identity ', / ' matrix, ',
916 $ / ' 2: A and B are upper triangular matrices, ',
917 $ / ' 3: A and B are as type 2, but each second diagonal ',
918 $ 'block in A_11 and ', /
919 $ ' each third diagonal block in A_22 are 2x2 blocks,',
920 $ / ' 4: A and B are block diagonal matrices, ',
921 $ / ' 5: (A,B) has potentially close or common ',
922 $ 'eigenvalues.', / )
923*
924 9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
925 $ 'Q and Z are ', a, ',', / 19x,
926 $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
927 $ / ' 1 = | A - Q S Z', a,
928 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
929 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
930 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
931 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
932 $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
933 $ ' and diagonals of (S,T)', /
934 $ ' 7 = 1/ULP if SDIM is not the correct number of ',
935 $ 'selected eigenvalues', /
936 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
937 $ 'DIFTRU/DIFEST > 10*THRESH',
938 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
939 $ 'when reordering fails', /
940 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
941 $ 'PLTRU/PLEST > THRESH', /
942 $ ' ( Test 10 is only for input examples )', / )
943 9992 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
944 $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
945 9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
946 $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, e10.3 )
947 9990 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
948 $ ' result ', i2, ' is', 0p, f8.2 )
949 9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
950 $ ' result ', i2, ' is', 1p, e10.3 )
951*
952* End of CDRGSX
953*
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cget51(itype, n, a, lda, b, ldb, u, ldu, v, ldv, work, rwork, result)
CGET51
Definition cget51.f:155
subroutine clakf2(m, n, a, lda, b, d, e, z, ldz)
CLAKF2
Definition clakf2.f:105
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
logical function clctsx(alpha, beta)
CLCTSX
Definition clctsx.f:57
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 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
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
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 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
Here is the call graph for this function:
Here is the caller graph for this function: