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

◆ ddrgev()

subroutine ddrgev ( 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 
)

DDRGEV

Purpose:
 DDRGEV checks the nonsymmetric generalized eigenvalue problem driver
 routine DGGEV.

 DGGEV 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 DDRGEV 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 DGGEV:

 (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,
          DDRGES 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, DDRGES
          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 DDRGES 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 DGGES.  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 DGGES.
[out]Q
          Q is DOUBLE PRECISION array,
                                 dimension (LDQ, max(NN))
          The (left) eigenvectors matrix computed by DGGEV.
[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 DGGES.
[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 DGGEV.
          ( 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 DGGEV 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 ddrgev.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 = 26 )
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, dggev, 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 /
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 /
471 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
472 $ 1, 1, -4, 2, -4, 8*8, 0 /
473 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
474 $ 4*5, 4*3, 1 /
475 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
476 $ 4*6, 4*4, 1 /
477 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
478 $ 2, 1 /
479 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
480 $ 2, 1 /
481 DATA ktrian / 16*0, 10*1 /
482 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
483 $ 5*2, 0 /
484 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*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( 'DDRGEV', -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.
592* KATYPE: the "type" to be passed to DLATM4 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* IASIGN: 1 if the diagonal elements of A are to be
601* multiplied by a random magnitude 1 number, =2 if
602* randomly chosen diagonal blocks are to be rotated
603* to form 2x2 blocks.
604* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
605* KTRIAN: =0: don't fill in the upper triangle, =1: do.
606* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
607* RMAGN: used to implement KAMAGN and KBMAGN.
608*
609 IF( mtypes.GT.maxtyp )
610 $ GO TO 100
611 ierr = 0
612 IF( kclass( jtype ).LT.3 ) THEN
613*
614* Generate A (w/o rotation)
615*
616 IF( abs( katype( jtype ) ).EQ.3 ) THEN
617 in = 2*( ( n-1 ) / 2 ) + 1
618 IF( in.NE.n )
619 $ CALL dlaset( 'Full', n, n, zero, zero, a, lda )
620 ELSE
621 in = n
622 END IF
623 CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
624 $ kz2( kazero( jtype ) ), iasign( jtype ),
625 $ rmagn( kamagn( jtype ) ), ulp,
626 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
627 $ iseed, a, lda )
628 iadd = kadd( kazero( jtype ) )
629 IF( iadd.GT.0 .AND. iadd.LE.n )
630 $ a( iadd, iadd ) = one
631*
632* Generate B (w/o rotation)
633*
634 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
635 in = 2*( ( n-1 ) / 2 ) + 1
636 IF( in.NE.n )
637 $ CALL dlaset( 'Full', n, n, zero, zero, b, lda )
638 ELSE
639 in = n
640 END IF
641 CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
642 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
643 $ rmagn( kbmagn( jtype ) ), one,
644 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
645 $ iseed, b, lda )
646 iadd = kadd( kbzero( jtype ) )
647 IF( iadd.NE.0 .AND. iadd.LE.n )
648 $ b( iadd, iadd ) = one
649*
650 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
651*
652* Include rotations
653*
654* Generate Q, Z as Householder transformations times
655* a diagonal matrix.
656*
657 DO 40 jc = 1, n - 1
658 DO 30 jr = jc, n
659 q( jr, jc ) = dlarnd( 3, iseed )
660 z( jr, jc ) = dlarnd( 3, iseed )
661 30 CONTINUE
662 CALL dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
663 $ work( jc ) )
664 work( 2*n+jc ) = sign( one, q( jc, jc ) )
665 q( jc, jc ) = one
666 CALL dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
667 $ work( n+jc ) )
668 work( 3*n+jc ) = sign( one, z( jc, jc ) )
669 z( jc, jc ) = one
670 40 CONTINUE
671 q( n, n ) = one
672 work( n ) = zero
673 work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
674 z( n, n ) = one
675 work( 2*n ) = zero
676 work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
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 )*work( 3*n+jc )*
683 $ a( jr, jc )
684 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
685 $ b( jr, jc )
686 50 CONTINUE
687 60 CONTINUE
688 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
689 $ lda, work( 2*n+1 ), ierr )
690 IF( ierr.NE.0 )
691 $ GO TO 90
692 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
693 $ a, lda, work( 2*n+1 ), ierr )
694 IF( ierr.NE.0 )
695 $ GO TO 90
696 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
697 $ lda, work( 2*n+1 ), ierr )
698 IF( ierr.NE.0 )
699 $ GO TO 90
700 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
701 $ b, lda, work( 2*n+1 ), ierr )
702 IF( ierr.NE.0 )
703 $ GO TO 90
704 END IF
705 ELSE
706*
707* Random matrices
708*
709 DO 80 jc = 1, n
710 DO 70 jr = 1, n
711 a( jr, jc ) = rmagn( kamagn( jtype ) )*
712 $ dlarnd( 2, iseed )
713 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
714 $ dlarnd( 2, iseed )
715 70 CONTINUE
716 80 CONTINUE
717 END IF
718*
719 90 CONTINUE
720*
721 IF( ierr.NE.0 ) THEN
722 WRITE( nounit, fmt = 9999 )'Generator', ierr, n, jtype,
723 $ ioldsd
724 info = abs( ierr )
725 RETURN
726 END IF
727*
728 100 CONTINUE
729*
730 DO 110 i = 1, 7
731 result( i ) = -one
732 110 CONTINUE
733*
734* Call DGGEV to compute eigenvalues and eigenvectors.
735*
736 CALL dlacpy( ' ', n, n, a, lda, s, lda )
737 CALL dlacpy( ' ', n, n, b, lda, t, lda )
738 CALL dggev( 'V', 'V', n, s, lda, t, lda, alphar, alphai,
739 $ beta, q, ldq, z, ldq, work, lwork, ierr )
740 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
741 result( 1 ) = ulpinv
742 WRITE( nounit, fmt = 9999 )'DGGEV1', ierr, n, jtype,
743 $ ioldsd
744 info = abs( ierr )
745 GO TO 190
746 END IF
747*
748* Do the tests (1) and (2)
749*
750 CALL dget52( .true., n, a, lda, b, lda, q, ldq, alphar,
751 $ alphai, beta, work, result( 1 ) )
752 IF( result( 2 ).GT.thresh ) THEN
753 WRITE( nounit, fmt = 9998 )'Left', 'DGGEV1',
754 $ result( 2 ), n, jtype, ioldsd
755 END IF
756*
757* Do the tests (3) and (4)
758*
759 CALL dget52( .false., n, a, lda, b, lda, z, ldq, alphar,
760 $ alphai, beta, work, result( 3 ) )
761 IF( result( 4 ).GT.thresh ) THEN
762 WRITE( nounit, fmt = 9998 )'Right', 'DGGEV1',
763 $ result( 4 ), n, jtype, ioldsd
764 END IF
765*
766* Do the test (5)
767*
768 CALL dlacpy( ' ', n, n, a, lda, s, lda )
769 CALL dlacpy( ' ', n, n, b, lda, t, lda )
770 CALL dggev( 'N', 'N', n, s, lda, t, lda, alphr1, alphi1,
771 $ beta1, q, ldq, z, ldq, work, lwork, ierr )
772 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
773 result( 1 ) = ulpinv
774 WRITE( nounit, fmt = 9999 )'DGGEV2', ierr, n, jtype,
775 $ ioldsd
776 info = abs( ierr )
777 GO TO 190
778 END IF
779*
780 DO 120 j = 1, n
781 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
782 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 5 )
783 $ = ulpinv
784 120 CONTINUE
785*
786* Do the test (6): Compute eigenvalues and left eigenvectors,
787* and test them
788*
789 CALL dlacpy( ' ', n, n, a, lda, s, lda )
790 CALL dlacpy( ' ', n, n, b, lda, t, lda )
791 CALL dggev( 'V', 'N', n, s, lda, t, lda, alphr1, alphi1,
792 $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
793 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
794 result( 1 ) = ulpinv
795 WRITE( nounit, fmt = 9999 )'DGGEV3', ierr, n, jtype,
796 $ ioldsd
797 info = abs( ierr )
798 GO TO 190
799 END IF
800*
801 DO 130 j = 1, n
802 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
803 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 6 )
804 $ = ulpinv
805 130 CONTINUE
806*
807 DO 150 j = 1, n
808 DO 140 jc = 1, n
809 IF( q( j, jc ).NE.qe( j, jc ) )
810 $ result( 6 ) = ulpinv
811 140 CONTINUE
812 150 CONTINUE
813*
814* DO the test (7): Compute eigenvalues and right eigenvectors,
815* and test them
816*
817 CALL dlacpy( ' ', n, n, a, lda, s, lda )
818 CALL dlacpy( ' ', n, n, b, lda, t, lda )
819 CALL dggev( 'N', 'V', n, s, lda, t, lda, alphr1, alphi1,
820 $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
821 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
822 result( 1 ) = ulpinv
823 WRITE( nounit, fmt = 9999 )'DGGEV4', ierr, n, jtype,
824 $ ioldsd
825 info = abs( ierr )
826 GO TO 190
827 END IF
828*
829 DO 160 j = 1, n
830 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
831 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 7 )
832 $ = ulpinv
833 160 CONTINUE
834*
835 DO 180 j = 1, n
836 DO 170 jc = 1, n
837 IF( z( j, jc ).NE.qe( j, jc ) )
838 $ result( 7 ) = ulpinv
839 170 CONTINUE
840 180 CONTINUE
841*
842* End of Loop -- Check for RESULT(j) > THRESH
843*
844 190 CONTINUE
845*
846 ntestt = ntestt + 7
847*
848* Print out tests which fail.
849*
850 DO 200 jr = 1, 7
851 IF( result( jr ).GE.thresh ) THEN
852*
853* If this is the first test to fail,
854* print a header to the data file.
855*
856 IF( nerrs.EQ.0 ) THEN
857 WRITE( nounit, fmt = 9997 )'DGV'
858*
859* Matrix types
860*
861 WRITE( nounit, fmt = 9996 )
862 WRITE( nounit, fmt = 9995 )
863 WRITE( nounit, fmt = 9994 )'Orthogonal'
864*
865* Tests performed
866*
867 WRITE( nounit, fmt = 9993 )
868*
869 END IF
870 nerrs = nerrs + 1
871 IF( result( jr ).LT.10000.0d0 ) THEN
872 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
873 $ result( jr )
874 ELSE
875 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
876 $ result( jr )
877 END IF
878 END IF
879 200 CONTINUE
880*
881 210 CONTINUE
882 220 CONTINUE
883*
884* Summary
885*
886 CALL alasvm( 'DGV', nounit, nerrs, ntestt, 0 )
887*
888 work( 1 ) = maxwrk
889*
890 RETURN
891*
892 9999 FORMAT( ' DDRGEV: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
893 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
894*
895 9998 FORMAT( ' DDRGEV: ', a, ' Eigenvectors from ', a, ' incorrectly ',
896 $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 3x,
897 $ 'N=', i4, ', JTYPE=', i3, ', ISEED=(', 4( i4, ',' ), i5,
898 $ ')' )
899*
900 9997 FORMAT( / 1x, a3, ' -- Real Generalized eigenvalue problem driver'
901 $ )
902*
903 9996 FORMAT( ' Matrix types (see DDRGEV for details): ' )
904*
905 9995 FORMAT( ' Special Matrices:', 23x,
906 $ '(J''=transposed Jordan block)',
907 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
908 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
909 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
910 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
911 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
912 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
913 9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
914 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
915 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
916 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
917 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
918 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
919 $ '23=(small,large) 24=(small,small) 25=(large,large)',
920 $ / ' 26=random O(1) matrices.' )
921*
922 9993 FORMAT( / ' Tests performed: ',
923 $ / ' 1 = max | ( b A - a B )''*l | / const.,',
924 $ / ' 2 = | |VR(i)| - 1 | / ulp,',
925 $ / ' 3 = max | ( b A - a B )*r | / const.',
926 $ / ' 4 = | |VL(i)| - 1 | / ulp,',
927 $ / ' 5 = 0 if W same no matter if r or l computed,',
928 $ / ' 6 = 0 if l same no matter if l computed,',
929 $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
930 9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
931 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
932 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
933 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
934*
935* End of DDRGEV
936*
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 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 dggev(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
DGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: dggev.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: