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


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