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

◆ zdrges3()

subroutine zdrges3 ( integer  nsizes,
integer, dimension( * )  nn,
integer  ntypes,
logical, dimension( * )  dotype,
integer, dimension( 4 )  iseed,
double precision  thresh,
integer  nounit,
complex*16, dimension( lda, * )  a,
integer  lda,
complex*16, dimension( lda, * )  b,
complex*16, dimension( lda, * )  s,
complex*16, dimension( lda, * )  t,
complex*16, dimension( ldq, * )  q,
integer  ldq,
complex*16, dimension( ldq, * )  z,
complex*16, dimension( * )  alpha,
complex*16, dimension( * )  beta,
complex*16, dimension( * )  work,
integer  lwork,
double precision, dimension( * )  rwork,
double precision, dimension( 13 )  result,
logical, dimension( * )  bwork,
integer  info 
)

ZDRGES3

Purpose:
 ZDRGES3 checks the nonsymmetric generalized eigenvalue (Schur form)
 problem driver ZGGES3.

 ZGGES3 factors A and B as Q*S*Z'  and Q*T*Z' , where ' means conjugate
 transpose, S and T are  upper triangular (i.e., in generalized Schur
 form), and Q and Z are unitary. 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 ZDRGES3 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. triangular form) (no sorting of
       eigenvalues)

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

                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
           D(j) = ------------------------ + -----------------------
                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)

       (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 elements of the Schur form (S, T),
       i.e. test the maximum over j of D(j)  where:

                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
           D(j) = ------------------------ + -----------------------
                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)

       (with sorting of eigenvalues).

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

 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 COMPLEX*16 array, dimension(LDA, max(NN))
          Used to hold the original A matrix.  Used as input only
          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
          DOTYPE(MAXTYP+1)=.TRUE.
[in]LDA
          LDA is INTEGER
          The leading dimension of A, B, S, and T.
          It must be at least 1 and at least max( NN ).
[in,out]B
          B is COMPLEX*16 array, dimension(LDA, max(NN))
          Used to hold the original B matrix.  Used as input only
          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
          DOTYPE(MAXTYP+1)=.TRUE.
[out]S
          S is COMPLEX*16 array, dimension (LDA, max(NN))
          The Schur form matrix computed from A by ZGGES3.  On exit, S
          contains the Schur form matrix corresponding to the matrix
          in A.
[out]T
          T is COMPLEX*16 array, dimension (LDA, max(NN))
          The upper triangular matrix computed from B by ZGGES3.
[out]Q
          Q is COMPLEX*16 array, dimension (LDQ, max(NN))
          The (left) orthogonal matrix computed by ZGGES3.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of Q and Z. It must
          be at least 1 and at least max( NN ).
[out]Z
          Z is COMPLEX*16 array, dimension( LDQ, max(NN) )
          The (right) orthogonal matrix computed by ZGGES3.
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (max(NN))
[out]BETA
          BETA is COMPLEX*16 array, dimension (max(NN))

          The generalized eigenvalues of (A,B) computed by ZGGES3.
          ALPHA(k) / BETA(k) is the k-th generalized eigenvalue of A
          and B.
[out]WORK
          WORK is COMPLEX*16 array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= 3*N*N.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension ( 8*N )
          Real workspace.
[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 378 of file zdrges3.f.

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