schkhs.c
Go to the documentation of this file.
00001 /* schkhs.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static real c_b18 = 0.f;
00019 static integer c__0 = 0;
00020 static real c_b32 = 1.f;
00021 static integer c__4 = 4;
00022 static integer c__6 = 6;
00023 static integer c__1 = 1;
00024 
00025 /* Subroutine */ int schkhs_(integer *nsizes, integer *nn, integer *ntypes, 
00026         logical *dotype, integer *iseed, real *thresh, integer *nounit, real *
00027         a, integer *lda, real *h__, real *t1, real *t2, real *u, integer *ldu, 
00028          real *z__, real *uz, real *wr1, real *wi1, real *wr3, real *wi3, 
00029         real *evectl, real *evectr, real *evecty, real *evectx, real *uu, 
00030         real *tau, real *work, integer *nwork, integer *iwork, logical *
00031         select, real *result, integer *info)
00032 {
00033     /* Initialized data */
00034 
00035     static integer ktype[21] = { 1,2,3,4,4,4,4,4,6,6,6,6,6,6,6,6,6,6,9,9,9 };
00036     static integer kmagn[21] = { 1,1,1,1,1,1,2,3,1,1,1,1,1,1,1,1,2,3,1,2,3 };
00037     static integer kmode[21] = { 0,0,0,4,3,1,4,4,4,3,1,5,4,3,1,5,5,5,4,3,1 };
00038     static integer kconds[21] = { 0,0,0,0,0,0,0,0,1,1,1,1,2,2,2,2,2,2,0,0,0 };
00039 
00040     /* Format strings */
00041     static char fmt_9999[] = "(\002 SCHKHS: \002,a,\002 returned INFO=\002,i"
00042             "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, ISEED="
00043             "(\002,3(i5,\002,\002),i5,\002)\002)";
00044     static char fmt_9998[] = "(\002 SCHKHS: \002,a,\002 Eigenvectors from"
00045             " \002,a,\002 incorrectly \002,\002normalized.\002,/\002 Bits of "
00046             "error=\002,0p,g10.3,\002,\002,9x,\002N=\002,i6,\002, JTYPE=\002,"
00047             "i6,\002, ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
00048     static char fmt_9997[] = "(\002 SCHKHS: Selected \002,a,\002 Eigenvector"
00049             "s from \002,a,\002 do not match other eigenvectors \002,9x,\002N="
00050             "\002,i6,\002, JTYPE=\002,i6,\002, ISEED=(\002,3(i5,\002,\002),i5,"
00051             "\002)\002)";
00052 
00053     /* System generated locals */
00054     integer a_dim1, a_offset, evectl_dim1, evectl_offset, evectr_dim1, 
00055             evectr_offset, evectx_dim1, evectx_offset, evecty_dim1, 
00056             evecty_offset, h_dim1, h_offset, t1_dim1, t1_offset, t2_dim1, 
00057             t2_offset, u_dim1, u_offset, uu_dim1, uu_offset, uz_dim1, 
00058             uz_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4;
00059     real r__1, r__2, r__3, r__4, r__5, r__6;
00060 
00061     /* Builtin functions */
00062     double sqrt(doublereal);
00063     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00064 
00065     /* Local variables */
00066     integer i__, j, k, n, n1, jj, in, ihi, ilo;
00067     real ulp, cond;
00068     integer jcol, nmax;
00069     real unfl, ovfl, temp1, temp2;
00070     logical badnn, match;
00071     integer imode;
00072     real dumma[6];
00073     integer iinfo, nselc;
00074     real conds;
00075     extern /* Subroutine */ int sget10_(integer *, integer *, real *, integer 
00076             *, real *, integer *, real *, real *), sgemm_(char *, char *, 
00077             integer *, integer *, integer *, real *, real *, integer *, real *
00078 , integer *, real *, real *, integer *), sget22_(
00079             char *, char *, char *, integer *, real *, integer *, real *, 
00080             integer *, real *, real *, real *, real *)
00081             ;
00082     real aninv, anorm;
00083     integer nmats, nselr, jsize;
00084     extern /* Subroutine */ int shst01_(integer *, integer *, integer *, real 
00085             *, integer *, real *, integer *, real *, integer *, real *, 
00086             integer *, real *);
00087     integer nerrs, itype, jtype, ntest;
00088     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
00089             integer *);
00090     real rtulp;
00091     extern /* Subroutine */ int slabad_(real *, real *);
00092     char adumma[1*1];
00093     extern doublereal slamch_(char *);
00094     integer idumma[1];
00095     extern /* Subroutine */ int sgehrd_(integer *, integer *, integer *, real 
00096             *, integer *, real *, real *, integer *, integer *);
00097     integer ioldsd[4];
00098     extern /* Subroutine */ int xerbla_(char *, integer *), slatme_(
00099             integer *, char *, integer *, real *, integer *, real *, real *, 
00100             char *, char *, char *, char *, real *, integer *, real *, 
00101             integer *, integer *, real *, real *, integer *, real *, integer *
00102 ), shsein_(char *, char *, 
00103              char *, logical *, integer *, real *, integer *, real *, real *, 
00104             real *, integer *, real *, integer *, integer *, integer *, real *
00105 , integer *, integer *, integer *), 
00106             slacpy_(char *, integer *, integer *, real *, integer *, real *, 
00107             integer *), slafts_(char *, integer *, integer *, integer 
00108             *, integer *, real *, integer *, real *, integer *, integer *), slaset_(char *, integer *, integer *, real *, real *, 
00109             real *, integer *), slatmr_(integer *, integer *, char *, 
00110             integer *, char *, real *, integer *, real *, real *, char *, 
00111             char *, real *, integer *, real *, real *, integer *, real *, 
00112             char *, integer *, integer *, integer *, real *, real *, char *, 
00113             real *, integer *, integer *, integer *), slatms_(integer *, integer *, char *, 
00114             integer *, char *, real *, integer *, real *, real *, integer *, 
00115             integer *, char *, real *, integer *, real *, integer *), shseqr_(char *, char *, integer *, integer *, 
00116             integer *, real *, integer *, real *, real *, real *, integer *, 
00117             real *, integer *, integer *), slasum_(char *, 
00118             integer *, integer *, integer *), sorghr_(integer *, 
00119             integer *, integer *, real *, integer *, real *, real *, integer *
00120 , integer *), strevc_(char *, char *, logical *, integer *, real *
00121 , integer *, real *, integer *, real *, integer *, integer *, 
00122             integer *, real *, integer *);
00123     real rtunfl, rtovfl, rtulpi, ulpinv;
00124     integer mtypes, ntestt;
00125     extern /* Subroutine */ int sormhr_(char *, char *, integer *, integer *, 
00126             integer *, integer *, real *, integer *, real *, real *, integer *
00127 , real *, integer *, integer *);
00128 
00129     /* Fortran I/O blocks */
00130     static cilist io___36 = { 0, 0, 0, fmt_9999, 0 };
00131     static cilist io___39 = { 0, 0, 0, fmt_9999, 0 };
00132     static cilist io___41 = { 0, 0, 0, fmt_9999, 0 };
00133     static cilist io___42 = { 0, 0, 0, fmt_9999, 0 };
00134     static cilist io___43 = { 0, 0, 0, fmt_9999, 0 };
00135     static cilist io___50 = { 0, 0, 0, fmt_9999, 0 };
00136     static cilist io___51 = { 0, 0, 0, fmt_9998, 0 };
00137     static cilist io___52 = { 0, 0, 0, fmt_9999, 0 };
00138     static cilist io___56 = { 0, 0, 0, fmt_9997, 0 };
00139     static cilist io___57 = { 0, 0, 0, fmt_9999, 0 };
00140     static cilist io___58 = { 0, 0, 0, fmt_9998, 0 };
00141     static cilist io___59 = { 0, 0, 0, fmt_9999, 0 };
00142     static cilist io___60 = { 0, 0, 0, fmt_9997, 0 };
00143     static cilist io___61 = { 0, 0, 0, fmt_9999, 0 };
00144     static cilist io___62 = { 0, 0, 0, fmt_9998, 0 };
00145     static cilist io___63 = { 0, 0, 0, fmt_9999, 0 };
00146     static cilist io___64 = { 0, 0, 0, fmt_9998, 0 };
00147     static cilist io___65 = { 0, 0, 0, fmt_9999, 0 };
00148     static cilist io___66 = { 0, 0, 0, fmt_9999, 0 };
00149 
00150 
00151 
00152 /*  -- LAPACK test routine (version 3.1.1) -- */
00153 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00154 /*     February 2007 */
00155 
00156 /*     .. Scalar Arguments .. */
00157 /*     .. */
00158 /*     .. Array Arguments .. */
00159 /*     .. */
00160 
00161 /*  Purpose */
00162 /*  ======= */
00163 
00164 /*     SCHKHS  checks the nonsymmetric eigenvalue problem routines. */
00165 
00166 /*             SGEHRD factors A as  U H U' , where ' means transpose, */
00167 /*             H is hessenberg, and U is an orthogonal matrix. */
00168 
00169 /*             SORGHR generates the orthogonal matrix U. */
00170 
00171 /*             SORMHR multiplies a matrix by the orthogonal matrix U. */
00172 
00173 /*             SHSEQR factors H as  Z T Z' , where Z is orthogonal and */
00174 /*             T is "quasi-triangular", and the eigenvalue vector W. */
00175 
00176 /*             STREVC computes the left and right eigenvector matrices */
00177 /*             L and R for T. */
00178 
00179 /*             SHSEIN computes the left and right eigenvector matrices */
00180 /*             Y and X for H, using inverse iteration. */
00181 
00182 /*     When SCHKHS is called, a number of matrix "sizes" ("n's") and a */
00183 /*     number of matrix "types" are specified.  For each size ("n") */
00184 /*     and each type of matrix, one matrix will be generated and used */
00185 /*     to test the nonsymmetric eigenroutines.  For each matrix, 14 */
00186 /*     tests will be performed: */
00187 
00188 /*     (1)     | A - U H U**T | / ( |A| n ulp ) */
00189 
00190 /*     (2)     | I - UU**T | / ( n ulp ) */
00191 
00192 /*     (3)     | H - Z T Z**T | / ( |H| n ulp ) */
00193 
00194 /*     (4)     | I - ZZ**T | / ( n ulp ) */
00195 
00196 /*     (5)     | A - UZ H (UZ)**T | / ( |A| n ulp ) */
00197 
00198 /*     (6)     | I - UZ (UZ)**T | / ( n ulp ) */
00199 
00200 /*     (7)     | T(Z computed) - T(Z not computed) | / ( |T| ulp ) */
00201 
00202 /*     (8)     | W(Z computed) - W(Z not computed) | / ( |W| ulp ) */
00203 
00204 /*     (9)     | TR - RW | / ( |T| |R| ulp ) */
00205 
00206 /*     (10)    | L**H T - W**H L | / ( |T| |L| ulp ) */
00207 
00208 /*     (11)    | HX - XW | / ( |H| |X| ulp ) */
00209 
00210 /*     (12)    | Y**H H - W**H Y | / ( |H| |Y| ulp ) */
00211 
00212 /*     (13)    | AX - XW | / ( |A| |X| ulp ) */
00213 
00214 /*     (14)    | Y**H A - W**H Y | / ( |A| |Y| ulp ) */
00215 
00216 /*     The "sizes" are specified by an array NN(1:NSIZES); the value of */
00217 /*     each element NN(j) specifies one size. */
00218 /*     The "types" are specified by a logical array DOTYPE( 1:NTYPES ); */
00219 /*     if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. */
00220 /*     Currently, the list of possible types is: */
00221 
00222 /*     (1)  The zero matrix. */
00223 /*     (2)  The identity matrix. */
00224 /*     (3)  A (transposed) Jordan block, with 1's on the diagonal. */
00225 
00226 /*     (4)  A diagonal matrix with evenly spaced entries */
00227 /*          1, ..., ULP  and random signs. */
00228 /*          (ULP = (first number larger than 1) - 1 ) */
00229 /*     (5)  A diagonal matrix with geometrically spaced entries */
00230 /*          1, ..., ULP  and random signs. */
00231 /*     (6)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP */
00232 /*          and random signs. */
00233 
00234 /*     (7)  Same as (4), but multiplied by SQRT( overflow threshold ) */
00235 /*     (8)  Same as (4), but multiplied by SQRT( underflow threshold ) */
00236 
00237 /*     (9)  A matrix of the form  U' T U, where U is orthogonal and */
00238 /*          T has evenly spaced entries 1, ..., ULP with random signs */
00239 /*          on the diagonal and random O(1) entries in the upper */
00240 /*          triangle. */
00241 
00242 /*     (10) A matrix of the form  U' T U, where U is orthogonal and */
00243 /*          T has geometrically spaced entries 1, ..., ULP with random */
00244 /*          signs on the diagonal and random O(1) entries in the upper */
00245 /*          triangle. */
00246 
00247 /*     (11) A matrix of the form  U' T U, where U is orthogonal and */
00248 /*          T has "clustered" entries 1, ULP,..., ULP with random */
00249 /*          signs on the diagonal and random O(1) entries in the upper */
00250 /*          triangle. */
00251 
00252 /*     (12) A matrix of the form  U' T U, where U is orthogonal and */
00253 /*          T has real or complex conjugate paired eigenvalues randomly */
00254 /*          chosen from ( ULP, 1 ) and random O(1) entries in the upper */
00255 /*          triangle. */
00256 
00257 /*     (13) A matrix of the form  X' T X, where X has condition */
00258 /*          SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP */
00259 /*          with random signs on the diagonal and random O(1) entries */
00260 /*          in the upper triangle. */
00261 
00262 /*     (14) A matrix of the form  X' T X, where X has condition */
00263 /*          SQRT( ULP ) and T has geometrically spaced entries */
00264 /*          1, ..., ULP with random signs on the diagonal and random */
00265 /*          O(1) entries in the upper triangle. */
00266 
00267 /*     (15) A matrix of the form  X' T X, where X has condition */
00268 /*          SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP */
00269 /*          with random signs on the diagonal and random O(1) entries */
00270 /*          in the upper triangle. */
00271 
00272 /*     (16) A matrix of the form  X' T X, where X has condition */
00273 /*          SQRT( ULP ) and T has real or complex conjugate paired */
00274 /*          eigenvalues randomly chosen from ( ULP, 1 ) and random */
00275 /*          O(1) entries in the upper triangle. */
00276 
00277 /*     (17) Same as (16), but multiplied by SQRT( overflow threshold ) */
00278 /*     (18) Same as (16), but multiplied by SQRT( underflow threshold ) */
00279 
00280 /*     (19) Nonsymmetric matrix with random entries chosen from (-1,1). */
00281 /*     (20) Same as (19), but multiplied by SQRT( overflow threshold ) */
00282 /*     (21) Same as (19), but multiplied by SQRT( underflow threshold ) */
00283 
00284 /*  Arguments */
00285 /*  ========== */
00286 
00287 /*  NSIZES - INTEGER */
00288 /*           The number of sizes of matrices to use.  If it is zero, */
00289 /*           SCHKHS does nothing.  It must be at least zero. */
00290 /*           Not modified. */
00291 
00292 /*  NN     - INTEGER array, dimension (NSIZES) */
00293 /*           An array containing the sizes to be used for the matrices. */
00294 /*           Zero values will be skipped.  The values must be at least */
00295 /*           zero. */
00296 /*           Not modified. */
00297 
00298 /*  NTYPES - INTEGER */
00299 /*           The number of elements in DOTYPE.   If it is zero, SCHKHS */
00300 /*           does nothing.  It must be at least zero.  If it is MAXTYP+1 */
00301 /*           and NSIZES is 1, then an additional type, MAXTYP+1 is */
00302 /*           defined, which is to use whatever matrix is in A.  This */
00303 /*           is only useful if DOTYPE(1:MAXTYP) is .FALSE. and */
00304 /*           DOTYPE(MAXTYP+1) is .TRUE. . */
00305 /*           Not modified. */
00306 
00307 /*  DOTYPE - LOGICAL array, dimension (NTYPES) */
00308 /*           If DOTYPE(j) is .TRUE., then for each size in NN a */
00309 /*           matrix of that size and of type j will be generated. */
00310 /*           If NTYPES is smaller than the maximum number of types */
00311 /*           defined (PARAMETER MAXTYP), then types NTYPES+1 through */
00312 /*           MAXTYP will not be generated.  If NTYPES is larger */
00313 /*           than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) */
00314 /*           will be ignored. */
00315 /*           Not modified. */
00316 
00317 /*  ISEED  - INTEGER array, dimension (4) */
00318 /*           On entry ISEED specifies the seed of the random number */
00319 /*           generator. The array elements should be between 0 and 4095; */
00320 /*           if not they will be reduced mod 4096.  Also, ISEED(4) must */
00321 /*           be odd.  The random number generator uses a linear */
00322 /*           congruential sequence limited to small integers, and so */
00323 /*           should produce machine independent random numbers. The */
00324 /*           values of ISEED are changed on exit, and can be used in the */
00325 /*           next call to SCHKHS to continue the same random number */
00326 /*           sequence. */
00327 /*           Modified. */
00328 
00329 /*  THRESH - REAL */
00330 /*           A test will count as "failed" if the "error", computed as */
00331 /*           described above, exceeds THRESH.  Note that the error */
00332 /*           is scaled to be O(1), so THRESH should be a reasonably */
00333 /*           small multiple of 1, e.g., 10 or 100.  In particular, */
00334 /*           it should not depend on the precision (single vs. double) */
00335 /*           or the size of the matrix.  It must be at least zero. */
00336 /*           Not modified. */
00337 
00338 /*  NOUNIT - INTEGER */
00339 /*           The FORTRAN unit number for printing out error messages */
00340 /*           (e.g., if a routine returns IINFO not equal to 0.) */
00341 /*           Not modified. */
00342 
00343 /*  A      - REAL array, dimension (LDA,max(NN)) */
00344 /*           Used to hold the matrix whose eigenvalues are to be */
00345 /*           computed.  On exit, A contains the last matrix actually */
00346 /*           used. */
00347 /*           Modified. */
00348 
00349 /*  LDA    - INTEGER */
00350 /*           The leading dimension of A, H, T1 and T2.  It must be at */
00351 /*           least 1 and at least max( NN ). */
00352 /*           Not modified. */
00353 
00354 /*  H      - REAL array, dimension (LDA,max(NN)) */
00355 /*           The upper hessenberg matrix computed by SGEHRD.  On exit, */
00356 /*           H contains the Hessenberg form of the matrix in A. */
00357 /*           Modified. */
00358 
00359 /*  T1     - REAL array, dimension (LDA,max(NN)) */
00360 /*           The Schur (="quasi-triangular") matrix computed by SHSEQR */
00361 /*           if Z is computed.  On exit, T1 contains the Schur form of */
00362 /*           the matrix in A. */
00363 /*           Modified. */
00364 
00365 /*  T2     - REAL array, dimension (LDA,max(NN)) */
00366 /*           The Schur matrix computed by SHSEQR when Z is not computed. */
00367 /*           This should be identical to T1. */
00368 /*           Modified. */
00369 
00370 /*  LDU    - INTEGER */
00371 /*           The leading dimension of U, Z, UZ and UU.  It must be at */
00372 /*           least 1 and at least max( NN ). */
00373 /*           Not modified. */
00374 
00375 /*  U      - REAL array, dimension (LDU,max(NN)) */
00376 /*           The orthogonal matrix computed by SGEHRD. */
00377 /*           Modified. */
00378 
00379 /*  Z      - REAL array, dimension (LDU,max(NN)) */
00380 /*           The orthogonal matrix computed by SHSEQR. */
00381 /*           Modified. */
00382 
00383 /*  UZ     - REAL array, dimension (LDU,max(NN)) */
00384 /*           The product of U times Z. */
00385 /*           Modified. */
00386 
00387 /*  WR1    - REAL array, dimension (max(NN)) */
00388 /*  WI1    - REAL array, dimension (max(NN)) */
00389 /*           The real and imaginary parts of the eigenvalues of A, */
00390 /*           as computed when Z is computed. */
00391 /*           On exit, WR1 + WI1*i are the eigenvalues of the matrix in A. */
00392 /*           Modified. */
00393 
00394 /*  WR3    - REAL array, dimension (max(NN)) */
00395 /*  WI3    - REAL array, dimension (max(NN)) */
00396 /*           Like WR1, WI1, these arrays contain the eigenvalues of A, */
00397 /*           but those computed when SHSEQR only computes the */
00398 /*           eigenvalues, i.e., not the Schur vectors and no more of the */
00399 /*           Schur form than is necessary for computing the */
00400 /*           eigenvalues. */
00401 /*           Modified. */
00402 
00403 /*  EVECTL - REAL array, dimension (LDU,max(NN)) */
00404 /*           The (upper triangular) left eigenvector matrix for the */
00405 /*           matrix in T1.  For complex conjugate pairs, the real part */
00406 /*           is stored in one row and the imaginary part in the next. */
00407 /*           Modified. */
00408 
00409 /*  EVECTR - REAL array, dimension (LDU,max(NN)) */
00410 /*           The (upper triangular) right eigenvector matrix for the */
00411 /*           matrix in T1.  For complex conjugate pairs, the real part */
00412 /*           is stored in one column and the imaginary part in the next. */
00413 /*           Modified. */
00414 
00415 /*  EVECTY - REAL array, dimension (LDU,max(NN)) */
00416 /*           The left eigenvector matrix for the */
00417 /*           matrix in H.  For complex conjugate pairs, the real part */
00418 /*           is stored in one row and the imaginary part in the next. */
00419 /*           Modified. */
00420 
00421 /*  EVECTX - REAL array, dimension (LDU,max(NN)) */
00422 /*           The right eigenvector matrix for the */
00423 /*           matrix in H.  For complex conjugate pairs, the real part */
00424 /*           is stored in one column and the imaginary part in the next. */
00425 /*           Modified. */
00426 
00427 /*  UU     - REAL array, dimension (LDU,max(NN)) */
00428 /*           Details of the orthogonal matrix computed by SGEHRD. */
00429 /*           Modified. */
00430 
00431 /*  TAU    - REAL array, dimension(max(NN)) */
00432 /*           Further details of the orthogonal matrix computed by SGEHRD. */
00433 /*           Modified. */
00434 
00435 /*  WORK   - REAL array, dimension (NWORK) */
00436 /*           Workspace. */
00437 /*           Modified. */
00438 
00439 /*  NWORK  - INTEGER */
00440 /*           The number of entries in WORK.  NWORK >= 4*NN(j)*NN(j) + 2. */
00441 
00442 /*  IWORK  - INTEGER array, dimension (max(NN)) */
00443 /*           Workspace. */
00444 /*           Modified. */
00445 
00446 /*  SELECT - LOGICAL array, dimension (max(NN)) */
00447 /*           Workspace. */
00448 /*           Modified. */
00449 
00450 /*  RESULT - REAL array, dimension (14) */
00451 /*           The values computed by the fourteen tests described above. */
00452 /*           The values are currently limited to 1/ulp, to avoid */
00453 /*           overflow. */
00454 /*           Modified. */
00455 
00456 /*  INFO   - INTEGER */
00457 /*           If 0, then everything ran OK. */
00458 /*            -1: NSIZES < 0 */
00459 /*            -2: Some NN(j) < 0 */
00460 /*            -3: NTYPES < 0 */
00461 /*            -6: THRESH < 0 */
00462 /*            -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). */
00463 /*           -14: LDU < 1 or LDU < NMAX. */
00464 /*           -28: NWORK too small. */
00465 /*           If  SLATMR, SLATMS, or SLATME returns an error code, the */
00466 /*               absolute value of it is returned. */
00467 /*           If 1, then SHSEQR could not find all the shifts. */
00468 /*           If 2, then the EISPACK code (for small blocks) failed. */
00469 /*           If >2, then 30*N iterations were not enough to find an */
00470 /*               eigenvalue or to decompose the problem. */
00471 /*           Modified. */
00472 
00473 /* ----------------------------------------------------------------------- */
00474 
00475 /*     Some Local Variables and Parameters: */
00476 /*     ---- ----- --------- --- ---------- */
00477 
00478 /*     ZERO, ONE       Real 0 and 1. */
00479 /*     MAXTYP          The number of types defined. */
00480 /*     MTEST           The number of tests defined: care must be taken */
00481 /*                     that (1) the size of RESULT, (2) the number of */
00482 /*                     tests actually performed, and (3) MTEST agree. */
00483 /*     NTEST           The number of tests performed on this matrix */
00484 /*                     so far.  This should be less than MTEST, and */
00485 /*                     equal to it by the last test.  It will be less */
00486 /*                     if any of the routines being tested indicates */
00487 /*                     that it could not compute the matrices that */
00488 /*                     would be tested. */
00489 /*     NMAX            Largest value in NN. */
00490 /*     NMATS           The number of matrices generated so far. */
00491 /*     NERRS           The number of tests which have exceeded THRESH */
00492 /*                     so far (computed by SLAFTS). */
00493 /*     COND, CONDS, */
00494 /*     IMODE           Values to be passed to the matrix generators. */
00495 /*     ANORM           Norm of A; passed to matrix generators. */
00496 
00497 /*     OVFL, UNFL      Overflow and underflow thresholds. */
00498 /*     ULP, ULPINV     Finest relative precision and its inverse. */
00499 /*     RTOVFL, RTUNFL, */
00500 /*     RTULP, RTULPI   Square roots of the previous 4 values. */
00501 
00502 /*             The following four arrays decode JTYPE: */
00503 /*     KTYPE(j)        The general type (1-10) for type "j". */
00504 /*     KMODE(j)        The MODE value to be passed to the matrix */
00505 /*                     generator for type "j". */
00506 /*     KMAGN(j)        The order of magnitude ( O(1), */
00507 /*                     O(overflow^(1/2) ), O(underflow^(1/2) ) */
00508 /*     KCONDS(j)       Selects whether CONDS is to be 1 or */
00509 /*                     1/sqrt(ulp).  (0 means irrelevant.) */
00510 
00511 /*  ===================================================================== */
00512 
00513 /*     .. Parameters .. */
00514 /*     .. */
00515 /*     .. Local Scalars .. */
00516 /*     .. */
00517 /*     .. Local Arrays .. */
00518 /*     .. */
00519 /*     .. External Functions .. */
00520 /*     .. */
00521 /*     .. External Subroutines .. */
00522 /*     .. */
00523 /*     .. Intrinsic Functions .. */
00524 /*     .. */
00525 /*     .. Data statements .. */
00526     /* Parameter adjustments */
00527     --nn;
00528     --dotype;
00529     --iseed;
00530     t2_dim1 = *lda;
00531     t2_offset = 1 + t2_dim1;
00532     t2 -= t2_offset;
00533     t1_dim1 = *lda;
00534     t1_offset = 1 + t1_dim1;
00535     t1 -= t1_offset;
00536     h_dim1 = *lda;
00537     h_offset = 1 + h_dim1;
00538     h__ -= h_offset;
00539     a_dim1 = *lda;
00540     a_offset = 1 + a_dim1;
00541     a -= a_offset;
00542     uu_dim1 = *ldu;
00543     uu_offset = 1 + uu_dim1;
00544     uu -= uu_offset;
00545     evectx_dim1 = *ldu;
00546     evectx_offset = 1 + evectx_dim1;
00547     evectx -= evectx_offset;
00548     evecty_dim1 = *ldu;
00549     evecty_offset = 1 + evecty_dim1;
00550     evecty -= evecty_offset;
00551     evectr_dim1 = *ldu;
00552     evectr_offset = 1 + evectr_dim1;
00553     evectr -= evectr_offset;
00554     evectl_dim1 = *ldu;
00555     evectl_offset = 1 + evectl_dim1;
00556     evectl -= evectl_offset;
00557     uz_dim1 = *ldu;
00558     uz_offset = 1 + uz_dim1;
00559     uz -= uz_offset;
00560     z_dim1 = *ldu;
00561     z_offset = 1 + z_dim1;
00562     z__ -= z_offset;
00563     u_dim1 = *ldu;
00564     u_offset = 1 + u_dim1;
00565     u -= u_offset;
00566     --wr1;
00567     --wi1;
00568     --wr3;
00569     --wi3;
00570     --tau;
00571     --work;
00572     --iwork;
00573     --select;
00574     --result;
00575 
00576     /* Function Body */
00577 /*     .. */
00578 /*     .. Executable Statements .. */
00579 
00580 /*     Check for errors */
00581 
00582     ntestt = 0;
00583     *info = 0;
00584 
00585     badnn = FALSE_;
00586     nmax = 0;
00587     i__1 = *nsizes;
00588     for (j = 1; j <= i__1; ++j) {
00589 /* Computing MAX */
00590         i__2 = nmax, i__3 = nn[j];
00591         nmax = max(i__2,i__3);
00592         if (nn[j] < 0) {
00593             badnn = TRUE_;
00594         }
00595 /* L10: */
00596     }
00597 
00598 /*     Check for errors */
00599 
00600     if (*nsizes < 0) {
00601         *info = -1;
00602     } else if (badnn) {
00603         *info = -2;
00604     } else if (*ntypes < 0) {
00605         *info = -3;
00606     } else if (*thresh < 0.f) {
00607         *info = -6;
00608     } else if (*lda <= 1 || *lda < nmax) {
00609         *info = -9;
00610     } else if (*ldu <= 1 || *ldu < nmax) {
00611         *info = -14;
00612     } else if ((nmax << 2) * nmax + 2 > *nwork) {
00613         *info = -28;
00614     }
00615 
00616     if (*info != 0) {
00617         i__1 = -(*info);
00618         xerbla_("SCHKHS", &i__1);
00619         return 0;
00620     }
00621 
00622 /*     Quick return if possible */
00623 
00624     if (*nsizes == 0 || *ntypes == 0) {
00625         return 0;
00626     }
00627 
00628 /*     More important constants */
00629 
00630     unfl = slamch_("Safe minimum");
00631     ovfl = slamch_("Overflow");
00632     slabad_(&unfl, &ovfl);
00633     ulp = slamch_("Epsilon") * slamch_("Base");
00634     ulpinv = 1.f / ulp;
00635     rtunfl = sqrt(unfl);
00636     rtovfl = sqrt(ovfl);
00637     rtulp = sqrt(ulp);
00638     rtulpi = 1.f / rtulp;
00639 
00640 /*     Loop over sizes, types */
00641 
00642     nerrs = 0;
00643     nmats = 0;
00644 
00645     i__1 = *nsizes;
00646     for (jsize = 1; jsize <= i__1; ++jsize) {
00647         n = nn[jsize];
00648         if (n == 0) {
00649             goto L270;
00650         }
00651         n1 = max(1,n);
00652         aninv = 1.f / (real) n1;
00653 
00654         if (*nsizes != 1) {
00655             mtypes = min(21,*ntypes);
00656         } else {
00657             mtypes = min(22,*ntypes);
00658         }
00659 
00660         i__2 = mtypes;
00661         for (jtype = 1; jtype <= i__2; ++jtype) {
00662             if (! dotype[jtype]) {
00663                 goto L260;
00664             }
00665             ++nmats;
00666             ntest = 0;
00667 
00668 /*           Save ISEED in case of an error. */
00669 
00670             for (j = 1; j <= 4; ++j) {
00671                 ioldsd[j - 1] = iseed[j];
00672 /* L20: */
00673             }
00674 
00675 /*           Initialize RESULT */
00676 
00677             for (j = 1; j <= 14; ++j) {
00678                 result[j] = 0.f;
00679 /* L30: */
00680             }
00681 
00682 /*           Compute "A" */
00683 
00684 /*           Control parameters: */
00685 
00686 /*           KMAGN  KCONDS  KMODE        KTYPE */
00687 /*       =1  O(1)   1       clustered 1  zero */
00688 /*       =2  large  large   clustered 2  identity */
00689 /*       =3  small          exponential  Jordan */
00690 /*       =4                 arithmetic   diagonal, (w/ eigenvalues) */
00691 /*       =5                 random log   symmetric, w/ eigenvalues */
00692 /*       =6                 random       general, w/ eigenvalues */
00693 /*       =7                              random diagonal */
00694 /*       =8                              random symmetric */
00695 /*       =9                              random general */
00696 /*       =10                             random triangular */
00697 
00698             if (mtypes > 21) {
00699                 goto L100;
00700             }
00701 
00702             itype = ktype[jtype - 1];
00703             imode = kmode[jtype - 1];
00704 
00705 /*           Compute norm */
00706 
00707             switch (kmagn[jtype - 1]) {
00708                 case 1:  goto L40;
00709                 case 2:  goto L50;
00710                 case 3:  goto L60;
00711             }
00712 
00713 L40:
00714             anorm = 1.f;
00715             goto L70;
00716 
00717 L50:
00718             anorm = rtovfl * ulp * aninv;
00719             goto L70;
00720 
00721 L60:
00722             anorm = rtunfl * n * ulpinv;
00723             goto L70;
00724 
00725 L70:
00726 
00727             slaset_("Full", lda, &n, &c_b18, &c_b18, &a[a_offset], lda);
00728             iinfo = 0;
00729             cond = ulpinv;
00730 
00731 /*           Special Matrices */
00732 
00733             if (itype == 1) {
00734 
00735 /*              Zero */
00736 
00737                 iinfo = 0;
00738 
00739             } else if (itype == 2) {
00740 
00741 /*              Identity */
00742 
00743                 i__3 = n;
00744                 for (jcol = 1; jcol <= i__3; ++jcol) {
00745                     a[jcol + jcol * a_dim1] = anorm;
00746 /* L80: */
00747                 }
00748 
00749             } else if (itype == 3) {
00750 
00751 /*              Jordan Block */
00752 
00753                 i__3 = n;
00754                 for (jcol = 1; jcol <= i__3; ++jcol) {
00755                     a[jcol + jcol * a_dim1] = anorm;
00756                     if (jcol > 1) {
00757                         a[jcol + (jcol - 1) * a_dim1] = 1.f;
00758                     }
00759 /* L90: */
00760                 }
00761 
00762             } else if (itype == 4) {
00763 
00764 /*              Diagonal Matrix, [Eigen]values Specified */
00765 
00766                 slatms_(&n, &n, "S", &iseed[1], "S", &work[1], &imode, &cond, 
00767                         &anorm, &c__0, &c__0, "N", &a[a_offset], lda, &work[n 
00768                         + 1], &iinfo);
00769 
00770             } else if (itype == 5) {
00771 
00772 /*              Symmetric, eigenvalues specified */
00773 
00774                 slatms_(&n, &n, "S", &iseed[1], "S", &work[1], &imode, &cond, 
00775                         &anorm, &n, &n, "N", &a[a_offset], lda, &work[n + 1], 
00776                         &iinfo);
00777 
00778             } else if (itype == 6) {
00779 
00780 /*              General, eigenvalues specified */
00781 
00782                 if (kconds[jtype - 1] == 1) {
00783                     conds = 1.f;
00784                 } else if (kconds[jtype - 1] == 2) {
00785                     conds = rtulpi;
00786                 } else {
00787                     conds = 0.f;
00788                 }
00789 
00790                 *(unsigned char *)&adumma[0] = ' ';
00791                 slatme_(&n, "S", &iseed[1], &work[1], &imode, &cond, &c_b32, 
00792                         adumma, "T", "T", "T", &work[n + 1], &c__4, &conds, &
00793                         n, &n, &anorm, &a[a_offset], lda, &work[(n << 1) + 1], 
00794                          &iinfo);
00795 
00796             } else if (itype == 7) {
00797 
00798 /*              Diagonal, random eigenvalues */
00799 
00800                 slatmr_(&n, &n, "S", &iseed[1], "S", &work[1], &c__6, &c_b32, 
00801                         &c_b32, "T", "N", &work[n + 1], &c__1, &c_b32, &work[(
00802                         n << 1) + 1], &c__1, &c_b32, "N", idumma, &c__0, &
00803                         c__0, &c_b18, &anorm, "NO", &a[a_offset], lda, &iwork[
00804                         1], &iinfo);
00805 
00806             } else if (itype == 8) {
00807 
00808 /*              Symmetric, random eigenvalues */
00809 
00810                 slatmr_(&n, &n, "S", &iseed[1], "S", &work[1], &c__6, &c_b32, 
00811                         &c_b32, "T", "N", &work[n + 1], &c__1, &c_b32, &work[(
00812                         n << 1) + 1], &c__1, &c_b32, "N", idumma, &n, &n, &
00813                         c_b18, &anorm, "NO", &a[a_offset], lda, &iwork[1], &
00814                         iinfo);
00815 
00816             } else if (itype == 9) {
00817 
00818 /*              General, random eigenvalues */
00819 
00820                 slatmr_(&n, &n, "S", &iseed[1], "N", &work[1], &c__6, &c_b32, 
00821                         &c_b32, "T", "N", &work[n + 1], &c__1, &c_b32, &work[(
00822                         n << 1) + 1], &c__1, &c_b32, "N", idumma, &n, &n, &
00823                         c_b18, &anorm, "NO", &a[a_offset], lda, &iwork[1], &
00824                         iinfo);
00825 
00826             } else if (itype == 10) {
00827 
00828 /*              Triangular, random eigenvalues */
00829 
00830                 slatmr_(&n, &n, "S", &iseed[1], "N", &work[1], &c__6, &c_b32, 
00831                         &c_b32, "T", "N", &work[n + 1], &c__1, &c_b32, &work[(
00832                         n << 1) + 1], &c__1, &c_b32, "N", idumma, &n, &c__0, &
00833                         c_b18, &anorm, "NO", &a[a_offset], lda, &iwork[1], &
00834                         iinfo);
00835 
00836             } else {
00837 
00838                 iinfo = 1;
00839             }
00840 
00841             if (iinfo != 0) {
00842                 io___36.ciunit = *nounit;
00843                 s_wsfe(&io___36);
00844                 do_fio(&c__1, "Generator", (ftnlen)9);
00845                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00846                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00847                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00848                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00849                 e_wsfe();
00850                 *info = abs(iinfo);
00851                 return 0;
00852             }
00853 
00854 L100:
00855 
00856 /*           Call SGEHRD to compute H and U, do tests. */
00857 
00858             slacpy_(" ", &n, &n, &a[a_offset], lda, &h__[h_offset], lda);
00859 
00860             ntest = 1;
00861 
00862             ilo = 1;
00863             ihi = n;
00864 
00865             i__3 = *nwork - n;
00866             sgehrd_(&n, &ilo, &ihi, &h__[h_offset], lda, &work[1], &work[n + 
00867                     1], &i__3, &iinfo);
00868 
00869             if (iinfo != 0) {
00870                 result[1] = ulpinv;
00871                 io___39.ciunit = *nounit;
00872                 s_wsfe(&io___39);
00873                 do_fio(&c__1, "SGEHRD", (ftnlen)6);
00874                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00875                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00876                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00877                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00878                 e_wsfe();
00879                 *info = abs(iinfo);
00880                 goto L250;
00881             }
00882 
00883             i__3 = n - 1;
00884             for (j = 1; j <= i__3; ++j) {
00885                 uu[j + 1 + j * uu_dim1] = 0.f;
00886                 i__4 = n;
00887                 for (i__ = j + 2; i__ <= i__4; ++i__) {
00888                     u[i__ + j * u_dim1] = h__[i__ + j * h_dim1];
00889                     uu[i__ + j * uu_dim1] = h__[i__ + j * h_dim1];
00890                     h__[i__ + j * h_dim1] = 0.f;
00891 /* L110: */
00892                 }
00893 /* L120: */
00894             }
00895             i__3 = n - 1;
00896             scopy_(&i__3, &work[1], &c__1, &tau[1], &c__1);
00897             i__3 = *nwork - n;
00898             sorghr_(&n, &ilo, &ihi, &u[u_offset], ldu, &work[1], &work[n + 1], 
00899                      &i__3, &iinfo);
00900             ntest = 2;
00901 
00902             shst01_(&n, &ilo, &ihi, &a[a_offset], lda, &h__[h_offset], lda, &
00903                     u[u_offset], ldu, &work[1], nwork, &result[1]);
00904 
00905 /*           Call SHSEQR to compute T1, T2 and Z, do tests. */
00906 
00907 /*           Eigenvalues only (WR3,WI3) */
00908 
00909             slacpy_(" ", &n, &n, &h__[h_offset], lda, &t2[t2_offset], lda);
00910             ntest = 3;
00911             result[3] = ulpinv;
00912 
00913             shseqr_("E", "N", &n, &ilo, &ihi, &t2[t2_offset], lda, &wr3[1], &
00914                     wi3[1], &uz[uz_offset], ldu, &work[1], nwork, &iinfo);
00915             if (iinfo != 0) {
00916                 io___41.ciunit = *nounit;
00917                 s_wsfe(&io___41);
00918                 do_fio(&c__1, "SHSEQR(E)", (ftnlen)9);
00919                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00920                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00921                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00922                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00923                 e_wsfe();
00924                 if (iinfo <= n + 2) {
00925                     *info = abs(iinfo);
00926                     goto L250;
00927                 }
00928             }
00929 
00930 /*           Eigenvalues (WR1,WI1) and Full Schur Form (T2) */
00931 
00932             slacpy_(" ", &n, &n, &h__[h_offset], lda, &t2[t2_offset], lda);
00933 
00934             shseqr_("S", "N", &n, &ilo, &ihi, &t2[t2_offset], lda, &wr1[1], &
00935                     wi1[1], &uz[uz_offset], ldu, &work[1], nwork, &iinfo);
00936             if (iinfo != 0 && iinfo <= n + 2) {
00937                 io___42.ciunit = *nounit;
00938                 s_wsfe(&io___42);
00939                 do_fio(&c__1, "SHSEQR(S)", (ftnlen)9);
00940                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00941                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00942                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00943                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00944                 e_wsfe();
00945                 *info = abs(iinfo);
00946                 goto L250;
00947             }
00948 
00949 /*           Eigenvalues (WR1,WI1), Schur Form (T1), and Schur vectors */
00950 /*           (UZ) */
00951 
00952             slacpy_(" ", &n, &n, &h__[h_offset], lda, &t1[t1_offset], lda);
00953             slacpy_(" ", &n, &n, &u[u_offset], ldu, &uz[uz_offset], lda);
00954 
00955             shseqr_("S", "V", &n, &ilo, &ihi, &t1[t1_offset], lda, &wr1[1], &
00956                     wi1[1], &uz[uz_offset], ldu, &work[1], nwork, &iinfo);
00957             if (iinfo != 0 && iinfo <= n + 2) {
00958                 io___43.ciunit = *nounit;
00959                 s_wsfe(&io___43);
00960                 do_fio(&c__1, "SHSEQR(V)", (ftnlen)9);
00961                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00962                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00963                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00964                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00965                 e_wsfe();
00966                 *info = abs(iinfo);
00967                 goto L250;
00968             }
00969 
00970 /*           Compute Z = U' UZ */
00971 
00972             sgemm_("T", "N", &n, &n, &n, &c_b32, &u[u_offset], ldu, &uz[
00973                     uz_offset], ldu, &c_b18, &z__[z_offset], ldu);
00974             ntest = 8;
00975 
00976 /*           Do Tests 3: | H - Z T Z' | / ( |H| n ulp ) */
00977 /*                and 4: | I - Z Z' | / ( n ulp ) */
00978 
00979             shst01_(&n, &ilo, &ihi, &h__[h_offset], lda, &t1[t1_offset], lda, 
00980                     &z__[z_offset], ldu, &work[1], nwork, &result[3]);
00981 
00982 /*           Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp ) */
00983 /*                and 6: | I - UZ (UZ)' | / ( n ulp ) */
00984 
00985             shst01_(&n, &ilo, &ihi, &a[a_offset], lda, &t1[t1_offset], lda, &
00986                     uz[uz_offset], ldu, &work[1], nwork, &result[5]);
00987 
00988 /*           Do Test 7: | T2 - T1 | / ( |T| n ulp ) */
00989 
00990             sget10_(&n, &n, &t2[t2_offset], lda, &t1[t1_offset], lda, &work[1]
00991 , &result[7]);
00992 
00993 /*           Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp ) */
00994 
00995             temp1 = 0.f;
00996             temp2 = 0.f;
00997             i__3 = n;
00998             for (j = 1; j <= i__3; ++j) {
00999 /* Computing MAX */
01000                 r__5 = temp1, r__6 = (r__1 = wr1[j], dabs(r__1)) + (r__2 = 
01001                         wi1[j], dabs(r__2)), r__5 = max(r__5,r__6), r__6 = (
01002                         r__3 = wr3[j], dabs(r__3)) + (r__4 = wi3[j], dabs(
01003                         r__4));
01004                 temp1 = dmax(r__5,r__6);
01005 /* Computing MAX */
01006                 r__3 = temp2, r__4 = (r__1 = wr1[j] - wr3[j], dabs(r__1)) + (
01007                         r__2 = wr1[j] - wr3[j], dabs(r__2));
01008                 temp2 = dmax(r__3,r__4);
01009 /* L130: */
01010             }
01011 
01012 /* Computing MAX */
01013             r__1 = unfl, r__2 = ulp * dmax(temp1,temp2);
01014             result[8] = temp2 / dmax(r__1,r__2);
01015 
01016 /*           Compute the Left and Right Eigenvectors of T */
01017 
01018 /*           Compute the Right eigenvector Matrix: */
01019 
01020             ntest = 9;
01021             result[9] = ulpinv;
01022 
01023 /*           Select last max(N/4,1) real, max(N/4,1) complex eigenvectors */
01024 
01025             nselc = 0;
01026             nselr = 0;
01027             j = n;
01028 L140:
01029             if (wi1[j] == 0.f) {
01030 /* Computing MAX */
01031                 i__3 = n / 4;
01032                 if (nselr < max(i__3,1)) {
01033                     ++nselr;
01034                     select[j] = TRUE_;
01035                 } else {
01036                     select[j] = FALSE_;
01037                 }
01038                 --j;
01039             } else {
01040 /* Computing MAX */
01041                 i__3 = n / 4;
01042                 if (nselc < max(i__3,1)) {
01043                     ++nselc;
01044                     select[j] = TRUE_;
01045                     select[j - 1] = FALSE_;
01046                 } else {
01047                     select[j] = FALSE_;
01048                     select[j - 1] = FALSE_;
01049                 }
01050                 j += -2;
01051             }
01052             if (j > 0) {
01053                 goto L140;
01054             }
01055 
01056             strevc_("Right", "All", &select[1], &n, &t1[t1_offset], lda, 
01057                     dumma, ldu, &evectr[evectr_offset], ldu, &n, &in, &work[1]
01058 , &iinfo);
01059             if (iinfo != 0) {
01060                 io___50.ciunit = *nounit;
01061                 s_wsfe(&io___50);
01062                 do_fio(&c__1, "STREVC(R,A)", (ftnlen)11);
01063                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01064                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01065                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01066                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01067                 e_wsfe();
01068                 *info = abs(iinfo);
01069                 goto L250;
01070             }
01071 
01072 /*           Test 9:  | TR - RW | / ( |T| |R| ulp ) */
01073 
01074             sget22_("N", "N", "N", &n, &t1[t1_offset], lda, &evectr[
01075                     evectr_offset], ldu, &wr1[1], &wi1[1], &work[1], dumma);
01076             result[9] = dumma[0];
01077             if (dumma[1] > *thresh) {
01078                 io___51.ciunit = *nounit;
01079                 s_wsfe(&io___51);
01080                 do_fio(&c__1, "Right", (ftnlen)5);
01081                 do_fio(&c__1, "STREVC", (ftnlen)6);
01082                 do_fio(&c__1, (char *)&dumma[1], (ftnlen)sizeof(real));
01083                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01084                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01085                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01086                 e_wsfe();
01087             }
01088 
01089 /*           Compute selected right eigenvectors and confirm that */
01090 /*           they agree with previous right eigenvectors */
01091 
01092             strevc_("Right", "Some", &select[1], &n, &t1[t1_offset], lda, 
01093                     dumma, ldu, &evectl[evectl_offset], ldu, &n, &in, &work[1]
01094 , &iinfo);
01095             if (iinfo != 0) {
01096                 io___52.ciunit = *nounit;
01097                 s_wsfe(&io___52);
01098                 do_fio(&c__1, "STREVC(R,S)", (ftnlen)11);
01099                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01100                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01101                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01102                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01103                 e_wsfe();
01104                 *info = abs(iinfo);
01105                 goto L250;
01106             }
01107 
01108             k = 1;
01109             match = TRUE_;
01110             i__3 = n;
01111             for (j = 1; j <= i__3; ++j) {
01112                 if (select[j] && wi1[j] == 0.f) {
01113                     i__4 = n;
01114                     for (jj = 1; jj <= i__4; ++jj) {
01115                         if (evectr[jj + j * evectr_dim1] != evectl[jj + k * 
01116                                 evectl_dim1]) {
01117                             match = FALSE_;
01118                             goto L180;
01119                         }
01120 /* L150: */
01121                     }
01122                     ++k;
01123                 } else if (select[j] && wi1[j] != 0.f) {
01124                     i__4 = n;
01125                     for (jj = 1; jj <= i__4; ++jj) {
01126                         if (evectr[jj + j * evectr_dim1] != evectl[jj + k * 
01127                                 evectl_dim1] || evectr[jj + (j + 1) * 
01128                                 evectr_dim1] != evectl[jj + (k + 1) * 
01129                                 evectl_dim1]) {
01130                             match = FALSE_;
01131                             goto L180;
01132                         }
01133 /* L160: */
01134                     }
01135                     k += 2;
01136                 }
01137 /* L170: */
01138             }
01139 L180:
01140             if (! match) {
01141                 io___56.ciunit = *nounit;
01142                 s_wsfe(&io___56);
01143                 do_fio(&c__1, "Right", (ftnlen)5);
01144                 do_fio(&c__1, "STREVC", (ftnlen)6);
01145                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01146                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01147                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01148                 e_wsfe();
01149             }
01150 
01151 /*           Compute the Left eigenvector Matrix: */
01152 
01153             ntest = 10;
01154             result[10] = ulpinv;
01155             strevc_("Left", "All", &select[1], &n, &t1[t1_offset], lda, &
01156                     evectl[evectl_offset], ldu, dumma, ldu, &n, &in, &work[1], 
01157                      &iinfo);
01158             if (iinfo != 0) {
01159                 io___57.ciunit = *nounit;
01160                 s_wsfe(&io___57);
01161                 do_fio(&c__1, "STREVC(L,A)", (ftnlen)11);
01162                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01163                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01164                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01165                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01166                 e_wsfe();
01167                 *info = abs(iinfo);
01168                 goto L250;
01169             }
01170 
01171 /*           Test 10:  | LT - WL | / ( |T| |L| ulp ) */
01172 
01173             sget22_("Trans", "N", "Conj", &n, &t1[t1_offset], lda, &evectl[
01174                     evectl_offset], ldu, &wr1[1], &wi1[1], &work[1], &dumma[2]
01175 );
01176             result[10] = dumma[2];
01177             if (dumma[3] > *thresh) {
01178                 io___58.ciunit = *nounit;
01179                 s_wsfe(&io___58);
01180                 do_fio(&c__1, "Left", (ftnlen)4);
01181                 do_fio(&c__1, "STREVC", (ftnlen)6);
01182                 do_fio(&c__1, (char *)&dumma[3], (ftnlen)sizeof(real));
01183                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01184                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01185                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01186                 e_wsfe();
01187             }
01188 
01189 /*           Compute selected left eigenvectors and confirm that */
01190 /*           they agree with previous left eigenvectors */
01191 
01192             strevc_("Left", "Some", &select[1], &n, &t1[t1_offset], lda, &
01193                     evectr[evectr_offset], ldu, dumma, ldu, &n, &in, &work[1], 
01194                      &iinfo);
01195             if (iinfo != 0) {
01196                 io___59.ciunit = *nounit;
01197                 s_wsfe(&io___59);
01198                 do_fio(&c__1, "STREVC(L,S)", (ftnlen)11);
01199                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01200                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01201                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01202                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01203                 e_wsfe();
01204                 *info = abs(iinfo);
01205                 goto L250;
01206             }
01207 
01208             k = 1;
01209             match = TRUE_;
01210             i__3 = n;
01211             for (j = 1; j <= i__3; ++j) {
01212                 if (select[j] && wi1[j] == 0.f) {
01213                     i__4 = n;
01214                     for (jj = 1; jj <= i__4; ++jj) {
01215                         if (evectl[jj + j * evectl_dim1] != evectr[jj + k * 
01216                                 evectr_dim1]) {
01217                             match = FALSE_;
01218                             goto L220;
01219                         }
01220 /* L190: */
01221                     }
01222                     ++k;
01223                 } else if (select[j] && wi1[j] != 0.f) {
01224                     i__4 = n;
01225                     for (jj = 1; jj <= i__4; ++jj) {
01226                         if (evectl[jj + j * evectl_dim1] != evectr[jj + k * 
01227                                 evectr_dim1] || evectl[jj + (j + 1) * 
01228                                 evectl_dim1] != evectr[jj + (k + 1) * 
01229                                 evectr_dim1]) {
01230                             match = FALSE_;
01231                             goto L220;
01232                         }
01233 /* L200: */
01234                     }
01235                     k += 2;
01236                 }
01237 /* L210: */
01238             }
01239 L220:
01240             if (! match) {
01241                 io___60.ciunit = *nounit;
01242                 s_wsfe(&io___60);
01243                 do_fio(&c__1, "Left", (ftnlen)4);
01244                 do_fio(&c__1, "STREVC", (ftnlen)6);
01245                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01246                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01247                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01248                 e_wsfe();
01249             }
01250 
01251 /*           Call SHSEIN for Right eigenvectors of H, do test 11 */
01252 
01253             ntest = 11;
01254             result[11] = ulpinv;
01255             i__3 = n;
01256             for (j = 1; j <= i__3; ++j) {
01257                 select[j] = TRUE_;
01258 /* L230: */
01259             }
01260 
01261             shsein_("Right", "Qr", "Ninitv", &select[1], &n, &h__[h_offset], 
01262                     lda, &wr3[1], &wi3[1], dumma, ldu, &evectx[evectx_offset], 
01263                      ldu, &n1, &in, &work[1], &iwork[1], &iwork[1], &iinfo);
01264             if (iinfo != 0) {
01265                 io___61.ciunit = *nounit;
01266                 s_wsfe(&io___61);
01267                 do_fio(&c__1, "SHSEIN(R)", (ftnlen)9);
01268                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01269                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01270                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01271                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01272                 e_wsfe();
01273                 *info = abs(iinfo);
01274                 if (iinfo < 0) {
01275                     goto L250;
01276                 }
01277             } else {
01278 
01279 /*              Test 11:  | HX - XW | / ( |H| |X| ulp ) */
01280 
01281 /*                        (from inverse iteration) */
01282 
01283                 sget22_("N", "N", "N", &n, &h__[h_offset], lda, &evectx[
01284                         evectx_offset], ldu, &wr3[1], &wi3[1], &work[1], 
01285                         dumma);
01286                 if (dumma[0] < ulpinv) {
01287                     result[11] = dumma[0] * aninv;
01288                 }
01289                 if (dumma[1] > *thresh) {
01290                     io___62.ciunit = *nounit;
01291                     s_wsfe(&io___62);
01292                     do_fio(&c__1, "Right", (ftnlen)5);
01293                     do_fio(&c__1, "SHSEIN", (ftnlen)6);
01294                     do_fio(&c__1, (char *)&dumma[1], (ftnlen)sizeof(real));
01295                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01296                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01297                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
01298                             ;
01299                     e_wsfe();
01300                 }
01301             }
01302 
01303 /*           Call SHSEIN for Left eigenvectors of H, do test 12 */
01304 
01305             ntest = 12;
01306             result[12] = ulpinv;
01307             i__3 = n;
01308             for (j = 1; j <= i__3; ++j) {
01309                 select[j] = TRUE_;
01310 /* L240: */
01311             }
01312 
01313             shsein_("Left", "Qr", "Ninitv", &select[1], &n, &h__[h_offset], 
01314                     lda, &wr3[1], &wi3[1], &evecty[evecty_offset], ldu, dumma, 
01315                      ldu, &n1, &in, &work[1], &iwork[1], &iwork[1], &iinfo);
01316             if (iinfo != 0) {
01317                 io___63.ciunit = *nounit;
01318                 s_wsfe(&io___63);
01319                 do_fio(&c__1, "SHSEIN(L)", (ftnlen)9);
01320                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01321                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01322                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01323                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01324                 e_wsfe();
01325                 *info = abs(iinfo);
01326                 if (iinfo < 0) {
01327                     goto L250;
01328                 }
01329             } else {
01330 
01331 /*              Test 12:  | YH - WY | / ( |H| |Y| ulp ) */
01332 
01333 /*                        (from inverse iteration) */
01334 
01335                 sget22_("C", "N", "C", &n, &h__[h_offset], lda, &evecty[
01336                         evecty_offset], ldu, &wr3[1], &wi3[1], &work[1], &
01337                         dumma[2]);
01338                 if (dumma[2] < ulpinv) {
01339                     result[12] = dumma[2] * aninv;
01340                 }
01341                 if (dumma[3] > *thresh) {
01342                     io___64.ciunit = *nounit;
01343                     s_wsfe(&io___64);
01344                     do_fio(&c__1, "Left", (ftnlen)4);
01345                     do_fio(&c__1, "SHSEIN", (ftnlen)6);
01346                     do_fio(&c__1, (char *)&dumma[3], (ftnlen)sizeof(real));
01347                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01348                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01349                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
01350                             ;
01351                     e_wsfe();
01352                 }
01353             }
01354 
01355 /*           Call SORMHR for Right eigenvectors of A, do test 13 */
01356 
01357             ntest = 13;
01358             result[13] = ulpinv;
01359 
01360             sormhr_("Left", "No transpose", &n, &n, &ilo, &ihi, &uu[uu_offset]
01361 , ldu, &tau[1], &evectx[evectx_offset], ldu, &work[1], 
01362                     nwork, &iinfo);
01363             if (iinfo != 0) {
01364                 io___65.ciunit = *nounit;
01365                 s_wsfe(&io___65);
01366                 do_fio(&c__1, "SORMHR(R)", (ftnlen)9);
01367                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01368                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01369                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01370                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01371                 e_wsfe();
01372                 *info = abs(iinfo);
01373                 if (iinfo < 0) {
01374                     goto L250;
01375                 }
01376             } else {
01377 
01378 /*              Test 13:  | AX - XW | / ( |A| |X| ulp ) */
01379 
01380 /*                        (from inverse iteration) */
01381 
01382                 sget22_("N", "N", "N", &n, &a[a_offset], lda, &evectx[
01383                         evectx_offset], ldu, &wr3[1], &wi3[1], &work[1], 
01384                         dumma);
01385                 if (dumma[0] < ulpinv) {
01386                     result[13] = dumma[0] * aninv;
01387                 }
01388             }
01389 
01390 /*           Call SORMHR for Left eigenvectors of A, do test 14 */
01391 
01392             ntest = 14;
01393             result[14] = ulpinv;
01394 
01395             sormhr_("Left", "No transpose", &n, &n, &ilo, &ihi, &uu[uu_offset]
01396 , ldu, &tau[1], &evecty[evecty_offset], ldu, &work[1], 
01397                     nwork, &iinfo);
01398             if (iinfo != 0) {
01399                 io___66.ciunit = *nounit;
01400                 s_wsfe(&io___66);
01401                 do_fio(&c__1, "SORMHR(L)", (ftnlen)9);
01402                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01403                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01404                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01405                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01406                 e_wsfe();
01407                 *info = abs(iinfo);
01408                 if (iinfo < 0) {
01409                     goto L250;
01410                 }
01411             } else {
01412 
01413 /*              Test 14:  | YA - WY | / ( |A| |Y| ulp ) */
01414 
01415 /*                        (from inverse iteration) */
01416 
01417                 sget22_("C", "N", "C", &n, &a[a_offset], lda, &evecty[
01418                         evecty_offset], ldu, &wr3[1], &wi3[1], &work[1], &
01419                         dumma[2]);
01420                 if (dumma[2] < ulpinv) {
01421                     result[14] = dumma[2] * aninv;
01422                 }
01423             }
01424 
01425 /*           End of Loop -- Check for RESULT(j) > THRESH */
01426 
01427 L250:
01428 
01429             ntestt += ntest;
01430             slafts_("SHS", &n, &n, &jtype, &ntest, &result[1], ioldsd, thresh, 
01431                      nounit, &nerrs);
01432 
01433 L260:
01434             ;
01435         }
01436 L270:
01437         ;
01438     }
01439 
01440 /*     Summary */
01441 
01442     slasum_("SHS", nounit, &nerrs, &ntestt);
01443 
01444     return 0;
01445 
01446 
01447 /*     End of SCHKHS */
01448 
01449 } /* schkhs_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:59