LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ ddrgev3()

subroutine ddrgev3 ( integer  NSIZES,
integer, dimension( * )  NN,
integer  NTYPES,
logical, dimension( * )  DOTYPE,
integer, dimension( 4 )  ISEED,
double precision  THRESH,
integer  NOUNIT,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( lda, * )  B,
double precision, dimension( lda, * )  S,
double precision, dimension( lda, * )  T,
double precision, dimension( ldq, * )  Q,
integer  LDQ,
double precision, dimension( ldq, * )  Z,
double precision, dimension( ldqe, * )  QE,
integer  LDQE,
double precision, dimension( * )  ALPHAR,
double precision, dimension( * )  ALPHAI,
double precision, dimension( * )  BETA,
double precision, dimension( * )  ALPHR1,
double precision, dimension( * )  ALPHI1,
double precision, dimension( * )  BETA1,
double precision, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RESULT,
integer  INFO 
)

DDRGEV3

Purpose:
 DDRGEV3 checks the nonsymmetric generalized eigenvalue problem driver
 routine DGGEV3.

 DGGEV3 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 DDRGEV3 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 DGGEV3:

 (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,
          DDRGEV3 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, DDRGEV3
          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 DDRGEV3 to continue the same random number
          sequence.
[in]THRESH
          THRESH is DOUBLE PRECISION
          A test will count as "failed" if the "error", computed as
          described above, exceeds THRESH.  Note that the error is
          scaled to be O(1), so THRESH should be a reasonably small
          multiple of 1, e.g., 10 or 100.  In particular, it should
          not depend on the precision (single vs. double) or the size
          of the matrix.  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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array,
                                 dimension (LDA, max(NN))
          The Schur form matrix computed from A by DGGEV3.  On exit, S
          contains the Schur form matrix corresponding to the matrix
          in A.
[out]T
          T is DOUBLE PRECISION array,
                                 dimension (LDA, max(NN))
          The upper triangular matrix computed from B by DGGEV3.
[out]Q
          Q is DOUBLE PRECISION array,
                                 dimension (LDQ, max(NN))
          The (left) eigenvectors matrix computed by DGGEV3.
[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 DOUBLE PRECISION array, dimension( LDQ, max(NN) )
          The (right) orthogonal matrix computed by DGGEV3.
[out]QE
          QE is DOUBLE PRECISION 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]ALPHAR
          ALPHAR is DOUBLE PRECISION array, dimension (max(NN))
[out]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (max(NN))
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (max(NN))

          The generalized eigenvalues of (A,B) computed by DGGEV3.
          ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
          generalized eigenvalue of A and B.
[out]ALPHR1
          ALPHR1 is DOUBLE PRECISION array, dimension (max(NN))
[out]ALPHI1
          ALPHI1 is DOUBLE PRECISION array, dimension (max(NN))
[out]BETA1
          BETA1 is DOUBLE PRECISION array, dimension (max(NN))

          Like ALPHAR, ALPHAI, BETA, these arrays contain the
          eigenvalues of A and B, but those computed when DGGEV3 only
          computes a partial eigendecomposition, i.e. not the
          eigenvalues and left and right eigenvectors.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK.  LWORK >= MAX( 8*N, N*(N+1) ).
[out]RESULT
          RESULT is DOUBLE PRECISION 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 404 of file ddrgev3.f.

408 *
409 * -- LAPACK test routine --
410 * -- LAPACK is a software package provided by Univ. of Tennessee, --
411 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
412 *
413 * .. Scalar Arguments ..
414  INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
415  $ NTYPES
416  DOUBLE PRECISION THRESH
417 * ..
418 * .. Array Arguments ..
419  LOGICAL DOTYPE( * )
420  INTEGER ISEED( 4 ), NN( * )
421  DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
422  $ ALPHI1( * ), ALPHR1( * ), B( LDA, * ),
423  $ BETA( * ), BETA1( * ), Q( LDQ, * ),
424  $ QE( LDQE, * ), RESULT( * ), S( LDA, * ),
425  $ T( LDA, * ), WORK( * ), Z( LDQ, * )
426 * ..
427 *
428 * =====================================================================
429 *
430 * .. Parameters ..
431  DOUBLE PRECISION ZERO, ONE
432  parameter( zero = 0.0d+0, one = 1.0d+0 )
433  INTEGER MAXTYP
434  parameter( maxtyp = 27 )
435 * ..
436 * .. Local Scalars ..
437  LOGICAL BADNN
438  INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
439  $ MAXWRK, MINWRK, MTYPES, N, N1, NERRS, NMATS,
440  $ NMAX, NTESTT
441  DOUBLE PRECISION SAFMAX, SAFMIN, ULP, ULPINV
442 * ..
443 * .. Local Arrays ..
444  INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
445  $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
446  $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
447  $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
448  $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
449  $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
450  DOUBLE PRECISION RMAGN( 0: 3 )
451 * ..
452 * .. External Functions ..
453  INTEGER ILAENV
454  DOUBLE PRECISION DLAMCH, DLARND
455  EXTERNAL ilaenv, dlamch, dlarnd
456 * ..
457 * .. External Subroutines ..
458  EXTERNAL alasvm, dget52, dggev3, dlabad, dlacpy, dlarfg,
460 * ..
461 * .. Intrinsic Functions ..
462  INTRINSIC abs, dble, max, min, sign
463 * ..
464 * .. Data statements ..
465  DATA kclass / 15*1, 10*2, 1*3, 1*4 /
466  DATA kz1 / 0, 1, 2, 1, 3, 3 /
467  DATA kz2 / 0, 0, 1, 2, 1, 1 /
468  DATA kadd / 0, 0, 0, 0, 3, 2 /
469  DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
470  $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0, 0 /
471  DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
472  $ 1, 1, -4, 2, -4, 8*8, 0, 0 /
473  DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
474  $ 4*5, 4*3, 1, 1 /
475  DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
476  $ 4*6, 4*4, 1, 1 /
477  DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
478  $ 2, 1, 3 /
479  DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
480  $ 2, 1, 3 /
481  DATA ktrian / 16*0, 11*1 /
482  DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
483  $ 5*2, 2*0 /
484  DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 10*0 /
485 * ..
486 * .. Executable Statements ..
487 *
488 * Check for errors
489 *
490  info = 0
491 *
492  badnn = .false.
493  nmax = 1
494  DO 10 j = 1, nsizes
495  nmax = max( nmax, nn( j ) )
496  IF( nn( j ).LT.0 )
497  $ badnn = .true.
498  10 CONTINUE
499 *
500  IF( nsizes.LT.0 ) THEN
501  info = -1
502  ELSE IF( badnn ) THEN
503  info = -2
504  ELSE IF( ntypes.LT.0 ) THEN
505  info = -3
506  ELSE IF( thresh.LT.zero ) THEN
507  info = -6
508  ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
509  info = -9
510  ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
511  info = -14
512  ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax ) THEN
513  info = -17
514  END IF
515 *
516 * Compute workspace
517 * (Note: Comments in the code beginning "Workspace:" describe the
518 * minimal amount of workspace needed at that point in the code,
519 * as well as the preferred amount for good performance.
520 * NB refers to the optimal block size for the immediately
521 * following subroutine, as returned by ILAENV.
522 *
523  minwrk = 1
524  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
525  minwrk = max( 1, 8*nmax, nmax*( nmax+1 ) )
526  maxwrk = 7*nmax + nmax*ilaenv( 1, 'DGEQRF', ' ', nmax, 1, nmax,
527  $ 0 )
528  maxwrk = max( maxwrk, nmax*( nmax+1 ) )
529  work( 1 ) = maxwrk
530  END IF
531 *
532  IF( lwork.LT.minwrk )
533  $ info = -25
534 *
535  IF( info.NE.0 ) THEN
536  CALL xerbla( 'DDRGEV3', -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  safmin = dlamch( 'Safe minimum' )
546  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
547  safmin = safmin / ulp
548  safmax = one / safmin
549  CALL dlabad( 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 / dble( 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 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
591 * =3 means random, =4 means random generalized
592 * upper Hessenberg.
593 * KATYPE: the "type" to be passed to DLATM4 for computing A.
594 * KAZERO: the pattern of zeros on the diagonal for A:
595 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
596 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
597 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
598 * non-zero entries.)
599 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
600 * =2: large, =3: small.
601 * IASIGN: 1 if the diagonal elements of A are to be
602 * multiplied by a random magnitude 1 number, =2 if
603 * randomly chosen diagonal blocks are to be rotated
604 * to form 2x2 blocks.
605 * KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
606 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
607 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
608 * RMAGN: used to implement KAMAGN and KBMAGN.
609 *
610  IF( mtypes.GT.maxtyp )
611  $ GO TO 100
612  ierr = 0
613  IF( kclass( jtype ).LT.3 ) THEN
614 *
615 * Generate A (w/o rotation)
616 *
617  IF( abs( katype( jtype ) ).EQ.3 ) THEN
618  in = 2*( ( n-1 ) / 2 ) + 1
619  IF( in.NE.n )
620  $ CALL dlaset( 'Full', n, n, zero, zero, a, lda )
621  ELSE
622  in = n
623  END IF
624  CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
625  $ kz2( kazero( jtype ) ), iasign( jtype ),
626  $ rmagn( kamagn( jtype ) ), ulp,
627  $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
628  $ iseed, a, lda )
629  iadd = kadd( kazero( jtype ) )
630  IF( iadd.GT.0 .AND. iadd.LE.n )
631  $ a( iadd, iadd ) = one
632 *
633 * Generate B (w/o rotation)
634 *
635  IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
636  in = 2*( ( n-1 ) / 2 ) + 1
637  IF( in.NE.n )
638  $ CALL dlaset( 'Full', n, n, zero, zero, b, lda )
639  ELSE
640  in = n
641  END IF
642  CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
643  $ kz2( kbzero( jtype ) ), ibsign( jtype ),
644  $ rmagn( kbmagn( jtype ) ), one,
645  $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
646  $ iseed, b, lda )
647  iadd = kadd( kbzero( jtype ) )
648  IF( iadd.NE.0 .AND. iadd.LE.n )
649  $ b( iadd, iadd ) = one
650 *
651  IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
652 *
653 * Include rotations
654 *
655 * Generate Q, Z as Householder transformations times
656 * a diagonal matrix.
657 *
658  DO 40 jc = 1, n - 1
659  DO 30 jr = jc, n
660  q( jr, jc ) = dlarnd( 3, iseed )
661  z( jr, jc ) = dlarnd( 3, iseed )
662  30 CONTINUE
663  CALL dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
664  $ work( jc ) )
665  work( 2*n+jc ) = sign( one, q( jc, jc ) )
666  q( jc, jc ) = one
667  CALL dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
668  $ work( n+jc ) )
669  work( 3*n+jc ) = sign( one, z( jc, jc ) )
670  z( jc, jc ) = one
671  40 CONTINUE
672  q( n, n ) = one
673  work( n ) = zero
674  work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
675  z( n, n ) = one
676  work( 2*n ) = zero
677  work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
678 *
679 * Apply the diagonal matrices
680 *
681  DO 60 jc = 1, n
682  DO 50 jr = 1, n
683  a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
684  $ a( jr, jc )
685  b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
686  $ b( jr, jc )
687  50 CONTINUE
688  60 CONTINUE
689  CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
690  $ lda, work( 2*n+1 ), ierr )
691  IF( ierr.NE.0 )
692  $ GO TO 90
693  CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
694  $ a, lda, work( 2*n+1 ), ierr )
695  IF( ierr.NE.0 )
696  $ GO TO 90
697  CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
698  $ lda, work( 2*n+1 ), ierr )
699  IF( ierr.NE.0 )
700  $ GO TO 90
701  CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
702  $ b, lda, work( 2*n+1 ), ierr )
703  IF( ierr.NE.0 )
704  $ GO TO 90
705  END IF
706  ELSE IF (kclass( jtype ).EQ.3) THEN
707 *
708 * Random matrices
709 *
710  DO 80 jc = 1, n
711  DO 70 jr = 1, n
712  a( jr, jc ) = rmagn( kamagn( jtype ) )*
713  $ dlarnd( 2, iseed )
714  b( jr, jc ) = rmagn( kbmagn( jtype ) )*
715  $ dlarnd( 2, iseed )
716  70 CONTINUE
717  80 CONTINUE
718  ELSE
719 *
720 * Random upper Hessenberg pencil with singular B
721 *
722  DO 81 jc = 1, n
723  DO 71 jr = 1, min( jc + 1, n)
724  a( jr, jc ) = rmagn( kamagn( jtype ) )*
725  $ dlarnd( 2, iseed )
726  71 CONTINUE
727  DO 72 jr = jc + 2, n
728  a( jr, jc ) = zero
729  72 CONTINUE
730  81 CONTINUE
731  DO 82 jc = 1, n
732  DO 73 jr = 1, jc
733  b( jr, jc ) = rmagn( kamagn( jtype ) )*
734  $ dlarnd( 2, iseed )
735  73 CONTINUE
736  DO 74 jr = jc + 1, n
737  b( jr, jc ) = zero
738  74 CONTINUE
739  82 CONTINUE
740  DO 83 jc = 1, n, 4
741  b( jc, jc ) = zero
742  83 CONTINUE
743 
744  END IF
745 *
746  90 CONTINUE
747 *
748  IF( ierr.NE.0 ) THEN
749  WRITE( nounit, fmt = 9999 )'Generator', ierr, n, jtype,
750  $ ioldsd
751  info = abs( ierr )
752  RETURN
753  END IF
754 *
755  100 CONTINUE
756 *
757  DO 110 i = 1, 7
758  result( i ) = -one
759  110 CONTINUE
760 *
761 * Call XLAENV to set the parameters used in DLAQZ0
762 *
763  CALL xlaenv( 12, 10 )
764  CALL xlaenv( 13, 12 )
765  CALL xlaenv( 14, 13 )
766  CALL xlaenv( 15, 2 )
767  CALL xlaenv( 17, 10 )
768 *
769 * Call DGGEV3 to compute eigenvalues and eigenvectors.
770 *
771  CALL dlacpy( ' ', n, n, a, lda, s, lda )
772  CALL dlacpy( ' ', n, n, b, lda, t, lda )
773  CALL dggev3( 'V', 'V', n, s, lda, t, lda, alphar, alphai,
774  $ beta, q, ldq, z, ldq, work, lwork, ierr )
775  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
776  result( 1 ) = ulpinv
777  WRITE( nounit, fmt = 9999 )'DGGEV31', ierr, n, jtype,
778  $ ioldsd
779  info = abs( ierr )
780  GO TO 190
781  END IF
782 *
783 * Do the tests (1) and (2)
784 *
785  CALL dget52( .true., n, a, lda, b, lda, q, ldq, alphar,
786  $ alphai, beta, work, result( 1 ) )
787  IF( result( 2 ).GT.thresh ) THEN
788  WRITE( nounit, fmt = 9998 )'Left', 'DGGEV31',
789  $ result( 2 ), n, jtype, ioldsd
790  END IF
791 *
792 * Do the tests (3) and (4)
793 *
794  CALL dget52( .false., n, a, lda, b, lda, z, ldq, alphar,
795  $ alphai, beta, work, result( 3 ) )
796  IF( result( 4 ).GT.thresh ) THEN
797  WRITE( nounit, fmt = 9998 )'Right', 'DGGEV31',
798  $ result( 4 ), n, jtype, ioldsd
799  END IF
800 *
801 * Do the test (5)
802 *
803  CALL dlacpy( ' ', n, n, a, lda, s, lda )
804  CALL dlacpy( ' ', n, n, b, lda, t, lda )
805  CALL dggev3( 'N', 'N', n, s, lda, t, lda, alphr1, alphi1,
806  $ beta1, q, ldq, z, ldq, work, lwork, ierr )
807  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
808  result( 1 ) = ulpinv
809  WRITE( nounit, fmt = 9999 )'DGGEV32', ierr, n, jtype,
810  $ ioldsd
811  info = abs( ierr )
812  GO TO 190
813  END IF
814 *
815  DO 120 j = 1, n
816  IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
817  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 5 )
818  $ = ulpinv
819  120 CONTINUE
820 *
821 * Do the test (6): Compute eigenvalues and left eigenvectors,
822 * and test them
823 *
824  CALL dlacpy( ' ', n, n, a, lda, s, lda )
825  CALL dlacpy( ' ', n, n, b, lda, t, lda )
826  CALL dggev3( 'V', 'N', n, s, lda, t, lda, alphr1, alphi1,
827  $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
828  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
829  result( 1 ) = ulpinv
830  WRITE( nounit, fmt = 9999 )'DGGEV33', ierr, n, jtype,
831  $ ioldsd
832  info = abs( ierr )
833  GO TO 190
834  END IF
835 *
836  DO 130 j = 1, n
837  IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
838  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 6 )
839  $ = ulpinv
840  130 CONTINUE
841 *
842  DO 150 j = 1, n
843  DO 140 jc = 1, n
844  IF( q( j, jc ).NE.qe( j, jc ) )
845  $ result( 6 ) = ulpinv
846  140 CONTINUE
847  150 CONTINUE
848 *
849 * DO the test (7): Compute eigenvalues and right eigenvectors,
850 * and test them
851 *
852  CALL dlacpy( ' ', n, n, a, lda, s, lda )
853  CALL dlacpy( ' ', n, n, b, lda, t, lda )
854  CALL dggev3( 'N', 'V', n, s, lda, t, lda, alphr1, alphi1,
855  $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
856  IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
857  result( 1 ) = ulpinv
858  WRITE( nounit, fmt = 9999 )'DGGEV34', ierr, n, jtype,
859  $ ioldsd
860  info = abs( ierr )
861  GO TO 190
862  END IF
863 *
864  DO 160 j = 1, n
865  IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
866  $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 7 )
867  $ = ulpinv
868  160 CONTINUE
869 *
870  DO 180 j = 1, n
871  DO 170 jc = 1, n
872  IF( z( j, jc ).NE.qe( j, jc ) )
873  $ result( 7 ) = ulpinv
874  170 CONTINUE
875  180 CONTINUE
876 *
877 * End of Loop -- Check for RESULT(j) > THRESH
878 *
879  190 CONTINUE
880 *
881  ntestt = ntestt + 7
882 *
883 * Print out tests which fail.
884 *
885  DO 200 jr = 1, 7
886  IF( result( jr ).GE.thresh ) THEN
887 *
888 * If this is the first test to fail,
889 * print a header to the data file.
890 *
891  IF( nerrs.EQ.0 ) THEN
892  WRITE( nounit, fmt = 9997 )'DGV'
893 *
894 * Matrix types
895 *
896  WRITE( nounit, fmt = 9996 )
897  WRITE( nounit, fmt = 9995 )
898  WRITE( nounit, fmt = 9994 )'Orthogonal'
899 *
900 * Tests performed
901 *
902  WRITE( nounit, fmt = 9993 )
903 *
904  END IF
905  nerrs = nerrs + 1
906  IF( result( jr ).LT.10000.0d0 ) THEN
907  WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
908  $ result( jr )
909  ELSE
910  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
911  $ result( jr )
912  END IF
913  END IF
914  200 CONTINUE
915 *
916  210 CONTINUE
917  220 CONTINUE
918 *
919 * Summary
920 *
921  CALL alasvm( 'DGV', nounit, nerrs, ntestt, 0 )
922 *
923  work( 1 ) = maxwrk
924 *
925  RETURN
926 *
927  9999 FORMAT( ' DDRGEV3: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
928  $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
929 *
930  9998 FORMAT( ' DDRGEV3: ', a, ' Eigenvectors from ', a,
931  $ ' incorrectly normalized.', / ' Bits of error=', 0p, g10.3,
932  $ ',', 3x, 'N=', i4, ', JTYPE=', i3, ', ISEED=(',
933  $ 4( i4, ',' ), i5, ')' )
934 *
935  9997 FORMAT( / 1x, a3, ' -- Real Generalized eigenvalue problem driver'
936  $ )
937 *
938  9996 FORMAT( ' Matrix types (see DDRGEV3 for details): ' )
939 *
940  9995 FORMAT( ' Special Matrices:', 23x,
941  $ '(J''=transposed Jordan block)',
942  $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
943  $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
944  $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
945  $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
946  $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
947  $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
948  9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
949  $ / ' 16=Transposed Jordan Blocks 19=geometric ',
950  $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
951  $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
952  $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
953  $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
954  $ '23=(small,large) 24=(small,small) 25=(large,large)',
955  $ / ' 26=random O(1) matrices.' )
956 *
957  9993 FORMAT( / ' Tests performed: ',
958  $ / ' 1 = max | ( b A - a B )''*l | / const.,',
959  $ / ' 2 = | |VR(i)| - 1 | / ulp,',
960  $ / ' 3 = max | ( b A - a B )*r | / const.',
961  $ / ' 4 = | |VL(i)| - 1 | / ulp,',
962  $ / ' 5 = 0 if W same no matter if r or l computed,',
963  $ / ' 6 = 0 if l same no matter if l computed,',
964  $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
965  9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
966  $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
967  9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
968  $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
969 *
970 * End of DDRGEV3
971 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:73
subroutine xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:81
subroutine dlatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
DLATM4
Definition: dlatm4.f:175
subroutine dget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
DGET52
Definition: dget52.f:199
double precision function dlarnd(IDIST, ISEED)
DLARND
Definition: dlarnd.f:73
subroutine dggev3(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
DGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (...
Definition: dggev3.f:226
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
Definition: dlarfg.f:106
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: dorm2r.f:159
Here is the call graph for this function:
Here is the caller graph for this function: