 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ ddrges3()

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

DDRGES3

Purpose:
DDRGES3 checks the nonsymmetric generalized eigenvalue (Schur form)
problem driver DGGES3.

DGGES3 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 DDRGES3 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, DDRGES3 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, DDRGES3 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 DDRGES3 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 DGGES3. 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 DGGES3. [out] Q Q is DOUBLE PRECISION array, dimension (LDQ, max(NN)) The (left) orthogonal matrix computed by DGGES3. [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 DGGES3. [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 DGGES3. ( 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.

Definition at line 399 of file ddrges3.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 ..
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, dgges3, 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 *
489  nmax = 1
490  DO 10 j = 1, nsizes
491  nmax = max( nmax, nn( j ) )
492  IF( nn( j ).LT.0 )
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( 'DDRGES3', -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 )
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 )
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 XLAENV to set the parameters used in DLAQZ0
750 *
751  CALL xlaenv( 12, 10 )
752  CALL xlaenv( 13, 12 )
753  CALL xlaenv( 14, 13 )
754  CALL xlaenv( 15, 2 )
755  CALL xlaenv( 17, 10 )
756 *
757 * Call DGGES3 to compute H, T, Q, Z, alpha, and beta.
758 *
759  CALL dlacpy( 'Full', n, n, a, lda, s, lda )
760  CALL dlacpy( 'Full', n, n, b, lda, t, lda )
761  ntest = 1 + rsub + isort
762  result( 1+rsub+isort ) = ulpinv
763  CALL dgges3( 'V', 'V', sort, dlctes, n, s, lda, t, lda,
764  \$ sdim, alphar, alphai, beta, q, ldq, z, ldq,
765  \$ work, lwork, bwork, iinfo )
766  IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
767  result( 1+rsub+isort ) = ulpinv
768  WRITE( nounit, fmt = 9999 )'DGGES3', iinfo, n, jtype,
769  \$ ioldsd
770  info = abs( iinfo )
771  GO TO 160
772  END IF
773 *
774  ntest = 4 + rsub
775 *
776 * Do tests 1--4 (or tests 7--9 when reordering )
777 *
778  IF( isort.EQ.0 ) THEN
779  CALL dget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
780  \$ work, result( 1 ) )
781  CALL dget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
782  \$ work, result( 2 ) )
783  ELSE
784  CALL dget54( n, a, lda, b, lda, s, lda, t, lda, q,
785  \$ ldq, z, ldq, work, result( 7 ) )
786  END IF
787  CALL dget51( 3, n, a, lda, t, lda, q, ldq, q, ldq, work,
788  \$ result( 3+rsub ) )
789  CALL dget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
790  \$ result( 4+rsub ) )
791 *
792 * Do test 5 and 6 (or Tests 10 and 11 when reordering):
793 * check Schur form of A and compare eigenvalues with
794 * diagonals.
795 *
796  ntest = 6 + rsub
797  temp1 = zero
798 *
799  DO 130 j = 1, n
801  IF( alphai( j ).EQ.zero ) THEN
802  temp2 = ( abs( alphar( j )-s( j, j ) ) /
803  \$ max( safmin, abs( alphar( j ) ), abs( s( j,
804  \$ j ) ) )+abs( beta( j )-t( j, j ) ) /
805  \$ max( safmin, abs( beta( j ) ), abs( t( j,
806  \$ j ) ) ) ) / ulp
807 *
808  IF( j.LT.n ) THEN
809  IF( s( j+1, j ).NE.zero ) THEN
811  result( 5+rsub ) = ulpinv
812  END IF
813  END IF
814  IF( j.GT.1 ) THEN
815  IF( s( j, j-1 ).NE.zero ) THEN
817  result( 5+rsub ) = ulpinv
818  END IF
819  END IF
820 *
821  ELSE
822  IF( alphai( j ).GT.zero ) THEN
823  i1 = j
824  ELSE
825  i1 = j - 1
826  END IF
827  IF( i1.LE.0 .OR. i1.GE.n ) THEN
829  ELSE IF( i1.LT.n-1 ) THEN
830  IF( s( i1+2, i1+1 ).NE.zero ) THEN
832  result( 5+rsub ) = ulpinv
833  END IF
834  ELSE IF( i1.GT.1 ) THEN
835  IF( s( i1, i1-1 ).NE.zero ) THEN
837  result( 5+rsub ) = ulpinv
838  END IF
839  END IF
841  CALL dget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
842  \$ beta( j ), alphar( j ),
843  \$ alphai( j ), temp2, ierr )
844  IF( ierr.GE.3 ) THEN
845  WRITE( nounit, fmt = 9998 )ierr, j, n,
846  \$ jtype, ioldsd
847  info = abs( ierr )
848  END IF
849  ELSE
850  temp2 = ulpinv
851  END IF
852 *
853  END IF
854  temp1 = max( temp1, temp2 )
856  WRITE( nounit, fmt = 9997 )j, n, jtype, ioldsd
857  END IF
858  130 CONTINUE
859  result( 6+rsub ) = temp1
860 *
861  IF( isort.GE.1 ) THEN
862 *
863 * Do test 12
864 *
865  ntest = 12
866  result( 12 ) = zero
867  knteig = 0
868  DO 140 i = 1, n
869  IF( dlctes( alphar( i ), alphai( i ),
870  \$ beta( i ) ) .OR. dlctes( alphar( i ),
871  \$ -alphai( i ), beta( i ) ) ) THEN
872  knteig = knteig + 1
873  END IF
874  IF( i.LT.n ) THEN
875  IF( ( dlctes( alphar( i+1 ), alphai( i+1 ),
876  \$ beta( i+1 ) ) .OR. dlctes( alphar( i+1 ),
877  \$ -alphai( i+1 ), beta( i+1 ) ) ) .AND.
878  \$ ( .NOT.( dlctes( alphar( i ), alphai( i ),
879  \$ beta( i ) ) .OR. dlctes( alphar( i ),
880  \$ -alphai( i ), beta( i ) ) ) ) .AND.
881  \$ iinfo.NE.n+2 ) THEN
882  result( 12 ) = ulpinv
883  END IF
884  END IF
885  140 CONTINUE
886  IF( sdim.NE.knteig ) THEN
887  result( 12 ) = ulpinv
888  END IF
889  END IF
890 *
891  150 CONTINUE
892 *
893 * End of Loop -- Check for RESULT(j) > THRESH
894 *
895  160 CONTINUE
896 *
897  ntestt = ntestt + ntest
898 *
899 * Print out tests which fail.
900 *
901  DO 170 jr = 1, ntest
902  IF( result( jr ).GE.thresh ) THEN
903 *
904 * If this is the first test to fail,
905 * print a header to the data file.
906 *
907  IF( nerrs.EQ.0 ) THEN
908  WRITE( nounit, fmt = 9996 )'DGS'
909 *
910 * Matrix types
911 *
912  WRITE( nounit, fmt = 9995 )
913  WRITE( nounit, fmt = 9994 )
914  WRITE( nounit, fmt = 9993 )'Orthogonal'
915 *
916 * Tests performed
917 *
918  WRITE( nounit, fmt = 9992 )'orthogonal', '''',
919  \$ 'transpose', ( '''', j = 1, 8 )
920 *
921  END IF
922  nerrs = nerrs + 1
923  IF( result( jr ).LT.10000.0d0 ) THEN
924  WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
925  \$ result( jr )
926  ELSE
927  WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
928  \$ result( jr )
929  END IF
930  END IF
931  170 CONTINUE
932 *
933  180 CONTINUE
934  190 CONTINUE
935 *
936 * Summary
937 *
938  CALL alasvm( 'DGS', nounit, nerrs, ntestt, 0 )
939 *
940  work( 1 ) = maxwrk
941 *
942  RETURN
943 *
944  9999 FORMAT( ' DDRGES3: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
945  \$ i6, ', JTYPE=', i6, ', ISEED=(', 4( i4, ',' ), i5, ')' )
946 *
947  9998 FORMAT( ' DDRGES3: DGET53 returned INFO=', i1, ' for eigenvalue ',
948  \$ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(',
949  \$ 4( i4, ',' ), i5, ')' )
950 *
951  9997 FORMAT( ' DDRGES3: S not in Schur form at eigenvalue ', i6, '.',
952  \$ / 9x, 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
953  \$ i5, ')' )
954 *
955  9996 FORMAT( / 1x, a3, ' -- Real Generalized Schur form driver' )
956 *
957  9995 FORMAT( ' Matrix types (see DDRGES3 for details): ' )
958 *
959  9994 FORMAT( ' Special Matrices:', 23x,
960  \$ '(J''=transposed Jordan block)',
961  \$ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
962  \$ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
963  \$ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
964  \$ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
965  \$ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
966  \$ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
967  9993 FORMAT( ' Matrices Rotated by Random ', a, ' Matrices U, V:',
968  \$ / ' 16=Transposed Jordan Blocks 19=geometric ',
969  \$ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
970  \$ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
971  \$ 'alpha, beta=0,1 21=random alpha, beta=0,1',
972  \$ / ' Large & Small Matrices:', / ' 22=(large, small) ',
973  \$ '23=(small,large) 24=(small,small) 25=(large,large)',
974  \$ / ' 26=random O(1) matrices.' )
975 *
976  9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
977  \$ 'Q and Z are ', a, ',', / 19x,
978  \$ 'l and r are the appropriate left and right', / 19x,
979  \$ 'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
980  \$ ' means ', a, '.)', / ' Without ordering: ',
981  \$ / ' 1 = | A - Q S Z', a,
982  \$ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
983  \$ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
984  \$ ' | / ( n ulp ) 4 = | I - ZZ', a,
985  \$ ' | / ( n ulp )', / ' 5 = A is in Schur form S',
986  \$ / ' 6 = difference between (alpha,beta)',
987  \$ ' and diagonals of (S,T)', / ' With ordering: ',
988  \$ / ' 7 = | (A,B) - Q (S,T) Z', a,
989  \$ ' | / ( |(A,B)| n ulp ) ', / ' 8 = | I - QQ', a,
990  \$ ' | / ( n ulp ) 9 = | I - ZZ', a,
991  \$ ' | / ( n ulp )', / ' 10 = A is in Schur form S',
992  \$ / ' 11 = difference between (alpha,beta) and diagonals',
993  \$ ' of (S,T)', / ' 12 = SDIM is the correct number of ',
994  \$ 'selected eigenvalues', / )
995  9991 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
996  \$ 4( i4, ',' ), ' result ', i2, ' is', 0p, f8.2 )
997  9990 FORMAT( ' Matrix order=', i5, ', type=', i2, ', seed=',
998  \$ 4( i4, ',' ), ' result ', i2, ' is', 1p, d10.3 )
999 *
1000 * End of DDRGES3
1001 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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 xlaenv(ISPEC, NVALUE)
XLAENV
Definition: xlaenv.f:81
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 dgges3(JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, BWORK, INFO)
DGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: dgges3.f:282
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: