zdrves.c
Go to the documentation of this file.
00001 /* zdrves.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 /* Common Block Declarations */
00017 
00018 struct {
00019     integer selopt, seldim;
00020     logical selval[20];
00021     doublereal selwr[20], selwi[20];
00022 } sslct_;
00023 
00024 #define sslct_1 sslct_
00025 
00026 /* Table of constant values */
00027 
00028 static doublecomplex c_b1 = {0.,0.};
00029 static doublecomplex c_b2 = {1.,0.};
00030 static integer c__0 = 0;
00031 static integer c__4 = 4;
00032 static integer c__6 = 6;
00033 static doublereal c_b38 = 1.;
00034 static integer c__1 = 1;
00035 static doublereal c_b48 = 0.;
00036 static integer c__2 = 2;
00037 
00038 /* Subroutine */ int zdrves_(integer *nsizes, integer *nn, integer *ntypes, 
00039         logical *dotype, integer *iseed, doublereal *thresh, integer *nounit, 
00040         doublecomplex *a, integer *lda, doublecomplex *h__, doublecomplex *ht, 
00041          doublecomplex *w, doublecomplex *wt, doublecomplex *vs, integer *
00042         ldvs, doublereal *result, doublecomplex *work, integer *nwork, 
00043         doublereal *rwork, integer *iwork, logical *bwork, integer *info)
00044 {
00045     /* Initialized data */
00046 
00047     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 };
00048     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 };
00049     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 };
00050     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 };
00051 
00052     /* Format strings */
00053     static char fmt_9992[] = "(\002 ZDRVES: \002,a,\002 returned INFO=\002,i"
00054             "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, ISEED="
00055             "(\002,3(i5,\002,\002),i5,\002)\002)";
00056     static char fmt_9999[] = "(/1x,a3,\002 -- Complex Schur Form Decompositi"
00057             "on Driver\002,/\002 Matrix types (see ZDRVES for details): \002)";
00058     static char fmt_9998[] = "(/\002 Special Matrices:\002,/\002  1=Zero mat"
00059             "rix.             \002,\002           \002,\002  5=Diagonal: geom"
00060             "etr. spaced entries.\002,/\002  2=Identity matrix.              "
00061             "      \002,\002  6=Diagona\002,\002l: clustered entries.\002,"
00062             "/\002  3=Transposed Jordan block.  \002,\002          \002,\002 "
00063             " 7=Diagonal: large, evenly spaced.\002,/\002  \002,\0024=Diagona"
00064             "l: evenly spaced entries.    \002,\002  8=Diagonal: s\002,\002ma"
00065             "ll, evenly spaced.\002)";
00066     static char fmt_9997[] = "(\002 Dense, Non-Symmetric Matrices:\002,/\002"
00067             "  9=Well-cond., ev\002,\002enly spaced eigenvals.\002,\002 14=Il"
00068             "l-cond., geomet. spaced e\002,\002igenals.\002,/\002 10=Well-con"
00069             "d., geom. spaced eigenvals. \002,\002 15=Ill-conditioned, cluste"
00070             "red e.vals.\002,/\002 11=Well-cond\002,\002itioned, clustered e."
00071             "vals. \002,\002 16=Ill-cond., random comp\002,\002lex \002,a6,"
00072             "/\002 12=Well-cond., random complex \002,a6,\002   \002,\002 17="
00073             "Ill-cond., large rand. complx \002,a4,/\002 13=Ill-condi\002,"
00074             "\002tioned, evenly spaced.     \002,\002 18=Ill-cond., small ran"
00075             "d.\002,\002 complx \002,a4)";
00076     static char fmt_9996[] = "(\002 19=Matrix with random O(1) entries.   "
00077             " \002,\002 21=Matrix \002,\002with small random entries.\002,"
00078             "/\002 20=Matrix with large ran\002,\002dom entries.   \002,/)";
00079     static char fmt_9995[] = "(\002 Tests performed with test threshold ="
00080             "\002,f8.2,/\002 ( A denotes A on input and T denotes A on output)"
00081             "\002,//\002 1 = 0 if T in Schur form (no sort), \002,\002  1/ulp"
00082             " otherwise\002,/\002 2 = | A - VS T transpose(VS) | / ( n |A| ul"
00083             "p ) (no sort)\002,/\002 3 = | I - VS transpose(VS) | / ( n ulp )"
00084             " (no sort) \002,/\002 4 = 0 if W are eigenvalues of T (no sort)"
00085             ",\002,\002  1/ulp otherwise\002,/\002 5 = 0 if T same no matter "
00086             "if VS computed (no sort),\002,\002  1/ulp otherwise\002,/\002 6 "
00087             "= 0 if W same no matter if VS computed (no sort)\002,\002,  1/ul"
00088             "p otherwise\002)";
00089     static char fmt_9994[] = "(\002 7 = 0 if T in Schur form (sort), \002"
00090             ",\002  1/ulp otherwise\002,/\002 8 = | A - VS T transpose(VS) | "
00091             "/ ( n |A| ulp ) (sort)\002,/\002 9 = | I - VS transpose(VS) | / "
00092             "( n ulp ) (sort) \002,/\002 10 = 0 if W are eigenvalues of T (so"
00093             "rt),\002,\002  1/ulp otherwise\002,/\002 11 = 0 if T same no mat"
00094             "ter if VS computed (sort),\002,\002  1/ulp otherwise\002,/\002 1"
00095             "2 = 0 if W same no matter if VS computed (sort),\002,\002  1/ulp"
00096             " otherwise\002,/\002 13 = 0 if sorting succesful, 1/ulp otherwise"
00097             "\002,/)";
00098     static char fmt_9993[] = "(\002 N=\002,i5,\002, IWK=\002,i2,\002, seed"
00099             "=\002,4(i4,\002,\002),\002 type \002,i2,\002, test(\002,i2,\002)="
00100             "\002,g10.3)";
00101 
00102     /* System generated locals */
00103     integer a_dim1, a_offset, h_dim1, h_offset, ht_dim1, ht_offset, vs_dim1, 
00104             vs_offset, i__1, i__2, i__3, i__4, i__5, i__6;
00105     doublecomplex z__1;
00106 
00107     /* Builtin functions */
00108     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00109     double sqrt(doublereal);
00110     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00111 
00112     /* Local variables */
00113     integer i__, j, n;
00114     doublereal res[2];
00115     integer iwk;
00116     doublereal ulp, cond;
00117     integer jcol;
00118     char path[3];
00119     integer sdim, nmax;
00120     doublereal unfl, ovfl;
00121     integer rsub;
00122     char sort[1];
00123     logical badnn;
00124     integer nfail, imode, iinfo;
00125     doublereal conds, anorm;
00126     extern /* Subroutine */ int zgees_(char *, char *, L_fp, integer *, 
00127             doublecomplex *, integer *, integer *, doublecomplex *, 
00128             doublecomplex *, integer *, doublecomplex *, integer *, 
00129             doublereal *, logical *, integer *);
00130     integer jsize, nerrs, itype, jtype, ntest, lwork, isort;
00131     extern /* Subroutine */ int zhst01_(integer *, integer *, integer *, 
00132             doublecomplex *, integer *, doublecomplex *, integer *, 
00133             doublecomplex *, integer *, doublecomplex *, integer *, 
00134             doublereal *, doublereal *);
00135     doublereal rtulp;
00136     extern /* Subroutine */ int dlabad_(doublereal *, doublereal *);
00137     extern doublereal dlamch_(char *);
00138     integer idumma[1], ioldsd[4];
00139     extern /* Subroutine */ int xerbla_(char *, integer *);
00140     integer knteig;
00141     extern /* Subroutine */ int dlasum_(char *, integer *, integer *, integer 
00142             *), zlatme_(integer *, char *, integer *, doublecomplex *, 
00143              integer *, doublereal *, doublecomplex *, char *, char *, char *, 
00144              char *, doublereal *, integer *, doublereal *, integer *, 
00145             integer *, doublereal *, doublecomplex *, integer *, 
00146             doublecomplex *, integer *), zlacpy_(char *, integer *, integer *, doublecomplex *, 
00147             integer *, doublecomplex *, integer *);
00148     integer ntestf;
00149     extern logical zslect_(doublecomplex *);
00150     extern /* Subroutine */ int zlaset_(char *, integer *, integer *, 
00151             doublecomplex *, doublecomplex *, doublecomplex *, integer *), zlatmr_(integer *, integer *, char *, integer *, char *, 
00152             doublecomplex *, integer *, doublereal *, doublecomplex *, char *, 
00153              char *, doublecomplex *, integer *, doublereal *, doublecomplex *
00154 , integer *, doublereal *, char *, integer *, integer *, integer *
00155 , doublereal *, doublereal *, char *, doublecomplex *, integer *, 
00156             integer *, integer *), zlatms_(integer *, integer *, char *, integer *, char *, 
00157             doublereal *, integer *, doublereal *, doublereal *, integer *, 
00158             integer *, char *, doublecomplex *, integer *, doublecomplex *, 
00159             integer *);
00160     integer nnwork;
00161     doublereal rtulpi;
00162     integer mtypes, ntestt;
00163     doublereal ulpinv;
00164 
00165     /* Fortran I/O blocks */
00166     static cilist io___31 = { 0, 0, 0, fmt_9992, 0 };
00167     static cilist io___38 = { 0, 0, 0, fmt_9992, 0 };
00168     static cilist io___42 = { 0, 0, 0, fmt_9992, 0 };
00169     static cilist io___46 = { 0, 0, 0, fmt_9999, 0 };
00170     static cilist io___47 = { 0, 0, 0, fmt_9998, 0 };
00171     static cilist io___48 = { 0, 0, 0, fmt_9997, 0 };
00172     static cilist io___49 = { 0, 0, 0, fmt_9996, 0 };
00173     static cilist io___50 = { 0, 0, 0, fmt_9995, 0 };
00174     static cilist io___51 = { 0, 0, 0, fmt_9994, 0 };
00175     static cilist io___52 = { 0, 0, 0, fmt_9993, 0 };
00176 
00177 
00178 
00179 /*  -- LAPACK test routine (version 3.1) -- */
00180 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00181 /*     November 2006 */
00182 
00183 /*     .. Scalar Arguments .. */
00184 /*     .. */
00185 /*     .. Array Arguments .. */
00186 /*     .. */
00187 
00188 /*  Purpose */
00189 /*  ======= */
00190 
00191 /*     ZDRVES checks the nonsymmetric eigenvalue (Schur form) problem */
00192 /*     driver ZGEES. */
00193 
00194 /*     When ZDRVES is called, a number of matrix "sizes" ("n's") and a */
00195 /*     number of matrix "types" are specified.  For each size ("n") */
00196 /*     and each type of matrix, one matrix will be generated and used */
00197 /*     to test the nonsymmetric eigenroutines.  For each matrix, 13 */
00198 /*     tests will be performed: */
00199 
00200 /*     (1)     0 if T is in Schur form, 1/ulp otherwise */
00201 /*            (no sorting of eigenvalues) */
00202 
00203 /*     (2)     | A - VS T VS' | / ( n |A| ulp ) */
00204 
00205 /*       Here VS is the matrix of Schur eigenvectors, and T is in Schur */
00206 /*       form  (no sorting of eigenvalues). */
00207 
00208 /*     (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues). */
00209 
00210 /*     (4)     0     if W are eigenvalues of T */
00211 /*             1/ulp otherwise */
00212 /*             (no sorting of eigenvalues) */
00213 
00214 /*     (5)     0     if T(with VS) = T(without VS), */
00215 /*             1/ulp otherwise */
00216 /*             (no sorting of eigenvalues) */
00217 
00218 /*     (6)     0     if eigenvalues(with VS) = eigenvalues(without VS), */
00219 /*             1/ulp otherwise */
00220 /*             (no sorting of eigenvalues) */
00221 
00222 /*     (7)     0 if T is in Schur form, 1/ulp otherwise */
00223 /*             (with sorting of eigenvalues) */
00224 
00225 /*     (8)     | A - VS T VS' | / ( n |A| ulp ) */
00226 
00227 /*       Here VS is the matrix of Schur eigenvectors, and T is in Schur */
00228 /*       form  (with sorting of eigenvalues). */
00229 
00230 /*     (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues). */
00231 
00232 /*     (10)    0     if W are eigenvalues of T */
00233 /*             1/ulp otherwise */
00234 /*             (with sorting of eigenvalues) */
00235 
00236 /*     (11)    0     if T(with VS) = T(without VS), */
00237 /*             1/ulp otherwise */
00238 /*             (with sorting of eigenvalues) */
00239 
00240 /*     (12)    0     if eigenvalues(with VS) = eigenvalues(without VS), */
00241 /*             1/ulp otherwise */
00242 /*             (with sorting of eigenvalues) */
00243 
00244 /*     (13)    if sorting worked and SDIM is the number of */
00245 /*             eigenvalues which were SELECTed */
00246 
00247 /*     The "sizes" are specified by an array NN(1:NSIZES); the value of */
00248 /*     each element NN(j) specifies one size. */
00249 /*     The "types" are specified by a logical array DOTYPE( 1:NTYPES ); */
00250 /*     if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. */
00251 /*     Currently, the list of possible types is: */
00252 
00253 /*     (1)  The zero matrix. */
00254 /*     (2)  The identity matrix. */
00255 /*     (3)  A (transposed) Jordan block, with 1's on the diagonal. */
00256 
00257 /*     (4)  A diagonal matrix with evenly spaced entries */
00258 /*          1, ..., ULP  and random complex angles. */
00259 /*          (ULP = (first number larger than 1) - 1 ) */
00260 /*     (5)  A diagonal matrix with geometrically spaced entries */
00261 /*          1, ..., ULP  and random complex angles. */
00262 /*     (6)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP */
00263 /*          and random complex angles. */
00264 
00265 /*     (7)  Same as (4), but multiplied by a constant near */
00266 /*          the overflow threshold */
00267 /*     (8)  Same as (4), but multiplied by a constant near */
00268 /*          the underflow threshold */
00269 
00270 /*     (9)  A matrix of the form  U' T U, where U is unitary and */
00271 /*          T has evenly spaced entries 1, ..., ULP with random */
00272 /*          complex angles on the diagonal and random O(1) entries in */
00273 /*          the upper triangle. */
00274 
00275 /*     (10) A matrix of the form  U' T U, where U is unitary and */
00276 /*          T has geometrically spaced entries 1, ..., ULP with random */
00277 /*          complex angles on the diagonal and random O(1) entries in */
00278 /*          the upper triangle. */
00279 
00280 /*     (11) A matrix of the form  U' T U, where U is orthogonal and */
00281 /*          T has "clustered" entries 1, ULP,..., ULP with random */
00282 /*          complex angles on the diagonal and random O(1) entries in */
00283 /*          the upper triangle. */
00284 
00285 /*     (12) A matrix of the form  U' T U, where U is unitary and */
00286 /*          T has complex eigenvalues randomly chosen from */
00287 /*          ULP < |z| < 1   and random O(1) entries in the upper */
00288 /*          triangle. */
00289 
00290 /*     (13) A matrix of the form  X' T X, where X has condition */
00291 /*          SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP */
00292 /*          with random complex angles on the diagonal and random O(1) */
00293 /*          entries in the upper triangle. */
00294 
00295 /*     (14) A matrix of the form  X' T X, where X has condition */
00296 /*          SQRT( ULP ) and T has geometrically spaced entries */
00297 /*          1, ..., ULP with random complex angles on the diagonal */
00298 /*          and random O(1) entries in the upper triangle. */
00299 
00300 /*     (15) A matrix of the form  X' T X, where X has condition */
00301 /*          SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP */
00302 /*          with random complex angles on the diagonal and random O(1) */
00303 /*          entries in the upper triangle. */
00304 
00305 /*     (16) A matrix of the form  X' T X, where X has condition */
00306 /*          SQRT( ULP ) and T has complex eigenvalues randomly chosen */
00307 /*          from ULP < |z| < 1 and random O(1) entries in the upper */
00308 /*          triangle. */
00309 
00310 /*     (17) Same as (16), but multiplied by a constant */
00311 /*          near the overflow threshold */
00312 /*     (18) Same as (16), but multiplied by a constant */
00313 /*          near the underflow threshold */
00314 
00315 /*     (19) Nonsymmetric matrix with random entries chosen from (-1,1). */
00316 /*          If N is at least 4, all entries in first two rows and last */
00317 /*          row, and first column and last two columns are zero. */
00318 /*     (20) Same as (19), but multiplied by a constant */
00319 /*          near the overflow threshold */
00320 /*     (21) Same as (19), but multiplied by a constant */
00321 /*          near the underflow threshold */
00322 
00323 /*  Arguments */
00324 /*  ========= */
00325 
00326 /*  NSIZES  (input) INTEGER */
00327 /*          The number of sizes of matrices to use.  If it is zero, */
00328 /*          ZDRVES does nothing.  It must be at least zero. */
00329 
00330 /*  NN      (input) INTEGER array, dimension (NSIZES) */
00331 /*          An array containing the sizes to be used for the matrices. */
00332 /*          Zero values will be skipped.  The values must be at least */
00333 /*          zero. */
00334 
00335 /*  NTYPES  (input) INTEGER */
00336 /*          The number of elements in DOTYPE.   If it is zero, ZDRVES */
00337 /*          does nothing.  It must be at least zero.  If it is MAXTYP+1 */
00338 /*          and NSIZES is 1, then an additional type, MAXTYP+1 is */
00339 /*          defined, which is to use whatever matrix is in A.  This */
00340 /*          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and */
00341 /*          DOTYPE(MAXTYP+1) is .TRUE. . */
00342 
00343 /*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) */
00344 /*          If DOTYPE(j) is .TRUE., then for each size in NN a */
00345 /*          matrix of that size and of type j will be generated. */
00346 /*          If NTYPES is smaller than the maximum number of types */
00347 /*          defined (PARAMETER MAXTYP), then types NTYPES+1 through */
00348 /*          MAXTYP will not be generated.  If NTYPES is larger */
00349 /*          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) */
00350 /*          will be ignored. */
00351 
00352 /*  ISEED   (input/output) INTEGER array, dimension (4) */
00353 /*          On entry ISEED specifies the seed of the random number */
00354 /*          generator. The array elements should be between 0 and 4095; */
00355 /*          if not they will be reduced mod 4096.  Also, ISEED(4) must */
00356 /*          be odd.  The random number generator uses a linear */
00357 /*          congruential sequence limited to small integers, and so */
00358 /*          should produce machine independent random numbers. The */
00359 /*          values of ISEED are changed on exit, and can be used in the */
00360 /*          next call to ZDRVES to continue the same random number */
00361 /*          sequence. */
00362 
00363 /*  THRESH  (input) DOUBLE PRECISION */
00364 /*          A test will count as "failed" if the "error", computed as */
00365 /*          described above, exceeds THRESH.  Note that the error */
00366 /*          is scaled to be O(1), so THRESH should be a reasonably */
00367 /*          small multiple of 1, e.g., 10 or 100.  In particular, */
00368 /*          it should not depend on the precision (single vs. double) */
00369 /*          or the size of the matrix.  It must be at least zero. */
00370 
00371 /*  NOUNIT  (input) INTEGER */
00372 /*          The FORTRAN unit number for printing out error messages */
00373 /*          (e.g., if a routine returns INFO not equal to 0.) */
00374 
00375 /*  A       (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) */
00376 /*          Used to hold the matrix whose eigenvalues are to be */
00377 /*          computed.  On exit, A contains the last matrix actually used. */
00378 
00379 /*  LDA     (input) INTEGER */
00380 /*          The leading dimension of A, and H. LDA must be at */
00381 /*          least 1 and at least max( NN ). */
00382 
00383 /*  H       (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) */
00384 /*          Another copy of the test matrix A, modified by ZGEES. */
00385 
00386 /*  HT      (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) */
00387 /*          Yet another copy of the test matrix A, modified by ZGEES. */
00388 
00389 /*  W       (workspace) COMPLEX*16 array, dimension (max(NN)) */
00390 /*          The computed eigenvalues of A. */
00391 
00392 /*  WT      (workspace) COMPLEX*16 array, dimension (max(NN)) */
00393 /*          Like W, this array contains the eigenvalues of A, */
00394 /*          but those computed when ZGEES only computes a partial */
00395 /*          eigendecomposition, i.e. not Schur vectors */
00396 
00397 /*  VS      (workspace) COMPLEX*16 array, dimension (LDVS, max(NN)) */
00398 /*          VS holds the computed Schur vectors. */
00399 
00400 /*  LDVS    (input) INTEGER */
00401 /*          Leading dimension of VS. Must be at least max(1,max(NN)). */
00402 
00403 /*  RESULT  (output) DOUBLE PRECISION array, dimension (13) */
00404 /*          The values computed by the 13 tests described above. */
00405 /*          The values are currently limited to 1/ulp, to avoid overflow. */
00406 
00407 /*  WORK    (workspace) COMPLEX*16 array, dimension (NWORK) */
00408 
00409 /*  NWORK   (input) INTEGER */
00410 /*          The number of entries in WORK.  This must be at least */
00411 /*          5*NN(j)+2*NN(j)**2 for all j. */
00412 
00413 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (max(NN)) */
00414 
00415 /*  IWORK   (workspace) INTEGER array, dimension (max(NN)) */
00416 
00417 /*  INFO    (output) INTEGER */
00418 /*          If 0, then everything ran OK. */
00419 /*           -1: NSIZES < 0 */
00420 /*           -2: Some NN(j) < 0 */
00421 /*           -3: NTYPES < 0 */
00422 /*           -6: THRESH < 0 */
00423 /*           -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). */
00424 /*          -15: LDVS < 1 or LDVS < NMAX, where NMAX is max( NN(j) ). */
00425 /*          -18: NWORK too small. */
00426 /*          If  ZLATMR, CLATMS, CLATME or ZGEES returns an error code, */
00427 /*              the absolute value of it is returned. */
00428 
00429 /* ----------------------------------------------------------------------- */
00430 
00431 /*     Some Local Variables and Parameters: */
00432 /*     ---- ----- --------- --- ---------- */
00433 /*     ZERO, ONE       Real 0 and 1. */
00434 /*     MAXTYP          The number of types defined. */
00435 /*     NMAX            Largest value in NN. */
00436 /*     NERRS           The number of tests which have exceeded THRESH */
00437 /*     COND, CONDS, */
00438 /*     IMODE           Values to be passed to the matrix generators. */
00439 /*     ANORM           Norm of A; passed to matrix generators. */
00440 
00441 /*     OVFL, UNFL      Overflow and underflow thresholds. */
00442 /*     ULP, ULPINV     Finest relative precision and its inverse. */
00443 /*     RTULP, RTULPI   Square roots of the previous 4 values. */
00444 /*             The following four arrays decode JTYPE: */
00445 /*     KTYPE(j)        The general type (1-10) for type "j". */
00446 /*     KMODE(j)        The MODE value to be passed to the matrix */
00447 /*                     generator for type "j". */
00448 /*     KMAGN(j)        The order of magnitude ( O(1), */
00449 /*                     O(overflow^(1/2) ), O(underflow^(1/2) ) */
00450 /*     KCONDS(j)       Select whether CONDS is to be 1 or */
00451 /*                     1/sqrt(ulp).  (0 means irrelevant.) */
00452 
00453 /*  ===================================================================== */
00454 
00455 /*     .. Parameters .. */
00456 /*     .. */
00457 /*     .. Local Scalars .. */
00458 /*     .. */
00459 /*     .. Local Arrays .. */
00460 /*     .. */
00461 /*     .. Arrays in Common .. */
00462 /*     .. */
00463 /*     .. Scalars in Common .. */
00464 /*     .. */
00465 /*     .. Common blocks .. */
00466 /*     .. */
00467 /*     .. External Functions .. */
00468 /*     .. */
00469 /*     .. External Subroutines .. */
00470 /*     .. */
00471 /*     .. Intrinsic Functions .. */
00472 /*     .. */
00473 /*     .. Data statements .. */
00474     /* Parameter adjustments */
00475     --nn;
00476     --dotype;
00477     --iseed;
00478     ht_dim1 = *lda;
00479     ht_offset = 1 + ht_dim1;
00480     ht -= ht_offset;
00481     h_dim1 = *lda;
00482     h_offset = 1 + h_dim1;
00483     h__ -= h_offset;
00484     a_dim1 = *lda;
00485     a_offset = 1 + a_dim1;
00486     a -= a_offset;
00487     --w;
00488     --wt;
00489     vs_dim1 = *ldvs;
00490     vs_offset = 1 + vs_dim1;
00491     vs -= vs_offset;
00492     --result;
00493     --work;
00494     --rwork;
00495     --iwork;
00496     --bwork;
00497 
00498     /* Function Body */
00499 /*     .. */
00500 /*     .. Executable Statements .. */
00501 
00502     s_copy(path, "Zomplex precision", (ftnlen)1, (ftnlen)17);
00503     s_copy(path + 1, "ES", (ftnlen)2, (ftnlen)2);
00504 
00505 /*     Check for errors */
00506 
00507     ntestt = 0;
00508     ntestf = 0;
00509     *info = 0;
00510     sslct_1.selopt = 0;
00511 
00512 /*     Important constants */
00513 
00514     badnn = FALSE_;
00515     nmax = 0;
00516     i__1 = *nsizes;
00517     for (j = 1; j <= i__1; ++j) {
00518 /* Computing MAX */
00519         i__2 = nmax, i__3 = nn[j];
00520         nmax = max(i__2,i__3);
00521         if (nn[j] < 0) {
00522             badnn = TRUE_;
00523         }
00524 /* L10: */
00525     }
00526 
00527 /*     Check for errors */
00528 
00529     if (*nsizes < 0) {
00530         *info = -1;
00531     } else if (badnn) {
00532         *info = -2;
00533     } else if (*ntypes < 0) {
00534         *info = -3;
00535     } else if (*thresh < 0.) {
00536         *info = -6;
00537     } else if (*nounit <= 0) {
00538         *info = -7;
00539     } else if (*lda < 1 || *lda < nmax) {
00540         *info = -9;
00541     } else if (*ldvs < 1 || *ldvs < nmax) {
00542         *info = -15;
00543     } else /* if(complicated condition) */ {
00544 /* Computing 2nd power */
00545         i__1 = nmax;
00546         if (nmax * 5 + (i__1 * i__1 << 1) > *nwork) {
00547             *info = -18;
00548         }
00549     }
00550 
00551     if (*info != 0) {
00552         i__1 = -(*info);
00553         xerbla_("ZDRVES", &i__1);
00554         return 0;
00555     }
00556 
00557 /*     Quick return if nothing to do */
00558 
00559     if (*nsizes == 0 || *ntypes == 0) {
00560         return 0;
00561     }
00562 
00563 /*     More Important constants */
00564 
00565     unfl = dlamch_("Safe minimum");
00566     ovfl = 1. / unfl;
00567     dlabad_(&unfl, &ovfl);
00568     ulp = dlamch_("Precision");
00569     ulpinv = 1. / ulp;
00570     rtulp = sqrt(ulp);
00571     rtulpi = 1. / rtulp;
00572 
00573 /*     Loop over sizes, types */
00574 
00575     nerrs = 0;
00576 
00577     i__1 = *nsizes;
00578     for (jsize = 1; jsize <= i__1; ++jsize) {
00579         n = nn[jsize];
00580         if (*nsizes != 1) {
00581             mtypes = min(21,*ntypes);
00582         } else {
00583             mtypes = min(22,*ntypes);
00584         }
00585 
00586         i__2 = mtypes;
00587         for (jtype = 1; jtype <= i__2; ++jtype) {
00588             if (! dotype[jtype]) {
00589                 goto L230;
00590             }
00591 
00592 /*           Save ISEED in case of an error. */
00593 
00594             for (j = 1; j <= 4; ++j) {
00595                 ioldsd[j - 1] = iseed[j];
00596 /* L20: */
00597             }
00598 
00599 /*           Compute "A" */
00600 
00601 /*           Control parameters: */
00602 
00603 /*           KMAGN  KCONDS  KMODE        KTYPE */
00604 /*       =1  O(1)   1       clustered 1  zero */
00605 /*       =2  large  large   clustered 2  identity */
00606 /*       =3  small          exponential  Jordan */
00607 /*       =4                 arithmetic   diagonal, (w/ eigenvalues) */
00608 /*       =5                 random log   symmetric, w/ eigenvalues */
00609 /*       =6                 random       general, w/ eigenvalues */
00610 /*       =7                              random diagonal */
00611 /*       =8                              random symmetric */
00612 /*       =9                              random general */
00613 /*       =10                             random triangular */
00614 
00615             if (mtypes > 21) {
00616                 goto L90;
00617             }
00618 
00619             itype = ktype[jtype - 1];
00620             imode = kmode[jtype - 1];
00621 
00622 /*           Compute norm */
00623 
00624             switch (kmagn[jtype - 1]) {
00625                 case 1:  goto L30;
00626                 case 2:  goto L40;
00627                 case 3:  goto L50;
00628             }
00629 
00630 L30:
00631             anorm = 1.;
00632             goto L60;
00633 
00634 L40:
00635             anorm = ovfl * ulp;
00636             goto L60;
00637 
00638 L50:
00639             anorm = unfl * ulpinv;
00640             goto L60;
00641 
00642 L60:
00643 
00644             zlaset_("Full", lda, &n, &c_b1, &c_b1, &a[a_offset], lda);
00645             iinfo = 0;
00646             cond = ulpinv;
00647 
00648 /*           Special Matrices -- Identity & Jordan block */
00649 
00650             if (itype == 1) {
00651 
00652 /*              Zero */
00653 
00654                 iinfo = 0;
00655 
00656             } else if (itype == 2) {
00657 
00658 /*              Identity */
00659 
00660                 i__3 = n;
00661                 for (jcol = 1; jcol <= i__3; ++jcol) {
00662                     i__4 = jcol + jcol * a_dim1;
00663                     z__1.r = anorm, z__1.i = 0.;
00664                     a[i__4].r = z__1.r, a[i__4].i = z__1.i;
00665 /* L70: */
00666                 }
00667 
00668             } else if (itype == 3) {
00669 
00670 /*              Jordan Block */
00671 
00672                 i__3 = n;
00673                 for (jcol = 1; jcol <= i__3; ++jcol) {
00674                     i__4 = jcol + jcol * a_dim1;
00675                     z__1.r = anorm, z__1.i = 0.;
00676                     a[i__4].r = z__1.r, a[i__4].i = z__1.i;
00677                     if (jcol > 1) {
00678                         i__4 = jcol + (jcol - 1) * a_dim1;
00679                         a[i__4].r = 1., a[i__4].i = 0.;
00680                     }
00681 /* L80: */
00682                 }
00683 
00684             } else if (itype == 4) {
00685 
00686 /*              Diagonal Matrix, [Eigen]values Specified */
00687 
00688                 zlatms_(&n, &n, "S", &iseed[1], "H", &rwork[1], &imode, &cond, 
00689                          &anorm, &c__0, &c__0, "N", &a[a_offset], lda, &work[
00690                         n + 1], &iinfo);
00691 
00692             } else if (itype == 5) {
00693 
00694 /*              Symmetric, eigenvalues specified */
00695 
00696                 zlatms_(&n, &n, "S", &iseed[1], "H", &rwork[1], &imode, &cond, 
00697                          &anorm, &n, &n, "N", &a[a_offset], lda, &work[n + 1], 
00698                          &iinfo);
00699 
00700             } else if (itype == 6) {
00701 
00702 /*              General, eigenvalues specified */
00703 
00704                 if (kconds[jtype - 1] == 1) {
00705                     conds = 1.;
00706                 } else if (kconds[jtype - 1] == 2) {
00707                     conds = rtulpi;
00708                 } else {
00709                     conds = 0.;
00710                 }
00711 
00712                 zlatme_(&n, "D", &iseed[1], &work[1], &imode, &cond, &c_b2, 
00713                         " ", "T", "T", "T", &rwork[1], &c__4, &conds, &n, &n, 
00714                         &anorm, &a[a_offset], lda, &work[(n << 1) + 1], &
00715                         iinfo);
00716 
00717             } else if (itype == 7) {
00718 
00719 /*              Diagonal, random eigenvalues */
00720 
00721                 zlatmr_(&n, &n, "D", &iseed[1], "N", &work[1], &c__6, &c_b38, 
00722                         &c_b2, "T", "N", &work[n + 1], &c__1, &c_b38, &work[(
00723                         n << 1) + 1], &c__1, &c_b38, "N", idumma, &c__0, &
00724                         c__0, &c_b48, &anorm, "NO", &a[a_offset], lda, &iwork[
00725                         1], &iinfo);
00726 
00727             } else if (itype == 8) {
00728 
00729 /*              Symmetric, random eigenvalues */
00730 
00731                 zlatmr_(&n, &n, "D", &iseed[1], "H", &work[1], &c__6, &c_b38, 
00732                         &c_b2, "T", "N", &work[n + 1], &c__1, &c_b38, &work[(
00733                         n << 1) + 1], &c__1, &c_b38, "N", idumma, &n, &n, &
00734                         c_b48, &anorm, "NO", &a[a_offset], lda, &iwork[1], &
00735                         iinfo);
00736 
00737             } else if (itype == 9) {
00738 
00739 /*              General, random eigenvalues */
00740 
00741                 zlatmr_(&n, &n, "D", &iseed[1], "N", &work[1], &c__6, &c_b38, 
00742                         &c_b2, "T", "N", &work[n + 1], &c__1, &c_b38, &work[(
00743                         n << 1) + 1], &c__1, &c_b38, "N", idumma, &n, &n, &
00744                         c_b48, &anorm, "NO", &a[a_offset], lda, &iwork[1], &
00745                         iinfo);
00746                 if (n >= 4) {
00747                     zlaset_("Full", &c__2, &n, &c_b1, &c_b1, &a[a_offset], 
00748                             lda);
00749                     i__3 = n - 3;
00750                     zlaset_("Full", &i__3, &c__1, &c_b1, &c_b1, &a[a_dim1 + 3]
00751 , lda);
00752                     i__3 = n - 3;
00753                     zlaset_("Full", &i__3, &c__2, &c_b1, &c_b1, &a[(n - 1) * 
00754                             a_dim1 + 3], lda);
00755                     zlaset_("Full", &c__1, &n, &c_b1, &c_b1, &a[n + a_dim1], 
00756                             lda);
00757                 }
00758 
00759             } else if (itype == 10) {
00760 
00761 /*              Triangular, random eigenvalues */
00762 
00763                 zlatmr_(&n, &n, "D", &iseed[1], "N", &work[1], &c__6, &c_b38, 
00764                         &c_b2, "T", "N", &work[n + 1], &c__1, &c_b38, &work[(
00765                         n << 1) + 1], &c__1, &c_b38, "N", idumma, &n, &c__0, &
00766                         c_b48, &anorm, "NO", &a[a_offset], lda, &iwork[1], &
00767                         iinfo);
00768 
00769             } else {
00770 
00771                 iinfo = 1;
00772             }
00773 
00774             if (iinfo != 0) {
00775                 io___31.ciunit = *nounit;
00776                 s_wsfe(&io___31);
00777                 do_fio(&c__1, "Generator", (ftnlen)9);
00778                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00779                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00780                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00781                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00782                 e_wsfe();
00783                 *info = abs(iinfo);
00784                 return 0;
00785             }
00786 
00787 L90:
00788 
00789 /*           Test for minimal and generous workspace */
00790 
00791             for (iwk = 1; iwk <= 2; ++iwk) {
00792                 if (iwk == 1) {
00793                     nnwork = n * 3;
00794                 } else {
00795 /* Computing 2nd power */
00796                     i__3 = n;
00797                     nnwork = n * 5 + (i__3 * i__3 << 1);
00798                 }
00799                 nnwork = max(nnwork,1);
00800 
00801 /*              Initialize RESULT */
00802 
00803                 for (j = 1; j <= 13; ++j) {
00804                     result[j] = -1.;
00805 /* L100: */
00806                 }
00807 
00808 /*              Test with and without sorting of eigenvalues */
00809 
00810                 for (isort = 0; isort <= 1; ++isort) {
00811                     if (isort == 0) {
00812                         *(unsigned char *)sort = 'N';
00813                         rsub = 0;
00814                     } else {
00815                         *(unsigned char *)sort = 'S';
00816                         rsub = 6;
00817                     }
00818 
00819 /*                 Compute Schur form and Schur vectors, and test them */
00820 
00821                     zlacpy_("F", &n, &n, &a[a_offset], lda, &h__[h_offset], 
00822                             lda);
00823                     zgees_("V", sort, (L_fp)zslect_, &n, &h__[h_offset], lda, 
00824                             &sdim, &w[1], &vs[vs_offset], ldvs, &work[1], &
00825                             nnwork, &rwork[1], &bwork[1], &iinfo);
00826                     if (iinfo != 0) {
00827                         result[rsub + 1] = ulpinv;
00828                         io___38.ciunit = *nounit;
00829                         s_wsfe(&io___38);
00830                         do_fio(&c__1, "ZGEES1", (ftnlen)6);
00831                         do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
00832                                 ;
00833                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00834                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
00835                                 ;
00836                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
00837                                 integer));
00838                         e_wsfe();
00839                         *info = abs(iinfo);
00840                         goto L190;
00841                     }
00842 
00843 /*                 Do Test (1) or Test (7) */
00844 
00845                     result[rsub + 1] = 0.;
00846                     i__3 = n - 1;
00847                     for (j = 1; j <= i__3; ++j) {
00848                         i__4 = n;
00849                         for (i__ = j + 1; i__ <= i__4; ++i__) {
00850                             i__5 = i__ + j * h_dim1;
00851                             if (h__[i__5].r != 0. || h__[i__5].i != 0.) {
00852                                 result[rsub + 1] = ulpinv;
00853                             }
00854 /* L110: */
00855                         }
00856 /* L120: */
00857                     }
00858 
00859 /*                 Do Tests (2) and (3) or Tests (8) and (9) */
00860 
00861 /* Computing MAX */
00862                     i__3 = 1, i__4 = (n << 1) * n;
00863                     lwork = max(i__3,i__4);
00864                     zhst01_(&n, &c__1, &n, &a[a_offset], lda, &h__[h_offset], 
00865                             lda, &vs[vs_offset], ldvs, &work[1], &lwork, &
00866                             rwork[1], res);
00867                     result[rsub + 2] = res[0];
00868                     result[rsub + 3] = res[1];
00869 
00870 /*                 Do Test (4) or Test (10) */
00871 
00872                     result[rsub + 4] = 0.;
00873                     i__3 = n;
00874                     for (i__ = 1; i__ <= i__3; ++i__) {
00875                         i__4 = i__ + i__ * h_dim1;
00876                         i__5 = i__;
00877                         if (h__[i__4].r != w[i__5].r || h__[i__4].i != w[i__5]
00878                                 .i) {
00879                             result[rsub + 4] = ulpinv;
00880                         }
00881 /* L130: */
00882                     }
00883 
00884 /*                 Do Test (5) or Test (11) */
00885 
00886                     zlacpy_("F", &n, &n, &a[a_offset], lda, &ht[ht_offset], 
00887                             lda);
00888                     zgees_("N", sort, (L_fp)zslect_, &n, &ht[ht_offset], lda, 
00889                             &sdim, &wt[1], &vs[vs_offset], ldvs, &work[1], &
00890                             nnwork, &rwork[1], &bwork[1], &iinfo);
00891                     if (iinfo != 0) {
00892                         result[rsub + 5] = ulpinv;
00893                         io___42.ciunit = *nounit;
00894                         s_wsfe(&io___42);
00895                         do_fio(&c__1, "ZGEES2", (ftnlen)6);
00896                         do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
00897                                 ;
00898                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00899                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
00900                                 ;
00901                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
00902                                 integer));
00903                         e_wsfe();
00904                         *info = abs(iinfo);
00905                         goto L190;
00906                     }
00907 
00908                     result[rsub + 5] = 0.;
00909                     i__3 = n;
00910                     for (j = 1; j <= i__3; ++j) {
00911                         i__4 = n;
00912                         for (i__ = 1; i__ <= i__4; ++i__) {
00913                             i__5 = i__ + j * h_dim1;
00914                             i__6 = i__ + j * ht_dim1;
00915                             if (h__[i__5].r != ht[i__6].r || h__[i__5].i != 
00916                                     ht[i__6].i) {
00917                                 result[rsub + 5] = ulpinv;
00918                             }
00919 /* L140: */
00920                         }
00921 /* L150: */
00922                     }
00923 
00924 /*                 Do Test (6) or Test (12) */
00925 
00926                     result[rsub + 6] = 0.;
00927                     i__3 = n;
00928                     for (i__ = 1; i__ <= i__3; ++i__) {
00929                         i__4 = i__;
00930                         i__5 = i__;
00931                         if (w[i__4].r != wt[i__5].r || w[i__4].i != wt[i__5]
00932                                 .i) {
00933                             result[rsub + 6] = ulpinv;
00934                         }
00935 /* L160: */
00936                     }
00937 
00938 /*                 Do Test (13) */
00939 
00940                     if (isort == 1) {
00941                         result[13] = 0.;
00942                         knteig = 0;
00943                         i__3 = n;
00944                         for (i__ = 1; i__ <= i__3; ++i__) {
00945                             if (zslect_(&w[i__])) {
00946                                 ++knteig;
00947                             }
00948                             if (i__ < n) {
00949                                 if (zslect_(&w[i__ + 1]) && ! zslect_(&w[i__])
00950                                         ) {
00951                                     result[13] = ulpinv;
00952                                 }
00953                             }
00954 /* L170: */
00955                         }
00956                         if (sdim != knteig) {
00957                             result[13] = ulpinv;
00958                         }
00959                     }
00960 
00961 /* L180: */
00962                 }
00963 
00964 /*              End of Loop -- Check for RESULT(j) > THRESH */
00965 
00966 L190:
00967 
00968                 ntest = 0;
00969                 nfail = 0;
00970                 for (j = 1; j <= 13; ++j) {
00971                     if (result[j] >= 0.) {
00972                         ++ntest;
00973                     }
00974                     if (result[j] >= *thresh) {
00975                         ++nfail;
00976                     }
00977 /* L200: */
00978                 }
00979 
00980                 if (nfail > 0) {
00981                     ++ntestf;
00982                 }
00983                 if (ntestf == 1) {
00984                     io___46.ciunit = *nounit;
00985                     s_wsfe(&io___46);
00986                     do_fio(&c__1, path, (ftnlen)3);
00987                     e_wsfe();
00988                     io___47.ciunit = *nounit;
00989                     s_wsfe(&io___47);
00990                     e_wsfe();
00991                     io___48.ciunit = *nounit;
00992                     s_wsfe(&io___48);
00993                     e_wsfe();
00994                     io___49.ciunit = *nounit;
00995                     s_wsfe(&io___49);
00996                     e_wsfe();
00997                     io___50.ciunit = *nounit;
00998                     s_wsfe(&io___50);
00999                     do_fio(&c__1, (char *)&(*thresh), (ftnlen)sizeof(
01000                             doublereal));
01001                     e_wsfe();
01002                     io___51.ciunit = *nounit;
01003                     s_wsfe(&io___51);
01004                     e_wsfe();
01005                     ntestf = 2;
01006                 }
01007 
01008                 for (j = 1; j <= 13; ++j) {
01009                     if (result[j] >= *thresh) {
01010                         io___52.ciunit = *nounit;
01011                         s_wsfe(&io___52);
01012                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01013                         do_fio(&c__1, (char *)&iwk, (ftnlen)sizeof(integer));
01014                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
01015                                 integer));
01016                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
01017                                 ;
01018                         do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01019                         do_fio(&c__1, (char *)&result[j], (ftnlen)sizeof(
01020                                 doublereal));
01021                         e_wsfe();
01022                     }
01023 /* L210: */
01024                 }
01025 
01026                 nerrs += nfail;
01027                 ntestt += ntest;
01028 
01029 /* L220: */
01030             }
01031 L230:
01032             ;
01033         }
01034 /* L240: */
01035     }
01036 
01037 /*     Summary */
01038 
01039     dlasum_(path, nounit, &nerrs, &ntestt);
01040 
01041 
01042 
01043     return 0;
01044 
01045 /*     End of ZDRVES */
01046 
01047 } /* zdrves_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:21