LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cdrgev3()

subroutine cdrgev3 ( 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( ldqe, * )  QE,
integer  LDQE,
complex, dimension( * )  ALPHA,
complex, dimension( * )  BETA,
complex, dimension( * )  ALPHA1,
complex, dimension( * )  BETA1,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
real, dimension( * )  RESULT,
integer  INFO 
)

CDRGEV3

Purpose:
 CDRGEV3 checks the nonsymmetric generalized eigenvalue problem driver
 routine CGGEV3.

 CGGEV3 computes for a pair of n-by-n nonsymmetric matrices (A,B) the
 generalized eigenvalues and, optionally, the left and right
 eigenvectors.

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
 or a ratio  alpha/beta = w, such that A - w*B is singular.  It is
 usually represented as the pair (alpha,beta), as there is reasonable
 interpretation for beta=0, and even for both being zero.

 A right generalized eigenvector corresponding to a generalized
 eigenvalue  w  for a pair of matrices (A,B) is a vector r  such that
 (A - wB) * r = 0.  A left generalized eigenvector is a vector l such
 that l**H * (A - wB) = 0, where l**H is the conjugate-transpose of l.

 When CDRGEV3 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 tests
 will be performed and compared with the threshold THRESH.

 Results from CGGEV3:

 (1)  max over all left eigenvalue/-vector pairs (alpha/beta,l) of

      | VL**H * (beta A - alpha B) |/( ulp max(|beta A|, |alpha B|) )

      where VL**H is the conjugate-transpose of VL.

 (2)  | |VL(i)| - 1 | / ulp and whether largest component real

      VL(i) denotes the i-th column of VL.

 (3)  max over all left eigenvalue/-vector pairs (alpha/beta,r) of

      | (beta A - alpha B) * VR | / ( ulp max(|beta A|, |alpha B|) )

 (4)  | |VR(i)| - 1 | / ulp and whether largest component real

      VR(i) denotes the i-th column of VR.

 (5)  W(full) = W(partial)
      W(full) denotes the eigenvalues computed when both l and r
      are also computed, and W(partial) denotes the eigenvalues
      computed when only W, only W and r, or only W and l are
      computed.

 (6)  VL(full) = VL(partial)
      VL(full) denotes the left eigenvectors computed when both l
      and r are computed, and VL(partial) denotes the result
      when only l is computed.

 (7)  VR(full) = VR(partial)
      VR(full) denotes the right eigenvectors computed when both l
      and r are also computed, and VR(partial) denotes the result
      when only l is computed.


 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,
          CDRGEV3 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, CDRGEV3
          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.  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 CDRGEV3 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.  It must be at least zero.
[in]NOUNIT
          NOUNIT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns IERR 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 CGGEV3.  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 CGGEV3.
[out]Q
          Q is COMPLEX array, dimension (LDQ, max(NN))
          The (left) eigenvectors matrix computed by CGGEV3.
[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 CGGEV3.
[out]QE
          QE is COMPLEX array, dimension( LDQ, max(NN) )
          QE holds the computed right or left eigenvectors.
[in]LDQE
          LDQE is INTEGER
          The leading dimension of QE. LDQE >= max(1,max(NN)).
[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 CGGEV3.
          ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
          generalized eigenvalue of A and B.
[out]ALPHA1
          ALPHA1 is COMPLEX array, dimension (max(NN))
[out]BETA1
          BETA1 is COMPLEX array, dimension (max(NN))

          Like ALPHAR, ALPHAI, BETA, these arrays contain the
          eigenvalues of A and B, but those computed when CGGEV3 only
          computes a partial eigendecomposition, i.e. not the
          eigenvalues and left and right eigenvectors.
[out]WORK
          WORK is COMPLEX array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK.  LWORK >= N*(N+1)
[out]RWORK
          RWORK is REAL array, dimension (8*N)
          Real workspace.
[out]RESULT
          RESULT is REAL array, dimension (2)
          The values computed by the tests described above.
          The values are currently limited to 1/ulp, to avoid overflow.
[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 395 of file cdrgev3.f.

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