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

## ◆ dchkhs()

 subroutine dchkhs ( 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, * ) h, double precision, dimension( lda, * ) t1, double precision, dimension( lda, * ) t2, double precision, dimension( ldu, * ) u, integer ldu, double precision, dimension( ldu, * ) z, double precision, dimension( ldu, * ) uz, double precision, dimension( * ) wr1, double precision, dimension( * ) wi1, double precision, dimension( * ) wr2, double precision, dimension( * ) wi2, double precision, dimension( * ) wr3, double precision, dimension( * ) wi3, double precision, dimension( ldu, * ) evectl, double precision, dimension( ldu, * ) evectr, double precision, dimension( ldu, * ) evecty, double precision, dimension( ldu, * ) evectx, double precision, dimension( ldu, * ) uu, double precision, dimension( * ) tau, double precision, dimension( * ) work, integer nwork, integer, dimension( * ) iwork, logical, dimension( * ) select, double precision, dimension( 16 ) result, integer info )

DCHKHS

Purpose:
```    DCHKHS  checks the nonsymmetric eigenvalue problem routines.

DGEHRD factors A as  U H U' , where ' means transpose,
H is hessenberg, and U is an orthogonal matrix.

DORGHR generates the orthogonal matrix U.

DORMHR multiplies a matrix by the orthogonal matrix U.

DHSEQR factors H as  Z T Z' , where Z is orthogonal and
T is "quasi-triangular", and the eigenvalue vector W.

DTREVC computes the left and right eigenvector matrices
L and R for T. L is lower quasi-triangular, and R is
upper quasi-triangular.

DHSEIN computes the left and right eigenvector matrices
Y and X for H, using inverse iteration.

DTREVC3 computes left and right eigenvector matrices
from a Schur matrix T and backtransforms them with Z
to eigenvector matrices L and R for A. L and R are
GE matrices.

When DCHKHS 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, one matrix will be generated and used
to test the nonsymmetric eigenroutines.  For each matrix, 16
tests will be performed:

(1)     | A - U H U**T | / ( |A| n ulp )

(2)     | I - UU**T | / ( n ulp )

(3)     | H - Z T Z**T | / ( |H| n ulp )

(4)     | I - ZZ**T | / ( n ulp )

(5)     | A - UZ H (UZ)**T | / ( |A| n ulp )

(6)     | I - UZ (UZ)**T | / ( n ulp )

(7)     | T(Z computed) - T(Z not computed) | / ( |T| ulp )

(8)     | W(Z computed) - W(Z not computed) | / ( |W| ulp )

(9)     | TR - RW | / ( |T| |R| ulp )

(10)    | L**H T - W**H L | / ( |T| |L| ulp )

(11)    | HX - XW | / ( |H| |X| ulp )

(12)    | Y**H H - W**H Y | / ( |H| |Y| ulp )

(13)    | AX - XW | / ( |A| |X| ulp )

(14)    | Y**H A - W**H Y | / ( |A| |Y| ulp )

(15)    | AR - RW | / ( |A| |R| ulp )

(16)    | LA - WL | / ( |A| |L| ulp )

The "sizes" 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)  The zero matrix.
(2)  The identity matrix.
(3)  A (transposed) Jordan block, with 1's on the diagonal.

(4)  A diagonal matrix with evenly spaced entries
1, ..., ULP  and random signs.
(ULP = (first number larger than 1) - 1 )
(5)  A diagonal matrix with geometrically spaced entries
1, ..., ULP  and random signs.
(6)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
and random signs.

(7)  Same as (4), but multiplied by SQRT( overflow threshold )
(8)  Same as (4), but multiplied by SQRT( underflow threshold )

(9)  A matrix of the form  U' T U, where U is orthogonal and
T has evenly spaced entries 1, ..., ULP with random signs
on the diagonal and random O(1) entries in the upper
triangle.

(10) A matrix of the form  U' T U, where U is orthogonal and
T has geometrically spaced entries 1, ..., ULP with random
signs on the diagonal and random O(1) entries in the upper
triangle.

(11) A matrix of the form  U' T U, where U is orthogonal and
T has "clustered" entries 1, ULP,..., ULP with random
signs on the diagonal and random O(1) entries in the upper
triangle.

(12) A matrix of the form  U' T U, where U is orthogonal and
T has real or complex conjugate paired eigenvalues randomly
chosen from ( ULP, 1 ) and random O(1) entries in the upper
triangle.

(13) A matrix of the form  X' T X, where X has condition
SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
with random signs on the diagonal and random O(1) entries
in the upper triangle.

(14) A matrix of the form  X' T X, where X has condition
SQRT( ULP ) and T has geometrically spaced entries
1, ..., ULP with random signs on the diagonal and random
O(1) entries in the upper triangle.

(15) A matrix of the form  X' T X, where X has condition
SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
with random signs on the diagonal and random O(1) entries
in the upper triangle.

(16) A matrix of the form  X' T X, where X has condition
SQRT( ULP ) and T has real or complex conjugate paired
eigenvalues randomly chosen from ( ULP, 1 ) and random
O(1) entries in the upper triangle.

(17) Same as (16), but multiplied by SQRT( overflow threshold )
(18) Same as (16), but multiplied by SQRT( underflow threshold )

(19) Nonsymmetric matrix with random entries chosen from (-1,1).
(20) Same as (19), but multiplied by SQRT( overflow threshold )
(21) Same as (19), but multiplied by SQRT( underflow threshold )```
```  NSIZES - INTEGER
The number of sizes of matrices to use.  If it is zero,
DCHKHS does nothing.  It must be at least zero.
Not modified.

NN     - INTEGER array, dimension (NSIZES)
An array containing the sizes to be used for the matrices.
Zero values will be skipped.  The values must be at least
zero.
Not modified.

NTYPES - INTEGER
The number of elements in DOTYPE.   If it is zero, DCHKHS
does nothing.  It must be at least zero.  If it is MAXTYP+1
and NSIZES is 1, then an additional type, MAXTYP+1 is
defined, which is to use whatever matrix is in A.  This
is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
DOTYPE(MAXTYP+1) is .TRUE. .
Not modified.

DOTYPE - 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.
Not modified.

ISEED  - 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 DCHKHS to continue the same random number
sequence.
Modified.

THRESH - 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.  It must be at least zero.
Not modified.

NOUNIT - INTEGER
The FORTRAN unit number for printing out error messages
(e.g., if a routine returns IINFO not equal to 0.)
Not modified.

A      - DOUBLE PRECISION array, dimension (LDA,max(NN))
Used to hold the matrix whose eigenvalues are to be
computed.  On exit, A contains the last matrix actually
used.
Modified.

LDA    - INTEGER
The leading dimension of A, H, T1 and T2.  It must be at
least 1 and at least max( NN ).
Not modified.

H      - DOUBLE PRECISION array, dimension (LDA,max(NN))
The upper hessenberg matrix computed by DGEHRD.  On exit,
H contains the Hessenberg form of the matrix in A.
Modified.

T1     - DOUBLE PRECISION array, dimension (LDA,max(NN))
The Schur (="quasi-triangular") matrix computed by DHSEQR
if Z is computed.  On exit, T1 contains the Schur form of
the matrix in A.
Modified.

T2     - DOUBLE PRECISION array, dimension (LDA,max(NN))
The Schur matrix computed by DHSEQR when Z is not computed.
This should be identical to T1.
Modified.

LDU    - INTEGER
The leading dimension of U, Z, UZ and UU.  It must be at
least 1 and at least max( NN ).
Not modified.

U      - DOUBLE PRECISION array, dimension (LDU,max(NN))
The orthogonal matrix computed by DGEHRD.
Modified.

Z      - DOUBLE PRECISION array, dimension (LDU,max(NN))
The orthogonal matrix computed by DHSEQR.
Modified.

UZ     - DOUBLE PRECISION array, dimension (LDU,max(NN))
The product of U times Z.
Modified.

WR1    - DOUBLE PRECISION array, dimension (max(NN))
WI1    - DOUBLE PRECISION array, dimension (max(NN))
The real and imaginary parts of the eigenvalues of A,
as computed when Z is computed.
On exit, WR1 + WI1*i are the eigenvalues of the matrix in A.
Modified.

WR2    - DOUBLE PRECISION array, dimension (max(NN))
WI2    - DOUBLE PRECISION array, dimension (max(NN))
The real and imaginary parts of the eigenvalues of A,
as computed when T is computed but not Z.
On exit, WR2 + WI2*i are the eigenvalues of the matrix in A.
Modified.

WR3    - DOUBLE PRECISION array, dimension (max(NN))
WI3    - DOUBLE PRECISION array, dimension (max(NN))
Like WR1, WI1, these arrays contain the eigenvalues of A,
but those computed when DHSEQR only computes the
eigenvalues, i.e., not the Schur vectors and no more of the
Schur form than is necessary for computing the
eigenvalues.
Modified.

EVECTL - DOUBLE PRECISION array, dimension (LDU,max(NN))
The (upper triangular) left eigenvector matrix for the
matrix in T1.  For complex conjugate pairs, the real part
is stored in one row and the imaginary part in the next.
Modified.

EVEZTR - DOUBLE PRECISION array, dimension (LDU,max(NN))
The (upper triangular) right eigenvector matrix for the
matrix in T1.  For complex conjugate pairs, the real part
is stored in one column and the imaginary part in the next.
Modified.

EVECTY - DOUBLE PRECISION array, dimension (LDU,max(NN))
The left eigenvector matrix for the
matrix in H.  For complex conjugate pairs, the real part
is stored in one row and the imaginary part in the next.
Modified.

EVECTX - DOUBLE PRECISION array, dimension (LDU,max(NN))
The right eigenvector matrix for the
matrix in H.  For complex conjugate pairs, the real part
is stored in one column and the imaginary part in the next.
Modified.

UU     - DOUBLE PRECISION array, dimension (LDU,max(NN))
Details of the orthogonal matrix computed by DGEHRD.
Modified.

TAU    - DOUBLE PRECISION array, dimension(max(NN))
Further details of the orthogonal matrix computed by DGEHRD.
Modified.

WORK   - DOUBLE PRECISION array, dimension (NWORK)
Workspace.
Modified.

NWORK  - INTEGER
The number of entries in WORK.  NWORK >= 4*NN(j)*NN(j) + 2.

IWORK  - INTEGER array, dimension (max(NN))
Workspace.
Modified.

SELECT - LOGICAL array, dimension (max(NN))
Workspace.
Modified.

RESULT - DOUBLE PRECISION array, dimension (16)
The values computed by the fourteen tests described above.
The values are currently limited to 1/ulp, to avoid
overflow.
Modified.

INFO   - INTEGER
If 0, then everything ran OK.
-1: NSIZES < 0
-2: Some NN(j) < 0
-3: NTYPES < 0
-6: THRESH < 0
-9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
-14: LDU < 1 or LDU < NMAX.
-28: NWORK too small.
If  DLATMR, SLATMS, or SLATME returns an error code, the
absolute value of it is returned.
If 1, then DHSEQR could not find all the shifts.
If 2, then the EISPACK code (for small blocks) failed.
If >2, then 30*N iterations were not enough to find an
eigenvalue or to decompose the problem.
Modified.

-----------------------------------------------------------------------

Some Local Variables and Parameters:
---- ----- --------- --- ----------

ZERO, ONE       Real 0 and 1.
MAXTYP          The number of types defined.
MTEST           The number of tests defined: care must be taken
that (1) the size of RESULT, (2) the number of
tests actually performed, and (3) MTEST agree.
NTEST           The number of tests performed on this matrix
so far.  This should be less than MTEST, and
equal to it by the last test.  It will be less
if any of the routines being tested indicates
that it could not compute the matrices that
would be tested.
NMAX            Largest value in NN.
NMATS           The number of matrices generated so far.
NERRS           The number of tests which have exceeded THRESH
so far (computed by DLAFTS).
COND, CONDS,
IMODE           Values to be passed to the matrix generators.
ANORM           Norm of A; passed to matrix generators.

OVFL, UNFL      Overflow and underflow thresholds.
ULP, ULPINV     Finest relative precision and its inverse.
RTOVFL, RTUNFL,
RTULP, RTULPI   Square roots of the previous 4 values.

The following four arrays decode JTYPE:
KTYPE(j)        The general type (1-10) for type "j".
KMODE(j)        The MODE value to be passed to the matrix
generator for type "j".
KMAGN(j)        The order of magnitude ( O(1),
O(overflow^(1/2) ), O(underflow^(1/2) )
KCONDS(j)       Selects whether CONDS is to be 1 or
1/sqrt(ulp).  (0 means irrelevant.)```

Definition at line 417 of file dchkhs.f.

422*
423* -- LAPACK test routine --
424* -- LAPACK is a software package provided by Univ. of Tennessee, --
425* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
426*
427* .. Scalar Arguments ..
428 INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
429 DOUBLE PRECISION THRESH
430* ..
431* .. Array Arguments ..
432 LOGICAL DOTYPE( * ), SELECT( * )
433 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
434 DOUBLE PRECISION A( LDA, * ), EVECTL( LDU, * ),
435 \$ EVECTR( LDU, * ), EVECTX( LDU, * ),
436 \$ EVECTY( LDU, * ), H( LDA, * ), RESULT( 16 ),
437 \$ T1( LDA, * ), T2( LDA, * ), TAU( * ),
438 \$ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
439 \$ WI1( * ), WI2( * ), WI3( * ), WORK( * ),
440 \$ WR1( * ), WR2( * ), WR3( * ), Z( LDU, * )
441* ..
442*
443* =====================================================================
444*
445* .. Parameters ..
446 DOUBLE PRECISION ZERO, ONE
447 parameter( zero = 0.0d0, one = 1.0d0 )
448 INTEGER MAXTYP
449 parameter( maxtyp = 21 )
450* ..
451* .. Local Scalars ..
453 INTEGER I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL,
454 \$ JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS,
455 \$ NMATS, NMAX, NSELC, NSELR, NTEST, NTESTT
456 DOUBLE PRECISION ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP,
457 \$ RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL
458* ..
459* .. Local Arrays ..
461 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
462 \$ KMAGN( MAXTYP ), KMODE( MAXTYP ),
463 \$ KTYPE( MAXTYP )
464 DOUBLE PRECISION DUMMA( 6 )
465* ..
466* .. External Functions ..
467 DOUBLE PRECISION DLAMCH
468 EXTERNAL dlamch
469* ..
470* .. External Subroutines ..
471 EXTERNAL dcopy, dgehrd, dgemm, dget10, dget22, dhsein,
474 \$ dtrevc3, xerbla
475* ..
476* .. Intrinsic Functions ..
477 INTRINSIC abs, dble, max, min, sqrt
478* ..
479* .. Data statements ..
480 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
481 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
482 \$ 3, 1, 2, 3 /
483 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
484 \$ 1, 5, 5, 5, 4, 3, 1 /
485 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
486* ..
487* .. Executable Statements ..
488*
489* Check for errors
490*
491 ntestt = 0
492 info = 0
493*
495 nmax = 0
496 DO 10 j = 1, nsizes
497 nmax = max( nmax, nn( j ) )
498 IF( nn( j ).LT.0 )
500 10 CONTINUE
501*
502* Check for errors
503*
504 IF( nsizes.LT.0 ) THEN
505 info = -1
506 ELSE IF( badnn ) THEN
507 info = -2
508 ELSE IF( ntypes.LT.0 ) THEN
509 info = -3
510 ELSE IF( thresh.LT.zero ) THEN
511 info = -6
512 ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
513 info = -9
514 ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax ) THEN
515 info = -14
516 ELSE IF( 4*nmax*nmax+2.GT.nwork ) THEN
517 info = -28
518 END IF
519*
520 IF( info.NE.0 ) THEN
521 CALL xerbla( 'DCHKHS', -info )
522 RETURN
523 END IF
524*
525* Quick return if possible
526*
527 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
528 \$ RETURN
529*
530* More important constants
531*
532 unfl = dlamch( 'Safe minimum' )
533 ovfl = dlamch( 'Overflow' )
534 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
535 ulpinv = one / ulp
536 rtunfl = sqrt( unfl )
537 rtovfl = sqrt( ovfl )
538 rtulp = sqrt( ulp )
539 rtulpi = one / rtulp
540*
541* Loop over sizes, types
542*
543 nerrs = 0
544 nmats = 0
545*
546 DO 270 jsize = 1, nsizes
547 n = nn( jsize )
548 IF( n.EQ.0 )
549 \$ GO TO 270
550 n1 = max( 1, n )
551 aninv = one / dble( n1 )
552*
553 IF( nsizes.NE.1 ) THEN
554 mtypes = min( maxtyp, ntypes )
555 ELSE
556 mtypes = min( maxtyp+1, ntypes )
557 END IF
558*
559 DO 260 jtype = 1, mtypes
560 IF( .NOT.dotype( jtype ) )
561 \$ GO TO 260
562 nmats = nmats + 1
563 ntest = 0
564*
565* Save ISEED in case of an error.
566*
567 DO 20 j = 1, 4
568 ioldsd( j ) = iseed( j )
569 20 CONTINUE
570*
571* Initialize RESULT
572*
573 DO 30 j = 1, 16
574 result( j ) = zero
575 30 CONTINUE
576*
577* Compute "A"
578*
579* Control parameters:
580*
581* KMAGN KCONDS KMODE KTYPE
582* =1 O(1) 1 clustered 1 zero
583* =2 large large clustered 2 identity
584* =3 small exponential Jordan
585* =4 arithmetic diagonal, (w/ eigenvalues)
586* =5 random log symmetric, w/ eigenvalues
587* =6 random general, w/ eigenvalues
588* =7 random diagonal
589* =8 random symmetric
590* =9 random general
591* =10 random triangular
592*
593 IF( mtypes.GT.maxtyp )
594 \$ GO TO 100
595*
596 itype = ktype( jtype )
597 imode = kmode( jtype )
598*
599* Compute norm
600*
601 GO TO ( 40, 50, 60 )kmagn( jtype )
602*
603 40 CONTINUE
604 anorm = one
605 GO TO 70
606*
607 50 CONTINUE
608 anorm = ( rtovfl*ulp )*aninv
609 GO TO 70
610*
611 60 CONTINUE
612 anorm = rtunfl*n*ulpinv
613 GO TO 70
614*
615 70 CONTINUE
616*
617 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
618 iinfo = 0
619 cond = ulpinv
620*
621* Special Matrices
622*
623 IF( itype.EQ.1 ) THEN
624*
625* Zero
626*
627 iinfo = 0
628*
629 ELSE IF( itype.EQ.2 ) THEN
630*
631* Identity
632*
633 DO 80 jcol = 1, n
634 a( jcol, jcol ) = anorm
635 80 CONTINUE
636*
637 ELSE IF( itype.EQ.3 ) THEN
638*
639* Jordan Block
640*
641 DO 90 jcol = 1, n
642 a( jcol, jcol ) = anorm
643 IF( jcol.GT.1 )
644 \$ a( jcol, jcol-1 ) = one
645 90 CONTINUE
646*
647 ELSE IF( itype.EQ.4 ) THEN
648*
649* Diagonal Matrix, [Eigen]values Specified
650*
651 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
652 \$ anorm, 0, 0, 'N', a, lda, work( n+1 ),
653 \$ iinfo )
654*
655 ELSE IF( itype.EQ.5 ) THEN
656*
657* Symmetric, eigenvalues specified
658*
659 CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
660 \$ anorm, n, n, 'N', a, lda, work( n+1 ),
661 \$ iinfo )
662*
663 ELSE IF( itype.EQ.6 ) THEN
664*
665* General, eigenvalues specified
666*
667 IF( kconds( jtype ).EQ.1 ) THEN
668 conds = one
669 ELSE IF( kconds( jtype ).EQ.2 ) THEN
670 conds = rtulpi
671 ELSE
672 conds = zero
673 END IF
674*
675 adumma( 1 ) = ' '
676 CALL dlatme( n, 'S', iseed, work, imode, cond, one,
677 \$ adumma, 'T', 'T', 'T', work( n+1 ), 4,
678 \$ conds, n, n, anorm, a, lda, work( 2*n+1 ),
679 \$ iinfo )
680*
681 ELSE IF( itype.EQ.7 ) THEN
682*
683* Diagonal, random eigenvalues
684*
685 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
686 \$ 'T', 'N', work( n+1 ), 1, one,
687 \$ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
688 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
689*
690 ELSE IF( itype.EQ.8 ) THEN
691*
692* Symmetric, random eigenvalues
693*
694 CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
695 \$ 'T', 'N', work( n+1 ), 1, one,
696 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
697 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
698*
699 ELSE IF( itype.EQ.9 ) THEN
700*
701* General, random eigenvalues
702*
703 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
704 \$ 'T', 'N', work( n+1 ), 1, one,
705 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
706 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
707*
708 ELSE IF( itype.EQ.10 ) THEN
709*
710* Triangular, random eigenvalues
711*
712 CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
713 \$ 'T', 'N', work( n+1 ), 1, one,
714 \$ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
715 \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
716*
717 ELSE
718*
719 iinfo = 1
720 END IF
721*
722 IF( iinfo.NE.0 ) THEN
723 WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
724 \$ ioldsd
725 info = abs( iinfo )
726 RETURN
727 END IF
728*
729 100 CONTINUE
730*
731* Call DGEHRD to compute H and U, do tests.
732*
733 CALL dlacpy( ' ', n, n, a, lda, h, lda )
734*
735 ntest = 1
736*
737 ilo = 1
738 ihi = n
739*
740 CALL dgehrd( n, ilo, ihi, h, lda, work, work( n+1 ),
741 \$ nwork-n, iinfo )
742*
743 IF( iinfo.NE.0 ) THEN
744 result( 1 ) = ulpinv
745 WRITE( nounit, fmt = 9999 )'DGEHRD', iinfo, n, jtype,
746 \$ ioldsd
747 info = abs( iinfo )
748 GO TO 250
749 END IF
750*
751 DO 120 j = 1, n - 1
752 uu( j+1, j ) = zero
753 DO 110 i = j + 2, n
754 u( i, j ) = h( i, j )
755 uu( i, j ) = h( i, j )
756 h( i, j ) = zero
757 110 CONTINUE
758 120 CONTINUE
759 CALL dcopy( n-1, work, 1, tau, 1 )
760 CALL dorghr( n, ilo, ihi, u, ldu, work, work( n+1 ),
761 \$ nwork-n, iinfo )
762 ntest = 2
763*
764 CALL dhst01( n, ilo, ihi, a, lda, h, lda, u, ldu, work,
765 \$ nwork, result( 1 ) )
766*
767* Call DHSEQR to compute T1, T2 and Z, do tests.
768*
769* Eigenvalues only (WR3,WI3)
770*
771 CALL dlacpy( ' ', n, n, h, lda, t2, lda )
772 ntest = 3
773 result( 3 ) = ulpinv
774*
775 CALL dhseqr( 'E', 'N', n, ilo, ihi, t2, lda, wr3, wi3, uz,
776 \$ ldu, work, nwork, iinfo )
777 IF( iinfo.NE.0 ) THEN
778 WRITE( nounit, fmt = 9999 )'DHSEQR(E)', iinfo, n, jtype,
779 \$ ioldsd
780 IF( iinfo.LE.n+2 ) THEN
781 info = abs( iinfo )
782 GO TO 250
783 END IF
784 END IF
785*
786* Eigenvalues (WR2,WI2) and Full Schur Form (T2)
787*
788 CALL dlacpy( ' ', n, n, h, lda, t2, lda )
789*
790 CALL dhseqr( 'S', 'N', n, ilo, ihi, t2, lda, wr2, wi2, uz,
791 \$ ldu, work, nwork, iinfo )
792 IF( iinfo.NE.0 .AND. iinfo.LE.n+2 ) THEN
793 WRITE( nounit, fmt = 9999 )'DHSEQR(S)', iinfo, n, jtype,
794 \$ ioldsd
795 info = abs( iinfo )
796 GO TO 250
797 END IF
798*
799* Eigenvalues (WR1,WI1), Schur Form (T1), and Schur vectors
800* (UZ)
801*
802 CALL dlacpy( ' ', n, n, h, lda, t1, lda )
803 CALL dlacpy( ' ', n, n, u, ldu, uz, ldu )
804*
805 CALL dhseqr( 'S', 'V', n, ilo, ihi, t1, lda, wr1, wi1, uz,
806 \$ ldu, work, nwork, iinfo )
807 IF( iinfo.NE.0 .AND. iinfo.LE.n+2 ) THEN
808 WRITE( nounit, fmt = 9999 )'DHSEQR(V)', iinfo, n, jtype,
809 \$ ioldsd
810 info = abs( iinfo )
811 GO TO 250
812 END IF
813*
814* Compute Z = U' UZ
815*
816 CALL dgemm( 'T', 'N', n, n, n, one, u, ldu, uz, ldu, zero,
817 \$ z, ldu )
818 ntest = 8
819*
820* Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
821* and 4: | I - Z Z' | / ( n ulp )
822*
823 CALL dhst01( n, ilo, ihi, h, lda, t1, lda, z, ldu, work,
824 \$ nwork, result( 3 ) )
825*
826* Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
827* and 6: | I - UZ (UZ)' | / ( n ulp )
828*
829 CALL dhst01( n, ilo, ihi, a, lda, t1, lda, uz, ldu, work,
830 \$ nwork, result( 5 ) )
831*
832* Do Test 7: | T2 - T1 | / ( |T| n ulp )
833*
834 CALL dget10( n, n, t2, lda, t1, lda, work, result( 7 ) )
835*
836* Do Test 8: | W2 - W1 | / ( max(|W1|,|W2|) ulp )
837*
838 temp1 = zero
839 temp2 = zero
840 DO 130 j = 1, n
841 temp1 = max( temp1, abs( wr1( j ) )+abs( wi1( j ) ),
842 \$ abs( wr2( j ) )+abs( wi2( j ) ) )
843 temp2 = max( temp2, abs( wr1( j )-wr2( j ) )+
844 & abs( wi1( j )-wi2( j ) ) )
845 130 CONTINUE
846*
847 result( 8 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
848*
849* Compute the Left and Right Eigenvectors of T
850*
851* Compute the Right eigenvector Matrix:
852*
853 ntest = 9
854 result( 9 ) = ulpinv
855*
856* Select last max(N/4,1) real, max(N/4,1) complex eigenvectors
857*
858 nselc = 0
859 nselr = 0
860 j = n
861 140 CONTINUE
862 IF( wi1( j ).EQ.zero ) THEN
863 IF( nselr.LT.max( n / 4, 1 ) ) THEN
864 nselr = nselr + 1
865 SELECT( j ) = .true.
866 ELSE
867 SELECT( j ) = .false.
868 END IF
869 j = j - 1
870 ELSE
871 IF( nselc.LT.max( n / 4, 1 ) ) THEN
872 nselc = nselc + 1
873 SELECT( j ) = .true.
874 SELECT( j-1 ) = .false.
875 ELSE
876 SELECT( j ) = .false.
877 SELECT( j-1 ) = .false.
878 END IF
879 j = j - 2
880 END IF
881 IF( j.GT.0 )
882 \$ GO TO 140
883*
884 CALL dtrevc( 'Right', 'All', SELECT, n, t1, lda, dumma, ldu,
885 \$ evectr, ldu, n, in, work, iinfo )
886 IF( iinfo.NE.0 ) THEN
887 WRITE( nounit, fmt = 9999 )'DTREVC(R,A)', iinfo, n,
888 \$ jtype, ioldsd
889 info = abs( iinfo )
890 GO TO 250
891 END IF
892*
893* Test 9: | TR - RW | / ( |T| |R| ulp )
894*
895 CALL dget22( 'N', 'N', 'N', n, t1, lda, evectr, ldu, wr1,
896 \$ wi1, work, dumma( 1 ) )
897 result( 9 ) = dumma( 1 )
898 IF( dumma( 2 ).GT.thresh ) THEN
899 WRITE( nounit, fmt = 9998 )'Right', 'DTREVC',
900 \$ dumma( 2 ), n, jtype, ioldsd
901 END IF
902*
903* Compute selected right eigenvectors and confirm that
904* they agree with previous right eigenvectors
905*
906 CALL dtrevc( 'Right', 'Some', SELECT, n, t1, lda, dumma,
907 \$ ldu, evectl, ldu, n, in, work, iinfo )
908 IF( iinfo.NE.0 ) THEN
909 WRITE( nounit, fmt = 9999 )'DTREVC(R,S)', iinfo, n,
910 \$ jtype, ioldsd
911 info = abs( iinfo )
912 GO TO 250
913 END IF
914*
915 k = 1
916 match = .true.
917 DO 170 j = 1, n
918 IF( SELECT( j ) .AND. wi1( j ).EQ.zero ) THEN
919 DO 150 jj = 1, n
920 IF( evectr( jj, j ).NE.evectl( jj, k ) ) THEN
921 match = .false.
922 GO TO 180
923 END IF
924 150 CONTINUE
925 k = k + 1
926 ELSE IF( SELECT( j ) .AND. wi1( j ).NE.zero ) THEN
927 DO 160 jj = 1, n
928 IF( evectr( jj, j ).NE.evectl( jj, k ) .OR.
929 \$ evectr( jj, j+1 ).NE.evectl( jj, k+1 ) ) THEN
930 match = .false.
931 GO TO 180
932 END IF
933 160 CONTINUE
934 k = k + 2
935 END IF
936 170 CONTINUE
937 180 CONTINUE
938 IF( .NOT.match )
939 \$ WRITE( nounit, fmt = 9997 )'Right', 'DTREVC', n, jtype,
940 \$ ioldsd
941*
942* Compute the Left eigenvector Matrix:
943*
944 ntest = 10
945 result( 10 ) = ulpinv
946 CALL dtrevc( 'Left', 'All', SELECT, n, t1, lda, evectl, ldu,
947 \$ dumma, ldu, n, in, work, iinfo )
948 IF( iinfo.NE.0 ) THEN
949 WRITE( nounit, fmt = 9999 )'DTREVC(L,A)', iinfo, n,
950 \$ jtype, ioldsd
951 info = abs( iinfo )
952 GO TO 250
953 END IF
954*
955* Test 10: | LT - WL | / ( |T| |L| ulp )
956*
957 CALL dget22( 'Trans', 'N', 'Conj', n, t1, lda, evectl, ldu,
958 \$ wr1, wi1, work, dumma( 3 ) )
959 result( 10 ) = dumma( 3 )
960 IF( dumma( 4 ).GT.thresh ) THEN
961 WRITE( nounit, fmt = 9998 )'Left', 'DTREVC', dumma( 4 ),
962 \$ n, jtype, ioldsd
963 END IF
964*
965* Compute selected left eigenvectors and confirm that
966* they agree with previous left eigenvectors
967*
968 CALL dtrevc( 'Left', 'Some', SELECT, n, t1, lda, evectr,
969 \$ ldu, dumma, ldu, n, in, work, iinfo )
970 IF( iinfo.NE.0 ) THEN
971 WRITE( nounit, fmt = 9999 )'DTREVC(L,S)', iinfo, n,
972 \$ jtype, ioldsd
973 info = abs( iinfo )
974 GO TO 250
975 END IF
976*
977 k = 1
978 match = .true.
979 DO 210 j = 1, n
980 IF( SELECT( j ) .AND. wi1( j ).EQ.zero ) THEN
981 DO 190 jj = 1, n
982 IF( evectl( jj, j ).NE.evectr( jj, k ) ) THEN
983 match = .false.
984 GO TO 220
985 END IF
986 190 CONTINUE
987 k = k + 1
988 ELSE IF( SELECT( j ) .AND. wi1( j ).NE.zero ) THEN
989 DO 200 jj = 1, n
990 IF( evectl( jj, j ).NE.evectr( jj, k ) .OR.
991 \$ evectl( jj, j+1 ).NE.evectr( jj, k+1 ) ) THEN
992 match = .false.
993 GO TO 220
994 END IF
995 200 CONTINUE
996 k = k + 2
997 END IF
998 210 CONTINUE
999 220 CONTINUE
1000 IF( .NOT.match )
1001 \$ WRITE( nounit, fmt = 9997 )'Left', 'DTREVC', n, jtype,
1002 \$ ioldsd
1003*
1004* Call DHSEIN for Right eigenvectors of H, do test 11
1005*
1006 ntest = 11
1007 result( 11 ) = ulpinv
1008 DO 230 j = 1, n
1009 SELECT( j ) = .true.
1010 230 CONTINUE
1011*
1012 CALL dhsein( 'Right', 'Qr', 'Ninitv', SELECT, n, h, lda,
1013 \$ wr3, wi3, dumma, ldu, evectx, ldu, n1, in,
1014 \$ work, iwork, iwork, iinfo )
1015 IF( iinfo.NE.0 ) THEN
1016 WRITE( nounit, fmt = 9999 )'DHSEIN(R)', iinfo, n, jtype,
1017 \$ ioldsd
1018 info = abs( iinfo )
1019 IF( iinfo.LT.0 )
1020 \$ GO TO 250
1021 ELSE
1022*
1023* Test 11: | HX - XW | / ( |H| |X| ulp )
1024*
1025* (from inverse iteration)
1026*
1027 CALL dget22( 'N', 'N', 'N', n, h, lda, evectx, ldu, wr3,
1028 \$ wi3, work, dumma( 1 ) )
1029 IF( dumma( 1 ).LT.ulpinv )
1030 \$ result( 11 ) = dumma( 1 )*aninv
1031 IF( dumma( 2 ).GT.thresh ) THEN
1032 WRITE( nounit, fmt = 9998 )'Right', 'DHSEIN',
1033 \$ dumma( 2 ), n, jtype, ioldsd
1034 END IF
1035 END IF
1036*
1037* Call DHSEIN for Left eigenvectors of H, do test 12
1038*
1039 ntest = 12
1040 result( 12 ) = ulpinv
1041 DO 240 j = 1, n
1042 SELECT( j ) = .true.
1043 240 CONTINUE
1044*
1045 CALL dhsein( 'Left', 'Qr', 'Ninitv', SELECT, n, h, lda, wr3,
1046 \$ wi3, evecty, ldu, dumma, ldu, n1, in, work,
1047 \$ iwork, iwork, iinfo )
1048 IF( iinfo.NE.0 ) THEN
1049 WRITE( nounit, fmt = 9999 )'DHSEIN(L)', iinfo, n, jtype,
1050 \$ ioldsd
1051 info = abs( iinfo )
1052 IF( iinfo.LT.0 )
1053 \$ GO TO 250
1054 ELSE
1055*
1056* Test 12: | YH - WY | / ( |H| |Y| ulp )
1057*
1058* (from inverse iteration)
1059*
1060 CALL dget22( 'C', 'N', 'C', n, h, lda, evecty, ldu, wr3,
1061 \$ wi3, work, dumma( 3 ) )
1062 IF( dumma( 3 ).LT.ulpinv )
1063 \$ result( 12 ) = dumma( 3 )*aninv
1064 IF( dumma( 4 ).GT.thresh ) THEN
1065 WRITE( nounit, fmt = 9998 )'Left', 'DHSEIN',
1066 \$ dumma( 4 ), n, jtype, ioldsd
1067 END IF
1068 END IF
1069*
1070* Call DORMHR for Right eigenvectors of A, do test 13
1071*
1072 ntest = 13
1073 result( 13 ) = ulpinv
1074*
1075 CALL dormhr( 'Left', 'No transpose', n, n, ilo, ihi, uu,
1076 \$ ldu, tau, evectx, ldu, work, nwork, iinfo )
1077 IF( iinfo.NE.0 ) THEN
1078 WRITE( nounit, fmt = 9999 )'DORMHR(R)', iinfo, n, jtype,
1079 \$ ioldsd
1080 info = abs( iinfo )
1081 IF( iinfo.LT.0 )
1082 \$ GO TO 250
1083 ELSE
1084*
1085* Test 13: | AX - XW | / ( |A| |X| ulp )
1086*
1087* (from inverse iteration)
1088*
1089 CALL dget22( 'N', 'N', 'N', n, a, lda, evectx, ldu, wr3,
1090 \$ wi3, work, dumma( 1 ) )
1091 IF( dumma( 1 ).LT.ulpinv )
1092 \$ result( 13 ) = dumma( 1 )*aninv
1093 END IF
1094*
1095* Call DORMHR for Left eigenvectors of A, do test 14
1096*
1097 ntest = 14
1098 result( 14 ) = ulpinv
1099*
1100 CALL dormhr( 'Left', 'No transpose', n, n, ilo, ihi, uu,
1101 \$ ldu, tau, evecty, ldu, work, nwork, iinfo )
1102 IF( iinfo.NE.0 ) THEN
1103 WRITE( nounit, fmt = 9999 )'DORMHR(L)', iinfo, n, jtype,
1104 \$ ioldsd
1105 info = abs( iinfo )
1106 IF( iinfo.LT.0 )
1107 \$ GO TO 250
1108 ELSE
1109*
1110* Test 14: | YA - WY | / ( |A| |Y| ulp )
1111*
1112* (from inverse iteration)
1113*
1114 CALL dget22( 'C', 'N', 'C', n, a, lda, evecty, ldu, wr3,
1115 \$ wi3, work, dumma( 3 ) )
1116 IF( dumma( 3 ).LT.ulpinv )
1117 \$ result( 14 ) = dumma( 3 )*aninv
1118 END IF
1119*
1120* Compute Left and Right Eigenvectors of A
1121*
1122* Compute a Right eigenvector matrix:
1123*
1124 ntest = 15
1125 result( 15 ) = ulpinv
1126*
1127 CALL dlacpy( ' ', n, n, uz, ldu, evectr, ldu )
1128*
1129 CALL dtrevc3( 'Right', 'Back', SELECT, n, t1, lda, dumma,
1130 \$ ldu, evectr, ldu, n, in, work, nwork, iinfo )
1131 IF( iinfo.NE.0 ) THEN
1132 WRITE( nounit, fmt = 9999 )'DTREVC3(R,B)', iinfo, n,
1133 \$ jtype, ioldsd
1134 info = abs( iinfo )
1135 GO TO 250
1136 END IF
1137*
1138* Test 15: | AR - RW | / ( |A| |R| ulp )
1139*
1140* (from Schur decomposition)
1141*
1142 CALL dget22( 'N', 'N', 'N', n, a, lda, evectr, ldu, wr1,
1143 \$ wi1, work, dumma( 1 ) )
1144 result( 15 ) = dumma( 1 )
1145 IF( dumma( 2 ).GT.thresh ) THEN
1146 WRITE( nounit, fmt = 9998 )'Right', 'DTREVC3',
1147 \$ dumma( 2 ), n, jtype, ioldsd
1148 END IF
1149*
1150* Compute a Left eigenvector matrix:
1151*
1152 ntest = 16
1153 result( 16 ) = ulpinv
1154*
1155 CALL dlacpy( ' ', n, n, uz, ldu, evectl, ldu )
1156*
1157 CALL dtrevc3( 'Left', 'Back', SELECT, n, t1, lda, evectl,
1158 \$ ldu, dumma, ldu, n, in, work, nwork, iinfo )
1159 IF( iinfo.NE.0 ) THEN
1160 WRITE( nounit, fmt = 9999 )'DTREVC3(L,B)', iinfo, n,
1161 \$ jtype, ioldsd
1162 info = abs( iinfo )
1163 GO TO 250
1164 END IF
1165*
1166* Test 16: | LA - WL | / ( |A| |L| ulp )
1167*
1168* (from Schur decomposition)
1169*
1170 CALL dget22( 'Trans', 'N', 'Conj', n, a, lda, evectl, ldu,
1171 \$ wr1, wi1, work, dumma( 3 ) )
1172 result( 16 ) = dumma( 3 )
1173 IF( dumma( 4 ).GT.thresh ) THEN
1174 WRITE( nounit, fmt = 9998 )'Left', 'DTREVC3', dumma( 4 ),
1175 \$ n, jtype, ioldsd
1176 END IF
1177*
1178* End of Loop -- Check for RESULT(j) > THRESH
1179*
1180 250 CONTINUE
1181*
1182 ntestt = ntestt + ntest
1183 CALL dlafts( 'DHS', n, n, jtype, ntest, result, ioldsd,
1184 \$ thresh, nounit, nerrs )
1185*
1186 260 CONTINUE
1187 270 CONTINUE
1188*
1189* Summary
1190*
1191 CALL dlasum( 'DHS', nounit, nerrs, ntestt )
1192*
1193 RETURN
1194*
1195 9999 FORMAT( ' DCHKHS: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
1196 \$ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1197 9998 FORMAT( ' DCHKHS: ', a, ' Eigenvectors from ', a, ' incorrectly ',
1198 \$ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
1199 \$ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5,
1200 \$ ')' )
1201 9997 FORMAT( ' DCHKHS: Selected ', a, ' Eigenvectors from ', a,
1202 \$ ' do not match other eigenvectors ', 9x, 'N=', i6,
1203 \$ ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1204*
1205* End of DCHKHS
1206*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dget10(m, n, a, lda, b, ldb, work, result)
DGET10
Definition dget10.f:93
subroutine dget22(transa, transe, transw, n, a, lda, e, lde, wr, wi, work, result)
DGET22
Definition dget22.f:168
subroutine dhst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
DHST01
Definition dhst01.f:134
subroutine dlafts(type, m, n, imat, ntests, result, iseed, thresh, iounit, ie)
DLAFTS
Definition dlafts.f:99
subroutine dlasum(type, iounit, ie, nrun)
DLASUM
Definition dlasum.f:43
subroutine dlatme(n, dist, iseed, d, mode, cond, dmax, ei, rsign, upper, sim, ds, modes, conds, kl, ku, anorm, a, lda, work, info)
DLATME
Definition dlatme.f:332
subroutine dlatmr(m, n, dist, iseed, sym, d, mode, cond, dmax, rsign, grade, dl, model, condl, dr, moder, condr, pivtng, ipivot, kl, ku, sparse, anorm, pack, a, lda, iwork, info)
DLATMR
Definition dlatmr.f:471
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
Definition dlatms.f:321
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:167
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dhsein(side, eigsrc, initv, select, n, h, ldh, wr, wi, vl, ldvl, vr, ldvr, mm, m, work, ifaill, ifailr, info)
DHSEIN
Definition dhsein.f:263
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
Definition dhseqr.f:316
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 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 dtrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
DTREVC3
Definition dtrevc3.f:237
subroutine dtrevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, info)
DTREVC
Definition dtrevc.f:222
subroutine dorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
DORGHR
Definition dorghr.f:126
subroutine dormhr(side, trans, m, n, ilo, ihi, a, lda, tau, c, ldc, work, lwork, info)
DORMHR
Definition dormhr.f:178
Here is the call graph for this function:
Here is the caller graph for this function: