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


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