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

◆ 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 ulpinv = one / ulp
550*
551* The values RMAGN(2:3) depend on N, see below.
552*
553 rmagn( 0 ) = zero
554 rmagn( 1 ) = one
555*
556* Loop over sizes, types
557*
558 ntestt = 0
559 nerrs = 0
560 nmats = 0
561*
562 DO 220 jsize = 1, nsizes
563 n = nn( jsize )
564 n1 = max( 1, n )
565 rmagn( 2 ) = safmax*ulp / real( n1 )
566 rmagn( 3 ) = safmin*ulpinv*n1
567*
568 IF( nsizes.NE.1 ) THEN
569 mtypes = min( maxtyp, ntypes )
570 ELSE
571 mtypes = min( maxtyp+1, ntypes )
572 END IF
573*
574 DO 210 jtype = 1, mtypes
575 IF( .NOT.dotype( jtype ) )
576 $ GO TO 210
577 nmats = nmats + 1
578*
579* Save ISEED in case of an error.
580*
581 DO 20 j = 1, 4
582 ioldsd( j ) = iseed( j )
583 20 CONTINUE
584*
585* Generate test matrices A and B
586*
587* Description of control parameters:
588*
589* KCLASS: =1 means w/o rotation, =2 means w/ rotation,
590* =3 means random.
591* KATYPE: the "type" to be passed to CLATM4 for computing A.
592* KAZERO: the pattern of zeros on the diagonal for A:
593* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
594* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
595* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
596* non-zero entries.)
597* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
598* =2: large, =3: small.
599* LASIGN: .TRUE. if the diagonal elements of A are to be
600* multiplied by a random magnitude 1 number.
601* KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
602* KTRIAN: =0: don't fill in the upper triangle, =1: do.
603* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
604* RMAGN: used to implement KAMAGN and KBMAGN.
605*
606 IF( mtypes.GT.maxtyp )
607 $ GO TO 100
608 ierr = 0
609 IF( kclass( jtype ).LT.3 ) THEN
610*
611* Generate A (w/o rotation)
612*
613 IF( abs( katype( jtype ) ).EQ.3 ) THEN
614 in = 2*( ( n-1 ) / 2 ) + 1
615 IF( in.NE.n )
616 $ CALL claset( 'Full', n, n, czero, czero, a, lda )
617 ELSE
618 in = n
619 END IF
620 CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
621 $ kz2( kazero( jtype ) ), lasign( jtype ),
622 $ rmagn( kamagn( jtype ) ), ulp,
623 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
624 $ iseed, a, lda )
625 iadd = kadd( kazero( jtype ) )
626 IF( iadd.GT.0 .AND. iadd.LE.n )
627 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
628*
629* Generate B (w/o rotation)
630*
631 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
632 in = 2*( ( n-1 ) / 2 ) + 1
633 IF( in.NE.n )
634 $ CALL claset( 'Full', n, n, czero, czero, b, lda )
635 ELSE
636 in = n
637 END IF
638 CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
639 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
640 $ rmagn( kbmagn( jtype ) ), one,
641 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
642 $ iseed, b, lda )
643 iadd = kadd( kbzero( jtype ) )
644 IF( iadd.NE.0 .AND. iadd.LE.n )
645 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
646*
647 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
648*
649* Include rotations
650*
651* Generate Q, Z as Householder transformations times
652* a diagonal matrix.
653*
654 DO 40 jc = 1, n - 1
655 DO 30 jr = jc, n
656 q( jr, jc ) = clarnd( 3, iseed )
657 z( jr, jc ) = clarnd( 3, iseed )
658 30 CONTINUE
659 CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
660 $ work( jc ) )
661 work( 2*n+jc ) = sign( one, real( q( jc, jc ) ) )
662 q( jc, jc ) = cone
663 CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
664 $ work( n+jc ) )
665 work( 3*n+jc ) = sign( one, real( z( jc, jc ) ) )
666 z( jc, jc ) = cone
667 40 CONTINUE
668 ctemp = clarnd( 3, iseed )
669 q( n, n ) = cone
670 work( n ) = czero
671 work( 3*n ) = ctemp / abs( ctemp )
672 ctemp = clarnd( 3, iseed )
673 z( n, n ) = cone
674 work( 2*n ) = czero
675 work( 4*n ) = ctemp / abs( ctemp )
676*
677* Apply the diagonal matrices
678*
679 DO 60 jc = 1, n
680 DO 50 jr = 1, n
681 a( jr, jc ) = work( 2*n+jr )*
682 $ conjg( work( 3*n+jc ) )*
683 $ a( jr, jc )
684 b( jr, jc ) = work( 2*n+jr )*
685 $ conjg( work( 3*n+jc ) )*
686 $ b( jr, jc )
687 50 CONTINUE
688 60 CONTINUE
689 CALL cunm2r( '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 cunm2r( 'R', 'C', 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 cunm2r( '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 cunm2r( 'R', 'C', 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
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 $ clarnd( 4, iseed )
714 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
715 $ clarnd( 4, iseed )
716 70 CONTINUE
717 80 CONTINUE
718 END IF
719*
720 90 CONTINUE
721*
722 IF( ierr.NE.0 ) THEN
723 WRITE( nounit, fmt = 9999 )'Generator', ierr, n, jtype,
724 $ ioldsd
725 info = abs( ierr )
726 RETURN
727 END IF
728*
729 100 CONTINUE
730*
731 DO 110 i = 1, 7
732 result( i ) = -one
733 110 CONTINUE
734*
735* Call XLAENV to set the parameters used in CLAQZ0
736*
737 CALL xlaenv( 12, 10 )
738 CALL xlaenv( 13, 12 )
739 CALL xlaenv( 14, 13 )
740 CALL xlaenv( 15, 2 )
741 CALL xlaenv( 17, 10 )
742*
743* Call CGGEV3 to compute eigenvalues and eigenvectors.
744*
745 CALL clacpy( ' ', n, n, a, lda, s, lda )
746 CALL clacpy( ' ', n, n, b, lda, t, lda )
747 CALL cggev3( 'V', 'V', n, s, lda, t, lda, alpha, beta, q,
748 $ ldq, z, ldq, work, lwork, rwork, ierr )
749 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
750 result( 1 ) = ulpinv
751 WRITE( nounit, fmt = 9999 )'CGGEV31', ierr, n, jtype,
752 $ ioldsd
753 info = abs( ierr )
754 GO TO 190
755 END IF
756*
757* Do the tests (1) and (2)
758*
759 CALL cget52( .true., n, a, lda, b, lda, q, ldq, alpha, beta,
760 $ work, rwork, result( 1 ) )
761 IF( result( 2 ).GT.thresh ) THEN
762 WRITE( nounit, fmt = 9998 )'Left', 'CGGEV31',
763 $ result( 2 ), n, jtype, ioldsd
764 END IF
765*
766* Do the tests (3) and (4)
767*
768 CALL cget52( .false., n, a, lda, b, lda, z, ldq, alpha,
769 $ beta, work, rwork, result( 3 ) )
770 IF( result( 4 ).GT.thresh ) THEN
771 WRITE( nounit, fmt = 9998 )'Right', 'CGGEV31',
772 $ result( 4 ), n, jtype, ioldsd
773 END IF
774*
775* Do test (5)
776*
777 CALL clacpy( ' ', n, n, a, lda, s, lda )
778 CALL clacpy( ' ', n, n, b, lda, t, lda )
779 CALL cggev3( 'N', 'N', n, s, lda, t, lda, alpha1, beta1, q,
780 $ ldq, z, ldq, work, lwork, rwork, ierr )
781 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
782 result( 1 ) = ulpinv
783 WRITE( nounit, fmt = 9999 )'CGGEV32', ierr, n, jtype,
784 $ ioldsd
785 info = abs( ierr )
786 GO TO 190
787 END IF
788*
789 DO 120 j = 1, n
790 IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
791 $ beta1( j ) ) result( 5 ) = ulpinv
792 120 CONTINUE
793*
794* Do the test (6): Compute eigenvalues and left eigenvectors,
795* and test them
796*
797 CALL clacpy( ' ', n, n, a, lda, s, lda )
798 CALL clacpy( ' ', n, n, b, lda, t, lda )
799 CALL cggev3( 'V', 'N', n, s, lda, t, lda, alpha1, beta1, qe,
800 $ ldqe, z, ldq, work, lwork, rwork, ierr )
801 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
802 result( 1 ) = ulpinv
803 WRITE( nounit, fmt = 9999 )'CGGEV33', ierr, n, jtype,
804 $ ioldsd
805 info = abs( ierr )
806 GO TO 190
807 END IF
808
809*
810 DO 130 j = 1, n
811 IF( alpha( j ).NE.alpha1( j ) .OR.
812 $ beta( j ).NE.beta1( j ) ) THEN
813 result( 6 ) = ulpinv
814 ENDIF
815 130 CONTINUE
816*
817 DO 150 j = 1, n
818 DO 140 jc = 1, n
819 IF( q( j, jc ).NE.qe( j, jc ) ) THEN
820 result( 6 ) = ulpinv
821 END IF
822 140 CONTINUE
823 150 CONTINUE
824*
825* DO the test (7): Compute eigenvalues and right eigenvectors,
826* and test them
827*
828 CALL clacpy( ' ', n, n, a, lda, s, lda )
829 CALL clacpy( ' ', n, n, b, lda, t, lda )
830 CALL cggev3( 'N', 'V', n, s, lda, t, lda, alpha1, beta1, q,
831 $ ldq, qe, ldqe, work, lwork, rwork, ierr )
832 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
833 result( 1 ) = ulpinv
834 WRITE( nounit, fmt = 9999 )'CGGEV34', ierr, n, jtype,
835 $ ioldsd
836 info = abs( ierr )
837 GO TO 190
838 END IF
839*
840 DO 160 j = 1, n
841 IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
842 $ beta1( j ) )result( 7 ) = ulpinv
843 160 CONTINUE
844*
845 DO 180 j = 1, n
846 DO 170 jc = 1, n
847 IF( z( j, jc ).NE.qe( j, jc ) )
848 $ result( 7 ) = ulpinv
849 170 CONTINUE
850 180 CONTINUE
851*
852* End of Loop -- Check for RESULT(j) > THRESH
853*
854 190 CONTINUE
855*
856 ntestt = ntestt + 7
857*
858* Print out tests which fail.
859*
860 DO 200 jr = 1, 7
861 IF( result( jr ).GE.thresh ) THEN
862*
863* If this is the first test to fail,
864* print a header to the data file.
865*
866 IF( nerrs.EQ.0 ) THEN
867 WRITE( nounit, fmt = 9997 )'CGV'
868*
869* Matrix types
870*
871 WRITE( nounit, fmt = 9996 )
872 WRITE( nounit, fmt = 9995 )
873 WRITE( nounit, fmt = 9994 )'Orthogonal'
874*
875* Tests performed
876*
877 WRITE( nounit, fmt = 9993 )
878*
879 END IF
880 nerrs = nerrs + 1
881 IF( result( jr ).LT.10000.0 ) THEN
882 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
883 $ result( jr )
884 ELSE
885 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
886 $ result( jr )
887 END IF
888 END IF
889 200 CONTINUE
890*
891 210 CONTINUE
892 220 CONTINUE
893*
894* Summary
895*
896 CALL alasvm( 'CGV3', nounit, nerrs, ntestt, 0 )
897*
898 work( 1 ) = maxwrk
899*
900 RETURN
901*
902 9999 FORMAT( ' CDRGEV3: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
903 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
904*
905 9998 FORMAT( ' CDRGEV3: ', a, ' Eigenvectors from ', a,
906 $ ' incorrectly normalized.', / ' Bits of error=', 0p, g10.3,
907 $ ',', 3x, 'N=', i4, ', JTYPE=', i3, ', ISEED=(',
908 $ 3( i4, ',' ), i5, ')' )
909*
910 9997 FORMAT( / 1x, a3, ' -- Complex Generalized eigenvalue problem ',
911 $ 'driver' )
912*
913 9996 FORMAT( ' Matrix types (see CDRGEV3 for details): ' )
914*
915 9995 FORMAT( ' Special Matrices:', 23x,
916 $ '(J''=transposed Jordan block)',
917 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
918 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
919 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
920 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
921 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
922 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
923 9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
924 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
925 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
926 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
927 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
928 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
929 $ '23=(small,large) 24=(small,small) 25=(large,large)',
930 $ / ' 26=random O(1) matrices.' )
931*
932 9993 FORMAT( / ' Tests performed: ',
933 $ / ' 1 = max | ( b A - a B )''*l | / const.,',
934 $ / ' 2 = | |VR(i)| - 1 | / ulp,',
935 $ / ' 3 = max | ( b A - a B )*r | / const.',
936 $ / ' 4 = | |VL(i)| - 1 | / ulp,',
937 $ / ' 5 = 0 if W same no matter if r or l computed,',
938 $ / ' 6 = 0 if l same no matter if l computed,',
939 $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
940 9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
941 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
942 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
943 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
944*
945* End of CDRGEV3
946*
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 clatm4(itype, n, nz1, nz2, rsign, amagn, rcond, triang, idist, iseed, a, lda)
CLATM4
Definition clatm4.f:171
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
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine clarfg(n, alpha, x, incx, tau)
CLARFG generates an elementary reflector (Householder matrix).
Definition clarfg.f:106
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 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
Here is the call graph for this function:
Here is the caller graph for this function: