LAPACK 3.11.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, dlabad,
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 CALL dlabad( safmin, safmax )
545 ulpinv = one / ulp
546*
547* The values RMAGN(2:3) depend on N, see below.
548*
549 rmagn( 0 ) = zero
550 rmagn( 1 ) = one
551*
552* Loop over matrix sizes
553*
554 ntestt = 0
555 nerrs = 0
556 nmats = 0
557*
558 DO 190 jsize = 1, nsizes
559 n = nn( jsize )
560 n1 = max( 1, n )
561 rmagn( 2 ) = safmax*ulp / dble( n1 )
562 rmagn( 3 ) = safmin*ulpinv*dble( n1 )
563*
564 IF( nsizes.NE.1 ) THEN
565 mtypes = min( maxtyp, ntypes )
566 ELSE
567 mtypes = min( maxtyp+1, ntypes )
568 END IF
569*
570* Loop over matrix types
571*
572 DO 180 jtype = 1, mtypes
573 IF( .NOT.dotype( jtype ) )
574 $ GO TO 180
575 nmats = nmats + 1
576 ntest = 0
577*
578* Save ISEED in case of an error.
579*
580 DO 20 j = 1, 4
581 ioldsd( j ) = iseed( j )
582 20 CONTINUE
583*
584* Initialize RESULT
585*
586 DO 30 j = 1, 13
587 result( j ) = zero
588 30 CONTINUE
589*
590* Generate test matrices A and B
591*
592* Description of control parameters:
593*
594* KZLASS: =1 means w/o rotation, =2 means w/ rotation,
595* =3 means random.
596* KATYPE: the "type" to be passed to DLATM4 for computing A.
597* KAZERO: the pattern of zeros on the diagonal for A:
598* =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
599* =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
600* =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
601* non-zero entries.)
602* KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
603* =2: large, =3: small.
604* IASIGN: 1 if the diagonal elements of A are to be
605* multiplied by a random magnitude 1 number, =2 if
606* randomly chosen diagonal blocks are to be rotated
607* to form 2x2 blocks.
608* KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
609* KTRIAN: =0: don't fill in the upper triangle, =1: do.
610* KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
611* RMAGN: used to implement KAMAGN and KBMAGN.
612*
613 IF( mtypes.GT.maxtyp )
614 $ GO TO 110
615 iinfo = 0
616 IF( kclass( jtype ).LT.3 ) THEN
617*
618* Generate A (w/o rotation)
619*
620 IF( abs( katype( jtype ) ).EQ.3 ) THEN
621 in = 2*( ( n-1 ) / 2 ) + 1
622 IF( in.NE.n )
623 $ CALL dlaset( 'Full', n, n, zero, zero, a, lda )
624 ELSE
625 in = n
626 END IF
627 CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
628 $ kz2( kazero( jtype ) ), iasign( jtype ),
629 $ rmagn( kamagn( jtype ) ), ulp,
630 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
631 $ iseed, a, lda )
632 iadd = kadd( kazero( jtype ) )
633 IF( iadd.GT.0 .AND. iadd.LE.n )
634 $ a( iadd, iadd ) = one
635*
636* Generate B (w/o rotation)
637*
638 IF( abs( kbtype( jtype ) ).EQ.3 ) THEN
639 in = 2*( ( n-1 ) / 2 ) + 1
640 IF( in.NE.n )
641 $ CALL dlaset( 'Full', n, n, zero, zero, b, lda )
642 ELSE
643 in = n
644 END IF
645 CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
646 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
647 $ rmagn( kbmagn( jtype ) ), one,
648 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
649 $ iseed, b, lda )
650 iadd = kadd( kbzero( jtype ) )
651 IF( iadd.NE.0 .AND. iadd.LE.n )
652 $ b( iadd, iadd ) = one
653*
654 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 ) THEN
655*
656* Include rotations
657*
658* Generate Q, Z as Householder transformations times
659* a diagonal matrix.
660*
661 DO 50 jc = 1, n - 1
662 DO 40 jr = jc, n
663 q( jr, jc ) = dlarnd( 3, iseed )
664 z( jr, jc ) = dlarnd( 3, iseed )
665 40 CONTINUE
666 CALL dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
667 $ work( jc ) )
668 work( 2*n+jc ) = sign( one, q( jc, jc ) )
669 q( jc, jc ) = one
670 CALL dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
671 $ work( n+jc ) )
672 work( 3*n+jc ) = sign( one, z( jc, jc ) )
673 z( jc, jc ) = one
674 50 CONTINUE
675 q( n, n ) = one
676 work( n ) = zero
677 work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
678 z( n, n ) = one
679 work( 2*n ) = zero
680 work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
681*
682* Apply the diagonal matrices
683*
684 DO 70 jc = 1, n
685 DO 60 jr = 1, n
686 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
687 $ a( jr, jc )
688 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
689 $ b( jr, jc )
690 60 CONTINUE
691 70 CONTINUE
692 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, a,
693 $ lda, work( 2*n+1 ), iinfo )
694 IF( iinfo.NE.0 )
695 $ GO TO 100
696 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
697 $ a, lda, work( 2*n+1 ), iinfo )
698 IF( iinfo.NE.0 )
699 $ GO TO 100
700 CALL dorm2r( 'L', 'N', n, n, n-1, q, ldq, work, b,
701 $ lda, work( 2*n+1 ), iinfo )
702 IF( iinfo.NE.0 )
703 $ GO TO 100
704 CALL dorm2r( 'R', 'T', n, n, n-1, z, ldq, work( n+1 ),
705 $ b, lda, work( 2*n+1 ), iinfo )
706 IF( iinfo.NE.0 )
707 $ GO TO 100
708 END IF
709 ELSE
710*
711* Random matrices
712*
713 DO 90 jc = 1, n
714 DO 80 jr = 1, n
715 a( jr, jc ) = rmagn( kamagn( jtype ) )*
716 $ dlarnd( 2, iseed )
717 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
718 $ dlarnd( 2, iseed )
719 80 CONTINUE
720 90 CONTINUE
721 END IF
722*
723 100 CONTINUE
724*
725 IF( iinfo.NE.0 ) THEN
726 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
727 $ ioldsd
728 info = abs( iinfo )
729 RETURN
730 END IF
731*
732 110 CONTINUE
733*
734 DO 120 i = 1, 13
735 result( i ) = -one
736 120 CONTINUE
737*
738* Test with and without sorting of eigenvalues
739*
740 DO 150 isort = 0, 1
741 IF( isort.EQ.0 ) THEN
742 sort = 'N'
743 rsub = 0
744 ELSE
745 sort = 'S'
746 rsub = 5
747 END IF
748*
749* Call DGGES to compute H, T, Q, Z, alpha, and beta.
750*
751 CALL dlacpy( 'Full', n, n, a, lda, s, lda )
752 CALL dlacpy( 'Full', n, n, b, lda, t, lda )
753 ntest = 1 + rsub + isort
754 result( 1+rsub+isort ) = ulpinv
755 CALL dgges( 'V', 'V', sort, dlctes, n, s, lda, t, lda,
756 $ sdim, alphar, alphai, beta, q, ldq, z, ldq,
757 $ work, lwork, bwork, iinfo )
758 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
759 result( 1+rsub+isort ) = ulpinv
760 WRITE( nounit, fmt = 9999 )'DGGES', iinfo, n, jtype,
761 $ ioldsd
762 info = abs( iinfo )
763 GO TO 160
764 END IF
765*
766 ntest = 4 + rsub
767*
768* Do tests 1--4 (or tests 7--9 when reordering )
769*
770 IF( isort.EQ.0 ) THEN
771 CALL dget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
772 $ work, result( 1 ) )
773 CALL dget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
774 $ work, result( 2 ) )
775 ELSE
776 CALL dget54( n, a, lda, b, lda, s, lda, t, lda, q,
777 $ ldq, z, ldq, work, result( 7 ) )
778 END IF
779 CALL dget51( 3, n, a, lda, t, lda, q, ldq, q, ldq, work,
780 $ result( 3+rsub ) )
781 CALL dget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
782 $ result( 4+rsub ) )
783*
784* Do test 5 and 6 (or Tests 10 and 11 when reordering):
785* check Schur form of A and compare eigenvalues with
786* diagonals.
787*
788 ntest = 6 + rsub
789 temp1 = zero
790*
791 DO 130 j = 1, n
792 ilabad = .false.
793 IF( alphai( j ).EQ.zero ) THEN
794 temp2 = ( abs( alphar( j )-s( j, j ) ) /
795 $ max( safmin, abs( alphar( j ) ), abs( s( j,
796 $ j ) ) )+abs( beta( j )-t( j, j ) ) /
797 $ max( safmin, abs( beta( j ) ), abs( t( j,
798 $ j ) ) ) ) / ulp
799*
800 IF( j.LT.n ) THEN
801 IF( s( j+1, j ).NE.zero ) THEN
802 ilabad = .true.
803 result( 5+rsub ) = ulpinv
804 END IF
805 END IF
806 IF( j.GT.1 ) THEN
807 IF( s( j, j-1 ).NE.zero ) THEN
808 ilabad = .true.
809 result( 5+rsub ) = ulpinv
810 END IF
811 END IF
812*
813 ELSE
814 IF( alphai( j ).GT.zero ) THEN
815 i1 = j
816 ELSE
817 i1 = j - 1
818 END IF
819 IF( i1.LE.0 .OR. i1.GE.n ) THEN
820 ilabad = .true.
821 ELSE IF( i1.LT.n-1 ) THEN
822 IF( s( i1+2, i1+1 ).NE.zero ) THEN
823 ilabad = .true.
824 result( 5+rsub ) = ulpinv
825 END IF
826 ELSE IF( i1.GT.1 ) THEN
827 IF( s( i1, i1-1 ).NE.zero ) THEN
828 ilabad = .true.
829 result( 5+rsub ) = ulpinv
830 END IF
831 END IF
832 IF( .NOT.ilabad ) THEN
833 CALL dget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
834 $ beta( j ), alphar( j ),
835 $ alphai( j ), temp2, ierr )
836 IF( ierr.GE.3 ) THEN
837 WRITE( nounit, fmt = 9998 )ierr, j, n,
838 $ jtype, ioldsd
839 info = abs( ierr )
840 END IF
841 ELSE
842 temp2 = ulpinv
843 END IF
844*
845 END IF
846 temp1 = max( temp1, temp2 )
847 IF( ilabad ) THEN
848 WRITE( nounit, fmt = 9997 )j, n, jtype, ioldsd
849 END IF
850 130 CONTINUE
851 result( 6+rsub ) = temp1
852*
853 IF( isort.GE.1 ) THEN
854*
855* Do test 12
856*
857 ntest = 12
858 result( 12 ) = zero
859 knteig = 0
860 DO 140 i = 1, n
861 IF( dlctes( alphar( i ), alphai( i ),
862 $ beta( i ) ) .OR. dlctes( alphar( i ),
863 $ -alphai( i ), beta( i ) ) ) THEN
864 knteig = knteig + 1
865 END IF
866 IF( i.LT.n ) THEN
867 IF( ( dlctes( alphar( i+1 ), alphai( i+1 ),
868 $ beta( i+1 ) ) .OR. dlctes( alphar( i+1 ),
869 $ -alphai( i+1 ), beta( i+1 ) ) ) .AND.
870 $ ( .NOT.( dlctes( alphar( i ), alphai( i ),
871 $ beta( i ) ) .OR. dlctes( alphar( i ),
872 $ -alphai( i ), beta( i ) ) ) ) .AND.
873 $ iinfo.NE.n+2 ) THEN
874 result( 12 ) = ulpinv
875 END IF
876 END IF
877 140 CONTINUE
878 IF( sdim.NE.knteig ) THEN
879 result( 12 ) = ulpinv
880 END IF
881 END IF
882*
883 150 CONTINUE
884*
885* End of Loop -- Check for RESULT(j) > THRESH
886*
887 160 CONTINUE
888*
889 ntestt = ntestt + ntest
890*
891* Print out tests which fail.
892*
893 DO 170 jr = 1, ntest
894 IF( result( jr ).GE.thresh ) THEN
895*
896* If this is the first test to fail,
897* print a header to the data file.
898*
899 IF( nerrs.EQ.0 ) THEN
900 WRITE( nounit, fmt = 9996 )'DGS'
901*
902* Matrix types
903*
904 WRITE( nounit, fmt = 9995 )
905 WRITE( nounit, fmt = 9994 )
906 WRITE( nounit, fmt = 9993 )'Orthogonal'
907*
908* Tests performed
909*
910 WRITE( nounit, fmt = 9992 )'orthogonal', '''',
911 $ 'transpose', ( '''', j = 1, 8 )
912*
913 END IF
914 nerrs = nerrs + 1
915 IF( result( jr ).LT.10000.0d0 ) THEN
916 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
917 $ result( jr )
918 ELSE
919 WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
920 $ result( jr )
921 END IF
922 END IF
923 170 CONTINUE
924*
925 180 CONTINUE
926 190 CONTINUE
927*
928* Summary
929*
930 CALL alasvm( 'DGS', nounit, nerrs, ntestt, 0 )
931*
932 work( 1 ) = maxwrk
933*
934 RETURN
935*
936 9999 FORMAT( ' DDRGES: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
937 $ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
938*
939 9998 FORMAT( ' DDRGES: DGET53 returned INFO=', i1, ' for eigenvalue ',
940 $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(',
941 $ 4( i4, ',' ), i5, ')' )
942*
943 9997 FORMAT( ' DDRGES: S not in Schur form at eigenvalue ', i6, '.',
944 $ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
945 $ i5, ')' )
946*
947 9996 FORMAT( / 1x, a3, ' -- Real Generalized Schur form driver' )
948*
949 9995 FORMAT( ' Matrix types (see DDRGES for details): ' )
950*
951 9994 FORMAT( ' Special Matrices:', 23x,
952 $ '(J''=transposed Jordan block)',
953 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
954 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
955 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
956 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
957 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
958 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
959 9993 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
960 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
961 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
962 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
963 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
964 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
965 $ '23=(small,large) 24=(small,small) 25=(large,large)',
966 $ / ' 26=random O(1) matrices.' )
967*
968 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
969 $ 'Q and Z are ', a, ',', / 19x,
970 $ 'l and r are the appropriate left and right', / 19x,
971 $ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
972 $ ' means ', a, '.)', / ' Without ordering: ',
973 $ / ' 1 = | A - Q S Z', a,
974 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
975 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
976 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
977 $ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
978 $ / ' 6 = difference between (alpha,beta)',
979 $ ' and diagonals of (S,T)', / ' With ordering: ',
980 $ / ' 7 = | (A,B) - Q (S,T) Z', a,
981 $ ' | / ( |(A,B)| n ulp ) ', / ' 8 = | I - QQ', a,
982 $ ' | / ( n ulp ) 9 = | I - ZZ', a,
983 $ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
984 $ / ' 11 = difference between (alpha,beta) and diagonals',
985 $ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
986 $ 'selected eigenvalues', / )
987 9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
988 $ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
989 9990 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
990 $ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
991*
992* End of DDRGES
993*
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 dget54(N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V, LDV, WORK, RESULT)
DGET54
Definition: dget54.f:156
subroutine dget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
DGET53
Definition: dget53.f:126
logical function dlctes(ZR, ZI, D)
DLCTES
Definition: dlctes.f:68
subroutine dget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
DGET51
Definition: dget51.f:149
double precision function dlarnd(IDIST, ISEED)
DLARND
Definition: dlarnd.f:73
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
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: