LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cdrges3()

subroutine cdrges3 ( integer  NSIZES,
integer, dimension( * )  NN,
integer  NTYPES,
logical, dimension( * )  DOTYPE,
integer, dimension( 4 )  ISEED,
real  THRESH,
integer  NOUNIT,
complex, dimension( lda, * )  A,
integer  LDA,
complex, dimension( lda, * )  B,
complex, dimension( lda, * )  S,
complex, dimension( lda, * )  T,
complex, dimension( ldq, * )  Q,
integer  LDQ,
complex, dimension( ldq, * )  Z,
complex, dimension( * )  ALPHA,
complex, dimension( * )  BETA,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
real, dimension( 13 )  RESULT,
logical, dimension( * )  BWORK,
integer  INFO 
)

CDRGES3

Purpose:
 CDRGES3 checks the nonsymmetric generalized eigenvalue (Schur form)
 problem driver CGGES3.

 CGGES3 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 reorder the eigenvalues so that a selected
 cluster of eigenvalues appears in the leading diagonal block of the
 Schur forms.

 When CDRGES3 is called, a number of matrix "sizes" ("N's") and a
 number of matrix "TYPES" are specified.  For each size ("N")
 and each TYPE of matrix, a pair of matrices (A, B) will be generated
 and used for testing. For each matrix pair, the following 13 tests
 will be performed and compared with the threshold THRESH except
 the tests (5), (11) and (13).


 (1)   | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues)


 (2)   | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues)


 (3)   | I - QQ' | / ( n ulp ) (no sorting of eigenvalues)


 (4)   | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues)

 (5)   if A is in Schur form (i.e. triangular form) (no sorting of
       eigenvalues)

 (6)   if eigenvalues = diagonal elements of the Schur form (S, T),
       i.e., test the 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)|)

       (no sorting of eigenvalues)

 (7)   | (A,B) - Q (S,T) Z' | / ( |(A,B)| n ulp )
       (with sorting of eigenvalues).

 (8)   | I - QQ' | / ( n ulp ) (with sorting of eigenvalues).

 (9)   | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues).

 (10)  if A is in Schur form (i.e. quasi-triangular form)
       (with sorting of eigenvalues).

 (11)  if eigenvalues = diagonal elements of the Schur form (S, T),
       i.e. test the 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)|)

       (with sorting of eigenvalues).

 (12)  if sorting worked and SDIM is the number of eigenvalues
       which were CELECTed.

 Test Matrices
 =============

 The sizes of the test matrices are specified by an array
 NN(1:NSIZES); the value of each element NN(j) specifies one size.
 The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
 DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
 Currently, the list of possible types is:

 (1)  ( 0, 0 )         (a pair of zero matrices)

 (2)  ( I, 0 )         (an identity and a zero matrix)

 (3)  ( 0, I )         (an identity and a zero matrix)

 (4)  ( I, I )         (a pair of identity matrices)

         t   t
 (5)  ( J , J  )       (a pair of transposed Jordan blocks)

                                     t                ( I   0  )
 (6)  ( X, Y )         where  X = ( J   0  )  and Y = (      t )
                                  ( 0   I  )          ( 0   J  )
                       and I is a k x k identity and J a (k+1)x(k+1)
                       Jordan block; k=(N-1)/2

 (7)  ( D, I )         where D is diag( 0, 1,..., N-1 ) (a diagonal
                       matrix with those diagonal entries.)
 (8)  ( I, D )

 (9)  ( big*D, small*I ) where "big" is near overflow and small=1/big

 (10) ( small*D, big*I )

 (11) ( big*I, small*D )

 (12) ( small*I, big*D )

 (13) ( big*D, big*I )

 (14) ( small*D, small*I )

 (15) ( D1, D2 )        where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
                        D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
           t   t
 (16) Q ( J , J ) Z     where Q and Z are random orthogonal matrices.

 (17) Q ( T1, T2 ) Z    where T1 and T2 are upper triangular matrices
                        with random O(1) entries above the diagonal
                        and diagonal entries diag(T1) =
                        ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
                        ( 0, N-3, N-4,..., 1, 0, 0 )

 (18) Q ( T1, T2 ) Z    diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
                        diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
                        s = machine precision.

 (19) Q ( T1, T2 ) Z    diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )

                                                        N-5
 (20) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, 1, a, ..., a   =s, 0 )
                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )

 (21) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
                        where r1,..., r(N-4) are random.

 (22) Q ( big*T1, small*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )

 (23) Q ( small*T1, big*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )

 (24) Q ( small*T1, small*T2 ) Z  diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )

 (25) Q ( big*T1, big*T2 ) Z      diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )

 (26) Q ( T1, T2 ) Z     where T1 and T2 are random upper-triangular
                         matrices.
Parameters
[in]NSIZES
          NSIZES is INTEGER
          The number of sizes of matrices to use.  If it is zero,
          SDRGES3 does nothing.  NSIZES >= 0.
[in]NN
          NN is INTEGER array, dimension (NSIZES)
          An array containing the sizes to be used for the matrices.
          Zero values will be skipped.  NN >= 0.
[in]NTYPES
          NTYPES is INTEGER
          The number of elements in DOTYPE.   If it is zero, SDRGES3
          does nothing.  It must be at least zero.  If it is MAXTYP+1
          and NSIZES is 1, then an additional type, MAXTYP+1 is
          defined, which is to use whatever matrix is in A on input.
          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
          DOTYPE(MAXTYP+1) is .TRUE. .
[in]DOTYPE
          DOTYPE is LOGICAL array, dimension (NTYPES)
          If DOTYPE(j) is .TRUE., then for each size in NN a
          matrix of that size and of type j will be generated.
          If NTYPES is smaller than the maximum number of types
          defined (PARAMETER MAXTYP), then types NTYPES+1 through
          MAXTYP will not be generated. If NTYPES is larger
          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
          will be ignored.
[in,out]ISEED
          ISEED is INTEGER array, dimension (4)
          On entry ISEED specifies the seed of the random number
          generator. The array elements should be between 0 and 4095;
          if not they will be reduced mod 4096. Also, ISEED(4) must
          be odd.  The random number generator uses a linear
          congruential sequence limited to small integers, and so
          should produce machine independent random numbers. The
          values of ISEED are changed on exit, and can be used in the
          next call to SDRGES3 to continue the same random number
          sequence.
[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]NOUNIT
          NOUNIT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns IINFO not equal to 0.)
[in,out]A
          A is COMPLEX array, dimension(LDA, max(NN))
          Used to hold the original A matrix.  Used as input only
          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
          DOTYPE(MAXTYP+1)=.TRUE.
[in]LDA
          LDA is INTEGER
          The leading dimension of A, B, S, and T.
          It must be at least 1 and at least max( NN ).
[in,out]B
          B is COMPLEX array, dimension(LDA, max(NN))
          Used to hold the original B matrix.  Used as input only
          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
          DOTYPE(MAXTYP+1)=.TRUE.
[out]S
          S is COMPLEX array, dimension (LDA, max(NN))
          The Schur form matrix computed from A by CGGES3.  On exit, S
          contains the Schur form matrix corresponding to the matrix
          in A.
[out]T
          T is COMPLEX array, dimension (LDA, max(NN))
          The upper triangular matrix computed from B by CGGES3.
[out]Q
          Q is COMPLEX array, dimension (LDQ, max(NN))
          The (left) orthogonal matrix computed by CGGES3.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of Q and Z. It must
          be at least 1 and at least max( NN ).
[out]Z
          Z is COMPLEX array, dimension( LDQ, max(NN) )
          The (right) orthogonal matrix computed by CGGES3.
[out]ALPHA
          ALPHA is COMPLEX array, dimension (max(NN))
[out]BETA
          BETA is COMPLEX array, dimension (max(NN))

          The generalized eigenvalues of (A,B) computed by CGGES3.
          ALPHA(k) / BETA(k) is the k-th generalized eigenvalue of A
          and B.
[out]WORK
          WORK is COMPLEX array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= 3*N*N.
[out]RWORK
          RWORK is REAL array, dimension ( 8*N )
          Real workspace.
[out]RESULT
          RESULT is REAL array, dimension (15)
          The values computed by the tests described above.
          The values are currently limited to 1/ulp, to avoid overflow.
[out]BWORK
          BWORK is LOGICAL array, dimension (N)
[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.  INFO is the
                absolute value of the INFO value returned.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 378 of file cdrges3.f.

382 *
383 * -- LAPACK test routine --
384 * -- LAPACK is a software package provided by Univ. of Tennessee, --
385 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
386 *
387 * .. Scalar Arguments ..
388  INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
389  REAL THRESH
390 * ..
391 * .. Array Arguments ..
392  LOGICAL BWORK( * ), DOTYPE( * )
393  INTEGER ISEED( 4 ), NN( * )
394  REAL RESULT( 13 ), RWORK( * )
395  COMPLEX A( LDA, * ), ALPHA( * ), B( LDA, * ),
396  $ BETA( * ), Q( LDQ, * ), S( LDA, * ),
397  $ T( LDA, * ), WORK( * ), Z( LDQ, * )
398 * ..
399 *
400 * =====================================================================
401 *
402 * .. Parameters ..
403  REAL ZERO, ONE
404  parameter( zero = 0.0e+0, one = 1.0e+0 )
405  COMPLEX CZERO, CONE
406  parameter( czero = ( 0.0e+0, 0.0e+0 ),
407  $ cone = ( 1.0e+0, 0.0e+0 ) )
408  INTEGER MAXTYP
409  parameter( maxtyp = 26 )
410 * ..
411 * .. Local Scalars ..
412  LOGICAL BADNN, ILABAD
413  CHARACTER SORT
414  INTEGER I, IADD, IINFO, IN, ISORT, J, JC, JR, JSIZE,
415  $ JTYPE, KNTEIG, MAXWRK, MINWRK, MTYPES, N, N1,
416  $ NB, NERRS, NMATS, NMAX, NTEST, NTESTT, RSUB,
417  $ SDIM
418  REAL SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
419  COMPLEX CTEMP, X
420 * ..
421 * .. Local Arrays ..
422  LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
423  INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
424  $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
425  $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
426  $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
427  $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
428  REAL RMAGN( 0: 3 )
429 * ..
430 * .. External Functions ..
431  LOGICAL CLCTES
432  INTEGER ILAENV
433  REAL SLAMCH
434  COMPLEX CLARND
435  EXTERNAL clctes, ilaenv, slamch, clarnd
436 * ..
437 * .. External Subroutines ..
438  EXTERNAL alasvm, cget51, cget54, cgges3, clacpy, clarfg,
440 * ..
441 * .. Intrinsic Functions ..
442  INTRINSIC abs, aimag, conjg, max, min, real, sign
443 * ..
444 * .. Statement Functions ..
445  REAL ABS1
446 * ..
447 * .. Statement Function definitions ..
448  abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
449 * ..
450 * .. Data statements ..
451  DATA kclass / 15*1, 10*2, 1*3 /
452  DATA kz1 / 0, 1, 2, 1, 3, 3 /
453  DATA kz2 / 0, 0, 1, 2, 1, 1 /
454  DATA kadd / 0, 0, 0, 0, 3, 2 /
455  DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
456  $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
457  DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
458  $ 1, 1, -4, 2, -4, 8*8, 0 /
459  DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
460  $ 4*5, 4*3, 1 /
461  DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
462  $ 4*6, 4*4, 1 /
463  DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
464  $ 2, 1 /
465  DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
466  $ 2, 1 /
467  DATA ktrian / 16*0, 10*1 /
468  DATA lasign / 6*.false., .true., .false., 2*.true.,
469  $ 2*.false., 3*.true., .false., .true.,
470  $ 3*.false., 5*.true., .false. /
471  DATA lbsign / 7*.false., .true., 2*.false.,
472  $ 2*.true., 2*.false., .true., .false., .true.,
473  $ 9*.false. /
474 * ..
475 * .. Executable Statements ..
476 *
477 * Check for errors
478 *
479  info = 0
480 *
481  badnn = .false.
482  nmax = 1
483  DO 10 j = 1, nsizes
484  nmax = max( nmax, nn( j ) )
485  IF( nn( j ).LT.0 )
486  $ badnn = .true.
487  10 CONTINUE
488 *
489  IF( nsizes.LT.0 ) THEN
490  info = -1
491  ELSE IF( badnn ) THEN
492  info = -2
493  ELSE IF( ntypes.LT.0 ) THEN
494  info = -3
495  ELSE IF( thresh.LT.zero ) THEN
496  info = -6
497  ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
498  info = -9
499  ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
500  info = -14
501  END IF
502 *
503 * Compute workspace
504 * (Note: Comments in the code beginning "Workspace:" describe the
505 * minimal amount of workspace needed at that point in the code,
506 * as well as the preferred amount for good performance.
507 * NB refers to the optimal block size for the immediately
508 * following subroutine, as returned by ILAENV.
509 *
510  minwrk = 1
511  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
512  minwrk = 3*nmax*nmax
513  nb = max( 1, ilaenv( 1, 'CGEQRF', ' ', nmax, nmax, -1, -1 ),
514  $ ilaenv( 1, 'CUNMQR', 'LC', nmax, nmax, nmax, -1 ),
515  $ ilaenv( 1, 'CUNGQR', ' ', nmax, nmax, nmax, -1 ) )
516  maxwrk = max( nmax+nmax*nb, 3*nmax*nmax)
517  work( 1 ) = maxwrk
518  END IF
519 *
520  IF( lwork.LT.minwrk )
521  $ info = -19
522 *
523  IF( info.NE.0 ) THEN
524  CALL xerbla( 'CDRGES3', -info )
525  RETURN
526  END IF
527 *
528 * Quick return if possible
529 *
530  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
531  $ RETURN
532 *
533  ulp = slamch( 'Precision' )
534  safmin = slamch( 'Safe minimum' )
535  safmin = safmin / ulp
536  safmax = one / safmin
537  CALL slabad( safmin, safmax )
538  ulpinv = one / ulp
539 *
540 * The values RMAGN(2:3) depend on N, see below.
541 *
542  rmagn( 0 ) = zero
543  rmagn( 1 ) = one
544 *
545 * Loop over matrix sizes
546 *
547  ntestt = 0
548  nerrs = 0
549  nmats = 0
550 *
551  DO 190 jsize = 1, nsizes
552  n = nn( jsize )
553  n1 = max( 1, n )
554  rmagn( 2 ) = safmax*ulp / real( n1 )
555  rmagn( 3 ) = safmin*ulpinv*real( n1 )
556 *
557  IF( nsizes.NE.1 ) THEN
558  mtypes = min( maxtyp, ntypes )
559  ELSE
560  mtypes = min( maxtyp+1, ntypes )
561  END IF
562 *
563 * Loop over matrix types
564 *
565  DO 180 jtype = 1, mtypes
566  IF( .NOT.dotype( jtype ) )
567  $ GO TO 180
568  nmats = nmats + 1
569  ntest = 0
570 *
571 * Save ISEED in case of an error.
572 *
573  DO 20 j = 1, 4
574  ioldsd( j ) = iseed( j )
575  20 CONTINUE
576 *
577 * Initialize RESULT
578 *
579  DO 30 j = 1, 13
580  result( j ) = zero
581  30 CONTINUE
582 *
583 * Generate test matrices A and B
584 *
585 * Description of control parameters:
586 *
587 * KCLASS: =1 means w/o rotation, =2 means w/ rotation,
588 * =3 means random.
589 * KATYPE: the "type" to be passed to CLATM4 for computing A.
590 * KAZERO: the pattern of zeros on the diagonal for A:
591 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
592 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
593 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
594 * non-zero entries.)
595 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
596 * =2: large, =3: small.
597 * LASIGN: .TRUE. if the diagonal elements of A are to be
598 * multiplied by a random magnitude 1 number.
599 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
600 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
601 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
602 * RMAGN: used to implement KAMAGN and KBMAGN.
603 *
604  IF( mtypes.GT.maxtyp )
605  $ GO TO 110
606  iinfo = 0
607  IF( kclass( jtype ).LT.3 ) THEN
608 *
609 * Generate A (w/o rotation)
610 *
611  IF( abs( katype( jtype ) ).EQ.3 ) THEN
612  in = 2*( ( n-1 ) / 2 ) + 1
613  IF( in.NE.n )
614  $ CALL claset( 'Full', n, n, czero, czero, a, lda )
615  ELSE
616  in = n
617  END IF
618  CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
619  $ kz2( kazero( jtype ) ), lasign( jtype ),
620  $ rmagn( kamagn( jtype ) ), ulp,
621  $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
622  $ iseed, a, lda )
623  iadd = kadd( kazero( jtype ) )
624  IF( iadd.GT.0 .AND. iadd.LE.n )
625  $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
626 *
627 * Generate B (w/o rotation)
628 *
629  IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
630  in = 2*( ( n-1 ) / 2 ) + 1
631  IF( in.NE.n )
632  $ CALL claset( 'Full', n, n, czero, czero, b, lda )
633  ELSE
634  in = n
635  END IF
636  CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
637  $ kz2( kbzero( jtype ) ), lbsign( jtype ),
638  $ rmagn( kbmagn( jtype ) ), one,
639  $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
640  $ iseed, b, lda )
641  iadd = kadd( kbzero( jtype ) )
642  IF( iadd.NE.0 .AND. iadd.LE.n )
643  $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
644 *
645  IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
646 *
647 * Include rotations
648 *
649 * Generate Q, Z as Householder transformations times
650 * a diagonal matrix.
651 *
652  DO 50 jc = 1, n - 1
653  DO 40 jr = jc, n
654  q( jr, jc ) = clarnd( 3, iseed )
655  z( jr, jc ) = clarnd( 3, iseed )
656  40 CONTINUE
657  CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
658  $ work( jc ) )
659  work( 2*n+jc ) = sign( one, real( q( jc, jc ) ) )
660  q( jc, jc ) = cone
661  CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
662  $ work( n+jc ) )
663  work( 3*n+jc ) = sign( one, real( z( jc, jc ) ) )
664  z( jc, jc ) = cone
665  50 CONTINUE
666  ctemp = clarnd( 3, iseed )
667  q( n, n ) = cone
668  work( n ) = czero
669  work( 3*n ) = ctemp / abs( ctemp )
670  ctemp = clarnd( 3, iseed )
671  z( n, n ) = cone
672  work( 2*n ) = czero
673  work( 4*n ) = ctemp / abs( ctemp )
674 *
675 * Apply the diagonal matrices
676 *
677  DO 70 jc = 1, n
678  DO 60 jr = 1, n
679  a( jr, jc ) = work( 2*n+jr )*
680  $ conjg( work( 3*n+jc ) )*
681  $ a( jr, jc )
682  b( jr, jc ) = work( 2*n+jr )*
683  $ conjg( work( 3*n+jc ) )*
684  $ b( jr, jc )
685  60 CONTINUE
686  70 CONTINUE
687  CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
688  $ lda, work( 2*n+1 ), iinfo )
689  IF( iinfo.NE.0 )
690  $ GO TO 100
691  CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
692  $ a, lda, work( 2*n+1 ), iinfo )
693  IF( iinfo.NE.0 )
694  $ GO TO 100
695  CALL cunm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
696  $ lda, work( 2*n+1 ), iinfo )
697  IF( iinfo.NE.0 )
698  $ GO TO 100
699  CALL cunm2r( 'R', 'C', n, n, n-1, z, ldq, work( n+1 ),
700  $ b, lda, work( 2*n+1 ), iinfo )
701  IF( iinfo.NE.0 )
702  $ GO TO 100
703  END IF
704  ELSE
705 *
706 * Random matrices
707 *
708  DO 90 jc = 1, n
709  DO 80 jr = 1, n
710  a( jr, jc ) = rmagn( kamagn( jtype ) )*
711  $ clarnd( 4, iseed )
712  b( jr, jc ) = rmagn( kbmagn( jtype ) )*
713  $ clarnd( 4, iseed )
714  80 CONTINUE
715  90 CONTINUE
716  END IF
717 *
718  100 CONTINUE
719 *
720  IF( iinfo.NE.0 ) THEN
721  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
722  $ ioldsd
723  info = abs( iinfo )
724  RETURN
725  END IF
726 *
727  110 CONTINUE
728 *
729  DO 120 i = 1, 13
730  result( i ) = -one
731  120 CONTINUE
732 *
733 * Test with and without sorting of eigenvalues
734 *
735  DO 150 isort = 0, 1
736  IF( isort.EQ.0 ) THEN
737  sort = 'N'
738  rsub = 0
739  ELSE
740  sort = 'S'
741  rsub = 5
742  END IF
743 *
744 * Call XLAENV to set the parameters used in CLAQZ0
745 *
746  CALL xlaenv( 12, 10 )
747  CALL xlaenv( 13, 12 )
748  CALL xlaenv( 14, 13 )
749  CALL xlaenv( 15, 2 )
750  CALL xlaenv( 17, 10 )
751 *
752 * Call CGGES3 to compute H, T, Q, Z, alpha, and beta.
753 *
754  CALL clacpy( 'Full', n, n, a, lda, s, lda )
755  CALL clacpy( 'Full', n, n, b, lda, t, lda )
756  ntest = 1 + rsub + isort
757  result( 1+rsub+isort ) = ulpinv
758  CALL cgges3( 'V', 'V', sort, clctes, n, s, lda, t, lda,
759  $ sdim, alpha, beta, q, ldq, z, ldq, work,
760  $ lwork, rwork, bwork, iinfo )
761  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
762  result( 1+rsub+isort ) = ulpinv
763  WRITE( nounit, fmt = 9999 )'CGGES3', iinfo, n, jtype,
764  $ ioldsd
765  info = abs( iinfo )
766  GO TO 160
767  END IF
768 *
769  ntest = 4 + rsub
770 *
771 * Do tests 1--4 (or tests 7--9 when reordering )
772 *
773  IF( isort.EQ.0 ) THEN
774  CALL cget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
775  $ work, rwork, result( 1 ) )
776  CALL cget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
777  $ work, rwork, result( 2 ) )
778  ELSE
779  CALL cget54( n, a, lda, b, lda, s, lda, t, lda, q,
780  $ ldq, z, ldq, work, result( 2+rsub ) )
781  END IF
782 *
783  CALL cget51( 3, n, b, lda, t, lda, q, ldq, q, ldq, work,
784  $ rwork, result( 3+rsub ) )
785  CALL cget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
786  $ rwork, result( 4+rsub ) )
787 *
788 * Do test 5 and 6 (or Tests 10 and 11 when reordering):
789 * check Schur form of A and compare eigenvalues with
790 * diagonals.
791 *
792  ntest = 6 + rsub
793  temp1 = zero
794 *
795  DO 130 j = 1, n
796  ilabad = .false.
797  temp2 = ( abs1( alpha( j )-s( j, j ) ) /
798  $ max( safmin, abs1( alpha( j ) ), abs1( s( j,
799  $ j ) ) )+abs1( beta( j )-t( j, j ) ) /
800  $ max( safmin, abs1( beta( j ) ), abs1( t( j,
801  $ j ) ) ) ) / ulp
802 *
803  IF( j.LT.n ) THEN
804  IF( s( j+1, j ).NE.zero ) THEN
805  ilabad = .true.
806  result( 5+rsub ) = ulpinv
807  END IF
808  END IF
809  IF( j.GT.1 ) THEN
810  IF( s( j, j-1 ).NE.zero ) THEN
811  ilabad = .true.
812  result( 5+rsub ) = ulpinv
813  END IF
814  END IF
815  temp1 = max( temp1, temp2 )
816  IF( ilabad ) THEN
817  WRITE( nounit, fmt = 9998 )j, n, jtype, ioldsd
818  END IF
819  130 CONTINUE
820  result( 6+rsub ) = temp1
821 *
822  IF( isort.GE.1 ) THEN
823 *
824 * Do test 12
825 *
826  ntest = 12
827  result( 12 ) = zero
828  knteig = 0
829  DO 140 i = 1, n
830  IF( clctes( alpha( i ), beta( i ) ) )
831  $ knteig = knteig + 1
832  140 CONTINUE
833  IF( sdim.NE.knteig )
834  $ result( 13 ) = ulpinv
835  END IF
836 *
837  150 CONTINUE
838 *
839 * End of Loop -- Check for RESULT(j) > THRESH
840 *
841  160 CONTINUE
842 *
843  ntestt = ntestt + ntest
844 *
845 * Print out tests which fail.
846 *
847  DO 170 jr = 1, ntest
848  IF( result( jr ).GE.thresh ) THEN
849 *
850 * If this is the first test to fail,
851 * print a header to the data file.
852 *
853  IF( nerrs.EQ.0 ) THEN
854  WRITE( nounit, fmt = 9997 )'CGS'
855 *
856 * Matrix types
857 *
858  WRITE( nounit, fmt = 9996 )
859  WRITE( nounit, fmt = 9995 )
860  WRITE( nounit, fmt = 9994 )'Unitary'
861 *
862 * Tests performed
863 *
864  WRITE( nounit, fmt = 9993 )'unitary', '''',
865  $ 'transpose', ( '''', j = 1, 8 )
866 *
867  END IF
868  nerrs = nerrs + 1
869  IF( result( jr ).LT.10000.0 ) THEN
870  WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
871  $ result( jr )
872  ELSE
873  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
874  $ result( jr )
875  END IF
876  END IF
877  170 CONTINUE
878 *
879  180 CONTINUE
880  190 CONTINUE
881 *
882 * Summary
883 *
884  CALL alasvm( 'CGS', nounit, nerrs, ntestt, 0 )
885 *
886  work( 1 ) = maxwrk
887 *
888  RETURN
889 *
890  9999 FORMAT( ' CDRGES3: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
891  $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
892 *
893  9998 FORMAT( ' CDRGES3: S not in Schur form at eigenvalue ', i6, '.',
894  $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
895  $ i5, ')' )
896 *
897  9997 FORMAT( / 1x, a3, ' -- Complex Generalized Schur from problem ',
898  $ 'driver' )
899 *
900  9996 FORMAT( ' Matrix types (see CDRGES3 for details): ' )
901 *
902  9995 FORMAT( ' Special Matrices:', 23x,
903  $ '(J''=transposed Jordan block)',
904  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
905  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
906  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
907  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
908  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
909  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
910  9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
911  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
912  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
913  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
914  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
915  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
916  $ '23=(small,large) 24=(small,small) 25=(large,large)',
917  $ / ' 26=random O(1) matrices.' )
918 *
919  9993 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
920  $ 'Q and Z are ', a, ',', / 19x,
921  $ 'l and r are the appropriate left and right', / 19x,
922  $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
923  $ ' means ', a, '.)', / ' Without ordering: ',
924  $ / ' 1 = | A - Q S Z', a,
925  $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
926  $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
927  $ ' | / ( n ulp ) 4 = | I - ZZ', a,
928  $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
929  $ / ' 6 = difference between (alpha,beta)',
930  $ ' and diagonals of (S,T)', / ' With ordering: ',
931  $ / ' 7 = | (A,B) - Q (S,T) Z', a, ' | / ( |(A,B)| n ulp )',
932  $ / ' 8 = | I - QQ', a,
933  $ ' | / ( n ulp ) 9 = | I - ZZ', a,
934  $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
935  $ / ' 11 = difference between (alpha,beta) and diagonals',
936  $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
937  $ 'selected eigenvalues', / )
938  9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
939  $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
940  9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
941  $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
942 *
943 * End of CDRGES3
944 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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 xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:81
subroutine clatm4(ITYPE, N, NZ1, NZ2, RSIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
CLATM4
Definition: clatm4.f:171
logical function clctes(Z, D)
CLCTES
Definition: clctes.f:58
subroutine cget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RWORK, RESULT)
CGET51
Definition: cget51.f:155
subroutine cget54(N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V, LDV, WORK, RESULT)
CGET54
Definition: cget54.f:156
complex function clarnd(IDIST, ISEED)
CLARND
Definition: clarnd.f:75
subroutine cgges3(JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK, BWORK, INFO)
CGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: cgges3.f:269
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 clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
Definition: clarfg.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
subroutine cunm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
Definition: cunm2r.f:159
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: