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

◆ sdrgev3()

subroutine sdrgev3 ( integer  nsizes,
integer, dimension( * )  nn,
integer  ntypes,
logical, dimension( * )  dotype,
integer, dimension( 4 )  iseed,
real  thresh,
integer  nounit,
real, dimension( lda, * )  a,
integer  lda,
real, dimension( lda, * )  b,
real, dimension( lda, * )  s,
real, dimension( lda, * )  t,
real, dimension( ldq, * )  q,
integer  ldq,
real, dimension( ldq, * )  z,
real, dimension( ldqe, * )  qe,
integer  ldqe,
real, dimension( * )  alphar,
real, dimension( * )  alphai,
real, dimension( * )  beta,
real, dimension( * )  alphr1,
real, dimension( * )  alphi1,
real, dimension( * )  beta1,
real, dimension( * )  work,
integer  lwork,
real, dimension( * )  result,
integer  info 
)

SDRGEV3

Purpose:
 SDRGEV3 checks the nonsymmetric generalized eigenvalue problem driver
 routine SGGEV3.

 SGGEV3 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 SDRGEV3 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 SGGEV3:

 (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,
          SDRGEV3 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, SDRGEV3
          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 SDRGEV3 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 REAL 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 REAL 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 REAL array,
                                 dimension (LDA, max(NN))
          The Schur form matrix computed from A by SGGEV3.  On exit, S
          contains the Schur form matrix corresponding to the matrix
          in A.
[out]T
          T is REAL array,
                                 dimension (LDA, max(NN))
          The upper triangular matrix computed from B by SGGEV3.
[out]Q
          Q is REAL array,
                                 dimension (LDQ, max(NN))
          The (left) eigenvectors matrix computed by SGGEV3.
[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 REAL array, dimension( LDQ, max(NN) )
          The (right) orthogonal matrix computed by SGGEV3.
[out]QE
          QE is REAL 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 REAL array, dimension (max(NN))
[out]ALPHAI
          ALPHAI is REAL array, dimension (max(NN))
[out]BETA
          BETA is REAL array, dimension (max(NN))
 \verbatim
          The generalized eigenvalues of (A,B) computed by SGGEV3.
          ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
          generalized eigenvalue of A and B.
[out]ALPHR1
          ALPHR1 is REAL array, dimension (max(NN))
[out]ALPHI1
          ALPHI1 is REAL array, dimension (max(NN))
[out]BETA1
          BETA1 is REAL array, dimension (max(NN))

          Like ALPHAR, ALPHAI, BETA, these arrays contain the
          eigenvalues of A and B, but those computed when SGGEV3 only
          computes a partial eigendecomposition, i.e. not the
          eigenvalues and left and right eigenvectors.
[out]WORK
          WORK is REAL 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 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 404 of file sdrgev3.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 REAL THRESH
417* ..
418* .. Array Arguments ..
419 LOGICAL DOTYPE( * )
420 INTEGER ISEED( 4 ), NN( * )
421 REAL A( LDA, * ), ALPHAI( * ), ALPHI1( * ),
422 $ ALPHAR( * ), 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 REAL ZERO, ONE
432 parameter( zero = 0.0e+0, one = 1.0e+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 REAL 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 REAL RMAGN( 0: 3 )
451* ..
452* .. External Functions ..
453 INTEGER ILAENV
454 REAL SLAMCH, SLARND
455 EXTERNAL ilaenv, slamch, slarnd
456* ..
457* .. External Subroutines ..
458 EXTERNAL alasvm, sget52, sggev3, slacpy, slarfg, slaset,
460* ..
461* .. Intrinsic Functions ..
462 INTRINSIC abs, max, min, real, 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, 'SGEQRF', ' ', 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( 'SDRGEV3', -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 = slamch( 'Safe minimum' )
546 ulp = slamch( 'Epsilon' )*slamch( 'Base' )
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 SLATM4 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* IASIGN: 1 if the diagonal elements of A are to be
600* multiplied by a random magnitude 1 number, =2 if
601* randomly chosen diagonal blocks are to be rotated
602* to form 2x2 blocks.
603* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
604* KTRIAN: =0: don't fill in the upper triangle, =1: do.
605* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
606* RMAGN: used to implement KAMAGN and KBMAGN.
607*
608 IF( mtypes.GT.maxtyp )
609 $ GO TO 100
610 ierr = 0
611 IF( kclass( jtype ).LT.3 ) THEN
612*
613* Generate A (w/o rotation)
614*
615 IF( abs( katype( jtype ) ).EQ.3 ) THEN
616 in = 2*( ( n-1 ) / 2 ) + 1
617 IF( in.NE.n )
618 $ CALL slaset( 'Full', n, n, zero, zero, a, lda )
619 ELSE
620 in = n
621 END IF
622 CALL slatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
623 $ kz2( kazero( jtype ) ), iasign( jtype ),
624 $ rmagn( kamagn( jtype ) ), ulp,
625 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
626 $ iseed, a, lda )
627 iadd = kadd( kazero( jtype ) )
628 IF( iadd.GT.0 .AND. iadd.LE.n )
629 $ a( iadd, iadd ) = one
630*
631* Generate B (w/o rotation)
632*
633 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
634 in = 2*( ( n-1 ) / 2 ) + 1
635 IF( in.NE.n )
636 $ CALL slaset( 'Full', n, n, zero, zero, b, lda )
637 ELSE
638 in = n
639 END IF
640 CALL slatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
641 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
642 $ rmagn( kbmagn( jtype ) ), one,
643 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
644 $ iseed, b, lda )
645 iadd = kadd( kbzero( jtype ) )
646 IF( iadd.NE.0 .AND. iadd.LE.n )
647 $ b( iadd, iadd ) = one
648*
649 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
650*
651* Include rotations
652*
653* Generate Q, Z as Householder transformations times
654* a diagonal matrix.
655*
656 DO 40 jc = 1, n - 1
657 DO 30 jr = jc, n
658 q( jr, jc ) = slarnd( 3, iseed )
659 z( jr, jc ) = slarnd( 3, iseed )
660 30 CONTINUE
661 CALL slarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
662 $ work( jc ) )
663 work( 2*n+jc ) = sign( one, q( jc, jc ) )
664 q( jc, jc ) = one
665 CALL slarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
666 $ work( n+jc ) )
667 work( 3*n+jc ) = sign( one, z( jc, jc ) )
668 z( jc, jc ) = one
669 40 CONTINUE
670 q( n, n ) = one
671 work( n ) = zero
672 work( 3*n ) = sign( one, slarnd( 2, iseed ) )
673 z( n, n ) = one
674 work( 2*n ) = zero
675 work( 4*n ) = sign( one, slarnd( 2, iseed ) )
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 )*work( 3*n+jc )*
682 $ a( jr, jc )
683 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
684 $ b( jr, jc )
685 50 CONTINUE
686 60 CONTINUE
687 CALL sorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
688 $ lda, work( 2*n+1 ), ierr )
689 IF( ierr.NE.0 )
690 $ GO TO 90
691 CALL sorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
692 $ a, lda, work( 2*n+1 ), ierr )
693 IF( ierr.NE.0 )
694 $ GO TO 90
695 CALL sorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
696 $ lda, work( 2*n+1 ), ierr )
697 IF( ierr.NE.0 )
698 $ GO TO 90
699 CALL sorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
700 $ b, lda, work( 2*n+1 ), ierr )
701 IF( ierr.NE.0 )
702 $ GO TO 90
703 END IF
704 ELSE
705*
706* Random matrices
707*
708 DO 80 jc = 1, n
709 DO 70 jr = 1, n
710 a( jr, jc ) = rmagn( kamagn( jtype ) )*
711 $ slarnd( 2, iseed )
712 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
713 $ slarnd( 2, iseed )
714 70 CONTINUE
715 80 CONTINUE
716 END IF
717*
718 90 CONTINUE
719*
720 IF( ierr.NE.0 ) THEN
721 WRITE( nounit, fmt = 9999 )'Generator', ierr, n, jtype,
722 $ ioldsd
723 info = abs( ierr )
724 RETURN
725 END IF
726*
727 100 CONTINUE
728*
729 DO 110 i = 1, 7
730 result( i ) = -one
731 110 CONTINUE
732*
733* Call XLAENV to set the parameters used in SLAQZ0
734*
735 CALL xlaenv( 12, 10 )
736 CALL xlaenv( 13, 12 )
737 CALL xlaenv( 14, 13 )
738 CALL xlaenv( 15, 2 )
739 CALL xlaenv( 17, 10 )
740*
741* Call SGGEV3 to compute eigenvalues and eigenvectors.
742*
743 CALL slacpy( ' ', n, n, a, lda, s, lda )
744 CALL slacpy( ' ', n, n, b, lda, t, lda )
745 CALL sggev3( 'V', 'V', n, s, lda, t, lda, alphar, alphai,
746 $ beta, q, ldq, z, ldq, work, lwork, ierr )
747 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
748 result( 1 ) = ulpinv
749 WRITE( nounit, fmt = 9999 )'SGGEV31', ierr, n, jtype,
750 $ ioldsd
751 info = abs( ierr )
752 GO TO 190
753 END IF
754*
755* Do the tests (1) and (2)
756*
757 CALL sget52( .true., n, a, lda, b, lda, q, ldq, alphar,
758 $ alphai, beta, work, result( 1 ) )
759 IF( result( 2 ).GT.thresh ) THEN
760 WRITE( nounit, fmt = 9998 )'Left', 'SGGEV31',
761 $ result( 2 ), n, jtype, ioldsd
762 END IF
763*
764* Do the tests (3) and (4)
765*
766 CALL sget52( .false., n, a, lda, b, lda, z, ldq, alphar,
767 $ alphai, beta, work, result( 3 ) )
768 IF( result( 4 ).GT.thresh ) THEN
769 WRITE( nounit, fmt = 9998 )'Right', 'SGGEV31',
770 $ result( 4 ), n, jtype, ioldsd
771 END IF
772*
773* Do the test (5)
774*
775 CALL slacpy( ' ', n, n, a, lda, s, lda )
776 CALL slacpy( ' ', n, n, b, lda, t, lda )
777 CALL sggev3( 'N', 'N', n, s, lda, t, lda, alphr1, alphi1,
778 $ beta1, q, ldq, z, ldq, work, lwork, ierr )
779 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
780 result( 1 ) = ulpinv
781 WRITE( nounit, fmt = 9999 )'SGGEV32', ierr, n, jtype,
782 $ ioldsd
783 info = abs( ierr )
784 GO TO 190
785 END IF
786*
787 DO 120 j = 1, n
788 IF( alphar( j ).NE.alphr1( j ) .OR.
789 $ beta( j ).NE. beta1( j ) ) THEN
790 result( 5 ) = ulpinv
791 END IF
792 120 CONTINUE
793*
794* Do the test (6): Compute eigenvalues and left eigenvectors,
795* and test them
796*
797 CALL slacpy( ' ', n, n, a, lda, s, lda )
798 CALL slacpy( ' ', n, n, b, lda, t, lda )
799 CALL sggev3( 'V', 'N', n, s, lda, t, lda, alphr1, alphi1,
800 $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
801 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
802 result( 1 ) = ulpinv
803 WRITE( nounit, fmt = 9999 )'SGGEV33', ierr, n, jtype,
804 $ ioldsd
805 info = abs( ierr )
806 GO TO 190
807 END IF
808*
809 DO 130 j = 1, n
810 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
811 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
812 $ result( 6 ) = ulpinv
813 130 CONTINUE
814*
815 DO 150 j = 1, n
816 DO 140 jc = 1, n
817 IF( q( j, jc ).NE.qe( j, jc ) )
818 $ result( 6 ) = ulpinv
819 140 CONTINUE
820 150 CONTINUE
821*
822* DO the test (7): Compute eigenvalues and right eigenvectors,
823* and test them
824*
825 CALL slacpy( ' ', n, n, a, lda, s, lda )
826 CALL slacpy( ' ', n, n, b, lda, t, lda )
827 CALL sggev3( 'N', 'V', n, s, lda, t, lda, alphr1, alphi1,
828 $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
829 IF( ierr.NE.0 .AND. ierr.NE.n+1 ) THEN
830 result( 1 ) = ulpinv
831 WRITE( nounit, fmt = 9999 )'SGGEV34', ierr, n, jtype,
832 $ ioldsd
833 info = abs( ierr )
834 GO TO 190
835 END IF
836*
837 DO 160 j = 1, n
838 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
839 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
840 $ result( 7 ) = ulpinv
841 160 CONTINUE
842*
843 DO 180 j = 1, n
844 DO 170 jc = 1, n
845 IF( z( j, jc ).NE.qe( j, jc ) )
846 $ result( 7 ) = ulpinv
847 170 CONTINUE
848 180 CONTINUE
849*
850* End of Loop -- Check for RESULT(j) > THRESH
851*
852 190 CONTINUE
853*
854 ntestt = ntestt + 7
855*
856* Print out tests which fail.
857*
858 DO 200 jr = 1, 7
859 IF( result( jr ).GE.thresh ) THEN
860*
861* If this is the first test to fail,
862* print a header to the data file.
863*
864 IF( nerrs.EQ.0 ) THEN
865 WRITE( nounit, fmt = 9997 )'SGV'
866*
867* Matrix types
868*
869 WRITE( nounit, fmt = 9996 )
870 WRITE( nounit, fmt = 9995 )
871 WRITE( nounit, fmt = 9994 )'Orthogonal'
872*
873* Tests performed
874*
875 WRITE( nounit, fmt = 9993 )
876*
877 END IF
878 nerrs = nerrs + 1
879 IF( result( jr ).LT.10000.0 ) THEN
880 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
881 $ result( jr )
882 ELSE
883 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
884 $ result( jr )
885 END IF
886 END IF
887 200 CONTINUE
888*
889 210 CONTINUE
890 220 CONTINUE
891*
892* Summary
893*
894 CALL alasvm( 'SGV', nounit, nerrs, ntestt, 0 )
895*
896 work( 1 ) = maxwrk
897*
898 RETURN
899*
900 9999 FORMAT( ' SDRGEV3: ', a, ' returned INFO=', i6, '.', / 3x, 'N=',
901 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
902*
903 9998 FORMAT( ' SDRGEV3: ', a, ' Eigenvectors from ', a,
904 $ ' incorrectly normalized.', / ' Bits of error=', 0p, g10.3,
905 $ ',', 3x, 'N=', i4, ', JTYPE=', i3, ', ISEED=(',
906 $ 4( i4, ',' ), i5, ')' )
907*
908 9997 FORMAT( / 1x, a3, ' -- Real Generalized eigenvalue problem driver'
909 $ )
910*
911 9996 FORMAT( ' Matrix types (see SDRGEV3 for details): ' )
912*
913 9995 FORMAT( ' Special Matrices:', 23x,
914 $ '(J''=transposed Jordan block)',
915 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
916 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
917 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
918 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
919 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
920 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
921 9994 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
922 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
923 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
924 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
925 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
926 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
927 $ '23=(small,large) 24=(small,small) 25=(large,large)',
928 $ / ' 26=random O(1) matrices.' )
929*
930 9993 FORMAT( / ' Tests performed: ',
931 $ / ' 1 = max | ( b A - a B )''*l | / const.,',
932 $ / ' 2 = | |VR(i)| - 1 | / ulp,',
933 $ / ' 3 = max | ( b A - a B )*r | / const.',
934 $ / ' 4 = | |VL(i)| - 1 | / ulp,',
935 $ / ' 5 = 0 if W same no matter if r or l computed,',
936 $ / ' 6 = 0 if l same no matter if l computed,',
937 $ / ' 7 = 0 if r same no matter if r computed,', / 1x )
938 9992 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
939 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
940 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
941 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, e10.3 )
942*
943* End of SDRGEV3
944*
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 sggev3(jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
SGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (...
Definition sggev3.f:226
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
subroutine sorm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition sorm2r.f:159
subroutine sget52(left, n, a, lda, b, ldb, e, lde, alphar, alphai, beta, work, result)
SGET52
Definition sget52.f:199
real function slarnd(idist, iseed)
SLARND
Definition slarnd.f:73
subroutine slatm4(itype, n, nz1, nz2, isign, amagn, rcond, triang, idist, iseed, a, lda)
SLATM4
Definition slatm4.f:175
Here is the call graph for this function:
Here is the caller graph for this function: