LAPACK: TESTING/EIG/cchkhs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00002      $                   NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, W1,
00003      $                   W3, EVECTL, EVECTR, EVECTY, EVECTX, UU, TAU,
00004      $                   WORK, NWORK, RWORK, IWORK, SELECT, RESULT,
00005      $                   INFO )
00006 *
00007 *  -- LAPACK test routine (version 3.1.1) --
00008 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00009 *     February 2007
00010 *
00011 *     .. Scalar Arguments ..
00012       INTEGER            INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
00013       REAL               THRESH
00014 *     ..
00015 *     .. Array Arguments ..
00016       LOGICAL            DOTYPE( * ), SELECT( * )
00017       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
00018       REAL               RESULT( 14 ), RWORK( * )
00019       COMPLEX            A( LDA, * ), EVECTL( LDU, * ),
00020      $                   EVECTR( LDU, * ), EVECTX( LDU, * ),
00021      $                   EVECTY( LDU, * ), H( LDA, * ), T1( LDA, * ),
00022      $                   T2( LDA, * ), TAU( * ), U( LDU, * ),
00023      $                   UU( LDU, * ), UZ( LDU, * ), W1( * ), W3( * ),
00024      $                   WORK( * ), Z( LDU, * )
00025 *     ..
00026 *
00027 *  Purpose
00028 *  =======
00029 *
00030 *     CCHKHS  checks the nonsymmetric eigenvalue problem routines.
00031 *
00032 *             CGEHRD factors A as  U H U' , where ' means conjugate
00033 *             transpose, H is hessenberg, and U is unitary.
00034 *
00035 *             CUNGHR generates the unitary matrix U.
00036 *
00037 *             CUNMHR multiplies a matrix by the unitary matrix U.
00038 *
00039 *             CHSEQR factors H as  Z T Z' , where Z is unitary and T
00040 *             is upper triangular.  It also computes the eigenvalues,
00041 *             w(1), ..., w(n); we define a diagonal matrix W whose
00042 *             (diagonal) entries are the eigenvalues.
00043 *
00044 *             CTREVC computes the left eigenvector matrix L and the
00045 *             right eigenvector matrix R for the matrix T.  The
00046 *             columns of L are the complex conjugates of the left
00047 *             eigenvectors of T.  The columns of R are the right
00048 *             eigenvectors of T.  L is lower triangular, and R is
00049 *             upper triangular.
00050 *
00051 *             CHSEIN computes the left eigenvector matrix Y and the
00052 *             right eigenvector matrix X for the matrix H.  The
00053 *             columns of Y are the complex conjugates of the left
00054 *             eigenvectors of H.  The columns of X are the right
00055 *             eigenvectors of H.  Y is lower triangular, and X is
00056 *             upper triangular.
00057 *
00058 *     When CCHKHS is called, a number of matrix "sizes" ("n's") and a
00059 *     number of matrix "types" are specified.  For each size ("n")
00060 *     and each type of matrix, one matrix will be generated and used
00061 *     to test the nonsymmetric eigenroutines.  For each matrix, 14
00062 *     tests will be performed:
00063 *
00064 *     (1)     | A - U H U**H | / ( |A| n ulp )
00065 *
00066 *     (2)     | I - UU**H | / ( n ulp )
00067 *
00068 *     (3)     | H - Z T Z**H | / ( |H| n ulp )
00069 *
00070 *     (4)     | I - ZZ**H | / ( n ulp )
00071 *
00072 *     (5)     | A - UZ H (UZ)**H | / ( |A| n ulp )
00073 *
00074 *     (6)     | I - UZ (UZ)**H | / ( n ulp )
00075 *
00076 *     (7)     | T(Z computed) - T(Z not computed) | / ( |T| ulp )
00077 *
00078 *     (8)     | W(Z computed) - W(Z not computed) | / ( |W| ulp )
00079 *
00080 *     (9)     | TR - RW | / ( |T| |R| ulp )
00081 *
00082 *     (10)    | L**H T - W**H L | / ( |T| |L| ulp )
00083 *
00084 *     (11)    | HX - XW | / ( |H| |X| ulp )
00085 *
00086 *     (12)    | Y**H H - W**H Y | / ( |H| |Y| ulp )
00087 *
00088 *     (13)    | AX - XW | / ( |A| |X| ulp )
00089 *
00090 *     (14)    | Y**H A - W**H Y | / ( |A| |Y| ulp )
00091 *
00092 *     The "sizes" are specified by an array NN(1:NSIZES); the value of
00093 *     each element NN(j) specifies one size.
00094 *     The "types" are specified by a logical array DOTYPE( 1:NTYPES );
00095 *     if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
00096 *     Currently, the list of possible types is:
00097 *
00098 *     (1)  The zero matrix.
00099 *     (2)  The identity matrix.
00100 *     (3)  A (transposed) Jordan block, with 1's on the diagonal.
00101 *
00102 *     (4)  A diagonal matrix with evenly spaced entries
00103 *          1, ..., ULP  and random complex angles.
00104 *          (ULP = (first number larger than 1) - 1 )
00105 *     (5)  A diagonal matrix with geometrically spaced entries
00106 *          1, ..., ULP  and random complex angles.
00107 *     (6)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
00108 *          and random complex angles.
00109 *
00110 *     (7)  Same as (4), but multiplied by SQRT( overflow threshold )
00111 *     (8)  Same as (4), but multiplied by SQRT( underflow threshold )
00112 *
00113 *     (9)  A matrix of the form  U' T U, where U is unitary and
00114 *          T has evenly spaced entries 1, ..., ULP with random complex
00115 *          angles on the diagonal and random O(1) entries in the upper
00116 *          triangle.
00117 *
00118 *     (10) A matrix of the form  U' T U, where U is unitary and
00119 *          T has geometrically spaced entries 1, ..., ULP with random
00120 *          complex angles on the diagonal and random O(1) entries in
00121 *          the upper triangle.
00122 *
00123 *     (11) A matrix of the form  U' T U, where U is unitary and
00124 *          T has "clustered" entries 1, ULP,..., ULP with random
00125 *          complex angles on the diagonal and random O(1) entries in
00126 *          the upper triangle.
00127 *
00128 *     (12) A matrix of the form  U' T U, where U is unitary and
00129 *          T has complex eigenvalues randomly chosen from
00130 *          ULP < |z| < 1   and random O(1) entries in the upper
00131 *          triangle.
00132 *
00133 *     (13) A matrix of the form  X' T X, where X has condition
00134 *          SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
00135 *          with random complex angles on the diagonal and random O(1)
00136 *          entries in the upper triangle.
00137 *
00138 *     (14) A matrix of the form  X' T X, where X has condition
00139 *          SQRT( ULP ) and T has geometrically spaced entries
00140 *          1, ..., ULP with random complex angles on the diagonal
00141 *          and random O(1) entries in the upper triangle.
00142 *
00143 *     (15) A matrix of the form  X' T X, where X has condition
00144 *          SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
00145 *          with random complex angles on the diagonal and random O(1)
00146 *          entries in the upper triangle.
00147 *
00148 *     (16) A matrix of the form  X' T X, where X has condition
00149 *          SQRT( ULP ) and T has complex eigenvalues randomly chosen
00150 *          from   ULP < |z| < 1   and random O(1) entries in the upper
00151 *          triangle.
00152 *
00153 *     (17) Same as (16), but multiplied by SQRT( overflow threshold )
00154 *     (18) Same as (16), but multiplied by SQRT( underflow threshold )
00155 *
00156 *     (19) Nonsymmetric matrix with random entries chosen from |z| < 1
00157 *     (20) Same as (19), but multiplied by SQRT( overflow threshold )
00158 *     (21) Same as (19), but multiplied by SQRT( underflow threshold )
00159 *
00160 *  Arguments
00161 *  ==========
00162 *
00163 *  NSIZES - INTEGER
00164 *           The number of sizes of matrices to use.  If it is zero,
00165 *           CCHKHS does nothing.  It must be at least zero.
00166 *           Not modified.
00167 *
00168 *  NN     - INTEGER array, dimension (NSIZES)
00169 *           An array containing the sizes to be used for the matrices.
00170 *           Zero values will be skipped.  The values must be at least
00171 *           zero.
00172 *           Not modified.
00173 *
00174 *  NTYPES - INTEGER
00175 *           The number of elements in DOTYPE.   If it is zero, CCHKHS
00176 *           does nothing.  It must be at least zero.  If it is MAXTYP+1
00177 *           and NSIZES is 1, then an additional type, MAXTYP+1 is
00178 *           defined, which is to use whatever matrix is in A.  This
00179 *           is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
00180 *           DOTYPE(MAXTYP+1) is .TRUE. .
00181 *           Not modified.
00182 *
00183 *  DOTYPE - LOGICAL array, dimension (NTYPES)
00184 *           If DOTYPE(j) is .TRUE., then for each size in NN a
00185 *           matrix of that size and of type j will be generated.
00186 *           If NTYPES is smaller than the maximum number of types
00187 *           defined (PARAMETER MAXTYP), then types NTYPES+1 through
00188 *           MAXTYP will not be generated.  If NTYPES is larger
00189 *           than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
00190 *           will be ignored.
00191 *           Not modified.
00192 *
00193 *  ISEED  - INTEGER array, dimension (4)
00194 *           On entry ISEED specifies the seed of the random number
00195 *           generator. The array elements should be between 0 and 4095;
00196 *           if not they will be reduced mod 4096.  Also, ISEED(4) must
00197 *           be odd.  The random number generator uses a linear
00198 *           congruential sequence limited to small integers, and so
00199 *           should produce machine independent random numbers. The
00200 *           values of ISEED are changed on exit, and can be used in the
00201 *           next call to CCHKHS to continue the same random number
00202 *           sequence.
00203 *           Modified.
00204 *
00205 *  THRESH - REAL
00206 *           A test will count as "failed" if the "error", computed as
00207 *           described above, exceeds THRESH.  Note that the error
00208 *           is scaled to be O(1), so THRESH should be a reasonably
00209 *           small multiple of 1, e.g., 10 or 100.  In particular,
00210 *           it should not depend on the precision (single vs. double)
00211 *           or the size of the matrix.  It must be at least zero.
00212 *           Not modified.
00213 *
00214 *  NOUNIT - INTEGER
00215 *           The FORTRAN unit number for printing out error messages
00216 *           (e.g., if a routine returns IINFO not equal to 0.)
00217 *           Not modified.
00218 *
00219 *  A      - COMPLEX array, dimension (LDA,max(NN))
00220 *           Used to hold the matrix whose eigenvalues are to be
00221 *           computed.  On exit, A contains the last matrix actually
00222 *           used.
00223 *           Modified.
00224 *
00225 *  LDA    - INTEGER
00226 *           The leading dimension of A, H, T1 and T2.  It must be at
00227 *           least 1 and at least max( NN ).
00228 *           Not modified.
00229 *
00230 *  H      - COMPLEX array, dimension (LDA,max(NN))
00231 *           The upper hessenberg matrix computed by CGEHRD.  On exit,
00232 *           H contains the Hessenberg form of the matrix in A.
00233 *           Modified.
00234 *
00235 *  T1     - COMPLEX array, dimension (LDA,max(NN))
00236 *           The Schur (="quasi-triangular") matrix computed by CHSEQR
00237 *           if Z is computed.  On exit, T1 contains the Schur form of
00238 *           the matrix in A.
00239 *           Modified.
00240 *
00241 *  T2     - COMPLEX array, dimension (LDA,max(NN))
00242 *           The Schur matrix computed by CHSEQR when Z is not computed.
00243 *           This should be identical to T1.
00244 *           Modified.
00245 *
00246 *  LDU    - INTEGER
00247 *           The leading dimension of U, Z, UZ and UU.  It must be at
00248 *           least 1 and at least max( NN ).
00249 *           Not modified.
00250 *
00251 *  U      - COMPLEX array, dimension (LDU,max(NN))
00252 *           The unitary matrix computed by CGEHRD.
00253 *           Modified.
00254 *
00255 *  Z      - COMPLEX array, dimension (LDU,max(NN))
00256 *           The unitary matrix computed by CHSEQR.
00257 *           Modified.
00258 *
00259 *  UZ     - COMPLEX array, dimension (LDU,max(NN))
00260 *           The product of U times Z.
00261 *           Modified.
00262 *
00263 *  W1     - COMPLEX array, dimension (max(NN))
00264 *           The eigenvalues of A, as computed by a full Schur
00265 *           decomposition H = Z T Z'.  On exit, W1 contains the
00266 *           eigenvalues of the matrix in A.
00267 *           Modified.
00268 *
00269 *  W3     - COMPLEX array, dimension (max(NN))
00270 *           The eigenvalues of A, as computed by a partial Schur
00271 *           decomposition (Z not computed, T only computed as much
00272 *           as is necessary for determining eigenvalues).  On exit,
00273 *           W3 contains the eigenvalues of the matrix in A, possibly
00274 *           perturbed by CHSEIN.
00275 *           Modified.
00276 *
00277 *  EVECTL - COMPLEX array, dimension (LDU,max(NN))
00278 *           The conjugate transpose of the (upper triangular) left
00279 *           eigenvector matrix for the matrix in T1.
00280 *           Modified.
00281 *
00282 *  EVECTR - COMPLEX array, dimension (LDU,max(NN))
00283 *           The (upper triangular) right eigenvector matrix for the
00284 *           matrix in T1.
00285 *           Modified.
00286 *
00287 *  EVECTY - COMPLEX array, dimension (LDU,max(NN))
00288 *           The conjugate transpose of the left eigenvector matrix
00289 *           for the matrix in H.
00290 *           Modified.
00291 *
00292 *  EVECTX - COMPLEX array, dimension (LDU,max(NN))
00293 *           The right eigenvector matrix for the matrix in H.
00294 *           Modified.
00295 *
00296 *  UU     - COMPLEX array, dimension (LDU,max(NN))
00297 *           Details of the unitary matrix computed by CGEHRD.
00298 *           Modified.
00299 *
00300 *  TAU    - COMPLEX array, dimension (max(NN))
00301 *           Further details of the unitary matrix computed by CGEHRD.
00302 *           Modified.
00303 *
00304 *  WORK   - COMPLEX array, dimension (NWORK)
00305 *           Workspace.
00306 *           Modified.
00307 *
00308 *  NWORK  - INTEGER
00309 *           The number of entries in WORK.  NWORK >= 4*NN(j)*NN(j) + 2.
00310 *
00311 *  RWORK  - REAL array, dimension (max(NN))
00312 *           Workspace.  Could be equivalenced to IWORK, but not SELECT.
00313 *           Modified.
00314 *
00315 *  IWORK  - INTEGER array, dimension (max(NN))
00316 *           Workspace.
00317 *           Modified.
00318 *
00319 *  SELECT - LOGICAL array, dimension (max(NN))
00320 *           Workspace.  Could be equivalenced to IWORK, but not RWORK.
00321 *           Modified.
00322 *
00323 *  RESULT - REAL array, dimension (14)
00324 *           The values computed by the fourteen tests described above.
00325 *           The values are currently limited to 1/ulp, to avoid
00326 *           overflow.
00327 *           Modified.
00328 *
00329 *  INFO   - INTEGER
00330 *           If 0, then everything ran OK.
00331 *            -1: NSIZES < 0
00332 *            -2: Some NN(j) < 0
00333 *            -3: NTYPES < 0
00334 *            -6: THRESH < 0
00335 *            -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
00336 *           -14: LDU < 1 or LDU < NMAX.
00337 *           -26: NWORK too small.
00338 *           If  CLATMR, CLATMS, or CLATME returns an error code, the
00339 *               absolute value of it is returned.
00340 *           If 1, then CHSEQR could not find all the shifts.
00341 *           If 2, then the EISPACK code (for small blocks) failed.
00342 *           If >2, then 30*N iterations were not enough to find an
00343 *               eigenvalue or to decompose the problem.
00344 *           Modified.
00345 *
00346 *-----------------------------------------------------------------------
00347 *
00348 *     Some Local Variables and Parameters:
00349 *     ---- ----- --------- --- ----------
00350 *
00351 *     ZERO, ONE       Real 0 and 1.
00352 *     MAXTYP          The number of types defined.
00353 *     MTEST           The number of tests defined: care must be taken
00354 *                     that (1) the size of RESULT, (2) the number of
00355 *                     tests actually performed, and (3) MTEST agree.
00356 *     NTEST           The number of tests performed on this matrix
00357 *                     so far.  This should be less than MTEST, and
00358 *                     equal to it by the last test.  It will be less
00359 *                     if any of the routines being tested indicates
00360 *                     that it could not compute the matrices that
00361 *                     would be tested.
00362 *     NMAX            Largest value in NN.
00363 *     NMATS           The number of matrices generated so far.
00364 *     NERRS           The number of tests which have exceeded THRESH
00365 *                     so far (computed by SLAFTS).
00366 *     COND, CONDS,
00367 *     IMODE           Values to be passed to the matrix generators.
00368 *     ANORM           Norm of A; passed to matrix generators.
00369 *
00370 *     OVFL, UNFL      Overflow and underflow thresholds.
00371 *     ULP, ULPINV     Finest relative precision and its inverse.
00372 *     RTOVFL, RTUNFL,
00373 *     RTULP, RTULPI   Square roots of the previous 4 values.
00374 *
00375 *             The following four arrays decode JTYPE:
00376 *     KTYPE(j)        The general type (1-10) for type "j".
00377 *     KMODE(j)        The MODE value to be passed to the matrix
00378 *                     generator for type "j".
00379 *     KMAGN(j)        The order of magnitude ( O(1),
00380 *                     O(overflow^(1/2) ), O(underflow^(1/2) )
00381 *     KCONDS(j)       Selects whether CONDS is to be 1 or
00382 *                     1/sqrt(ulp).  (0 means irrelevant.)
00383 *
00384 *  =====================================================================
00385 *
00386 *     .. Parameters ..
00387       REAL               ZERO, ONE
00388       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00389       COMPLEX            CZERO, CONE
00390       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00391      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00392       INTEGER            MAXTYP
00393       PARAMETER          ( MAXTYP = 21 )
00394 *     ..
00395 *     .. Local Scalars ..
00396       LOGICAL            BADNN, MATCH
00397       INTEGER            I, IHI, IINFO, ILO, IMODE, IN, ITYPE, J, JCOL,
00398      $                   JJ, JSIZE, JTYPE, K, MTYPES, N, N1, NERRS,
00399      $                   NMATS, NMAX, NTEST, NTESTT
00400       REAL               ANINV, ANORM, COND, CONDS, OVFL, RTOVFL, RTULP,
00401      $                   RTULPI, RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL
00402 *     ..
00403 *     .. Local Arrays ..
00404       INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
00405      $                   KMAGN( MAXTYP ), KMODE( MAXTYP ),
00406      $                   KTYPE( MAXTYP )
00407       REAL               DUMMA( 4 )
00408       COMPLEX            CDUMMA( 4 )
00409 *     ..
00410 *     .. External Functions ..
00411       REAL               SLAMCH
00412       EXTERNAL           SLAMCH
00413 *     ..
00414 *     .. External Subroutines ..
00415       EXTERNAL           CCOPY, CGEHRD, CGEMM, CGET10, CGET22, CHSEIN,
00416      $                   CHSEQR, CHST01, CLACPY, CLASET, CLATME, CLATMR,
00417      $                   CLATMS, CTREVC, CUNGHR, CUNMHR, SLABAD, SLAFTS,
00418      $                   SLASUM, XERBLA
00419 *     ..
00420 *     .. Intrinsic Functions ..
00421       INTRINSIC          ABS, MAX, MIN, REAL, SQRT
00422 *     ..
00423 *     .. Data statements ..
00424       DATA               KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
00425       DATA               KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
00426      $                   3, 1, 2, 3 /
00427       DATA               KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
00428      $                   1, 5, 5, 5, 4, 3, 1 /
00429       DATA               KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
00430 *     ..
00431 *     .. Executable Statements ..
00432 *
00433 *     Check for errors
00434 *
00435       NTESTT = 0
00436       INFO = 0
00437 *
00438       BADNN = .FALSE.
00439       NMAX = 0
00440       DO 10 J = 1, NSIZES
00441          NMAX = MAX( NMAX, NN( J ) )
00442          IF( NN( J ).LT.0 )
00443      $      BADNN = .TRUE.
00444    10 CONTINUE
00445 *
00446 *     Check for errors
00447 *
00448       IF( NSIZES.LT.0 ) THEN
00449          INFO = -1
00450       ELSE IF( BADNN ) THEN
00451          INFO = -2
00452       ELSE IF( NTYPES.LT.0 ) THEN
00453          INFO = -3
00454       ELSE IF( THRESH.LT.ZERO ) THEN
00455          INFO = -6
00456       ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
00457          INFO = -9
00458       ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
00459          INFO = -14
00460       ELSE IF( 4*NMAX*NMAX+2.GT.NWORK ) THEN
00461          INFO = -26
00462       END IF
00463 *
00464       IF( INFO.NE.0 ) THEN
00465          CALL XERBLA( 'CCHKHS', -INFO )
00466          RETURN
00467       END IF
00468 *
00469 *     Quick return if possible
00470 *
00471       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
00472      $   RETURN
00473 *
00474 *     More important constants
00475 *
00476       UNFL = SLAMCH( 'Safe minimum' )
00477       OVFL = SLAMCH( 'Overflow' )
00478       CALL SLABAD( UNFL, OVFL )
00479       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00480       ULPINV = ONE / ULP
00481       RTUNFL = SQRT( UNFL )
00482       RTOVFL = SQRT( OVFL )
00483       RTULP = SQRT( ULP )
00484       RTULPI = ONE / RTULP
00485 *
00486 *     Loop over sizes, types
00487 *
00488       NERRS = 0
00489       NMATS = 0
00490 *
00491       DO 260 JSIZE = 1, NSIZES
00492          N = NN( JSIZE )
00493          N1 = MAX( 1, N )
00494          ANINV = ONE / REAL( N1 )
00495 *
00496          IF( NSIZES.NE.1 ) THEN
00497             MTYPES = MIN( MAXTYP, NTYPES )
00498          ELSE
00499             MTYPES = MIN( MAXTYP+1, NTYPES )
00500          END IF
00501 *
00502          DO 250 JTYPE = 1, MTYPES
00503             IF( .NOT.DOTYPE( JTYPE ) )
00504      $         GO TO 250
00505             NMATS = NMATS + 1
00506             NTEST = 0
00507 *
00508 *           Save ISEED in case of an error.
00509 *
00510             DO 20 J = 1, 4
00511                IOLDSD( J ) = ISEED( J )
00512    20       CONTINUE
00513 *
00514 *           Initialize RESULT
00515 *
00516             DO 30 J = 1, 14
00517                RESULT( J ) = ZERO
00518    30       CONTINUE
00519 *
00520 *           Compute "A"
00521 *
00522 *           Control parameters:
00523 *
00524 *           KMAGN  KCONDS  KMODE        KTYPE
00525 *       =1  O(1)   1       clustered 1  zero
00526 *       =2  large  large   clustered 2  identity
00527 *       =3  small          exponential  Jordan
00528 *       =4                 arithmetic   diagonal, (w/ eigenvalues)
00529 *       =5                 random log   hermitian, w/ eigenvalues
00530 *       =6                 random       general, w/ eigenvalues
00531 *       =7                              random diagonal
00532 *       =8                              random hermitian
00533 *       =9                              random general
00534 *       =10                             random triangular
00535 *
00536             IF( MTYPES.GT.MAXTYP )
00537      $         GO TO 100
00538 *
00539             ITYPE = KTYPE( JTYPE )
00540             IMODE = KMODE( JTYPE )
00541 *
00542 *           Compute norm
00543 *
00544             GO TO ( 40, 50, 60 )KMAGN( JTYPE )
00545 *
00546    40       CONTINUE
00547             ANORM = ONE
00548             GO TO 70
00549 *
00550    50       CONTINUE
00551             ANORM = ( RTOVFL*ULP )*ANINV
00552             GO TO 70
00553 *
00554    60       CONTINUE
00555             ANORM = RTUNFL*N*ULPINV
00556             GO TO 70
00557 *
00558    70       CONTINUE
00559 *
00560             CALL CLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
00561             IINFO = 0
00562             COND = ULPINV
00563 *
00564 *           Special Matrices
00565 *
00566             IF( ITYPE.EQ.1 ) THEN
00567 *
00568 *              Zero
00569 *
00570                IINFO = 0
00571             ELSE IF( ITYPE.EQ.2 ) THEN
00572 *
00573 *              Identity
00574 *
00575                DO 80 JCOL = 1, N
00576                   A( JCOL, JCOL ) = ANORM
00577    80          CONTINUE
00578 *
00579             ELSE IF( ITYPE.EQ.3 ) THEN
00580 *
00581 *              Jordan Block
00582 *
00583                DO 90 JCOL = 1, N
00584                   A( JCOL, JCOL ) = ANORM
00585                   IF( JCOL.GT.1 )
00586      $               A( JCOL, JCOL-1 ) = ONE
00587    90          CONTINUE
00588 *
00589             ELSE IF( ITYPE.EQ.4 ) THEN
00590 *
00591 *              Diagonal Matrix, [Eigen]values Specified
00592 *
00593                CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, IMODE, COND,
00594      $                      CONE, 'T', 'N', WORK( N+1 ), 1, ONE,
00595      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
00596      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00597 *
00598             ELSE IF( ITYPE.EQ.5 ) THEN
00599 *
00600 *              Hermitian, eigenvalues specified
00601 *
00602                CALL CLATMS( N, N, 'D', ISEED, 'H', RWORK, IMODE, COND,
00603      $                      ANORM, N, N, 'N', A, LDA, WORK, IINFO )
00604 *
00605             ELSE IF( ITYPE.EQ.6 ) THEN
00606 *
00607 *              General, eigenvalues specified
00608 *
00609                IF( KCONDS( JTYPE ).EQ.1 ) THEN
00610                   CONDS = ONE
00611                ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
00612                   CONDS = RTULPI
00613                ELSE
00614                   CONDS = ZERO
00615                END IF
00616 *
00617                CALL CLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ',
00618      $                      'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
00619      $                      A, LDA, WORK( N+1 ), IINFO )
00620 *
00621             ELSE IF( ITYPE.EQ.7 ) THEN
00622 *
00623 *              Diagonal, random eigenvalues
00624 *
00625                CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00626      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00627      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
00628      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00629 *
00630             ELSE IF( ITYPE.EQ.8 ) THEN
00631 *
00632 *              Hermitian, random eigenvalues
00633 *
00634                CALL CLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
00635      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00636      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00637      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00638 *
00639             ELSE IF( ITYPE.EQ.9 ) THEN
00640 *
00641 *              General, random eigenvalues
00642 *
00643                CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00644      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00645      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
00646      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00647 *
00648             ELSE IF( ITYPE.EQ.10 ) THEN
00649 *
00650 *              Triangular, random eigenvalues
00651 *
00652                CALL CLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
00653      $                      'T', 'N', WORK( N+1 ), 1, ONE,
00654      $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
00655      $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
00656 *
00657             ELSE
00658 *
00659                IINFO = 1
00660             END IF
00661 *
00662             IF( IINFO.NE.0 ) THEN
00663                WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
00664      $            IOLDSD
00665                INFO = ABS( IINFO )
00666                RETURN
00667             END IF
00668 *
00669   100       CONTINUE
00670 *
00671 *           Call CGEHRD to compute H and U, do tests.
00672 *
00673             CALL CLACPY( ' ', N, N, A, LDA, H, LDA )
00674             NTEST = 1
00675 *
00676             ILO = 1
00677             IHI = N
00678 *
00679             CALL CGEHRD( N, ILO, IHI, H, LDA, WORK, WORK( N+1 ),
00680      $                   NWORK-N, IINFO )
00681 *
00682             IF( IINFO.NE.0 ) THEN
00683                RESULT( 1 ) = ULPINV
00684                WRITE( NOUNIT, FMT = 9999 )'CGEHRD', IINFO, N, JTYPE,
00685      $            IOLDSD
00686                INFO = ABS( IINFO )
00687                GO TO 240
00688             END IF
00689 *
00690             DO 120 J = 1, N - 1
00691                UU( J+1, J ) = CZERO
00692                DO 110 I = J + 2, N
00693                   U( I, J ) = H( I, J )
00694                   UU( I, J ) = H( I, J )
00695                   H( I, J ) = CZERO
00696   110          CONTINUE
00697   120       CONTINUE
00698             CALL CCOPY( N-1, WORK, 1, TAU, 1 )
00699             CALL CUNGHR( N, ILO, IHI, U, LDU, WORK, WORK( N+1 ),
00700      $                   NWORK-N, IINFO )
00701             NTEST = 2
00702 *
00703             CALL CHST01( N, ILO, IHI, A, LDA, H, LDA, U, LDU, WORK,
00704      $                   NWORK, RWORK, RESULT( 1 ) )
00705 *
00706 *           Call CHSEQR to compute T1, T2 and Z, do tests.
00707 *
00708 *           Eigenvalues only (W3)
00709 *
00710             CALL CLACPY( ' ', N, N, H, LDA, T2, LDA )
00711             NTEST = 3
00712             RESULT( 3 ) = ULPINV
00713 *
00714             CALL CHSEQR( 'E', 'N', N, ILO, IHI, T2, LDA, W3, UZ, LDU,
00715      $                   WORK, NWORK, IINFO )
00716             IF( IINFO.NE.0 ) THEN
00717                WRITE( NOUNIT, FMT = 9999 )'CHSEQR(E)', IINFO, N, JTYPE,
00718      $            IOLDSD
00719                IF( IINFO.LE.N+2 ) THEN
00720                   INFO = ABS( IINFO )
00721                   GO TO 240
00722                END IF
00723             END IF
00724 *
00725 *           Eigenvalues (W1) and Full Schur Form (T2)
00726 *
00727             CALL CLACPY( ' ', N, N, H, LDA, T2, LDA )
00728 *
00729             CALL CHSEQR( 'S', 'N', N, ILO, IHI, T2, LDA, W1, UZ, LDU,
00730      $                   WORK, NWORK, IINFO )
00731             IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
00732                WRITE( NOUNIT, FMT = 9999 )'CHSEQR(S)', IINFO, N, JTYPE,
00733      $            IOLDSD
00734                INFO = ABS( IINFO )
00735                GO TO 240
00736             END IF
00737 *
00738 *           Eigenvalues (W1), Schur Form (T1), and Schur Vectors (UZ)
00739 *
00740             CALL CLACPY( ' ', N, N, H, LDA, T1, LDA )
00741             CALL CLACPY( ' ', N, N, U, LDU, UZ, LDU )
00742 *
00743             CALL CHSEQR( 'S', 'V', N, ILO, IHI, T1, LDA, W1, UZ, LDU,
00744      $                   WORK, NWORK, IINFO )
00745             IF( IINFO.NE.0 .AND. IINFO.LE.N+2 ) THEN
00746                WRITE( NOUNIT, FMT = 9999 )'CHSEQR(V)', IINFO, N, JTYPE,
00747      $            IOLDSD
00748                INFO = ABS( IINFO )
00749                GO TO 240
00750             END IF
00751 *
00752 *           Compute Z = U' UZ
00753 *
00754             CALL CGEMM( 'C', 'N', N, N, N, CONE, U, LDU, UZ, LDU, CZERO,
00755      $                  Z, LDU )
00756             NTEST = 8
00757 *
00758 *           Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
00759 *                and 4: | I - Z Z' | / ( n ulp )
00760 *
00761             CALL CHST01( N, ILO, IHI, H, LDA, T1, LDA, Z, LDU, WORK,
00762      $                   NWORK, RWORK, RESULT( 3 ) )
00763 *
00764 *           Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
00765 *                and 6: | I - UZ (UZ)' | / ( n ulp )
00766 *
00767             CALL CHST01( N, ILO, IHI, A, LDA, T1, LDA, UZ, LDU, WORK,
00768      $                   NWORK, RWORK, RESULT( 5 ) )
00769 *
00770 *           Do Test 7: | T2 - T1 | / ( |T| n ulp )
00771 *
00772             CALL CGET10( N, N, T2, LDA, T1, LDA, WORK, RWORK,
00773      $                   RESULT( 7 ) )
00774 *
00775 *           Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp )
00776 *
00777             TEMP1 = ZERO
00778             TEMP2 = ZERO
00779             DO 130 J = 1, N
00780                TEMP1 = MAX( TEMP1, ABS( W1( J ) ), ABS( W3( J ) ) )
00781                TEMP2 = MAX( TEMP2, ABS( W1( J )-W3( J ) ) )
00782   130       CONTINUE
00783 *
00784             RESULT( 8 ) = TEMP2 / MAX( UNFL, ULP*MAX( TEMP1, TEMP2 ) )
00785 *
00786 *           Compute the Left and Right Eigenvectors of T
00787 *
00788 *           Compute the Right eigenvector Matrix:
00789 *
00790             NTEST = 9
00791             RESULT( 9 ) = ULPINV
00792 *
00793 *           Select every other eigenvector
00794 *
00795             DO 140 J = 1, N
00796                SELECT( J ) = .FALSE.
00797   140       CONTINUE
00798             DO 150 J = 1, N, 2
00799                SELECT( J ) = .TRUE.
00800   150       CONTINUE
00801             CALL CTREVC( 'Right', 'All', SELECT, N, T1, LDA, CDUMMA,
00802      $                   LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
00803             IF( IINFO.NE.0 ) THEN
00804                WRITE( NOUNIT, FMT = 9999 )'CTREVC(R,A)', IINFO, N,
00805      $            JTYPE, IOLDSD
00806                INFO = ABS( IINFO )
00807                GO TO 240
00808             END IF
00809 *
00810 *           Test 9:  | TR - RW | / ( |T| |R| ulp )
00811 *
00812             CALL CGET22( 'N', 'N', 'N', N, T1, LDA, EVECTR, LDU, W1,
00813      $                   WORK, RWORK, DUMMA( 1 ) )
00814             RESULT( 9 ) = DUMMA( 1 )
00815             IF( DUMMA( 2 ).GT.THRESH ) THEN
00816                WRITE( NOUNIT, FMT = 9998 )'Right', 'CTREVC',
00817      $            DUMMA( 2 ), N, JTYPE, IOLDSD
00818             END IF
00819 *
00820 *           Compute selected right eigenvectors and confirm that
00821 *           they agree with previous right eigenvectors
00822 *
00823             CALL CTREVC( 'Right', 'Some', SELECT, N, T1, LDA, CDUMMA,
00824      $                   LDU, EVECTL, LDU, N, IN, WORK, RWORK, IINFO )
00825             IF( IINFO.NE.0 ) THEN
00826                WRITE( NOUNIT, FMT = 9999 )'CTREVC(R,S)', IINFO, N,
00827      $            JTYPE, IOLDSD
00828                INFO = ABS( IINFO )
00829                GO TO 240
00830             END IF
00831 *
00832             K = 1
00833             MATCH = .TRUE.
00834             DO 170 J = 1, N
00835                IF( SELECT( J ) ) THEN
00836                   DO 160 JJ = 1, N
00837                      IF( EVECTR( JJ, J ).NE.EVECTL( JJ, K ) ) THEN
00838                         MATCH = .FALSE.
00839                         GO TO 180
00840                      END IF
00841   160             CONTINUE
00842                   K = K + 1
00843                END IF
00844   170       CONTINUE
00845   180       CONTINUE
00846             IF( .NOT.MATCH )
00847      $         WRITE( NOUNIT, FMT = 9997 )'Right', 'CTREVC', N, JTYPE,
00848      $         IOLDSD
00849 *
00850 *           Compute the Left eigenvector Matrix:
00851 *
00852             NTEST = 10
00853             RESULT( 10 ) = ULPINV
00854             CALL CTREVC( 'Left', 'All', SELECT, N, T1, LDA, EVECTL, LDU,
00855      $                   CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
00856             IF( IINFO.NE.0 ) THEN
00857                WRITE( NOUNIT, FMT = 9999 )'CTREVC(L,A)', IINFO, N,
00858      $            JTYPE, IOLDSD
00859                INFO = ABS( IINFO )
00860                GO TO 240
00861             END IF
00862 *
00863 *           Test 10:  | LT - WL | / ( |T| |L| ulp )
00864 *
00865             CALL CGET22( 'C', 'N', 'C', N, T1, LDA, EVECTL, LDU, W1,
00866      $                   WORK, RWORK, DUMMA( 3 ) )
00867             RESULT( 10 ) = DUMMA( 3 )
00868             IF( DUMMA( 4 ).GT.THRESH ) THEN
00869                WRITE( NOUNIT, FMT = 9998 )'Left', 'CTREVC', DUMMA( 4 ),
00870      $            N, JTYPE, IOLDSD
00871             END IF
00872 *
00873 *           Compute selected left eigenvectors and confirm that
00874 *           they agree with previous left eigenvectors
00875 *
00876             CALL CTREVC( 'Left', 'Some', SELECT, N, T1, LDA, EVECTR,
00877      $                   LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
00878             IF( IINFO.NE.0 ) THEN
00879                WRITE( NOUNIT, FMT = 9999 )'CTREVC(L,S)', IINFO, N,
00880      $            JTYPE, IOLDSD
00881                INFO = ABS( IINFO )
00882                GO TO 240
00883             END IF
00884 *
00885             K = 1
00886             MATCH = .TRUE.
00887             DO 200 J = 1, N
00888                IF( SELECT( J ) ) THEN
00889                   DO 190 JJ = 1, N
00890                      IF( EVECTL( JJ, J ).NE.EVECTR( JJ, K ) ) THEN
00891                         MATCH = .FALSE.
00892                         GO TO 210
00893                      END IF
00894   190             CONTINUE
00895                   K = K + 1
00896                END IF
00897   200       CONTINUE
00898   210       CONTINUE
00899             IF( .NOT.MATCH )
00900      $         WRITE( NOUNIT, FMT = 9997 )'Left', 'CTREVC', N, JTYPE,
00901      $         IOLDSD
00902 *
00903 *           Call CHSEIN for Right eigenvectors of H, do test 11
00904 *
00905             NTEST = 11
00906             RESULT( 11 ) = ULPINV
00907             DO 220 J = 1, N
00908                SELECT( J ) = .TRUE.
00909   220       CONTINUE
00910 *
00911             CALL CHSEIN( 'Right', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3,
00912      $                   CDUMMA, LDU, EVECTX, LDU, N1, IN, WORK, RWORK,
00913      $                   IWORK, IWORK, IINFO )
00914             IF( IINFO.NE.0 ) THEN
00915                WRITE( NOUNIT, FMT = 9999 )'CHSEIN(R)', IINFO, N, JTYPE,
00916      $            IOLDSD
00917                INFO = ABS( IINFO )
00918                IF( IINFO.LT.0 )
00919      $            GO TO 240
00920             ELSE
00921 *
00922 *              Test 11:  | HX - XW | / ( |H| |X| ulp )
00923 *
00924 *                        (from inverse iteration)
00925 *
00926                CALL CGET22( 'N', 'N', 'N', N, H, LDA, EVECTX, LDU, W3,
00927      $                      WORK, RWORK, DUMMA( 1 ) )
00928                IF( DUMMA( 1 ).LT.ULPINV )
00929      $            RESULT( 11 ) = DUMMA( 1 )*ANINV
00930                IF( DUMMA( 2 ).GT.THRESH ) THEN
00931                   WRITE( NOUNIT, FMT = 9998 )'Right', 'CHSEIN',
00932      $               DUMMA( 2 ), N, JTYPE, IOLDSD
00933                END IF
00934             END IF
00935 *
00936 *           Call CHSEIN for Left eigenvectors of H, do test 12
00937 *
00938             NTEST = 12
00939             RESULT( 12 ) = ULPINV
00940             DO 230 J = 1, N
00941                SELECT( J ) = .TRUE.
00942   230       CONTINUE
00943 *
00944             CALL CHSEIN( 'Left', 'Qr', 'Ninitv', SELECT, N, H, LDA, W3,
00945      $                   EVECTY, LDU, CDUMMA, LDU, N1, IN, WORK, RWORK,
00946      $                   IWORK, IWORK, IINFO )
00947             IF( IINFO.NE.0 ) THEN
00948                WRITE( NOUNIT, FMT = 9999 )'CHSEIN(L)', IINFO, N, JTYPE,
00949      $            IOLDSD
00950                INFO = ABS( IINFO )
00951                IF( IINFO.LT.0 )
00952      $            GO TO 240
00953             ELSE
00954 *
00955 *              Test 12:  | YH - WY | / ( |H| |Y| ulp )
00956 *
00957 *                        (from inverse iteration)
00958 *
00959                CALL CGET22( 'C', 'N', 'C', N, H, LDA, EVECTY, LDU, W3,
00960      $                      WORK, RWORK, DUMMA( 3 ) )
00961                IF( DUMMA( 3 ).LT.ULPINV )
00962      $            RESULT( 12 ) = DUMMA( 3 )*ANINV
00963                IF( DUMMA( 4 ).GT.THRESH ) THEN
00964                   WRITE( NOUNIT, FMT = 9998 )'Left', 'CHSEIN',
00965      $               DUMMA( 4 ), N, JTYPE, IOLDSD
00966                END IF
00967             END IF
00968 *
00969 *           Call CUNMHR for Right eigenvectors of A, do test 13
00970 *
00971             NTEST = 13
00972             RESULT( 13 ) = ULPINV
00973 *
00974             CALL CUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
00975      $                   LDU, TAU, EVECTX, LDU, WORK, NWORK, IINFO )
00976             IF( IINFO.NE.0 ) THEN
00977                WRITE( NOUNIT, FMT = 9999 )'CUNMHR(L)', IINFO, N, JTYPE,
00978      $            IOLDSD
00979                INFO = ABS( IINFO )
00980                IF( IINFO.LT.0 )
00981      $            GO TO 240
00982             ELSE
00983 *
00984 *              Test 13:  | AX - XW | / ( |A| |X| ulp )
00985 *
00986 *                        (from inverse iteration)
00987 *
00988                CALL CGET22( 'N', 'N', 'N', N, A, LDA, EVECTX, LDU, W3,
00989      $                      WORK, RWORK, DUMMA( 1 ) )
00990                IF( DUMMA( 1 ).LT.ULPINV )
00991      $            RESULT( 13 ) = DUMMA( 1 )*ANINV
00992             END IF
00993 *
00994 *           Call CUNMHR for Left eigenvectors of A, do test 14
00995 *
00996             NTEST = 14
00997             RESULT( 14 ) = ULPINV
00998 *
00999             CALL CUNMHR( 'Left', 'No transpose', N, N, ILO, IHI, UU,
01000      $                   LDU, TAU, EVECTY, LDU, WORK, NWORK, IINFO )
01001             IF( IINFO.NE.0 ) THEN
01002                WRITE( NOUNIT, FMT = 9999 )'CUNMHR(L)', IINFO, N, JTYPE,
01003      $            IOLDSD
01004                INFO = ABS( IINFO )
01005                IF( IINFO.LT.0 )
01006      $            GO TO 240
01007             ELSE
01008 *
01009 *              Test 14:  | YA - WY | / ( |A| |Y| ulp )
01010 *
01011 *                        (from inverse iteration)
01012 *
01013                CALL CGET22( 'C', 'N', 'C', N, A, LDA, EVECTY, LDU, W3,
01014      $                      WORK, RWORK, DUMMA( 3 ) )
01015                IF( DUMMA( 3 ).LT.ULPINV )
01016      $            RESULT( 14 ) = DUMMA( 3 )*ANINV
01017             END IF
01018 *
01019 *           End of Loop -- Check for RESULT(j) > THRESH
01020 *
01021   240       CONTINUE
01022 *
01023             NTESTT = NTESTT + NTEST
01024             CALL SLAFTS( 'CHS', N, N, JTYPE, NTEST, RESULT, IOLDSD,
01025      $                   THRESH, NOUNIT, NERRS )
01026 *
01027   250    CONTINUE
01028   260 CONTINUE
01029 *
01030 *     Summary
01031 *
01032       CALL SLASUM( 'CHS', NOUNIT, NERRS, NTESTT )
01033 *
01034       RETURN
01035 *
01036  9999 FORMAT( ' CCHKHS: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
01037      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
01038  9998 FORMAT( ' CCHKHS: ', A, ' Eigenvectors from ', A, ' incorrectly ',
01039      $      'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
01040      $      'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
01041      $      ')' )
01042  9997 FORMAT( ' CCHKHS: Selected ', A, ' Eigenvectors from ', A,
01043      $      ' do not match other eigenvectors ', 9X, 'N=', I6,
01044      $      ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
01045 *
01046 *     End of CCHKHS
01047 *
01048       END