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

◆ ddrges()

subroutine ddrges ( 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( * )  alphar,
double precision, dimension( * )  alphai,
double precision, dimension( * )  beta,
double precision, dimension( * )  work,
integer  lwork,
double precision, dimension( 13 )  result,
logical, dimension( * )  bwork,
integer  info 
)

DDRGES

Purpose:
 DDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
 problem driver DGGES.

 DGGES factors A and B as Q S Z'  and Q T Z' , where ' means
 transpose, T is upper triangular, S is in generalized Schur form
 (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
 the 2x2 blocks corresponding to complex conjugate pairs of
 generalized eigenvalues), and Q and Z are orthogonal. It also
 computes the generalized eigenvalues (alpha(j),beta(j)), j=1,...,n,
 Thus, w(j) = alpha(j)/beta(j) is a root of the characteristic
 equation
                 det( A - w(j) B ) = 0
 Optionally it also reorder the eigenvalues so that a selected
 cluster of eigenvalues appears in the leading diagonal block of the
 Schur forms.

 When DDRGES 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 13 tests
 will be performed and compared with the threshold THRESH except
 the tests (5), (11) and (13).


 (1)   | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues)


 (2)   | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues)


 (3)   | I - QQ' | / ( n ulp ) (no sorting of eigenvalues)


 (4)   | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues)

 (5)   if A is in Schur form (i.e. quasi-triangular form)
       (no sorting of eigenvalues)

 (6)   if eigenvalues = diagonal blocks of the Schur form (S, T),
       i.e., test the maximum over j of D(j)  where:

       if alpha(j) is real:
                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
           D(j) = ------------------------ + -----------------------
                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)

       if alpha(j) is complex:
                                 | det( s S - w T ) |
           D(j) = ---------------------------------------------------
                  ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )

       and S and T are here the 2 x 2 diagonal blocks of S and T
       corresponding to the j-th and j+1-th eigenvalues.
       (no sorting of eigenvalues)

 (7)   | (A,B) - Q (S,T) Z' | / ( | (A,B) | n ulp )
            (with sorting of eigenvalues).

 (8)   | I - QQ' | / ( n ulp ) (with sorting of eigenvalues).

 (9)   | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues).

 (10)  if A is in Schur form (i.e. quasi-triangular form)
       (with sorting of eigenvalues).

 (11)  if eigenvalues = diagonal blocks of the Schur form (S, T),
       i.e. test the maximum over j of D(j)  where:

       if alpha(j) is real:
                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
           D(j) = ------------------------ + -----------------------
                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)

       if alpha(j) is complex:
                                 | det( s S - w T ) |
           D(j) = ---------------------------------------------------
                  ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )

       and S and T are here the 2 x 2 diagonal blocks of S and T
       corresponding to the j-th and j+1-th eigenvalues.
       (with sorting of eigenvalues).

 (12)  if sorting worked and SDIM is the number of eigenvalues
       which were SELECTed.

 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 on input.
          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.  THRESH >= 0.
[in]NOUNIT
          NOUNIT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns IINFO 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) orthogonal matrix computed by DGGES.
[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]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 DGGES.
          ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
          generalized eigenvalue of A and B.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >= MAX( 10*(N+1), 3*N*N ), where N is the largest
          matrix dimension.
[out]RESULT
          RESULT is DOUBLE PRECISION array, dimension (15)
          The values computed by the tests described above.
          The values are currently limited to 1/ulp, to avoid overflow.
[out]BWORK
          BWORK is LOGICAL array, dimension (N)
[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 399 of file ddrges.f.

403*
404* -- LAPACK test routine --
405* -- LAPACK is a software package provided by Univ. of Tennessee, --
406* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
407*
408* .. Scalar Arguments ..
409 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
410 DOUBLE PRECISION THRESH
411* ..
412* .. Array Arguments ..
413 LOGICAL BWORK( * ), DOTYPE( * )
414 INTEGER ISEED( 4 ), NN( * )
415 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
416 $ B( LDA, * ), BETA( * ), Q( LDQ, * ),
417 $ RESULT( 13 ), S( LDA, * ), T( LDA, * ),
418 $ WORK( * ), Z( LDQ, * )
419* ..
420*
421* =====================================================================
422*
423* .. Parameters ..
424 DOUBLE PRECISION ZERO, ONE
425 parameter( zero = 0.0d+0, one = 1.0d+0 )
426 INTEGER MAXTYP
427 parameter( maxtyp = 26 )
428* ..
429* .. Local Scalars ..
430 LOGICAL BADNN, ILABAD
431 CHARACTER SORT
432 INTEGER I, I1, IADD, IERR, IINFO, IN, ISORT, J, JC, JR,
433 $ JSIZE, JTYPE, KNTEIG, MAXWRK, MINWRK, MTYPES,
434 $ N, N1, NB, NERRS, NMATS, NMAX, NTEST, NTESTT,
435 $ RSUB, SDIM
436 DOUBLE PRECISION SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
437* ..
438* .. Local Arrays ..
439 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
440 $ 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 DOUBLE PRECISION RMAGN( 0: 3 )
446* ..
447* .. External Functions ..
448 LOGICAL DLCTES
449 INTEGER ILAENV
450 DOUBLE PRECISION DLAMCH, DLARND
451 EXTERNAL dlctes, ilaenv, dlamch, dlarnd
452* ..
453* .. External Subroutines ..
454 EXTERNAL alasvm, dget51, dget53, dget54, dgges, dlacpy,
456* ..
457* .. Intrinsic Functions ..
458 INTRINSIC abs, dble, max, min, 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 iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
479 $ 5*2, 0 /
480 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
481* ..
482* .. Executable Statements ..
483*
484* Check for errors
485*
486 info = 0
487*
488 badnn = .false.
489 nmax = 1
490 DO 10 j = 1, nsizes
491 nmax = max( nmax, nn( j ) )
492 IF( nn( j ).LT.0 )
493 $ badnn = .true.
494 10 CONTINUE
495*
496 IF( nsizes.LT.0 ) THEN
497 info = -1
498 ELSE IF( badnn ) THEN
499 info = -2
500 ELSE IF( ntypes.LT.0 ) THEN
501 info = -3
502 ELSE IF( thresh.LT.zero ) THEN
503 info = -6
504 ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
505 info = -9
506 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax ) THEN
507 info = -14
508 END IF
509*
510* Compute workspace
511* (Note: Comments in the code beginning "Workspace:" describe the
512* minimal amount of workspace needed at that point in the code,
513* as well as the preferred amount for good performance.
514* NB refers to the optimal block size for the immediately
515* following subroutine, as returned by ILAENV.
516*
517 minwrk = 1
518 IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
519 minwrk = max( 10*( nmax+1 ), 3*nmax*nmax )
520 nb = max( 1, ilaenv( 1, 'DGEQRF', ' ', nmax, nmax, -1, -1 ),
521 $ ilaenv( 1, 'DORMQR', 'LT', nmax, nmax, nmax, -1 ),
522 $ ilaenv( 1, 'DORGQR', ' ', nmax, nmax, nmax, -1 ) )
523 maxwrk = max( 10*( nmax+1 ), 2*nmax+nmax*nb, 3*nmax*nmax )
524 work( 1 ) = maxwrk
525 END IF
526*
527 IF( lwork.LT.minwrk )
528 $ info = -20
529*
530 IF( info.NE.0 ) THEN
531 CALL xerbla( 'DDRGES', -info )
532 RETURN
533 END IF
534*
535* Quick return if possible
536*
537 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
538 $ RETURN
539*
540 safmin = dlamch( 'Safe minimum' )
541 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
542 safmin = safmin / ulp
543 safmax = one / safmin
544 ulpinv = one / ulp
545*
546* The values RMAGN(2:3) depend on N, see below.
547*
548 rmagn( 0 ) = zero
549 rmagn( 1 ) = one
550*
551* Loop over matrix sizes
552*
553 ntestt = 0
554 nerrs = 0
555 nmats = 0
556*
557 DO 190 jsize = 1, nsizes
558 n = nn( jsize )
559 n1 = max( 1, n )
560 rmagn( 2 ) = safmax*ulp / dble( n1 )
561 rmagn( 3 ) = safmin*ulpinv*dble( n1 )
562*
563 IF( nsizes.NE.1 ) THEN
564 mtypes = min( maxtyp, ntypes )
565 ELSE
566 mtypes = min( maxtyp+1, ntypes )
567 END IF
568*
569* Loop over matrix types
570*
571 DO 180 jtype = 1, mtypes
572 IF( .NOT.dotype( jtype ) )
573 $ GO TO 180
574 nmats = nmats + 1
575 ntest = 0
576*
577* Save ISEED in case of an error.
578*
579 DO 20 j = 1, 4
580 ioldsd( j ) = iseed( j )
581 20 CONTINUE
582*
583* Initialize RESULT
584*
585 DO 30 j = 1, 13
586 result( j ) = zero
587 30 CONTINUE
588*
589* Generate test matrices A and B
590*
591* Description of control parameters:
592*
593* KZLASS: =1 means w/o rotation, =2 means w/ rotation,
594* =3 means random.
595* KATYPE: the "type" to be passed to DLATM4 for computing A.
596* KAZERO: the pattern of zeros on the diagonal for A:
597* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
598* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
599* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
600* non-zero entries.)
601* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
602* =2: large, =3: small.
603* IASIGN: 1 if the diagonal elements of A are to be
604* multiplied by a random magnitude 1 number, =2 if
605* randomly chosen diagonal blocks are to be rotated
606* to form 2x2 blocks.
607* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
608* KTRIAN: =0: don't fill in the upper triangle, =1: do.
609* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
610* RMAGN: used to implement KAMAGN and KBMAGN.
611*
612 IF( mtypes.GT.maxtyp )
613 $ GO TO 110
614 iinfo = 0
615 IF( kclass( jtype ).LT.3 ) THEN
616*
617* Generate A (w/o rotation)
618*
619 IF( abs( katype( jtype ) ).EQ.3 ) THEN
620 in = 2*( ( n-1 ) / 2 ) + 1
621 IF( in.NE.n )
622 $ CALL dlaset( 'Full', n, n, zero, zero, a, lda )
623 ELSE
624 in = n
625 END IF
626 CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
627 $ kz2( kazero( jtype ) ), iasign( jtype ),
628 $ rmagn( kamagn( jtype ) ), ulp,
629 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
630 $ iseed, a, lda )
631 iadd = kadd( kazero( jtype ) )
632 IF( iadd.GT.0 .AND. iadd.LE.n )
633 $ a( iadd, iadd ) = one
634*
635* Generate B (w/o rotation)
636*
637 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
638 in = 2*( ( n-1 ) / 2 ) + 1
639 IF( in.NE.n )
640 $ CALL dlaset( 'Full', n, n, zero, zero, b, lda )
641 ELSE
642 in = n
643 END IF
644 CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
645 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
646 $ rmagn( kbmagn( jtype ) ), one,
647 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
648 $ iseed, b, lda )
649 iadd = kadd( kbzero( jtype ) )
650 IF( iadd.NE.0 .AND. iadd.LE.n )
651 $ b( iadd, iadd ) = one
652*
653 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
654*
655* Include rotations
656*
657* Generate Q, Z as Householder transformations times
658* a diagonal matrix.
659*
660 DO 50 jc = 1, n - 1
661 DO 40 jr = jc, n
662 q( jr, jc ) = dlarnd( 3, iseed )
663 z( jr, jc ) = dlarnd( 3, iseed )
664 40 CONTINUE
665 CALL dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
666 $ work( jc ) )
667 work( 2*n+jc ) = sign( one, q( jc, jc ) )
668 q( jc, jc ) = one
669 CALL dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
670 $ work( n+jc ) )
671 work( 3*n+jc ) = sign( one, z( jc, jc ) )
672 z( jc, jc ) = one
673 50 CONTINUE
674 q( n, n ) = one
675 work( n ) = zero
676 work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
677 z( n, n ) = one
678 work( 2*n ) = zero
679 work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
680*
681* Apply the diagonal matrices
682*
683 DO 70 jc = 1, n
684 DO 60 jr = 1, n
685 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
686 $ a( jr, jc )
687 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
688 $ b( jr, jc )
689 60 CONTINUE
690 70 CONTINUE
691 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
692 $ lda, work( 2*n+1 ), iinfo )
693 IF( iinfo.NE.0 )
694 $ GO TO 100
695 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
696 $ a, lda, work( 2*n+1 ), iinfo )
697 IF( iinfo.NE.0 )
698 $ GO TO 100
699 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
700 $ lda, work( 2*n+1 ), iinfo )
701 IF( iinfo.NE.0 )
702 $ GO TO 100
703 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
704 $ b, lda, work( 2*n+1 ), iinfo )
705 IF( iinfo.NE.0 )
706 $ GO TO 100
707 END IF
708 ELSE
709*
710* Random matrices
711*
712 DO 90 jc = 1, n
713 DO 80 jr = 1, n
714 a( jr, jc ) = rmagn( kamagn( jtype ) )*
715 $ dlarnd( 2, iseed )
716 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
717 $ dlarnd( 2, iseed )
718 80 CONTINUE
719 90 CONTINUE
720 END IF
721*
722 100 CONTINUE
723*
724 IF( iinfo.NE.0 ) THEN
725 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
726 $ ioldsd
727 info = abs( iinfo )
728 RETURN
729 END IF
730*
731 110 CONTINUE
732*
733 DO 120 i = 1, 13
734 result( i ) = -one
735 120 CONTINUE
736*
737* Test with and without sorting of eigenvalues
738*
739 DO 150 isort = 0, 1
740 IF( isort.EQ.0 ) THEN
741 sort = 'N'
742 rsub = 0
743 ELSE
744 sort = 'S'
745 rsub = 5
746 END IF
747*
748* Call DGGES to compute H, T, Q, Z, alpha, and beta.
749*
750 CALL dlacpy( 'Full', n, n, a, lda, s, lda )
751 CALL dlacpy( 'Full', n, n, b, lda, t, lda )
752 ntest = 1 + rsub + isort
753 result( 1+rsub+isort ) = ulpinv
754 CALL dgges( 'V', 'V', sort, dlctes, n, s, lda, t, lda,
755 $ sdim, alphar, alphai, beta, q, ldq, z, ldq,
756 $ work, lwork, bwork, iinfo )
757 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
758 result( 1+rsub+isort ) = ulpinv
759 WRITE( nounit, fmt = 9999 )'DGGES', iinfo, n, jtype,
760 $ ioldsd
761 info = abs( iinfo )
762 GO TO 160
763 END IF
764*
765 ntest = 4 + rsub
766*
767* Do tests 1--4 (or tests 7--9 when reordering )
768*
769 IF( isort.EQ.0 ) THEN
770 CALL dget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
771 $ work, result( 1 ) )
772 CALL dget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
773 $ work, result( 2 ) )
774 ELSE
775 CALL dget54( n, a, lda, b, lda, s, lda, t, lda, q,
776 $ ldq, z, ldq, work, result( 7 ) )
777 END IF
778 CALL dget51( 3, n, a, lda, t, lda, q, ldq, q, ldq, work,
779 $ result( 3+rsub ) )
780 CALL dget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
781 $ result( 4+rsub ) )
782*
783* Do test 5 and 6 (or Tests 10 and 11 when reordering):
784* check Schur form of A and compare eigenvalues with
785* diagonals.
786*
787 ntest = 6 + rsub
788 temp1 = zero
789*
790 DO 130 j = 1, n
791 ilabad = .false.
792 IF( alphai( j ).EQ.zero ) THEN
793 temp2 = ( abs( alphar( j )-s( j, j ) ) /
794 $ max( safmin, abs( alphar( j ) ), abs( s( j,
795 $ j ) ) )+abs( beta( j )-t( j, j ) ) /
796 $ max( safmin, abs( beta( j ) ), abs( t( j,
797 $ j ) ) ) ) / ulp
798*
799 IF( j.LT.n ) THEN
800 IF( s( j+1, j ).NE.zero ) THEN
801 ilabad = .true.
802 result( 5+rsub ) = ulpinv
803 END IF
804 END IF
805 IF( j.GT.1 ) THEN
806 IF( s( j, j-1 ).NE.zero ) THEN
807 ilabad = .true.
808 result( 5+rsub ) = ulpinv
809 END IF
810 END IF
811*
812 ELSE
813 IF( alphai( j ).GT.zero ) THEN
814 i1 = j
815 ELSE
816 i1 = j - 1
817 END IF
818 IF( i1.LE.0 .OR. i1.GE.n ) THEN
819 ilabad = .true.
820 ELSE IF( i1.LT.n-1 ) THEN
821 IF( s( i1+2, i1+1 ).NE.zero ) THEN
822 ilabad = .true.
823 result( 5+rsub ) = ulpinv
824 END IF
825 ELSE IF( i1.GT.1 ) THEN
826 IF( s( i1, i1-1 ).NE.zero ) THEN
827 ilabad = .true.
828 result( 5+rsub ) = ulpinv
829 END IF
830 END IF
831 IF( .NOT.ilabad ) THEN
832 CALL dget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
833 $ beta( j ), alphar( j ),
834 $ alphai( j ), temp2, ierr )
835 IF( ierr.GE.3 ) THEN
836 WRITE( nounit, fmt = 9998 )ierr, j, n,
837 $ jtype, ioldsd
838 info = abs( ierr )
839 END IF
840 ELSE
841 temp2 = ulpinv
842 END IF
843*
844 END IF
845 temp1 = max( temp1, temp2 )
846 IF( ilabad ) THEN
847 WRITE( nounit, fmt = 9997 )j, n, jtype, ioldsd
848 END IF
849 130 CONTINUE
850 result( 6+rsub ) = temp1
851*
852 IF( isort.GE.1 ) THEN
853*
854* Do test 12
855*
856 ntest = 12
857 result( 12 ) = zero
858 knteig = 0
859 DO 140 i = 1, n
860 IF( dlctes( alphar( i ), alphai( i ),
861 $ beta( i ) ) .OR. dlctes( alphar( i ),
862 $ -alphai( i ), beta( i ) ) ) THEN
863 knteig = knteig + 1
864 END IF
865 IF( i.LT.n ) THEN
866 IF( ( dlctes( alphar( i+1 ), alphai( i+1 ),
867 $ beta( i+1 ) ) .OR. dlctes( alphar( i+1 ),
868 $ -alphai( i+1 ), beta( i+1 ) ) ) .AND.
869 $ ( .NOT.( dlctes( alphar( i ), alphai( i ),
870 $ beta( i ) ) .OR. dlctes( alphar( i ),
871 $ -alphai( i ), beta( i ) ) ) ) .AND.
872 $ iinfo.NE.n+2 ) THEN
873 result( 12 ) = ulpinv
874 END IF
875 END IF
876 140 CONTINUE
877 IF( sdim.NE.knteig ) THEN
878 result( 12 ) = ulpinv
879 END IF
880 END IF
881*
882 150 CONTINUE
883*
884* End of Loop -- Check for RESULT(j) > THRESH
885*
886 160 CONTINUE
887*
888 ntestt = ntestt + ntest
889*
890* Print out tests which fail.
891*
892 DO 170 jr = 1, ntest
893 IF( result( jr ).GE.thresh ) THEN
894*
895* If this is the first test to fail,
896* print a header to the data file.
897*
898 IF( nerrs.EQ.0 ) THEN
899 WRITE( nounit, fmt = 9996 )'DGS'
900*
901* Matrix types
902*
903 WRITE( nounit, fmt = 9995 )
904 WRITE( nounit, fmt = 9994 )
905 WRITE( nounit, fmt = 9993 )'Orthogonal'
906*
907* Tests performed
908*
909 WRITE( nounit, fmt = 9992 )'orthogonal', '''',
910 $ 'transpose', ( '''', j = 1, 8 )
911*
912 END IF
913 nerrs = nerrs + 1
914 IF( result( jr ).LT.10000.0d0 ) THEN
915 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
916 $ result( jr )
917 ELSE
918 WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
919 $ result( jr )
920 END IF
921 END IF
922 170 CONTINUE
923*
924 180 CONTINUE
925 190 CONTINUE
926*
927* Summary
928*
929 CALL alasvm( 'DGS', nounit, nerrs, ntestt, 0 )
930*
931 work( 1 ) = maxwrk
932*
933 RETURN
934*
935 9999 FORMAT( ' DDRGES: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
936 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
937*
938 9998 FORMAT( ' DDRGES: DGET53 returned INFO=', i1, ' for eigenvalue ',
939 $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(',
940 $ 4( i4, ',' ), i5, ')' )
941*
942 9997 FORMAT( ' DDRGES: S not in Schur form at eigenvalue ', i6, '.',
943 $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
944 $ i5, ')' )
945*
946 9996 FORMAT( / 1x, a3, ' -- Real Generalized Schur form driver' )
947*
948 9995 FORMAT( ' Matrix types (see DDRGES for details): ' )
949*
950 9994 FORMAT( ' Special Matrices:', 23x,
951 $ '(J''=transposed Jordan block)',
952 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
953 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
954 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
955 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
956 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
957 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
958 9993 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
959 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
960 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
961 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
962 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
963 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
964 $ '23=(small,large) 24=(small,small) 25=(large,large)',
965 $ / ' 26=random O(1) matrices.' )
966*
967 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
968 $ 'Q and Z are ', a, ',', / 19x,
969 $ 'l and r are the appropriate left and right', / 19x,
970 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
971 $ ' means ', a, '.)', / ' Without ordering: ',
972 $ / ' 1 = | A - Q S Z', a,
973 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
974 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
975 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
976 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
977 $ / ' 6 = difference between (alpha,beta)',
978 $ ' and diagonals of (S,T)', / ' With ordering: ',
979 $ / ' 7 = | (A,B) - Q (S,T) Z', a,
980 $ ' | / ( |(A,B)| n ulp ) ', / ' 8 = | I - QQ', a,
981 $ ' | / ( n ulp ) 9 = | I - ZZ', a,
982 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
983 $ / ' 11 = difference between (alpha,beta) and diagonals',
984 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
985 $ 'selected eigenvalues', / )
986 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
987 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
988 9990 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
989 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
990*
991* End of DDRGES
992*
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dget51(itype, n, a, lda, b, ldb, u, ldu, v, ldv, work, result)
DGET51
Definition dget51.f:149
subroutine dget53(a, lda, b, ldb, scale, wr, wi, result, info)
DGET53
Definition dget53.f:126
subroutine dget54(n, a, lda, b, ldb, s, lds, t, ldt, u, ldu, v, ldv, work, result)
DGET54
Definition dget54.f:156
double precision function dlarnd(idist, iseed)
DLARND
Definition dlarnd.f:73
subroutine dlatm4(itype, n, nz1, nz2, isign, amagn, rcond, triang, idist, iseed, a, lda)
DLATM4
Definition dlatm4.f:175
logical function dlctes(zr, zi, d)
DLCTES
Definition dlctes.f:68
subroutine dgges(jobvsl, jobvsr, sort, selctg, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info)
DGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE m...
Definition dgges.f:284
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
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
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: