schkst.c
Go to the documentation of this file.
00001 /* schkst.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 static integer c_n1 = -1;
00020 static integer c__2 = 2;
00021 static real c_b25 = 0.f;
00022 static integer c__0 = 0;
00023 static integer c__6 = 6;
00024 static real c_b39 = 1.f;
00025 static integer c__4 = 4;
00026 static integer c__3 = 3;
00027 static integer c__10 = 10;
00028 static integer c__11 = 11;
00029 
00030 /* Subroutine */ int schkst_(integer *nsizes, integer *nn, integer *ntypes, 
00031         logical *dotype, integer *iseed, real *thresh, integer *nounit, real *
00032         a, integer *lda, real *ap, real *sd, real *se, real *d1, real *d2, 
00033         real *d3, real *d4, real *d5, real *wa1, real *wa2, real *wa3, real *
00034         wr, real *u, integer *ldu, real *v, real *vp, real *tau, real *z__, 
00035         real *work, integer *lwork, integer *iwork, integer *liwork, real *
00036         result, integer *info)
00037 {
00038     /* Initialized data */
00039 
00040     static integer ktype[21] = { 1,2,4,4,4,4,4,5,5,5,5,5,8,8,8,9,9,9,9,9,10 };
00041     static integer kmagn[21] = { 1,1,1,1,1,2,3,1,1,1,2,3,1,2,3,1,1,1,2,3,1 };
00042     static integer kmode[21] = { 0,0,4,3,1,4,4,4,3,1,4,4,0,0,0,4,3,1,4,4,3 };
00043 
00044     /* Format strings */
00045     static char fmt_9999[] = "(\002 SCHKST: \002,a,\002 returned INFO=\002,i"
00046             "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, ISEED="
00047             "(\002,3(i5,\002,\002),i5,\002)\002)";
00048     static char fmt_9998[] = "(/1x,a3,\002 -- Real Symmetric eigenvalue prob"
00049             "lem\002)";
00050     static char fmt_9997[] = "(\002 Matrix types (see SCHKST for details):"
00051             " \002)";
00052     static char fmt_9996[] = "(/\002 Special Matrices:\002,/\002  1=Zero mat"
00053             "rix.                        \002,\002  5=Diagonal: clustered ent"
00054             "ries.\002,/\002  2=Identity matrix.                    \002,\002"
00055             "  6=Diagonal: large, evenly spaced.\002,/\002  3=Diagonal: evenl"
00056             "y spaced entries.    \002,\002  7=Diagonal: small, evenly spaced."
00057             "\002,/\002  4=Diagonal: geometr. spaced entries.\002)";
00058     static char fmt_9995[] = "(\002 Dense \002,a,\002 Matrices:\002,/\002  8"
00059             "=Evenly spaced eigenvals.            \002,\002 12=Small, evenly "
00060             "spaced eigenvals.\002,/\002  9=Geometrically spaced eigenvals.  "
00061             "   \002,\002 13=Matrix with random O(1) entries.\002,/\002 10=Cl"
00062             "ustered eigenvalues.              \002,\002 14=Matrix with large"
00063             " random entries.\002,/\002 11=Large, evenly spaced eigenvals.   "
00064             "  \002,\002 15=Matrix with small random entries.\002)";
00065     static char fmt_9994[] = "(\002 16=Positive definite, evenly spaced eige"
00066             "nvalues\002,/\002 17=Positive definite, geometrically spaced eig"
00067             "envlaues\002,/\002 18=Positive definite, clustered eigenvalue"
00068             "s\002,/\002 19=Positive definite, small evenly spaced eigenvalues"
00069             "\002,/\002 20=Positive definite, large evenly spaced eigenvalue"
00070             "s\002,/\002 21=Diagonally dominant tridiagonal, geometrically"
00071             "\002,\002 spaced eigenvalues\002)";
00072     static char fmt_9988[] = "(/\002Test performed:  see SCHKST for details"
00073             ".\002,/)";
00074     static char fmt_9990[] = "(\002 N=\002,i5,\002, seed=\002,4(i4,\002,\002"
00075             "),\002 type \002,i2,\002, test(\002,i2,\002)=\002,g10.3)";
00076 
00077     /* System generated locals */
00078     integer a_dim1, a_offset, u_dim1, u_offset, v_dim1, v_offset, z_dim1, 
00079             z_offset, i__1, i__2, i__3, i__4;
00080     real r__1, r__2, r__3, r__4;
00081 
00082     /* Builtin functions */
00083     double log(doublereal), sqrt(doublereal);
00084     integer pow_ii(integer *, integer *), s_wsfe(cilist *), do_fio(integer *, 
00085             char *, ftnlen), e_wsfe(void);
00086 
00087     /* Local variables */
00088     integer i__, j, m, n, m2, m3, jc, il, jr, iu;
00089     real vl, vu;
00090     integer nap, lgn;
00091     real ulp, cond;
00092     integer nmax;
00093     real unfl, ovfl, temp1, temp2, temp3, temp4;
00094     logical badnn;
00095     extern doublereal ssxt1_(integer *, real *, integer *, real *, integer *, 
00096             real *, real *, real *);
00097     integer imode, lwedc;
00098     real dumma[1];
00099     integer iinfo;
00100     real aninv, anorm;
00101     integer itemp, nmats, jsize, nerrs, itype, jtype, ntest;
00102     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
00103             integer *), sspt21_(integer *, char *, integer *, integer *, real 
00104             *, real *, real *, real *, integer *, real *, real *, real *, 
00105             real *), sstt21_(integer *, integer *, real *, real *, 
00106             real *, real *, real *, integer *, real *, real *), sstt22_(
00107             integer *, integer *, integer *, real *, real *, real *, real *, 
00108             real *, integer *, real *, integer *, real *), ssyt21_(integer *, 
00109             char *, integer *, integer *, real *, integer *, real *, real *, 
00110             real *, integer *, real *, integer *, real *, real *, real *);
00111     integer iseed2[4], log2ui;
00112     extern /* Subroutine */ int slabad_(real *, real *);
00113     integer liwedc, nblock;
00114     extern doublereal slamch_(char *);
00115     integer idumma[1];
00116     extern /* Subroutine */ int xerbla_(char *, integer *);
00117     integer ioldsd[4];
00118     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00119             integer *, integer *);
00120     extern doublereal slarnd_(integer *, integer *);
00121     real abstol;
00122     extern /* Subroutine */ int sstedc_(char *, integer *, real *, real *, 
00123             real *, integer *, real *, integer *, integer *, integer *, 
00124             integer *), sstech_(integer *, real *, real *, real *, 
00125             real *, real *, integer *), slacpy_(char *, integer *, integer *, 
00126             real *, integer *, real *, integer *), slaset_(char *, 
00127             integer *, integer *, real *, real *, real *, integer *), 
00128             slatmr_(integer *, integer *, char *, integer *, char *, real *, 
00129             integer *, real *, real *, char *, char *, real *, integer *, 
00130             real *, real *, integer *, real *, char *, integer *, integer *, 
00131             integer *, real *, real *, char *, real *, integer *, integer *, 
00132             integer *), 
00133             slatms_(integer *, integer *, char *, integer *, char *, real *, 
00134             integer *, real *, real *, integer *, integer *, char *, real *, 
00135             integer *, real *, integer *);
00136     logical tryrac;
00137     extern /* Subroutine */ int slasum_(char *, integer *, integer *, integer 
00138             *), sstein_(integer *, real *, real *, integer *, real *, 
00139             integer *, integer *, real *, integer *, real *, integer *, 
00140             integer *, integer *);
00141     integer nsplit;
00142     real rtunfl, rtovfl, ulpinv;
00143     extern /* Subroutine */ int sopgtr_(char *, integer *, real *, real *, 
00144             real *, integer *, real *, integer *);
00145     integer mtypes, ntestt;
00146     extern /* Subroutine */ int sorgtr_(char *, integer *, real *, integer *, 
00147             real *, real *, integer *, integer *), spteqr_(char *, 
00148             integer *, real *, real *, real *, integer *, real *, integer *), ssptrd_(char *, integer *, real *, real *, real *, real *
00149 , integer *), sstebz_(char *, char *, integer *, real *, 
00150             real *, integer *, integer *, real *, real *, real *, integer *, 
00151             integer *, real *, integer *, integer *, real *, integer *, 
00152             integer *), sstemr_(char *, char *, integer *, 
00153             real *, real *, real *, real *, integer *, integer *, integer *, 
00154             real *, real *, integer *, integer *, integer *, logical *, real *
00155 , integer *, integer *, integer *, integer *), 
00156             ssteqr_(char *, integer *, real *, real *, real *, integer *, 
00157             real *, integer *), ssterf_(integer *, real *, real *, 
00158             integer *), ssytrd_(char *, integer *, real *, integer *, real *, 
00159             real *, real *, real *, integer *, integer *);
00160 
00161     /* Fortran I/O blocks */
00162     static cilist io___40 = { 0, 0, 0, fmt_9999, 0 };
00163     static cilist io___41 = { 0, 0, 0, fmt_9999, 0 };
00164     static cilist io___42 = { 0, 0, 0, fmt_9999, 0 };
00165     static cilist io___43 = { 0, 0, 0, fmt_9999, 0 };
00166     static cilist io___44 = { 0, 0, 0, fmt_9999, 0 };
00167     static cilist io___46 = { 0, 0, 0, fmt_9999, 0 };
00168     static cilist io___47 = { 0, 0, 0, fmt_9999, 0 };
00169     static cilist io___48 = { 0, 0, 0, fmt_9999, 0 };
00170     static cilist io___49 = { 0, 0, 0, fmt_9999, 0 };
00171     static cilist io___50 = { 0, 0, 0, fmt_9999, 0 };
00172     static cilist io___51 = { 0, 0, 0, fmt_9999, 0 };
00173     static cilist io___52 = { 0, 0, 0, fmt_9999, 0 };
00174     static cilist io___57 = { 0, 0, 0, fmt_9999, 0 };
00175     static cilist io___58 = { 0, 0, 0, fmt_9999, 0 };
00176     static cilist io___66 = { 0, 0, 0, fmt_9999, 0 };
00177     static cilist io___67 = { 0, 0, 0, fmt_9999, 0 };
00178     static cilist io___70 = { 0, 0, 0, fmt_9999, 0 };
00179     static cilist io___72 = { 0, 0, 0, fmt_9999, 0 };
00180     static cilist io___73 = { 0, 0, 0, fmt_9999, 0 };
00181     static cilist io___74 = { 0, 0, 0, fmt_9999, 0 };
00182     static cilist io___75 = { 0, 0, 0, fmt_9999, 0 };
00183     static cilist io___76 = { 0, 0, 0, fmt_9999, 0 };
00184     static cilist io___77 = { 0, 0, 0, fmt_9999, 0 };
00185     static cilist io___78 = { 0, 0, 0, fmt_9999, 0 };
00186     static cilist io___79 = { 0, 0, 0, fmt_9999, 0 };
00187     static cilist io___80 = { 0, 0, 0, fmt_9999, 0 };
00188     static cilist io___81 = { 0, 0, 0, fmt_9999, 0 };
00189     static cilist io___82 = { 0, 0, 0, fmt_9999, 0 };
00190     static cilist io___83 = { 0, 0, 0, fmt_9999, 0 };
00191     static cilist io___84 = { 0, 0, 0, fmt_9999, 0 };
00192     static cilist io___85 = { 0, 0, 0, fmt_9999, 0 };
00193     static cilist io___86 = { 0, 0, 0, fmt_9998, 0 };
00194     static cilist io___87 = { 0, 0, 0, fmt_9997, 0 };
00195     static cilist io___88 = { 0, 0, 0, fmt_9996, 0 };
00196     static cilist io___89 = { 0, 0, 0, fmt_9995, 0 };
00197     static cilist io___90 = { 0, 0, 0, fmt_9994, 0 };
00198     static cilist io___91 = { 0, 0, 0, fmt_9988, 0 };
00199     static cilist io___92 = { 0, 0, 0, fmt_9990, 0 };
00200 
00201 
00202 
00203 /*  -- LAPACK test routine (version 3.1) -- */
00204 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00205 /*     November 2006 */
00206 
00207 /*     .. Scalar Arguments .. */
00208 /*     .. */
00209 /*     .. Array Arguments .. */
00210 /*     .. */
00211 
00212 /*  Purpose */
00213 /*  ======= */
00214 
00215 /*  SCHKST  checks the symmetric eigenvalue problem routines. */
00216 
00217 /*     SSYTRD factors A as  U S U' , where ' means transpose, */
00218 /*     S is symmetric tridiagonal, and U is orthogonal. */
00219 /*     SSYTRD can use either just the lower or just the upper triangle */
00220 /*     of A; SCHKST checks both cases. */
00221 /*     U is represented as a product of Householder */
00222 /*     transformations, whose vectors are stored in the first */
00223 /*     n-1 columns of V, and whose scale factors are in TAU. */
00224 
00225 /*     SSPTRD does the same as SSYTRD, except that A and V are stored */
00226 /*     in "packed" format. */
00227 
00228 /*     SORGTR constructs the matrix U from the contents of V and TAU. */
00229 
00230 /*     SOPGTR constructs the matrix U from the contents of VP and TAU. */
00231 
00232 /*     SSTEQR factors S as  Z D1 Z' , where Z is the orthogonal */
00233 /*     matrix of eigenvectors and D1 is a diagonal matrix with */
00234 /*     the eigenvalues on the diagonal.  D2 is the matrix of */
00235 /*     eigenvalues computed when Z is not computed. */
00236 
00237 /*     SSTERF computes D3, the matrix of eigenvalues, by the */
00238 /*     PWK method, which does not yield eigenvectors. */
00239 
00240 /*     SPTEQR factors S as  Z4 D4 Z4' , for a */
00241 /*     symmetric positive definite tridiagonal matrix. */
00242 /*     D5 is the matrix of eigenvalues computed when Z is not */
00243 /*     computed. */
00244 
00245 /*     SSTEBZ computes selected eigenvalues.  WA1, WA2, and */
00246 /*     WA3 will denote eigenvalues computed to high */
00247 /*     absolute accuracy, with different range options. */
00248 /*     WR will denote eigenvalues computed to high relative */
00249 /*     accuracy. */
00250 
00251 /*     SSTEIN computes Y, the eigenvectors of S, given the */
00252 /*     eigenvalues. */
00253 
00254 /*     SSTEDC factors S as Z D1 Z' , where Z is the orthogonal */
00255 /*     matrix of eigenvectors and D1 is a diagonal matrix with */
00256 /*     the eigenvalues on the diagonal ('I' option). It may also */
00257 /*     update an input orthogonal matrix, usually the output */
00258 /*     from SSYTRD/SORGTR or SSPTRD/SOPGTR ('V' option). It may */
00259 /*     also just compute eigenvalues ('N' option). */
00260 
00261 /*     SSTEMR factors S as Z D1 Z' , where Z is the orthogonal */
00262 /*     matrix of eigenvectors and D1 is a diagonal matrix with */
00263 /*     the eigenvalues on the diagonal ('I' option).  SSTEMR */
00264 /*     uses the Relatively Robust Representation whenever possible. */
00265 
00266 /*  When SCHKST is called, a number of matrix "sizes" ("n's") and a */
00267 /*  number of matrix "types" are specified.  For each size ("n") */
00268 /*  and each type of matrix, one matrix will be generated and used */
00269 /*  to test the symmetric eigenroutines.  For each matrix, a number */
00270 /*  of tests will be performed: */
00271 
00272 /*  (1)     | A - V S V' | / ( |A| n ulp ) SSYTRD( UPLO='U', ... ) */
00273 
00274 /*  (2)     | I - UV' | / ( n ulp )        SORGTR( UPLO='U', ... ) */
00275 
00276 /*  (3)     | A - V S V' | / ( |A| n ulp ) SSYTRD( UPLO='L', ... ) */
00277 
00278 /*  (4)     | I - UV' | / ( n ulp )        SORGTR( UPLO='L', ... ) */
00279 
00280 /*  (5-8)   Same as 1-4, but for SSPTRD and SOPGTR. */
00281 
00282 /*  (9)     | S - Z D Z' | / ( |S| n ulp ) SSTEQR('V',...) */
00283 
00284 /*  (10)    | I - ZZ' | / ( n ulp )        SSTEQR('V',...) */
00285 
00286 /*  (11)    | D1 - D2 | / ( |D1| ulp )        SSTEQR('N',...) */
00287 
00288 /*  (12)    | D1 - D3 | / ( |D1| ulp )        SSTERF */
00289 
00290 /*  (13)    0 if the true eigenvalues (computed by sturm count) */
00291 /*          of S are within THRESH of */
00292 /*          those in D1.  2*THRESH if they are not.  (Tested using */
00293 /*          SSTECH) */
00294 
00295 /*  For S positive definite, */
00296 
00297 /*  (14)    | S - Z4 D4 Z4' | / ( |S| n ulp ) SPTEQR('V',...) */
00298 
00299 /*  (15)    | I - Z4 Z4' | / ( n ulp )        SPTEQR('V',...) */
00300 
00301 /*  (16)    | D4 - D5 | / ( 100 |D4| ulp )       SPTEQR('N',...) */
00302 
00303 /*  When S is also diagonally dominant by the factor gamma < 1, */
00304 
00305 /*  (17)    max | D4(i) - WR(i) | / ( |D4(i)| omega ) , */
00306 /*           i */
00307 /*          omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 */
00308 /*                                               SSTEBZ( 'A', 'E', ...) */
00309 
00310 /*  (18)    | WA1 - D3 | / ( |D3| ulp )          SSTEBZ( 'A', 'E', ...) */
00311 
00312 /*  (19)    ( max { min | WA2(i)-WA3(j) | } + */
00313 /*             i     j */
00314 /*            max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) */
00315 /*             i     j */
00316 /*                                               SSTEBZ( 'I', 'E', ...) */
00317 
00318 /*  (20)    | S - Y WA1 Y' | / ( |S| n ulp )  SSTEBZ, SSTEIN */
00319 
00320 /*  (21)    | I - Y Y' | / ( n ulp )          SSTEBZ, SSTEIN */
00321 
00322 /*  (22)    | S - Z D Z' | / ( |S| n ulp )    SSTEDC('I') */
00323 
00324 /*  (23)    | I - ZZ' | / ( n ulp )           SSTEDC('I') */
00325 
00326 /*  (24)    | S - Z D Z' | / ( |S| n ulp )    SSTEDC('V') */
00327 
00328 /*  (25)    | I - ZZ' | / ( n ulp )           SSTEDC('V') */
00329 
00330 /*  (26)    | D1 - D2 | / ( |D1| ulp )           SSTEDC('V') and */
00331 /*                                               SSTEDC('N') */
00332 
00333 /*  Test 27 is disabled at the moment because SSTEMR does not */
00334 /*  guarantee high relatvie accuracy. */
00335 
00336 /*  (27)    max | D6(i) - WR(i) | / ( |D6(i)| omega ) , */
00337 /*           i */
00338 /*          omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 */
00339 /*                                               SSTEMR('V', 'A') */
00340 
00341 /*  (28)    max | D6(i) - WR(i) | / ( |D6(i)| omega ) , */
00342 /*           i */
00343 /*          omega = 2 (2n-1) ULP (1 + 8 gamma**2) / (1 - gamma)**4 */
00344 /*                                               SSTEMR('V', 'I') */
00345 
00346 /*  Tests 29 through 34 are disable at present because SSTEMR */
00347 /*  does not handle partial specturm requests. */
00348 
00349 /*  (29)    | S - Z D Z' | / ( |S| n ulp )    SSTEMR('V', 'I') */
00350 
00351 /*  (30)    | I - ZZ' | / ( n ulp )           SSTEMR('V', 'I') */
00352 
00353 /*  (31)    ( max { min | WA2(i)-WA3(j) | } + */
00354 /*             i     j */
00355 /*            max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) */
00356 /*             i     j */
00357 /*          SSTEMR('N', 'I') vs. SSTEMR('V', 'I') */
00358 
00359 /*  (32)    | S - Z D Z' | / ( |S| n ulp )    SSTEMR('V', 'V') */
00360 
00361 /*  (33)    | I - ZZ' | / ( n ulp )           SSTEMR('V', 'V') */
00362 
00363 /*  (34)    ( max { min | WA2(i)-WA3(j) | } + */
00364 /*             i     j */
00365 /*            max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) */
00366 /*             i     j */
00367 /*          SSTEMR('N', 'V') vs. SSTEMR('V', 'V') */
00368 
00369 /*  (35)    | S - Z D Z' | / ( |S| n ulp )    SSTEMR('V', 'A') */
00370 
00371 /*  (36)    | I - ZZ' | / ( n ulp )           SSTEMR('V', 'A') */
00372 
00373 /*  (37)    ( max { min | WA2(i)-WA3(j) | } + */
00374 /*             i     j */
00375 /*            max { min | WA3(i)-WA2(j) | } ) / ( |D3| ulp ) */
00376 /*             i     j */
00377 /*          SSTEMR('N', 'A') vs. SSTEMR('V', 'A') */
00378 
00379 /*  The "sizes" are specified by an array NN(1:NSIZES); the value of */
00380 /*  each element NN(j) specifies one size. */
00381 /*  The "types" are specified by a logical array DOTYPE( 1:NTYPES ); */
00382 /*  if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. */
00383 /*  Currently, the list of possible types is: */
00384 
00385 /*  (1)  The zero matrix. */
00386 /*  (2)  The identity matrix. */
00387 
00388 /*  (3)  A diagonal matrix with evenly spaced entries */
00389 /*       1, ..., ULP  and random signs. */
00390 /*       (ULP = (first number larger than 1) - 1 ) */
00391 /*  (4)  A diagonal matrix with geometrically spaced entries */
00392 /*       1, ..., ULP  and random signs. */
00393 /*  (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP */
00394 /*       and random signs. */
00395 
00396 /*  (6)  Same as (4), but multiplied by SQRT( overflow threshold ) */
00397 /*  (7)  Same as (4), but multiplied by SQRT( underflow threshold ) */
00398 
00399 /*  (8)  A matrix of the form  U' D U, where U is orthogonal and */
00400 /*       D has evenly spaced entries 1, ..., ULP with random signs */
00401 /*       on the diagonal. */
00402 
00403 /*  (9)  A matrix of the form  U' D U, where U is orthogonal and */
00404 /*       D has geometrically spaced entries 1, ..., ULP with random */
00405 /*       signs on the diagonal. */
00406 
00407 /*  (10) A matrix of the form  U' D U, where U is orthogonal and */
00408 /*       D has "clustered" entries 1, ULP,..., ULP with random */
00409 /*       signs on the diagonal. */
00410 
00411 /*  (11) Same as (8), but multiplied by SQRT( overflow threshold ) */
00412 /*  (12) Same as (8), but multiplied by SQRT( underflow threshold ) */
00413 
00414 /*  (13) Symmetric matrix with random entries chosen from (-1,1). */
00415 /*  (14) Same as (13), but multiplied by SQRT( overflow threshold ) */
00416 /*  (15) Same as (13), but multiplied by SQRT( underflow threshold ) */
00417 /*  (16) Same as (8), but diagonal elements are all positive. */
00418 /*  (17) Same as (9), but diagonal elements are all positive. */
00419 /*  (18) Same as (10), but diagonal elements are all positive. */
00420 /*  (19) Same as (16), but multiplied by SQRT( overflow threshold ) */
00421 /*  (20) Same as (16), but multiplied by SQRT( underflow threshold ) */
00422 /*  (21) A diagonally dominant tridiagonal matrix with geometrically */
00423 /*       spaced diagonal entries 1, ..., ULP. */
00424 
00425 /*  Arguments */
00426 /*  ========= */
00427 
00428 /*  NSIZES  (input) INTEGER */
00429 /*          The number of sizes of matrices to use.  If it is zero, */
00430 /*          SCHKST does nothing.  It must be at least zero. */
00431 
00432 /*  NN      (input) INTEGER array, dimension (NSIZES) */
00433 /*          An array containing the sizes to be used for the matrices. */
00434 /*          Zero values will be skipped.  The values must be at least */
00435 /*          zero. */
00436 
00437 /*  NTYPES  (input) INTEGER */
00438 /*          The number of elements in DOTYPE.   If it is zero, SCHKST */
00439 /*          does nothing.  It must be at least zero.  If it is MAXTYP+1 */
00440 /*          and NSIZES is 1, then an additional type, MAXTYP+1 is */
00441 /*          defined, which is to use whatever matrix is in A.  This */
00442 /*          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and */
00443 /*          DOTYPE(MAXTYP+1) is .TRUE. . */
00444 
00445 /*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) */
00446 /*          If DOTYPE(j) is .TRUE., then for each size in NN a */
00447 /*          matrix of that size and of type j will be generated. */
00448 /*          If NTYPES is smaller than the maximum number of types */
00449 /*          defined (PARAMETER MAXTYP), then types NTYPES+1 through */
00450 /*          MAXTYP will not be generated.  If NTYPES is larger */
00451 /*          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) */
00452 /*          will be ignored. */
00453 
00454 /*  ISEED   (input/output) INTEGER array, dimension (4) */
00455 /*          On entry ISEED specifies the seed of the random number */
00456 /*          generator. The array elements should be between 0 and 4095; */
00457 /*          if not they will be reduced mod 4096.  Also, ISEED(4) must */
00458 /*          be odd.  The random number generator uses a linear */
00459 /*          congruential sequence limited to small integers, and so */
00460 /*          should produce machine independent random numbers. The */
00461 /*          values of ISEED are changed on exit, and can be used in the */
00462 /*          next call to SCHKST to continue the same random number */
00463 /*          sequence. */
00464 
00465 /*  THRESH  (input) REAL */
00466 /*          A test will count as "failed" if the "error", computed as */
00467 /*          described above, exceeds THRESH.  Note that the error */
00468 /*          is scaled to be O(1), so THRESH should be a reasonably */
00469 /*          small multiple of 1, e.g., 10 or 100.  In particular, */
00470 /*          it should not depend on the precision (single vs. double) */
00471 /*          or the size of the matrix.  It must be at least zero. */
00472 
00473 /*  NOUNIT  (input) INTEGER */
00474 /*          The FORTRAN unit number for printing out error messages */
00475 /*          (e.g., if a routine returns IINFO not equal to 0.) */
00476 
00477 /*  A       (input/workspace/output) REAL array of */
00478 /*                                  dimension ( LDA , max(NN) ) */
00479 /*          Used to hold the matrix whose eigenvalues are to be */
00480 /*          computed.  On exit, A contains the last matrix actually */
00481 /*          used. */
00482 
00483 /*  LDA     (input) INTEGER */
00484 /*          The leading dimension of A.  It must be at */
00485 /*          least 1 and at least max( NN ). */
00486 
00487 /*  AP      (workspace) REAL array of */
00488 /*                      dimension( max(NN)*max(NN+1)/2 ) */
00489 /*          The matrix A stored in packed format. */
00490 
00491 /*  SD      (workspace/output) REAL array of */
00492 /*                             dimension( max(NN) ) */
00493 /*          The diagonal of the tridiagonal matrix computed by SSYTRD. */
00494 /*          On exit, SD and SE contain the tridiagonal form of the */
00495 /*          matrix in A. */
00496 
00497 /*  SE      (workspace/output) REAL array of */
00498 /*                             dimension( max(NN) ) */
00499 /*          The off-diagonal of the tridiagonal matrix computed by */
00500 /*          SSYTRD.  On exit, SD and SE contain the tridiagonal form of */
00501 /*          the matrix in A. */
00502 
00503 /*  D1      (workspace/output) REAL array of */
00504 /*                             dimension( max(NN) ) */
00505 /*          The eigenvalues of A, as computed by SSTEQR simlutaneously */
00506 /*          with Z.  On exit, the eigenvalues in D1 correspond with the */
00507 /*          matrix in A. */
00508 
00509 /*  D2      (workspace/output) REAL array of */
00510 /*                             dimension( max(NN) ) */
00511 /*          The eigenvalues of A, as computed by SSTEQR if Z is not */
00512 /*          computed.  On exit, the eigenvalues in D2 correspond with */
00513 /*          the matrix in A. */
00514 
00515 /*  D3      (workspace/output) REAL array of */
00516 /*                             dimension( max(NN) ) */
00517 /*          The eigenvalues of A, as computed by SSTERF.  On exit, the */
00518 /*          eigenvalues in D3 correspond with the matrix in A. */
00519 
00520 /*  U       (workspace/output) REAL array of */
00521 /*                             dimension( LDU, max(NN) ). */
00522 /*          The orthogonal matrix computed by SSYTRD + SORGTR. */
00523 
00524 /*  LDU     (input) INTEGER */
00525 /*          The leading dimension of U, Z, and V.  It must be at least 1 */
00526 /*          and at least max( NN ). */
00527 
00528 /*  V       (workspace/output) REAL array of */
00529 /*                             dimension( LDU, max(NN) ). */
00530 /*          The Housholder vectors computed by SSYTRD in reducing A to */
00531 /*          tridiagonal form.  The vectors computed with UPLO='U' are */
00532 /*          in the upper triangle, and the vectors computed with UPLO='L' */
00533 /*          are in the lower triangle.  (As described in SSYTRD, the */
00534 /*          sub- and superdiagonal are not set to 1, although the */
00535 /*          true Householder vector has a 1 in that position.  The */
00536 /*          routines that use V, such as SORGTR, set those entries to */
00537 /*          1 before using them, and then restore them later.) */
00538 
00539 /*  VP      (workspace) REAL array of */
00540 /*                      dimension( max(NN)*max(NN+1)/2 ) */
00541 /*          The matrix V stored in packed format. */
00542 
00543 /*  TAU     (workspace/output) REAL array of */
00544 /*                             dimension( max(NN) ) */
00545 /*          The Householder factors computed by SSYTRD in reducing A */
00546 /*          to tridiagonal form. */
00547 
00548 /*  Z       (workspace/output) REAL array of */
00549 /*                             dimension( LDU, max(NN) ). */
00550 /*          The orthogonal matrix of eigenvectors computed by SSTEQR, */
00551 /*          SPTEQR, and SSTEIN. */
00552 
00553 /*  WORK    (workspace/output) REAL array of */
00554 /*                      dimension( LWORK ) */
00555 
00556 /*  LWORK   (input) INTEGER */
00557 /*          The number of entries in WORK.  This must be at least */
00558 /*          1 + 4 * Nmax + 2 * Nmax * lg Nmax + 3 * Nmax**2 */
00559 /*          where Nmax = max( NN(j), 2 ) and lg = log base 2. */
00560 
00561 /*  IWORK   (workspace/output) INTEGER array, */
00562 /*             dimension (6 + 6*Nmax + 5 * Nmax * lg Nmax ) */
00563 /*          where Nmax = max( NN(j), 2 ) and lg = log base 2. */
00564 /*          Workspace. */
00565 
00566 /*  RESULT  (output) REAL array, dimension (26) */
00567 /*          The values computed by the tests described above. */
00568 /*          The values are currently limited to 1/ulp, to avoid */
00569 /*          overflow. */
00570 
00571 /*  INFO    (output) INTEGER */
00572 /*          If 0, then everything ran OK. */
00573 /*           -1: NSIZES < 0 */
00574 /*           -2: Some NN(j) < 0 */
00575 /*           -3: NTYPES < 0 */
00576 /*           -5: THRESH < 0 */
00577 /*           -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ). */
00578 /*          -23: LDU < 1 or LDU < NMAX. */
00579 /*          -29: LWORK too small. */
00580 /*          If  SLATMR, SLATMS, SSYTRD, SORGTR, SSTEQR, SSTERF, */
00581 /*              or SORMC2 returns an error code, the */
00582 /*              absolute value of it is returned. */
00583 
00584 /* ----------------------------------------------------------------------- */
00585 
00586 /*       Some Local Variables and Parameters: */
00587 /*       ---- ----- --------- --- ---------- */
00588 /*       ZERO, ONE       Real 0 and 1. */
00589 /*       MAXTYP          The number of types defined. */
00590 /*       NTEST           The number of tests performed, or which can */
00591 /*                       be performed so far, for the current matrix. */
00592 /*       NTESTT          The total number of tests performed so far. */
00593 /*       NBLOCK          Blocksize as returned by ENVIR. */
00594 /*       NMAX            Largest value in NN. */
00595 /*       NMATS           The number of matrices generated so far. */
00596 /*       NERRS           The number of tests which have exceeded THRESH */
00597 /*                       so far. */
00598 /*       COND, IMODE     Values to be passed to the matrix generators. */
00599 /*       ANORM           Norm of A; passed to matrix generators. */
00600 
00601 /*       OVFL, UNFL      Overflow and underflow thresholds. */
00602 /*       ULP, ULPINV     Finest relative precision and its inverse. */
00603 /*       RTOVFL, RTUNFL  Square roots of the previous 2 values. */
00604 /*               The following four arrays decode JTYPE: */
00605 /*       KTYPE(j)        The general type (1-10) for type "j". */
00606 /*       KMODE(j)        The MODE value to be passed to the matrix */
00607 /*                       generator for type "j". */
00608 /*       KMAGN(j)        The order of magnitude ( O(1), */
00609 /*                       O(overflow^(1/2) ), O(underflow^(1/2) ) */
00610 
00611 /*  ===================================================================== */
00612 
00613 /*     .. Parameters .. */
00614 /*     .. */
00615 /*     .. Local Scalars .. */
00616 /*     .. */
00617 /*     .. Local Arrays .. */
00618 /*     .. */
00619 /*     .. External Functions .. */
00620 /*     .. */
00621 /*     .. External Subroutines .. */
00622 /*     .. */
00623 /*     .. Intrinsic Functions .. */
00624 /*     .. */
00625 /*     .. Data statements .. */
00626     /* Parameter adjustments */
00627     --nn;
00628     --dotype;
00629     --iseed;
00630     a_dim1 = *lda;
00631     a_offset = 1 + a_dim1;
00632     a -= a_offset;
00633     --ap;
00634     --sd;
00635     --se;
00636     --d1;
00637     --d2;
00638     --d3;
00639     --d4;
00640     --d5;
00641     --wa1;
00642     --wa2;
00643     --wa3;
00644     --wr;
00645     z_dim1 = *ldu;
00646     z_offset = 1 + z_dim1;
00647     z__ -= z_offset;
00648     v_dim1 = *ldu;
00649     v_offset = 1 + v_dim1;
00650     v -= v_offset;
00651     u_dim1 = *ldu;
00652     u_offset = 1 + u_dim1;
00653     u -= u_offset;
00654     --vp;
00655     --tau;
00656     --work;
00657     --iwork;
00658     --result;
00659 
00660     /* Function Body */
00661 /*     .. */
00662 /*     .. Executable Statements .. */
00663 
00664 /*     Keep ftnchek happy */
00665     idumma[0] = 1;
00666 
00667 /*     Check for errors */
00668 
00669     ntestt = 0;
00670     *info = 0;
00671 
00672 /*     Important constants */
00673 
00674     badnn = FALSE_;
00675     tryrac = TRUE_;
00676     nmax = 1;
00677     i__1 = *nsizes;
00678     for (j = 1; j <= i__1; ++j) {
00679 /* Computing MAX */
00680         i__2 = nmax, i__3 = nn[j];
00681         nmax = max(i__2,i__3);
00682         if (nn[j] < 0) {
00683             badnn = TRUE_;
00684         }
00685 /* L10: */
00686     }
00687 
00688     nblock = ilaenv_(&c__1, "SSYTRD", "L", &nmax, &c_n1, &c_n1, &c_n1);
00689 /* Computing MIN */
00690     i__1 = nmax, i__2 = max(1,nblock);
00691     nblock = min(i__1,i__2);
00692 
00693 /*     Check for errors */
00694 
00695     if (*nsizes < 0) {
00696         *info = -1;
00697     } else if (badnn) {
00698         *info = -2;
00699     } else if (*ntypes < 0) {
00700         *info = -3;
00701     } else if (*lda < nmax) {
00702         *info = -9;
00703     } else if (*ldu < nmax) {
00704         *info = -23;
00705     } else /* if(complicated condition) */ {
00706 /* Computing 2nd power */
00707         i__1 = max(2,nmax);
00708         if (i__1 * i__1 << 1 > *lwork) {
00709             *info = -29;
00710         }
00711     }
00712 
00713     if (*info != 0) {
00714         i__1 = -(*info);
00715         xerbla_("SCHKST", &i__1);
00716         return 0;
00717     }
00718 
00719 /*     Quick return if possible */
00720 
00721     if (*nsizes == 0 || *ntypes == 0) {
00722         return 0;
00723     }
00724 
00725 /*     More Important constants */
00726 
00727     unfl = slamch_("Safe minimum");
00728     ovfl = 1.f / unfl;
00729     slabad_(&unfl, &ovfl);
00730     ulp = slamch_("Epsilon") * slamch_("Base");
00731     ulpinv = 1.f / ulp;
00732     log2ui = (integer) (log(ulpinv) / log(2.f));
00733     rtunfl = sqrt(unfl);
00734     rtovfl = sqrt(ovfl);
00735 
00736 /*     Loop over sizes, types */
00737 
00738     for (i__ = 1; i__ <= 4; ++i__) {
00739         iseed2[i__ - 1] = iseed[i__];
00740 /* L20: */
00741     }
00742     nerrs = 0;
00743     nmats = 0;
00744 
00745     i__1 = *nsizes;
00746     for (jsize = 1; jsize <= i__1; ++jsize) {
00747         n = nn[jsize];
00748         if (n > 0) {
00749             lgn = (integer) (log((real) n) / log(2.f));
00750             if (pow_ii(&c__2, &lgn) < n) {
00751                 ++lgn;
00752             }
00753             if (pow_ii(&c__2, &lgn) < n) {
00754                 ++lgn;
00755             }
00756 /* Computing 2nd power */
00757             i__2 = n;
00758             lwedc = (n << 2) + 1 + (n << 1) * lgn + i__2 * i__2 * 3;
00759             liwedc = n * 6 + 6 + n * 5 * lgn;
00760         } else {
00761             lwedc = 8;
00762             liwedc = 12;
00763         }
00764         nap = n * (n + 1) / 2;
00765         aninv = 1.f / (real) max(1,n);
00766 
00767         if (*nsizes != 1) {
00768             mtypes = min(21,*ntypes);
00769         } else {
00770             mtypes = min(22,*ntypes);
00771         }
00772 
00773         i__2 = mtypes;
00774         for (jtype = 1; jtype <= i__2; ++jtype) {
00775             if (! dotype[jtype]) {
00776                 goto L300;
00777             }
00778             ++nmats;
00779             ntest = 0;
00780 
00781             for (j = 1; j <= 4; ++j) {
00782                 ioldsd[j - 1] = iseed[j];
00783 /* L30: */
00784             }
00785 
00786 /*           Compute "A" */
00787 
00788 /*           Control parameters: */
00789 
00790 /*               KMAGN  KMODE        KTYPE */
00791 /*           =1  O(1)   clustered 1  zero */
00792 /*           =2  large  clustered 2  identity */
00793 /*           =3  small  exponential  (none) */
00794 /*           =4         arithmetic   diagonal, (w/ eigenvalues) */
00795 /*           =5         random log   symmetric, w/ eigenvalues */
00796 /*           =6         random       (none) */
00797 /*           =7                      random diagonal */
00798 /*           =8                      random symmetric */
00799 /*           =9                      positive definite */
00800 /*           =10                     diagonally dominant tridiagonal */
00801 
00802             if (mtypes > 21) {
00803                 goto L100;
00804             }
00805 
00806             itype = ktype[jtype - 1];
00807             imode = kmode[jtype - 1];
00808 
00809 /*           Compute norm */
00810 
00811             switch (kmagn[jtype - 1]) {
00812                 case 1:  goto L40;
00813                 case 2:  goto L50;
00814                 case 3:  goto L60;
00815             }
00816 
00817 L40:
00818             anorm = 1.f;
00819             goto L70;
00820 
00821 L50:
00822             anorm = rtovfl * ulp * aninv;
00823             goto L70;
00824 
00825 L60:
00826             anorm = rtunfl * n * ulpinv;
00827             goto L70;
00828 
00829 L70:
00830 
00831             slaset_("Full", lda, &n, &c_b25, &c_b25, &a[a_offset], lda);
00832             iinfo = 0;
00833             if (jtype <= 15) {
00834                 cond = ulpinv;
00835             } else {
00836                 cond = ulpinv * aninv / 10.f;
00837             }
00838 
00839 /*           Special Matrices -- Identity & Jordan block */
00840 
00841 /*              Zero */
00842 
00843             if (itype == 1) {
00844                 iinfo = 0;
00845 
00846             } else if (itype == 2) {
00847 
00848 /*              Identity */
00849 
00850                 i__3 = n;
00851                 for (jc = 1; jc <= i__3; ++jc) {
00852                     a[jc + jc * a_dim1] = anorm;
00853 /* L80: */
00854                 }
00855 
00856             } else if (itype == 4) {
00857 
00858 /*              Diagonal Matrix, [Eigen]values Specified */
00859 
00860                 slatms_(&n, &n, "S", &iseed[1], "S", &work[1], &imode, &cond, 
00861                         &anorm, &c__0, &c__0, "N", &a[a_offset], lda, &work[n 
00862                         + 1], &iinfo);
00863 
00864 
00865             } else if (itype == 5) {
00866 
00867 /*              Symmetric, eigenvalues specified */
00868 
00869                 slatms_(&n, &n, "S", &iseed[1], "S", &work[1], &imode, &cond, 
00870                         &anorm, &n, &n, "N", &a[a_offset], lda, &work[n + 1], 
00871                         &iinfo);
00872 
00873             } else if (itype == 7) {
00874 
00875 /*              Diagonal, random eigenvalues */
00876 
00877                 slatmr_(&n, &n, "S", &iseed[1], "S", &work[1], &c__6, &c_b39, 
00878                         &c_b39, "T", "N", &work[n + 1], &c__1, &c_b39, &work[(
00879                         n << 1) + 1], &c__1, &c_b39, "N", idumma, &c__0, &
00880                         c__0, &c_b25, &anorm, "NO", &a[a_offset], lda, &iwork[
00881                         1], &iinfo);
00882 
00883             } else if (itype == 8) {
00884 
00885 /*              Symmetric, random eigenvalues */
00886 
00887                 slatmr_(&n, &n, "S", &iseed[1], "S", &work[1], &c__6, &c_b39, 
00888                         &c_b39, "T", "N", &work[n + 1], &c__1, &c_b39, &work[(
00889                         n << 1) + 1], &c__1, &c_b39, "N", idumma, &n, &n, &
00890                         c_b25, &anorm, "NO", &a[a_offset], lda, &iwork[1], &
00891                         iinfo);
00892 
00893             } else if (itype == 9) {
00894 
00895 /*              Positive definite, eigenvalues specified. */
00896 
00897                 slatms_(&n, &n, "S", &iseed[1], "P", &work[1], &imode, &cond, 
00898                         &anorm, &n, &n, "N", &a[a_offset], lda, &work[n + 1], 
00899                         &iinfo);
00900 
00901             } else if (itype == 10) {
00902 
00903 /*              Positive definite tridiagonal, eigenvalues specified. */
00904 
00905                 slatms_(&n, &n, "S", &iseed[1], "P", &work[1], &imode, &cond, 
00906                         &anorm, &c__1, &c__1, "N", &a[a_offset], lda, &work[n 
00907                         + 1], &iinfo);
00908                 i__3 = n;
00909                 for (i__ = 2; i__ <= i__3; ++i__) {
00910                     temp1 = (r__1 = a[i__ - 1 + i__ * a_dim1], dabs(r__1)) / 
00911                             sqrt((r__2 = a[i__ - 1 + (i__ - 1) * a_dim1] * a[
00912                             i__ + i__ * a_dim1], dabs(r__2)));
00913                     if (temp1 > .5f) {
00914                         a[i__ - 1 + i__ * a_dim1] = sqrt((r__1 = a[i__ - 1 + (
00915                                 i__ - 1) * a_dim1] * a[i__ + i__ * a_dim1], 
00916                                 dabs(r__1))) * .5f;
00917                         a[i__ + (i__ - 1) * a_dim1] = a[i__ - 1 + i__ * 
00918                                 a_dim1];
00919                     }
00920 /* L90: */
00921                 }
00922 
00923             } else {
00924 
00925                 iinfo = 1;
00926             }
00927 
00928             if (iinfo != 0) {
00929                 io___40.ciunit = *nounit;
00930                 s_wsfe(&io___40);
00931                 do_fio(&c__1, "Generator", (ftnlen)9);
00932                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00933                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00934                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00935                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00936                 e_wsfe();
00937                 *info = abs(iinfo);
00938                 return 0;
00939             }
00940 
00941 L100:
00942 
00943 /*           Call SSYTRD and SORGTR to compute S and U from */
00944 /*           upper triangle. */
00945 
00946             slacpy_("U", &n, &n, &a[a_offset], lda, &v[v_offset], ldu);
00947 
00948             ntest = 1;
00949             ssytrd_("U", &n, &v[v_offset], ldu, &sd[1], &se[1], &tau[1], &
00950                     work[1], lwork, &iinfo);
00951 
00952             if (iinfo != 0) {
00953                 io___41.ciunit = *nounit;
00954                 s_wsfe(&io___41);
00955                 do_fio(&c__1, "SSYTRD(U)", (ftnlen)9);
00956                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00957                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00958                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00959                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00960                 e_wsfe();
00961                 *info = abs(iinfo);
00962                 if (iinfo < 0) {
00963                     return 0;
00964                 } else {
00965                     result[1] = ulpinv;
00966                     goto L280;
00967                 }
00968             }
00969 
00970             slacpy_("U", &n, &n, &v[v_offset], ldu, &u[u_offset], ldu);
00971 
00972             ntest = 2;
00973             sorgtr_("U", &n, &u[u_offset], ldu, &tau[1], &work[1], lwork, &
00974                     iinfo);
00975             if (iinfo != 0) {
00976                 io___42.ciunit = *nounit;
00977                 s_wsfe(&io___42);
00978                 do_fio(&c__1, "SORGTR(U)", (ftnlen)9);
00979                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00980                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00981                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00982                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00983                 e_wsfe();
00984                 *info = abs(iinfo);
00985                 if (iinfo < 0) {
00986                     return 0;
00987                 } else {
00988                     result[2] = ulpinv;
00989                     goto L280;
00990                 }
00991             }
00992 
00993 /*           Do tests 1 and 2 */
00994 
00995             ssyt21_(&c__2, "Upper", &n, &c__1, &a[a_offset], lda, &sd[1], &se[
00996                     1], &u[u_offset], ldu, &v[v_offset], ldu, &tau[1], &work[
00997                     1], &result[1]);
00998             ssyt21_(&c__3, "Upper", &n, &c__1, &a[a_offset], lda, &sd[1], &se[
00999                     1], &u[u_offset], ldu, &v[v_offset], ldu, &tau[1], &work[
01000                     1], &result[2]);
01001 
01002 /*           Call SSYTRD and SORGTR to compute S and U from */
01003 /*           lower triangle, do tests. */
01004 
01005             slacpy_("L", &n, &n, &a[a_offset], lda, &v[v_offset], ldu);
01006 
01007             ntest = 3;
01008             ssytrd_("L", &n, &v[v_offset], ldu, &sd[1], &se[1], &tau[1], &
01009                     work[1], lwork, &iinfo);
01010 
01011             if (iinfo != 0) {
01012                 io___43.ciunit = *nounit;
01013                 s_wsfe(&io___43);
01014                 do_fio(&c__1, "SSYTRD(L)", (ftnlen)9);
01015                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01016                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01017                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01018                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01019                 e_wsfe();
01020                 *info = abs(iinfo);
01021                 if (iinfo < 0) {
01022                     return 0;
01023                 } else {
01024                     result[3] = ulpinv;
01025                     goto L280;
01026                 }
01027             }
01028 
01029             slacpy_("L", &n, &n, &v[v_offset], ldu, &u[u_offset], ldu);
01030 
01031             ntest = 4;
01032             sorgtr_("L", &n, &u[u_offset], ldu, &tau[1], &work[1], lwork, &
01033                     iinfo);
01034             if (iinfo != 0) {
01035                 io___44.ciunit = *nounit;
01036                 s_wsfe(&io___44);
01037                 do_fio(&c__1, "SORGTR(L)", (ftnlen)9);
01038                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01039                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01040                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01041                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01042                 e_wsfe();
01043                 *info = abs(iinfo);
01044                 if (iinfo < 0) {
01045                     return 0;
01046                 } else {
01047                     result[4] = ulpinv;
01048                     goto L280;
01049                 }
01050             }
01051 
01052             ssyt21_(&c__2, "Lower", &n, &c__1, &a[a_offset], lda, &sd[1], &se[
01053                     1], &u[u_offset], ldu, &v[v_offset], ldu, &tau[1], &work[
01054                     1], &result[3]);
01055             ssyt21_(&c__3, "Lower", &n, &c__1, &a[a_offset], lda, &sd[1], &se[
01056                     1], &u[u_offset], ldu, &v[v_offset], ldu, &tau[1], &work[
01057                     1], &result[4]);
01058 
01059 /*           Store the upper triangle of A in AP */
01060 
01061             i__ = 0;
01062             i__3 = n;
01063             for (jc = 1; jc <= i__3; ++jc) {
01064                 i__4 = jc;
01065                 for (jr = 1; jr <= i__4; ++jr) {
01066                     ++i__;
01067                     ap[i__] = a[jr + jc * a_dim1];
01068 /* L110: */
01069                 }
01070 /* L120: */
01071             }
01072 
01073 /*           Call SSPTRD and SOPGTR to compute S and U from AP */
01074 
01075             scopy_(&nap, &ap[1], &c__1, &vp[1], &c__1);
01076 
01077             ntest = 5;
01078             ssptrd_("U", &n, &vp[1], &sd[1], &se[1], &tau[1], &iinfo);
01079 
01080             if (iinfo != 0) {
01081                 io___46.ciunit = *nounit;
01082                 s_wsfe(&io___46);
01083                 do_fio(&c__1, "SSPTRD(U)", (ftnlen)9);
01084                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01085                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01086                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01087                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01088                 e_wsfe();
01089                 *info = abs(iinfo);
01090                 if (iinfo < 0) {
01091                     return 0;
01092                 } else {
01093                     result[5] = ulpinv;
01094                     goto L280;
01095                 }
01096             }
01097 
01098             ntest = 6;
01099             sopgtr_("U", &n, &vp[1], &tau[1], &u[u_offset], ldu, &work[1], &
01100                     iinfo);
01101             if (iinfo != 0) {
01102                 io___47.ciunit = *nounit;
01103                 s_wsfe(&io___47);
01104                 do_fio(&c__1, "SOPGTR(U)", (ftnlen)9);
01105                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01106                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01107                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01108                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01109                 e_wsfe();
01110                 *info = abs(iinfo);
01111                 if (iinfo < 0) {
01112                     return 0;
01113                 } else {
01114                     result[6] = ulpinv;
01115                     goto L280;
01116                 }
01117             }
01118 
01119 /*           Do tests 5 and 6 */
01120 
01121             sspt21_(&c__2, "Upper", &n, &c__1, &ap[1], &sd[1], &se[1], &u[
01122                     u_offset], ldu, &vp[1], &tau[1], &work[1], &result[5]);
01123             sspt21_(&c__3, "Upper", &n, &c__1, &ap[1], &sd[1], &se[1], &u[
01124                     u_offset], ldu, &vp[1], &tau[1], &work[1], &result[6]);
01125 
01126 /*           Store the lower triangle of A in AP */
01127 
01128             i__ = 0;
01129             i__3 = n;
01130             for (jc = 1; jc <= i__3; ++jc) {
01131                 i__4 = n;
01132                 for (jr = jc; jr <= i__4; ++jr) {
01133                     ++i__;
01134                     ap[i__] = a[jr + jc * a_dim1];
01135 /* L130: */
01136                 }
01137 /* L140: */
01138             }
01139 
01140 /*           Call SSPTRD and SOPGTR to compute S and U from AP */
01141 
01142             scopy_(&nap, &ap[1], &c__1, &vp[1], &c__1);
01143 
01144             ntest = 7;
01145             ssptrd_("L", &n, &vp[1], &sd[1], &se[1], &tau[1], &iinfo);
01146 
01147             if (iinfo != 0) {
01148                 io___48.ciunit = *nounit;
01149                 s_wsfe(&io___48);
01150                 do_fio(&c__1, "SSPTRD(L)", (ftnlen)9);
01151                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01152                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01153                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01154                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01155                 e_wsfe();
01156                 *info = abs(iinfo);
01157                 if (iinfo < 0) {
01158                     return 0;
01159                 } else {
01160                     result[7] = ulpinv;
01161                     goto L280;
01162                 }
01163             }
01164 
01165             ntest = 8;
01166             sopgtr_("L", &n, &vp[1], &tau[1], &u[u_offset], ldu, &work[1], &
01167                     iinfo);
01168             if (iinfo != 0) {
01169                 io___49.ciunit = *nounit;
01170                 s_wsfe(&io___49);
01171                 do_fio(&c__1, "SOPGTR(L)", (ftnlen)9);
01172                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01173                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01174                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01175                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01176                 e_wsfe();
01177                 *info = abs(iinfo);
01178                 if (iinfo < 0) {
01179                     return 0;
01180                 } else {
01181                     result[8] = ulpinv;
01182                     goto L280;
01183                 }
01184             }
01185 
01186             sspt21_(&c__2, "Lower", &n, &c__1, &ap[1], &sd[1], &se[1], &u[
01187                     u_offset], ldu, &vp[1], &tau[1], &work[1], &result[7]);
01188             sspt21_(&c__3, "Lower", &n, &c__1, &ap[1], &sd[1], &se[1], &u[
01189                     u_offset], ldu, &vp[1], &tau[1], &work[1], &result[8]);
01190 
01191 /*           Call SSTEQR to compute D1, D2, and Z, do tests. */
01192 
01193 /*           Compute D1 and Z */
01194 
01195             scopy_(&n, &sd[1], &c__1, &d1[1], &c__1);
01196             if (n > 0) {
01197                 i__3 = n - 1;
01198                 scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01199             }
01200             slaset_("Full", &n, &n, &c_b25, &c_b39, &z__[z_offset], ldu);
01201 
01202             ntest = 9;
01203             ssteqr_("V", &n, &d1[1], &work[1], &z__[z_offset], ldu, &work[n + 
01204                     1], &iinfo);
01205             if (iinfo != 0) {
01206                 io___50.ciunit = *nounit;
01207                 s_wsfe(&io___50);
01208                 do_fio(&c__1, "SSTEQR(V)", (ftnlen)9);
01209                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01210                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01211                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01212                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01213                 e_wsfe();
01214                 *info = abs(iinfo);
01215                 if (iinfo < 0) {
01216                     return 0;
01217                 } else {
01218                     result[9] = ulpinv;
01219                     goto L280;
01220                 }
01221             }
01222 
01223 /*           Compute D2 */
01224 
01225             scopy_(&n, &sd[1], &c__1, &d2[1], &c__1);
01226             if (n > 0) {
01227                 i__3 = n - 1;
01228                 scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01229             }
01230 
01231             ntest = 11;
01232             ssteqr_("N", &n, &d2[1], &work[1], &work[n + 1], ldu, &work[n + 1]
01233 , &iinfo);
01234             if (iinfo != 0) {
01235                 io___51.ciunit = *nounit;
01236                 s_wsfe(&io___51);
01237                 do_fio(&c__1, "SSTEQR(N)", (ftnlen)9);
01238                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01239                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01240                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01241                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01242                 e_wsfe();
01243                 *info = abs(iinfo);
01244                 if (iinfo < 0) {
01245                     return 0;
01246                 } else {
01247                     result[11] = ulpinv;
01248                     goto L280;
01249                 }
01250             }
01251 
01252 /*           Compute D3 (using PWK method) */
01253 
01254             scopy_(&n, &sd[1], &c__1, &d3[1], &c__1);
01255             if (n > 0) {
01256                 i__3 = n - 1;
01257                 scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01258             }
01259 
01260             ntest = 12;
01261             ssterf_(&n, &d3[1], &work[1], &iinfo);
01262             if (iinfo != 0) {
01263                 io___52.ciunit = *nounit;
01264                 s_wsfe(&io___52);
01265                 do_fio(&c__1, "SSTERF", (ftnlen)6);
01266                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01267                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01268                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01269                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01270                 e_wsfe();
01271                 *info = abs(iinfo);
01272                 if (iinfo < 0) {
01273                     return 0;
01274                 } else {
01275                     result[12] = ulpinv;
01276                     goto L280;
01277                 }
01278             }
01279 
01280 /*           Do Tests 9 and 10 */
01281 
01282             sstt21_(&n, &c__0, &sd[1], &se[1], &d1[1], dumma, &z__[z_offset], 
01283                     ldu, &work[1], &result[9]);
01284 
01285 /*           Do Tests 11 and 12 */
01286 
01287             temp1 = 0.f;
01288             temp2 = 0.f;
01289             temp3 = 0.f;
01290             temp4 = 0.f;
01291 
01292             i__3 = n;
01293             for (j = 1; j <= i__3; ++j) {
01294 /* Computing MAX */
01295                 r__3 = temp1, r__4 = (r__1 = d1[j], dabs(r__1)), r__3 = max(
01296                         r__3,r__4), r__4 = (r__2 = d2[j], dabs(r__2));
01297                 temp1 = dmax(r__3,r__4);
01298 /* Computing MAX */
01299                 r__2 = temp2, r__3 = (r__1 = d1[j] - d2[j], dabs(r__1));
01300                 temp2 = dmax(r__2,r__3);
01301 /* Computing MAX */
01302                 r__3 = temp3, r__4 = (r__1 = d1[j], dabs(r__1)), r__3 = max(
01303                         r__3,r__4), r__4 = (r__2 = d3[j], dabs(r__2));
01304                 temp3 = dmax(r__3,r__4);
01305 /* Computing MAX */
01306                 r__2 = temp4, r__3 = (r__1 = d1[j] - d3[j], dabs(r__1));
01307                 temp4 = dmax(r__2,r__3);
01308 /* L150: */
01309             }
01310 
01311 /* Computing MAX */
01312             r__1 = unfl, r__2 = ulp * dmax(temp1,temp2);
01313             result[11] = temp2 / dmax(r__1,r__2);
01314 /* Computing MAX */
01315             r__1 = unfl, r__2 = ulp * dmax(temp3,temp4);
01316             result[12] = temp4 / dmax(r__1,r__2);
01317 
01318 /*           Do Test 13 -- Sturm Sequence Test of Eigenvalues */
01319 /*                         Go up by factors of two until it succeeds */
01320 
01321             ntest = 13;
01322             temp1 = *thresh * (.5f - ulp);
01323 
01324             i__3 = log2ui;
01325             for (j = 0; j <= i__3; ++j) {
01326                 sstech_(&n, &sd[1], &se[1], &d1[1], &temp1, &work[1], &iinfo);
01327                 if (iinfo == 0) {
01328                     goto L170;
01329                 }
01330                 temp1 *= 2.f;
01331 /* L160: */
01332             }
01333 
01334 L170:
01335             result[13] = temp1;
01336 
01337 /*           For positive definite matrices ( JTYPE.GT.15 ) call SPTEQR */
01338 /*           and do tests 14, 15, and 16 . */
01339 
01340             if (jtype > 15) {
01341 
01342 /*              Compute D4 and Z4 */
01343 
01344                 scopy_(&n, &sd[1], &c__1, &d4[1], &c__1);
01345                 if (n > 0) {
01346                     i__3 = n - 1;
01347                     scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01348                 }
01349                 slaset_("Full", &n, &n, &c_b25, &c_b39, &z__[z_offset], ldu);
01350 
01351                 ntest = 14;
01352                 spteqr_("V", &n, &d4[1], &work[1], &z__[z_offset], ldu, &work[
01353                         n + 1], &iinfo);
01354                 if (iinfo != 0) {
01355                     io___57.ciunit = *nounit;
01356                     s_wsfe(&io___57);
01357                     do_fio(&c__1, "SPTEQR(V)", (ftnlen)9);
01358                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01359                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01360                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01361                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
01362                             ;
01363                     e_wsfe();
01364                     *info = abs(iinfo);
01365                     if (iinfo < 0) {
01366                         return 0;
01367                     } else {
01368                         result[14] = ulpinv;
01369                         goto L280;
01370                     }
01371                 }
01372 
01373 /*              Do Tests 14 and 15 */
01374 
01375                 sstt21_(&n, &c__0, &sd[1], &se[1], &d4[1], dumma, &z__[
01376                         z_offset], ldu, &work[1], &result[14]);
01377 
01378 /*              Compute D5 */
01379 
01380                 scopy_(&n, &sd[1], &c__1, &d5[1], &c__1);
01381                 if (n > 0) {
01382                     i__3 = n - 1;
01383                     scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01384                 }
01385 
01386                 ntest = 16;
01387                 spteqr_("N", &n, &d5[1], &work[1], &z__[z_offset], ldu, &work[
01388                         n + 1], &iinfo);
01389                 if (iinfo != 0) {
01390                     io___58.ciunit = *nounit;
01391                     s_wsfe(&io___58);
01392                     do_fio(&c__1, "SPTEQR(N)", (ftnlen)9);
01393                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01394                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01395                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01396                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
01397                             ;
01398                     e_wsfe();
01399                     *info = abs(iinfo);
01400                     if (iinfo < 0) {
01401                         return 0;
01402                     } else {
01403                         result[16] = ulpinv;
01404                         goto L280;
01405                     }
01406                 }
01407 
01408 /*              Do Test 16 */
01409 
01410                 temp1 = 0.f;
01411                 temp2 = 0.f;
01412                 i__3 = n;
01413                 for (j = 1; j <= i__3; ++j) {
01414 /* Computing MAX */
01415                     r__3 = temp1, r__4 = (r__1 = d4[j], dabs(r__1)), r__3 = 
01416                             max(r__3,r__4), r__4 = (r__2 = d5[j], dabs(r__2));
01417                     temp1 = dmax(r__3,r__4);
01418 /* Computing MAX */
01419                     r__2 = temp2, r__3 = (r__1 = d4[j] - d5[j], dabs(r__1));
01420                     temp2 = dmax(r__2,r__3);
01421 /* L180: */
01422                 }
01423 
01424 /* Computing MAX */
01425                 r__1 = unfl, r__2 = ulp * 100.f * dmax(temp1,temp2);
01426                 result[16] = temp2 / dmax(r__1,r__2);
01427             } else {
01428                 result[14] = 0.f;
01429                 result[15] = 0.f;
01430                 result[16] = 0.f;
01431             }
01432 
01433 /*           Call SSTEBZ with different options and do tests 17-18. */
01434 
01435 /*              If S is positive definite and diagonally dominant, */
01436 /*              ask for all eigenvalues with high relative accuracy. */
01437 
01438             vl = 0.f;
01439             vu = 0.f;
01440             il = 0;
01441             iu = 0;
01442             if (jtype == 21) {
01443                 ntest = 17;
01444                 abstol = unfl + unfl;
01445                 sstebz_("A", "E", &n, &vl, &vu, &il, &iu, &abstol, &sd[1], &
01446                         se[1], &m, &nsplit, &wr[1], &iwork[1], &iwork[n + 1], 
01447                         &work[1], &iwork[(n << 1) + 1], &iinfo);
01448                 if (iinfo != 0) {
01449                     io___66.ciunit = *nounit;
01450                     s_wsfe(&io___66);
01451                     do_fio(&c__1, "SSTEBZ(A,rel)", (ftnlen)13);
01452                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01453                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01454                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01455                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
01456                             ;
01457                     e_wsfe();
01458                     *info = abs(iinfo);
01459                     if (iinfo < 0) {
01460                         return 0;
01461                     } else {
01462                         result[17] = ulpinv;
01463                         goto L280;
01464                     }
01465                 }
01466 
01467 /*              Do test 17 */
01468 
01469                 temp2 = (n * 2.f - 1.f) * 2.f * ulp * 3.f / .0625f;
01470 
01471                 temp1 = 0.f;
01472                 i__3 = n;
01473                 for (j = 1; j <= i__3; ++j) {
01474 /* Computing MAX */
01475                     r__3 = temp1, r__4 = (r__2 = d4[j] - wr[n - j + 1], dabs(
01476                             r__2)) / (abstol + (r__1 = d4[j], dabs(r__1)));
01477                     temp1 = dmax(r__3,r__4);
01478 /* L190: */
01479                 }
01480 
01481                 result[17] = temp1 / temp2;
01482             } else {
01483                 result[17] = 0.f;
01484             }
01485 
01486 /*           Now ask for all eigenvalues with high absolute accuracy. */
01487 
01488             ntest = 18;
01489             abstol = unfl + unfl;
01490             sstebz_("A", "E", &n, &vl, &vu, &il, &iu, &abstol, &sd[1], &se[1], 
01491                      &m, &nsplit, &wa1[1], &iwork[1], &iwork[n + 1], &work[1], 
01492                      &iwork[(n << 1) + 1], &iinfo);
01493             if (iinfo != 0) {
01494                 io___67.ciunit = *nounit;
01495                 s_wsfe(&io___67);
01496                 do_fio(&c__1, "SSTEBZ(A)", (ftnlen)9);
01497                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01498                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01499                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01500                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01501                 e_wsfe();
01502                 *info = abs(iinfo);
01503                 if (iinfo < 0) {
01504                     return 0;
01505                 } else {
01506                     result[18] = ulpinv;
01507                     goto L280;
01508                 }
01509             }
01510 
01511 /*           Do test 18 */
01512 
01513             temp1 = 0.f;
01514             temp2 = 0.f;
01515             i__3 = n;
01516             for (j = 1; j <= i__3; ++j) {
01517 /* Computing MAX */
01518                 r__3 = temp1, r__4 = (r__1 = d3[j], dabs(r__1)), r__3 = max(
01519                         r__3,r__4), r__4 = (r__2 = wa1[j], dabs(r__2));
01520                 temp1 = dmax(r__3,r__4);
01521 /* Computing MAX */
01522                 r__2 = temp2, r__3 = (r__1 = d3[j] - wa1[j], dabs(r__1));
01523                 temp2 = dmax(r__2,r__3);
01524 /* L200: */
01525             }
01526 
01527 /* Computing MAX */
01528             r__1 = unfl, r__2 = ulp * dmax(temp1,temp2);
01529             result[18] = temp2 / dmax(r__1,r__2);
01530 
01531 /*           Choose random values for IL and IU, and ask for the */
01532 /*           IL-th through IU-th eigenvalues. */
01533 
01534             ntest = 19;
01535             if (n <= 1) {
01536                 il = 1;
01537                 iu = n;
01538             } else {
01539                 il = (n - 1) * (integer) slarnd_(&c__1, iseed2) + 1;
01540                 iu = (n - 1) * (integer) slarnd_(&c__1, iseed2) + 1;
01541                 if (iu < il) {
01542                     itemp = iu;
01543                     iu = il;
01544                     il = itemp;
01545                 }
01546             }
01547 
01548             sstebz_("I", "E", &n, &vl, &vu, &il, &iu, &abstol, &sd[1], &se[1], 
01549                      &m2, &nsplit, &wa2[1], &iwork[1], &iwork[n + 1], &work[1]
01550 , &iwork[(n << 1) + 1], &iinfo);
01551             if (iinfo != 0) {
01552                 io___70.ciunit = *nounit;
01553                 s_wsfe(&io___70);
01554                 do_fio(&c__1, "SSTEBZ(I)", (ftnlen)9);
01555                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01556                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01557                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01558                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01559                 e_wsfe();
01560                 *info = abs(iinfo);
01561                 if (iinfo < 0) {
01562                     return 0;
01563                 } else {
01564                     result[19] = ulpinv;
01565                     goto L280;
01566                 }
01567             }
01568 
01569 /*           Determine the values VL and VU of the IL-th and IU-th */
01570 /*           eigenvalues and ask for all eigenvalues in this range. */
01571 
01572             if (n > 0) {
01573                 if (il != 1) {
01574 /* Computing MAX */
01575                     r__1 = (wa1[il] - wa1[il - 1]) * .5f, r__2 = ulp * anorm, 
01576                             r__1 = max(r__1,r__2), r__2 = rtunfl * 2.f;
01577                     vl = wa1[il] - dmax(r__1,r__2);
01578                 } else {
01579 /* Computing MAX */
01580                     r__1 = (wa1[n] - wa1[1]) * .5f, r__2 = ulp * anorm, r__1 =
01581                              max(r__1,r__2), r__2 = rtunfl * 2.f;
01582                     vl = wa1[1] - dmax(r__1,r__2);
01583                 }
01584                 if (iu != n) {
01585 /* Computing MAX */
01586                     r__1 = (wa1[iu + 1] - wa1[iu]) * .5f, r__2 = ulp * anorm, 
01587                             r__1 = max(r__1,r__2), r__2 = rtunfl * 2.f;
01588                     vu = wa1[iu] + dmax(r__1,r__2);
01589                 } else {
01590 /* Computing MAX */
01591                     r__1 = (wa1[n] - wa1[1]) * .5f, r__2 = ulp * anorm, r__1 =
01592                              max(r__1,r__2), r__2 = rtunfl * 2.f;
01593                     vu = wa1[n] + dmax(r__1,r__2);
01594                 }
01595             } else {
01596                 vl = 0.f;
01597                 vu = 1.f;
01598             }
01599 
01600             sstebz_("V", "E", &n, &vl, &vu, &il, &iu, &abstol, &sd[1], &se[1], 
01601                      &m3, &nsplit, &wa3[1], &iwork[1], &iwork[n + 1], &work[1]
01602 , &iwork[(n << 1) + 1], &iinfo);
01603             if (iinfo != 0) {
01604                 io___72.ciunit = *nounit;
01605                 s_wsfe(&io___72);
01606                 do_fio(&c__1, "SSTEBZ(V)", (ftnlen)9);
01607                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01608                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01609                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01610                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01611                 e_wsfe();
01612                 *info = abs(iinfo);
01613                 if (iinfo < 0) {
01614                     return 0;
01615                 } else {
01616                     result[19] = ulpinv;
01617                     goto L280;
01618                 }
01619             }
01620 
01621             if (m3 == 0 && n != 0) {
01622                 result[19] = ulpinv;
01623                 goto L280;
01624             }
01625 
01626 /*           Do test 19 */
01627 
01628             temp1 = ssxt1_(&c__1, &wa2[1], &m2, &wa3[1], &m3, &abstol, &ulp, &
01629                     unfl);
01630             temp2 = ssxt1_(&c__1, &wa3[1], &m3, &wa2[1], &m2, &abstol, &ulp, &
01631                     unfl);
01632             if (n > 0) {
01633 /* Computing MAX */
01634                 r__2 = (r__1 = wa1[n], dabs(r__1)), r__3 = dabs(wa1[1]);
01635                 temp3 = dmax(r__2,r__3);
01636             } else {
01637                 temp3 = 0.f;
01638             }
01639 
01640 /* Computing MAX */
01641             r__1 = unfl, r__2 = temp3 * ulp;
01642             result[19] = (temp1 + temp2) / dmax(r__1,r__2);
01643 
01644 /*           Call SSTEIN to compute eigenvectors corresponding to */
01645 /*           eigenvalues in WA1.  (First call SSTEBZ again, to make sure */
01646 /*           it returns these eigenvalues in the correct order.) */
01647 
01648             ntest = 21;
01649             sstebz_("A", "B", &n, &vl, &vu, &il, &iu, &abstol, &sd[1], &se[1], 
01650                      &m, &nsplit, &wa1[1], &iwork[1], &iwork[n + 1], &work[1], 
01651                      &iwork[(n << 1) + 1], &iinfo);
01652             if (iinfo != 0) {
01653                 io___73.ciunit = *nounit;
01654                 s_wsfe(&io___73);
01655                 do_fio(&c__1, "SSTEBZ(A,B)", (ftnlen)11);
01656                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01657                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01658                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01659                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01660                 e_wsfe();
01661                 *info = abs(iinfo);
01662                 if (iinfo < 0) {
01663                     return 0;
01664                 } else {
01665                     result[20] = ulpinv;
01666                     result[21] = ulpinv;
01667                     goto L280;
01668                 }
01669             }
01670 
01671             sstein_(&n, &sd[1], &se[1], &m, &wa1[1], &iwork[1], &iwork[n + 1], 
01672                      &z__[z_offset], ldu, &work[1], &iwork[(n << 1) + 1], &
01673                     iwork[n * 3 + 1], &iinfo);
01674             if (iinfo != 0) {
01675                 io___74.ciunit = *nounit;
01676                 s_wsfe(&io___74);
01677                 do_fio(&c__1, "SSTEIN", (ftnlen)6);
01678                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01679                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01680                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01681                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01682                 e_wsfe();
01683                 *info = abs(iinfo);
01684                 if (iinfo < 0) {
01685                     return 0;
01686                 } else {
01687                     result[20] = ulpinv;
01688                     result[21] = ulpinv;
01689                     goto L280;
01690                 }
01691             }
01692 
01693 /*           Do tests 20 and 21 */
01694 
01695             sstt21_(&n, &c__0, &sd[1], &se[1], &wa1[1], dumma, &z__[z_offset], 
01696                      ldu, &work[1], &result[20]);
01697 
01698 /*           Call SSTEDC(I) to compute D1 and Z, do tests. */
01699 
01700 /*           Compute D1 and Z */
01701 
01702             scopy_(&n, &sd[1], &c__1, &d1[1], &c__1);
01703             if (n > 0) {
01704                 i__3 = n - 1;
01705                 scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01706             }
01707             slaset_("Full", &n, &n, &c_b25, &c_b39, &z__[z_offset], ldu);
01708 
01709             ntest = 22;
01710             i__3 = lwedc - n;
01711             sstedc_("I", &n, &d1[1], &work[1], &z__[z_offset], ldu, &work[n + 
01712                     1], &i__3, &iwork[1], &liwedc, &iinfo);
01713             if (iinfo != 0) {
01714                 io___75.ciunit = *nounit;
01715                 s_wsfe(&io___75);
01716                 do_fio(&c__1, "SSTEDC(I)", (ftnlen)9);
01717                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01718                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01719                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01720                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01721                 e_wsfe();
01722                 *info = abs(iinfo);
01723                 if (iinfo < 0) {
01724                     return 0;
01725                 } else {
01726                     result[22] = ulpinv;
01727                     goto L280;
01728                 }
01729             }
01730 
01731 /*           Do Tests 22 and 23 */
01732 
01733             sstt21_(&n, &c__0, &sd[1], &se[1], &d1[1], dumma, &z__[z_offset], 
01734                     ldu, &work[1], &result[22]);
01735 
01736 /*           Call SSTEDC(V) to compute D1 and Z, do tests. */
01737 
01738 /*           Compute D1 and Z */
01739 
01740             scopy_(&n, &sd[1], &c__1, &d1[1], &c__1);
01741             if (n > 0) {
01742                 i__3 = n - 1;
01743                 scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01744             }
01745             slaset_("Full", &n, &n, &c_b25, &c_b39, &z__[z_offset], ldu);
01746 
01747             ntest = 24;
01748             i__3 = lwedc - n;
01749             sstedc_("V", &n, &d1[1], &work[1], &z__[z_offset], ldu, &work[n + 
01750                     1], &i__3, &iwork[1], &liwedc, &iinfo);
01751             if (iinfo != 0) {
01752                 io___76.ciunit = *nounit;
01753                 s_wsfe(&io___76);
01754                 do_fio(&c__1, "SSTEDC(V)", (ftnlen)9);
01755                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01756                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01757                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01758                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01759                 e_wsfe();
01760                 *info = abs(iinfo);
01761                 if (iinfo < 0) {
01762                     return 0;
01763                 } else {
01764                     result[24] = ulpinv;
01765                     goto L280;
01766                 }
01767             }
01768 
01769 /*           Do Tests 24 and 25 */
01770 
01771             sstt21_(&n, &c__0, &sd[1], &se[1], &d1[1], dumma, &z__[z_offset], 
01772                     ldu, &work[1], &result[24]);
01773 
01774 /*           Call SSTEDC(N) to compute D2, do tests. */
01775 
01776 /*           Compute D2 */
01777 
01778             scopy_(&n, &sd[1], &c__1, &d2[1], &c__1);
01779             if (n > 0) {
01780                 i__3 = n - 1;
01781                 scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01782             }
01783             slaset_("Full", &n, &n, &c_b25, &c_b39, &z__[z_offset], ldu);
01784 
01785             ntest = 26;
01786             i__3 = lwedc - n;
01787             sstedc_("N", &n, &d2[1], &work[1], &z__[z_offset], ldu, &work[n + 
01788                     1], &i__3, &iwork[1], &liwedc, &iinfo);
01789             if (iinfo != 0) {
01790                 io___77.ciunit = *nounit;
01791                 s_wsfe(&io___77);
01792                 do_fio(&c__1, "SSTEDC(N)", (ftnlen)9);
01793                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01794                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01795                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
01796                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
01797                 e_wsfe();
01798                 *info = abs(iinfo);
01799                 if (iinfo < 0) {
01800                     return 0;
01801                 } else {
01802                     result[26] = ulpinv;
01803                     goto L280;
01804                 }
01805             }
01806 
01807 /*           Do Test 26 */
01808 
01809             temp1 = 0.f;
01810             temp2 = 0.f;
01811 
01812             i__3 = n;
01813             for (j = 1; j <= i__3; ++j) {
01814 /* Computing MAX */
01815                 r__3 = temp1, r__4 = (r__1 = d1[j], dabs(r__1)), r__3 = max(
01816                         r__3,r__4), r__4 = (r__2 = d2[j], dabs(r__2));
01817                 temp1 = dmax(r__3,r__4);
01818 /* Computing MAX */
01819                 r__2 = temp2, r__3 = (r__1 = d1[j] - d2[j], dabs(r__1));
01820                 temp2 = dmax(r__2,r__3);
01821 /* L210: */
01822             }
01823 
01824 /* Computing MAX */
01825             r__1 = unfl, r__2 = ulp * dmax(temp1,temp2);
01826             result[26] = temp2 / dmax(r__1,r__2);
01827 
01828 /*           Only test SSTEMR if IEEE compliant */
01829 
01830             if (ilaenv_(&c__10, "SSTEMR", "VA", &c__1, &c__0, &c__0, &c__0) == 1 && ilaenv_(&c__11, "SSTEMR", 
01831                     "VA", &c__1, &c__0, &c__0, &c__0) ==
01832                      1) {
01833 
01834 /*           Call SSTEMR, do test 27 (relative eigenvalue accuracy) */
01835 
01836 /*              If S is positive definite and diagonally dominant, */
01837 /*              ask for all eigenvalues with high relative accuracy. */
01838 
01839                 vl = 0.f;
01840                 vu = 0.f;
01841                 il = 0;
01842                 iu = 0;
01843                 if (FALSE_) {
01844                     ntest = 27;
01845                     abstol = unfl + unfl;
01846                     i__3 = *lwork - (n << 1);
01847                     sstemr_("V", "A", &n, &sd[1], &se[1], &vl, &vu, &il, &iu, 
01848                             &m, &wr[1], &z__[z_offset], ldu, &n, &iwork[1], &
01849                             tryrac, &work[1], lwork, &iwork[(n << 1) + 1], &
01850                             i__3, &iinfo);
01851                     if (iinfo != 0) {
01852                         io___78.ciunit = *nounit;
01853                         s_wsfe(&io___78);
01854                         do_fio(&c__1, "SSTEMR(V,A,rel)", (ftnlen)15);
01855                         do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
01856                                 ;
01857                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01858                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
01859                                 ;
01860                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
01861                                 integer));
01862                         e_wsfe();
01863                         *info = abs(iinfo);
01864                         if (iinfo < 0) {
01865                             return 0;
01866                         } else {
01867                             result[27] = ulpinv;
01868                             goto L270;
01869                         }
01870                     }
01871 
01872 /*              Do test 27 */
01873 
01874                     temp2 = (n * 2.f - 1.f) * 2.f * ulp * 3.f / .0625f;
01875 
01876                     temp1 = 0.f;
01877                     i__3 = n;
01878                     for (j = 1; j <= i__3; ++j) {
01879 /* Computing MAX */
01880                         r__3 = temp1, r__4 = (r__2 = d4[j] - wr[n - j + 1], 
01881                                 dabs(r__2)) / (abstol + (r__1 = d4[j], dabs(
01882                                 r__1)));
01883                         temp1 = dmax(r__3,r__4);
01884 /* L220: */
01885                     }
01886 
01887                     result[27] = temp1 / temp2;
01888 
01889                     il = (n - 1) * (integer) slarnd_(&c__1, iseed2) + 1;
01890                     iu = (n - 1) * (integer) slarnd_(&c__1, iseed2) + 1;
01891                     if (iu < il) {
01892                         itemp = iu;
01893                         iu = il;
01894                         il = itemp;
01895                     }
01896 
01897                     if (FALSE_) {
01898                         ntest = 28;
01899                         abstol = unfl + unfl;
01900                         i__3 = *lwork - (n << 1);
01901                         sstemr_("V", "I", &n, &sd[1], &se[1], &vl, &vu, &il, &
01902                                 iu, &m, &wr[1], &z__[z_offset], ldu, &n, &
01903                                 iwork[1], &tryrac, &work[1], lwork, &iwork[(n 
01904                                 << 1) + 1], &i__3, &iinfo);
01905 
01906                         if (iinfo != 0) {
01907                             io___79.ciunit = *nounit;
01908                             s_wsfe(&io___79);
01909                             do_fio(&c__1, "SSTEMR(V,I,rel)", (ftnlen)15);
01910                             do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(
01911                                     integer));
01912                             do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer))
01913                                     ;
01914                             do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(
01915                                     integer));
01916                             do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
01917                                     integer));
01918                             e_wsfe();
01919                             *info = abs(iinfo);
01920                             if (iinfo < 0) {
01921                                 return 0;
01922                             } else {
01923                                 result[28] = ulpinv;
01924                                 goto L270;
01925                             }
01926                         }
01927 
01928 
01929 /*                 Do test 28 */
01930 
01931                         temp2 = (n * 2.f - 1.f) * 2.f * ulp * 3.f / .0625f;
01932 
01933                         temp1 = 0.f;
01934                         i__3 = iu;
01935                         for (j = il; j <= i__3; ++j) {
01936 /* Computing MAX */
01937                             r__3 = temp1, r__4 = (r__2 = wr[j - il + 1] - d4[
01938                                     n - j + 1], dabs(r__2)) / (abstol + (r__1 
01939                                     = wr[j - il + 1], dabs(r__1)));
01940                             temp1 = dmax(r__3,r__4);
01941 /* L230: */
01942                         }
01943 
01944                         result[28] = temp1 / temp2;
01945                     } else {
01946                         result[28] = 0.f;
01947                     }
01948                 } else {
01949                     result[27] = 0.f;
01950                     result[28] = 0.f;
01951                 }
01952 
01953 /*           Call SSTEMR(V,I) to compute D1 and Z, do tests. */
01954 
01955 /*           Compute D1 and Z */
01956 
01957                 scopy_(&n, &sd[1], &c__1, &d5[1], &c__1);
01958                 if (n > 0) {
01959                     i__3 = n - 1;
01960                     scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
01961                 }
01962                 slaset_("Full", &n, &n, &c_b25, &c_b39, &z__[z_offset], ldu);
01963 
01964                 if (FALSE_) {
01965                     ntest = 29;
01966                     il = (n - 1) * (integer) slarnd_(&c__1, iseed2) + 1;
01967                     iu = (n - 1) * (integer) slarnd_(&c__1, iseed2) + 1;
01968                     if (iu < il) {
01969                         itemp = iu;
01970                         iu = il;
01971                         il = itemp;
01972                     }
01973                     i__3 = *lwork - n;
01974                     i__4 = *liwork - (n << 1);
01975                     sstemr_("V", "I", &n, &d5[1], &work[1], &vl, &vu, &il, &
01976                             iu, &m, &d1[1], &z__[z_offset], ldu, &n, &iwork[1]
01977 , &tryrac, &work[n + 1], &i__3, &iwork[(n << 1) + 
01978                             1], &i__4, &iinfo);
01979                     if (iinfo != 0) {
01980                         io___80.ciunit = *nounit;
01981                         s_wsfe(&io___80);
01982                         do_fio(&c__1, "SSTEMR(V,I)", (ftnlen)11);
01983                         do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
01984                                 ;
01985                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01986                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
01987                                 ;
01988                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
01989                                 integer));
01990                         e_wsfe();
01991                         *info = abs(iinfo);
01992                         if (iinfo < 0) {
01993                             return 0;
01994                         } else {
01995                             result[29] = ulpinv;
01996                             goto L280;
01997                         }
01998                     }
01999 
02000 /*           Do Tests 29 and 30 */
02001 
02002                     sstt22_(&n, &m, &c__0, &sd[1], &se[1], &d1[1], dumma, &
02003                             z__[z_offset], ldu, &work[1], &m, &result[29]);
02004 
02005 /*           Call SSTEMR to compute D2, do tests. */
02006 
02007 /*           Compute D2 */
02008 
02009                     scopy_(&n, &sd[1], &c__1, &d5[1], &c__1);
02010                     if (n > 0) {
02011                         i__3 = n - 1;
02012                         scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
02013                     }
02014 
02015                     ntest = 31;
02016                     i__3 = *lwork - n;
02017                     i__4 = *liwork - (n << 1);
02018                     sstemr_("N", "I", &n, &d5[1], &work[1], &vl, &vu, &il, &
02019                             iu, &m, &d2[1], &z__[z_offset], ldu, &n, &iwork[1]
02020 , &tryrac, &work[n + 1], &i__3, &iwork[(n << 1) + 
02021                             1], &i__4, &iinfo);
02022                     if (iinfo != 0) {
02023                         io___81.ciunit = *nounit;
02024                         s_wsfe(&io___81);
02025                         do_fio(&c__1, "SSTEMR(N,I)", (ftnlen)11);
02026                         do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
02027                                 ;
02028                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
02029                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
02030                                 ;
02031                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
02032                                 integer));
02033                         e_wsfe();
02034                         *info = abs(iinfo);
02035                         if (iinfo < 0) {
02036                             return 0;
02037                         } else {
02038                             result[31] = ulpinv;
02039                             goto L280;
02040                         }
02041                     }
02042 
02043 /*           Do Test 31 */
02044 
02045                     temp1 = 0.f;
02046                     temp2 = 0.f;
02047 
02048                     i__3 = iu - il + 1;
02049                     for (j = 1; j <= i__3; ++j) {
02050 /* Computing MAX */
02051                         r__3 = temp1, r__4 = (r__1 = d1[j], dabs(r__1)), r__3 
02052                                 = max(r__3,r__4), r__4 = (r__2 = d2[j], dabs(
02053                                 r__2));
02054                         temp1 = dmax(r__3,r__4);
02055 /* Computing MAX */
02056                         r__2 = temp2, r__3 = (r__1 = d1[j] - d2[j], dabs(r__1)
02057                                 );
02058                         temp2 = dmax(r__2,r__3);
02059 /* L240: */
02060                     }
02061 
02062 /* Computing MAX */
02063                     r__1 = unfl, r__2 = ulp * dmax(temp1,temp2);
02064                     result[31] = temp2 / dmax(r__1,r__2);
02065 
02066 
02067 /*           Call SSTEMR(V,V) to compute D1 and Z, do tests. */
02068 
02069 /*           Compute D1 and Z */
02070 
02071                     scopy_(&n, &sd[1], &c__1, &d5[1], &c__1);
02072                     if (n > 0) {
02073                         i__3 = n - 1;
02074                         scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
02075                     }
02076                     slaset_("Full", &n, &n, &c_b25, &c_b39, &z__[z_offset], 
02077                             ldu);
02078 
02079                     ntest = 32;
02080 
02081                     if (n > 0) {
02082                         if (il != 1) {
02083 /* Computing MAX */
02084                             r__1 = (d2[il] - d2[il - 1]) * .5f, r__2 = ulp * 
02085                                     anorm, r__1 = max(r__1,r__2), r__2 = 
02086                                     rtunfl * 2.f;
02087                             vl = d2[il] - dmax(r__1,r__2);
02088                         } else {
02089 /* Computing MAX */
02090                             r__1 = (d2[n] - d2[1]) * .5f, r__2 = ulp * anorm, 
02091                                     r__1 = max(r__1,r__2), r__2 = rtunfl * 
02092                                     2.f;
02093                             vl = d2[1] - dmax(r__1,r__2);
02094                         }
02095                         if (iu != n) {
02096 /* Computing MAX */
02097                             r__1 = (d2[iu + 1] - d2[iu]) * .5f, r__2 = ulp * 
02098                                     anorm, r__1 = max(r__1,r__2), r__2 = 
02099                                     rtunfl * 2.f;
02100                             vu = d2[iu] + dmax(r__1,r__2);
02101                         } else {
02102 /* Computing MAX */
02103                             r__1 = (d2[n] - d2[1]) * .5f, r__2 = ulp * anorm, 
02104                                     r__1 = max(r__1,r__2), r__2 = rtunfl * 
02105                                     2.f;
02106                             vu = d2[n] + dmax(r__1,r__2);
02107                         }
02108                     } else {
02109                         vl = 0.f;
02110                         vu = 1.f;
02111                     }
02112 
02113                     i__3 = *lwork - n;
02114                     i__4 = *liwork - (n << 1);
02115                     sstemr_("V", "V", &n, &d5[1], &work[1], &vl, &vu, &il, &
02116                             iu, &m, &d1[1], &z__[z_offset], ldu, &n, &iwork[1]
02117 , &tryrac, &work[n + 1], &i__3, &iwork[(n << 1) + 
02118                             1], &i__4, &iinfo);
02119                     if (iinfo != 0) {
02120                         io___82.ciunit = *nounit;
02121                         s_wsfe(&io___82);
02122                         do_fio(&c__1, "SSTEMR(V,V)", (ftnlen)11);
02123                         do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
02124                                 ;
02125                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
02126                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
02127                                 ;
02128                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
02129                                 integer));
02130                         e_wsfe();
02131                         *info = abs(iinfo);
02132                         if (iinfo < 0) {
02133                             return 0;
02134                         } else {
02135                             result[32] = ulpinv;
02136                             goto L280;
02137                         }
02138                     }
02139 
02140 /*           Do Tests 32 and 33 */
02141 
02142                     sstt22_(&n, &m, &c__0, &sd[1], &se[1], &d1[1], dumma, &
02143                             z__[z_offset], ldu, &work[1], &m, &result[32]);
02144 
02145 /*           Call SSTEMR to compute D2, do tests. */
02146 
02147 /*           Compute D2 */
02148 
02149                     scopy_(&n, &sd[1], &c__1, &d5[1], &c__1);
02150                     if (n > 0) {
02151                         i__3 = n - 1;
02152                         scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
02153                     }
02154 
02155                     ntest = 34;
02156                     i__3 = *lwork - n;
02157                     i__4 = *liwork - (n << 1);
02158                     sstemr_("N", "V", &n, &d5[1], &work[1], &vl, &vu, &il, &
02159                             iu, &m, &d2[1], &z__[z_offset], ldu, &n, &iwork[1]
02160 , &tryrac, &work[n + 1], &i__3, &iwork[(n << 1) + 
02161                             1], &i__4, &iinfo);
02162                     if (iinfo != 0) {
02163                         io___83.ciunit = *nounit;
02164                         s_wsfe(&io___83);
02165                         do_fio(&c__1, "SSTEMR(N,V)", (ftnlen)11);
02166                         do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
02167                                 ;
02168                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
02169                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
02170                                 ;
02171                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
02172                                 integer));
02173                         e_wsfe();
02174                         *info = abs(iinfo);
02175                         if (iinfo < 0) {
02176                             return 0;
02177                         } else {
02178                             result[34] = ulpinv;
02179                             goto L280;
02180                         }
02181                     }
02182 
02183 /*           Do Test 34 */
02184 
02185                     temp1 = 0.f;
02186                     temp2 = 0.f;
02187 
02188                     i__3 = iu - il + 1;
02189                     for (j = 1; j <= i__3; ++j) {
02190 /* Computing MAX */
02191                         r__3 = temp1, r__4 = (r__1 = d1[j], dabs(r__1)), r__3 
02192                                 = max(r__3,r__4), r__4 = (r__2 = d2[j], dabs(
02193                                 r__2));
02194                         temp1 = dmax(r__3,r__4);
02195 /* Computing MAX */
02196                         r__2 = temp2, r__3 = (r__1 = d1[j] - d2[j], dabs(r__1)
02197                                 );
02198                         temp2 = dmax(r__2,r__3);
02199 /* L250: */
02200                     }
02201 
02202 /* Computing MAX */
02203                     r__1 = unfl, r__2 = ulp * dmax(temp1,temp2);
02204                     result[34] = temp2 / dmax(r__1,r__2);
02205                 } else {
02206                     result[29] = 0.f;
02207                     result[30] = 0.f;
02208                     result[31] = 0.f;
02209                     result[32] = 0.f;
02210                     result[33] = 0.f;
02211                     result[34] = 0.f;
02212                 }
02213 
02214 
02215 /*           Call SSTEMR(V,A) to compute D1 and Z, do tests. */
02216 
02217 /*           Compute D1 and Z */
02218 
02219                 scopy_(&n, &sd[1], &c__1, &d5[1], &c__1);
02220                 if (n > 0) {
02221                     i__3 = n - 1;
02222                     scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
02223                 }
02224 
02225                 ntest = 35;
02226 
02227                 i__3 = *lwork - n;
02228                 i__4 = *liwork - (n << 1);
02229                 sstemr_("V", "A", &n, &d5[1], &work[1], &vl, &vu, &il, &iu, &
02230                         m, &d1[1], &z__[z_offset], ldu, &n, &iwork[1], &
02231                         tryrac, &work[n + 1], &i__3, &iwork[(n << 1) + 1], &
02232                         i__4, &iinfo);
02233                 if (iinfo != 0) {
02234                     io___84.ciunit = *nounit;
02235                     s_wsfe(&io___84);
02236                     do_fio(&c__1, "SSTEMR(V,A)", (ftnlen)11);
02237                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
02238                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
02239                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
02240                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
02241                             ;
02242                     e_wsfe();
02243                     *info = abs(iinfo);
02244                     if (iinfo < 0) {
02245                         return 0;
02246                     } else {
02247                         result[35] = ulpinv;
02248                         goto L280;
02249                     }
02250                 }
02251 
02252 /*           Do Tests 35 and 36 */
02253 
02254                 sstt22_(&n, &m, &c__0, &sd[1], &se[1], &d1[1], dumma, &z__[
02255                         z_offset], ldu, &work[1], &m, &result[35]);
02256 
02257 /*           Call SSTEMR to compute D2, do tests. */
02258 
02259 /*           Compute D2 */
02260 
02261                 scopy_(&n, &sd[1], &c__1, &d5[1], &c__1);
02262                 if (n > 0) {
02263                     i__3 = n - 1;
02264                     scopy_(&i__3, &se[1], &c__1, &work[1], &c__1);
02265                 }
02266 
02267                 ntest = 37;
02268                 i__3 = *lwork - n;
02269                 i__4 = *liwork - (n << 1);
02270                 sstemr_("N", "A", &n, &d5[1], &work[1], &vl, &vu, &il, &iu, &
02271                         m, &d2[1], &z__[z_offset], ldu, &n, &iwork[1], &
02272                         tryrac, &work[n + 1], &i__3, &iwork[(n << 1) + 1], &
02273                         i__4, &iinfo);
02274                 if (iinfo != 0) {
02275                     io___85.ciunit = *nounit;
02276                     s_wsfe(&io___85);
02277                     do_fio(&c__1, "SSTEMR(N,A)", (ftnlen)11);
02278                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
02279                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
02280                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
02281                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
02282                             ;
02283                     e_wsfe();
02284                     *info = abs(iinfo);
02285                     if (iinfo < 0) {
02286                         return 0;
02287                     } else {
02288                         result[37] = ulpinv;
02289                         goto L280;
02290                     }
02291                 }
02292 
02293 /*           Do Test 34 */
02294 
02295                 temp1 = 0.f;
02296                 temp2 = 0.f;
02297 
02298                 i__3 = n;
02299                 for (j = 1; j <= i__3; ++j) {
02300 /* Computing MAX */
02301                     r__3 = temp1, r__4 = (r__1 = d1[j], dabs(r__1)), r__3 = 
02302                             max(r__3,r__4), r__4 = (r__2 = d2[j], dabs(r__2));
02303                     temp1 = dmax(r__3,r__4);
02304 /* Computing MAX */
02305                     r__2 = temp2, r__3 = (r__1 = d1[j] - d2[j], dabs(r__1));
02306                     temp2 = dmax(r__2,r__3);
02307 /* L260: */
02308                 }
02309 
02310 /* Computing MAX */
02311                 r__1 = unfl, r__2 = ulp * dmax(temp1,temp2);
02312                 result[37] = temp2 / dmax(r__1,r__2);
02313             }
02314 L270:
02315 L280:
02316             ntestt += ntest;
02317 
02318 /*           End of Loop -- Check for RESULT(j) > THRESH */
02319 
02320 
02321 /*           Print out tests which fail. */
02322 
02323             i__3 = ntest;
02324             for (jr = 1; jr <= i__3; ++jr) {
02325                 if (result[jr] >= *thresh) {
02326 
02327 /*                 If this is the first test to fail, */
02328 /*                 print a header to the data file. */
02329 
02330                     if (nerrs == 0) {
02331                         io___86.ciunit = *nounit;
02332                         s_wsfe(&io___86);
02333                         do_fio(&c__1, "SST", (ftnlen)3);
02334                         e_wsfe();
02335                         io___87.ciunit = *nounit;
02336                         s_wsfe(&io___87);
02337                         e_wsfe();
02338                         io___88.ciunit = *nounit;
02339                         s_wsfe(&io___88);
02340                         e_wsfe();
02341                         io___89.ciunit = *nounit;
02342                         s_wsfe(&io___89);
02343                         do_fio(&c__1, "Symmetric", (ftnlen)9);
02344                         e_wsfe();
02345                         io___90.ciunit = *nounit;
02346                         s_wsfe(&io___90);
02347                         e_wsfe();
02348 
02349 /*                    Tests performed */
02350 
02351                         io___91.ciunit = *nounit;
02352                         s_wsfe(&io___91);
02353                         e_wsfe();
02354                     }
02355                     ++nerrs;
02356                     io___92.ciunit = *nounit;
02357                     s_wsfe(&io___92);
02358                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
02359                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
02360                             ;
02361                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
02362                     do_fio(&c__1, (char *)&jr, (ftnlen)sizeof(integer));
02363                     do_fio(&c__1, (char *)&result[jr], (ftnlen)sizeof(real));
02364                     e_wsfe();
02365                 }
02366 /* L290: */
02367             }
02368 L300:
02369             ;
02370         }
02371 /* L310: */
02372     }
02373 
02374 /*     Summary */
02375 
02376     slasum_("SST", nounit, &nerrs, &ntestt);
02377     return 0;
02378 
02379 
02380 
02381 
02382 /* L9993: */
02383 /* L9992: */
02384 /* L9991: */
02385 /* L9989: */
02386 
02387 /*     End of SCHKST */
02388 
02389 } /* schkst_ */


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