sdrgsx.c
Go to the documentation of this file.
00001 /* sdrgsx.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Common Block Declarations */
00017 
00018 struct {
00019     integer m, n, mplusn, k;
00020     logical fs;
00021 } mn_;
00022 
00023 #define mn_1 mn_
00024 
00025 /* Table of constant values */
00026 
00027 static integer c__1 = 1;
00028 static integer c__0 = 0;
00029 static integer c_n1 = -1;
00030 static real c_b26 = 0.f;
00031 static integer c__3 = 3;
00032 static integer c__4 = 4;
00033 
00034 /* Subroutine */ int sdrgsx_(integer *nsize, integer *ncmax, real *thresh, 
00035         integer *nin, integer *nout, real *a, integer *lda, real *b, real *ai, 
00036          real *bi, real *z__, real *q, real *alphar, real *alphai, real *beta, 
00037          real *c__, integer *ldc, real *s, real *work, integer *lwork, 
00038         integer *iwork, integer *liwork, logical *bwork, integer *info)
00039 {
00040     /* Format strings */
00041     static char fmt_9999[] = "(\002 SDRGSX: \002,a,\002 returned INFO=\002,i"
00042             "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002)\002)";
00043     static char fmt_9997[] = "(\002 SDRGSX: SGET53 returned INFO=\002,i1,"
00044             "\002 for eigenvalue \002,i6,\002.\002,/9x,\002N=\002,i6,\002, JT"
00045             "YPE=\002,i6,\002)\002)";
00046     static char fmt_9996[] = "(\002 SDRGSX: S not in Schur form at eigenvalu"
00047             "e \002,i6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002"
00048             ")\002)";
00049     static char fmt_9995[] = "(/1x,a3,\002 -- Real Expert Generalized Schur "
00050             "form\002,\002 problem driver\002)";
00051     static char fmt_9993[] = "(\002 Matrix types: \002,/\002  1:  A is a blo"
00052             "ck diagonal matrix of Jordan blocks \002,\002and B is the identi"
00053             "ty \002,/\002      matrix, \002,/\002  2:  A and B are upper tri"
00054             "angular matrices, \002,/\002  3:  A and B are as type 2, but eac"
00055             "h second diagonal \002,\002block in A_11 and \002,/\002      eac"
00056             "h third diaongal block in A_22 are 2x2 blocks,\002,/\002  4:  A "
00057             "and B are block diagonal matrices, \002,/\002  5:  (A,B) has pot"
00058             "entially close or common \002,\002eigenvalues.\002,/)";
00059     static char fmt_9992[] = "(/\002 Tests performed:  (S is Schur, T is tri"
00060             "angular, \002,\002Q and Z are \002,a,\002,\002,/19x,\002 a is al"
00061             "pha, b is beta, and \002,a,\002 means \002,a,\002.)\002,/\002  1"
00062             " = | A - Q S Z\002,a,\002 | / ( |A| n ulp )      2 = | B - Q T "
00063             "Z\002,a,\002 | / ( |B| n ulp )\002,/\002  3 = | I - QQ\002,a,"
00064             "\002 | / ( n ulp )             4 = | I - ZZ\002,a,\002 | / ( n u"
00065             "lp )\002,/\002  5 = 1/ULP  if A is not in \002,\002Schur form "
00066             "S\002,/\002  6 = difference between (alpha,beta)\002,\002 and di"
00067             "agonals of (S,T)\002,/\002  7 = 1/ULP  if SDIM is not the correc"
00068             "t number of \002,\002selected eigenvalues\002,/\002  8 = 1/ULP  "
00069             "if DIFEST/DIFTRU > 10*THRESH or \002,\002DIFTRU/DIFEST > 10*THRE"
00070             "SH\002,/\002  9 = 1/ULP  if DIFEST <> 0 or DIFTRU > ULP*norm(A,B"
00071             ") \002,\002when reordering fails\002,/\002 10 = 1/ULP  if PLEST/"
00072             "PLTRU > THRESH or \002,\002PLTRU/PLEST > THRESH\002,/\002    ( T"
00073             "est 10 is only for input examples )\002,/)";
00074     static char fmt_9991[] = "(\002 Matrix order=\002,i2,\002, type=\002,i2"
00075             ",\002, a=\002,e10.4,\002, order(A_11)=\002,i2,\002, result \002,"
00076             "i2,\002 is \002,0p,f8.2)";
00077     static char fmt_9990[] = "(\002 Matrix order=\002,i2,\002, type=\002,i2"
00078             ",\002, a=\002,e10.4,\002, order(A_11)=\002,i2,\002, result \002,"
00079             "i2,\002 is \002,0p,e10.4)";
00080     static char fmt_9998[] = "(\002 SDRGSX: \002,a,\002 returned INFO=\002,i"
00081             "6,\002.\002,/9x,\002N=\002,i6,\002, Input Example #\002,i2,\002"
00082             ")\002)";
00083     static char fmt_9994[] = "(\002Input Example\002)";
00084     static char fmt_9989[] = "(\002 Input example #\002,i2,\002, matrix orde"
00085             "r=\002,i4,\002,\002,\002 result \002,i2,\002 is\002,0p,f8.2)";
00086     static char fmt_9988[] = "(\002 Input example #\002,i2,\002, matrix orde"
00087             "r=\002,i4,\002,\002,\002 result \002,i2,\002 is\002,1p,e10.3)";
00088 
00089     /* System generated locals */
00090     integer a_dim1, a_offset, ai_dim1, ai_offset, b_dim1, b_offset, bi_dim1, 
00091             bi_offset, c_dim1, c_offset, q_dim1, q_offset, z_dim1, z_offset, 
00092             i__1, i__2, i__3, i__4;
00093     real r__1, r__2, r__3, r__4, r__5, r__6, r__7, r__8, r__9, r__10;
00094 
00095     /* Builtin functions */
00096     double sqrt(doublereal);
00097     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
00098              s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
00099             e_rsle(void);
00100 
00101     /* Local variables */
00102     integer i__, j, i1, mm;
00103     real pl[2];
00104     integer mn2, qba, qbb;
00105     real ulp, temp1, temp2, abnrm;
00106     integer ifunc, iinfo, linfo;
00107     extern /* Subroutine */ int sget51_(integer *, integer *, real *, integer 
00108             *, real *, integer *, real *, integer *, real *, integer *, real *
00109 , real *), sget53_(real *, integer *, real *, integer *, real *, 
00110             real *, real *, real *, integer *);
00111     char sense[1];
00112     integer nerrs, ntest;
00113     real pltru;
00114     extern /* Subroutine */ int slakf2_(integer *, integer *, real *, integer 
00115             *, real *, real *, real *, real *, integer *), slatm5_(integer *, 
00116             integer *, integer *, real *, integer *, real *, integer *, real *
00117 , integer *, real *, integer *, real *, integer *, real *, 
00118             integer *, real *, integer *, real *, integer *, real *, integer *
00119 , integer *);
00120     real thrsh2;
00121     logical ilabad;
00122     extern /* Subroutine */ int slabad_(real *, real *);
00123     integer bdspac;
00124     extern doublereal slamch_(char *), slange_(char *, integer *, 
00125             integer *, real *, integer *, real *);
00126     extern /* Subroutine */ int xerbla_(char *, integer *);
00127     real difest[2];
00128     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00129             integer *, integer *);
00130     real bignum;
00131     extern /* Subroutine */ int alasvm_(char *, integer *, integer *, integer 
00132             *, integer *);
00133     real weight;
00134     extern /* Subroutine */ int sgesvd_(char *, char *, integer *, integer *, 
00135             real *, integer *, real *, real *, integer *, real *, integer *, 
00136             real *, integer *, integer *), slacpy_(char *, 
00137             integer *, integer *, real *, integer *, real *, integer *);
00138     real diftru;
00139     extern /* Subroutine */ int slaset_(char *, integer *, integer *, real *, 
00140             real *, real *, integer *), sggesx_(char *, char *, char *
00141 , L_fp, char *, integer *, real *, integer *, real *, integer *, 
00142             integer *, real *, real *, real *, real *, integer *, real *, 
00143             integer *, real *, real *, real *, integer *, integer *, integer *
00144 , logical *, integer *);
00145     integer minwrk, maxwrk;
00146     real smlnum, ulpinv;
00147     integer nptknt;
00148     real result[10];
00149     extern logical slctsx_();
00150     integer ntestt, prtype;
00151 
00152     /* Fortran I/O blocks */
00153     static cilist io___22 = { 0, 0, 0, fmt_9999, 0 };
00154     static cilist io___31 = { 0, 0, 0, fmt_9997, 0 };
00155     static cilist io___32 = { 0, 0, 0, fmt_9996, 0 };
00156     static cilist io___35 = { 0, 0, 0, fmt_9995, 0 };
00157     static cilist io___36 = { 0, 0, 0, fmt_9993, 0 };
00158     static cilist io___37 = { 0, 0, 0, fmt_9992, 0 };
00159     static cilist io___39 = { 0, 0, 0, fmt_9991, 0 };
00160     static cilist io___40 = { 0, 0, 0, fmt_9990, 0 };
00161     static cilist io___42 = { 0, 0, 1, 0, 0 };
00162     static cilist io___43 = { 0, 0, 1, 0, 0 };
00163     static cilist io___44 = { 0, 0, 0, 0, 0 };
00164     static cilist io___45 = { 0, 0, 0, 0, 0 };
00165     static cilist io___46 = { 0, 0, 0, 0, 0 };
00166     static cilist io___48 = { 0, 0, 0, fmt_9998, 0 };
00167     static cilist io___49 = { 0, 0, 0, fmt_9997, 0 };
00168     static cilist io___50 = { 0, 0, 0, fmt_9996, 0 };
00169     static cilist io___51 = { 0, 0, 0, fmt_9995, 0 };
00170     static cilist io___52 = { 0, 0, 0, fmt_9994, 0 };
00171     static cilist io___53 = { 0, 0, 0, fmt_9992, 0 };
00172     static cilist io___54 = { 0, 0, 0, fmt_9989, 0 };
00173     static cilist io___55 = { 0, 0, 0, fmt_9988, 0 };
00174 
00175 
00176 
00177 /*  -- LAPACK test routine (version 3.1) -- */
00178 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00179 /*     November 2006 */
00180 
00181 /*     .. Scalar Arguments .. */
00182 /*     .. */
00183 /*     .. Array Arguments .. */
00184 /*     .. */
00185 
00186 /*  Purpose */
00187 /*  ======= */
00188 
00189 /*  SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form) */
00190 /*  problem expert driver SGGESX. */
00191 
00192 /*  SGGESX factors A and B as Q S Z' and Q T Z', where ' means */
00193 /*  transpose, T is upper triangular, S is in generalized Schur form */
00194 /*  (block upper triangular, with 1x1 and 2x2 blocks on the diagonal, */
00195 /*  the 2x2 blocks corresponding to complex conjugate pairs of */
00196 /*  generalized eigenvalues), and Q and Z are orthogonal.  It also */
00197 /*  computes the generalized eigenvalues (alpha(1),beta(1)), ..., */
00198 /*  (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the */
00199 /*  characteristic equation */
00200 
00201 /*      det( A - w(j) B ) = 0 */
00202 
00203 /*  Optionally it also reorders the eigenvalues so that a selected */
00204 /*  cluster of eigenvalues appears in the leading diagonal block of the */
00205 /*  Schur forms; computes a reciprocal condition number for the average */
00206 /*  of the selected eigenvalues; and computes a reciprocal condition */
00207 /*  number for the right and left deflating subspaces corresponding to */
00208 /*  the selected eigenvalues. */
00209 
00210 /*  When SDRGSX is called with NSIZE > 0, five (5) types of built-in */
00211 /*  matrix pairs are used to test the routine SGGESX. */
00212 
00213 /*  When SDRGSX is called with NSIZE = 0, it reads in test matrix data */
00214 /*  to test SGGESX. */
00215 
00216 /*  For each matrix pair, the following tests will be performed and */
00217 /*  compared with the threshhold THRESH except for the tests (7) and (9): */
00218 
00219 /*  (1)   | A - Q S Z' | / ( |A| n ulp ) */
00220 
00221 /*  (2)   | B - Q T Z' | / ( |B| n ulp ) */
00222 
00223 /*  (3)   | I - QQ' | / ( n ulp ) */
00224 
00225 /*  (4)   | I - ZZ' | / ( n ulp ) */
00226 
00227 /*  (5)   if A is in Schur form (i.e. quasi-triangular form) */
00228 
00229 /*  (6)   maximum over j of D(j)  where: */
00230 
00231 /*        if alpha(j) is real: */
00232 /*                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)| */
00233 /*            D(j) = ------------------------ + ----------------------- */
00234 /*                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|) */
00235 
00236 /*        if alpha(j) is complex: */
00237 /*                                  | det( s S - w T ) | */
00238 /*            D(j) = --------------------------------------------------- */
00239 /*                   ulp max( s norm(S), |w| norm(T) )*norm( s S - w T ) */
00240 
00241 /*            and S and T are here the 2 x 2 diagonal blocks of S and T */
00242 /*            corresponding to the j-th and j+1-th eigenvalues. */
00243 
00244 /*  (7)   if sorting worked and SDIM is the number of eigenvalues */
00245 /*        which were selected. */
00246 
00247 /*  (8)   the estimated value DIF does not differ from the true values of */
00248 /*        Difu and Difl more than a factor 10*THRESH. If the estimate DIF */
00249 /*        equals zero the corresponding true values of Difu and Difl */
00250 /*        should be less than EPS*norm(A, B). If the true value of Difu */
00251 /*        and Difl equal zero, the estimate DIF should be less than */
00252 /*        EPS*norm(A, B). */
00253 
00254 /*  (9)   If INFO = N+3 is returned by SGGESX, the reordering "failed" */
00255 /*        and we check that DIF = PL = PR = 0 and that the true value of */
00256 /*        Difu and Difl is < EPS*norm(A, B). We count the events when */
00257 /*        INFO=N+3. */
00258 
00259 /*  For read-in test matrices, the above tests are run except that the */
00260 /*  exact value for DIF (and PL) is input data.  Additionally, there is */
00261 /*  one more test run for read-in test matrices: */
00262 
00263 /*  (10)  the estimated value PL does not differ from the true value of */
00264 /*        PLTRU more than a factor THRESH. If the estimate PL equals */
00265 /*        zero the corresponding true value of PLTRU should be less than */
00266 /*        EPS*norm(A, B). If the true value of PLTRU equal zero, the */
00267 /*        estimate PL should be less than EPS*norm(A, B). */
00268 
00269 /*  Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1) */
00270 /*  matrix pairs are generated and tested. NSIZE should be kept small. */
00271 
00272 /*  SVD (routine SGESVD) is used for computing the true value of DIF_u */
00273 /*  and DIF_l when testing the built-in test problems. */
00274 
00275 /*  Built-in Test Matrices */
00276 /*  ====================== */
00277 
00278 /*  All built-in test matrices are the 2 by 2 block of triangular */
00279 /*  matrices */
00280 
00281 /*           A = [ A11 A12 ]    and      B = [ B11 B12 ] */
00282 /*               [     A22 ]                 [     B22 ] */
00283 
00284 /*  where for different type of A11 and A22 are given as the following. */
00285 /*  A12 and B12 are chosen so that the generalized Sylvester equation */
00286 
00287 /*           A11*R - L*A22 = -A12 */
00288 /*           B11*R - L*B22 = -B12 */
00289 
00290 /*  have prescribed solution R and L. */
00291 
00292 /*  Type 1:  A11 = J_m(1,-1) and A_22 = J_k(1-a,1). */
00293 /*           B11 = I_m, B22 = I_k */
00294 /*           where J_k(a,b) is the k-by-k Jordan block with ``a'' on */
00295 /*           diagonal and ``b'' on superdiagonal. */
00296 
00297 /*  Type 2:  A11 = (a_ij) = ( 2(.5-sin(i)) ) and */
00298 /*           B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m */
00299 /*           A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and */
00300 /*           B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k */
00301 
00302 /*  Type 3:  A11, A22 and B11, B22 are chosen as for Type 2, but each */
00303 /*           second diagonal block in A_11 and each third diagonal block */
00304 /*           in A_22 are made as 2 by 2 blocks. */
00305 
00306 /*  Type 4:  A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) ) */
00307 /*              for i=1,...,m,  j=1,...,m and */
00308 /*           A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) ) */
00309 /*              for i=m+1,...,k,  j=m+1,...,k */
00310 
00311 /*  Type 5:  (A,B) and have potentially close or common eigenvalues and */
00312 /*           very large departure from block diagonality A_11 is chosen */
00313 /*           as the m x m leading submatrix of A_1: */
00314 /*                   |  1  b                            | */
00315 /*                   | -b  1                            | */
00316 /*                   |        1+d  b                    | */
00317 /*                   |         -b 1+d                   | */
00318 /*            A_1 =  |                  d  1            | */
00319 /*                   |                 -1  d            | */
00320 /*                   |                        -d  1     | */
00321 /*                   |                        -1 -d     | */
00322 /*                   |                               1  | */
00323 /*           and A_22 is chosen as the k x k leading submatrix of A_2: */
00324 /*                   | -1  b                            | */
00325 /*                   | -b -1                            | */
00326 /*                   |       1-d  b                     | */
00327 /*                   |       -b  1-d                    | */
00328 /*            A_2 =  |                 d 1+b            | */
00329 /*                   |               -1-b d             | */
00330 /*                   |                       -d  1+b    | */
00331 /*                   |                      -1+b  -d    | */
00332 /*                   |                              1-d | */
00333 /*           and matrix B are chosen as identity matrices (see SLATM5). */
00334 
00335 
00336 /*  Arguments */
00337 /*  ========= */
00338 
00339 /*  NSIZE   (input) INTEGER */
00340 /*          The maximum size of the matrices to use. NSIZE >= 0. */
00341 /*          If NSIZE = 0, no built-in tests matrices are used, but */
00342 /*          read-in test matrices are used to test SGGESX. */
00343 
00344 /*  NCMAX   (input) INTEGER */
00345 /*          Maximum allowable NMAX for generating Kroneker matrix */
00346 /*          in call to SLAKF2 */
00347 
00348 /*  THRESH  (input) REAL */
00349 /*          A test will count as "failed" if the "error", computed as */
00350 /*          described above, exceeds THRESH.  Note that the error */
00351 /*          is scaled to be O(1), so THRESH should be a reasonably */
00352 /*          small multiple of 1, e.g., 10 or 100.  In particular, */
00353 /*          it should not depend on the precision (single vs. double) */
00354 /*          or the size of the matrix.  THRESH >= 0. */
00355 
00356 /*  NIN     (input) INTEGER */
00357 /*          The FORTRAN unit number for reading in the data file of */
00358 /*          problems to solve. */
00359 
00360 /*  NOUT    (input) INTEGER */
00361 /*          The FORTRAN unit number for printing out error messages */
00362 /*          (e.g., if a routine returns IINFO not equal to 0.) */
00363 
00364 /*  A       (workspace) REAL array, dimension (LDA, NSIZE) */
00365 /*          Used to store the matrix whose eigenvalues are to be */
00366 /*          computed.  On exit, A contains the last matrix actually used. */
00367 
00368 /*  LDA     (input) INTEGER */
00369 /*          The leading dimension of A, B, AI, BI, Z and Q, */
00370 /*          LDA >= max( 1, NSIZE ). For the read-in test, */
00371 /*          LDA >= max( 1, N ), N is the size of the test matrices. */
00372 
00373 /*  B       (workspace) REAL array, dimension (LDA, NSIZE) */
00374 /*          Used to store the matrix whose eigenvalues are to be */
00375 /*          computed.  On exit, B contains the last matrix actually used. */
00376 
00377 /*  AI      (workspace) REAL array, dimension (LDA, NSIZE) */
00378 /*          Copy of A, modified by SGGESX. */
00379 
00380 /*  BI      (workspace) REAL array, dimension (LDA, NSIZE) */
00381 /*          Copy of B, modified by SGGESX. */
00382 
00383 /*  Z       (workspace) REAL array, dimension (LDA, NSIZE) */
00384 /*          Z holds the left Schur vectors computed by SGGESX. */
00385 
00386 /*  Q       (workspace) REAL array, dimension (LDA, NSIZE) */
00387 /*          Q holds the right Schur vectors computed by SGGESX. */
00388 
00389 /*  ALPHAR  (workspace) REAL array, dimension (NSIZE) */
00390 /*  ALPHAI  (workspace) REAL array, dimension (NSIZE) */
00391 /*  BETA    (workspace) REAL array, dimension (NSIZE) */
00392 /*          On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues. */
00393 
00394 /*  C       (workspace) REAL array, dimension (LDC, LDC) */
00395 /*          Store the matrix generated by subroutine SLAKF2, this is the */
00396 /*          matrix formed by Kronecker products used for estimating */
00397 /*          DIF. */
00398 
00399 /*  LDC     (input) INTEGER */
00400 /*          The leading dimension of C. LDC >= max(1, LDA*LDA/2 ). */
00401 
00402 /*  S       (workspace) REAL array, dimension (LDC) */
00403 /*          Singular values of C */
00404 
00405 /*  WORK    (workspace) REAL array, dimension (LWORK) */
00406 
00407 /*  LWORK   (input) INTEGER */
00408 /*          The dimension of the array WORK. */
00409 /*          LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) ) */
00410 
00411 /*  IWORK   (workspace) INTEGER array, dimension (LIWORK) */
00412 
00413 /*  LIWORK  (input) INTEGER */
00414 /*          The dimension of the array IWORK. LIWORK >= NSIZE + 6. */
00415 
00416 /*  BWORK   (workspace) LOGICAL array, dimension (LDA) */
00417 
00418 /*  INFO    (output) INTEGER */
00419 /*          = 0:  successful exit */
00420 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00421 /*          > 0:  A routine returned an error code. */
00422 
00423 /*  ===================================================================== */
00424 
00425 /*     .. Parameters .. */
00426 /*     .. */
00427 /*     .. Local Scalars .. */
00428 /*     .. */
00429 /*     .. Local Arrays .. */
00430 /*     .. */
00431 /*     .. External Functions .. */
00432 /*     .. */
00433 /*     .. External Subroutines .. */
00434 /*     .. */
00435 /*     .. Intrinsic Functions .. */
00436 /*     .. */
00437 /*     .. Scalars in Common .. */
00438 /*     .. */
00439 /*     .. Common blocks .. */
00440 /*     .. */
00441 /*     .. Executable Statements .. */
00442 
00443 /*     Check for errors */
00444 
00445     /* Parameter adjustments */
00446     q_dim1 = *lda;
00447     q_offset = 1 + q_dim1;
00448     q -= q_offset;
00449     z_dim1 = *lda;
00450     z_offset = 1 + z_dim1;
00451     z__ -= z_offset;
00452     bi_dim1 = *lda;
00453     bi_offset = 1 + bi_dim1;
00454     bi -= bi_offset;
00455     ai_dim1 = *lda;
00456     ai_offset = 1 + ai_dim1;
00457     ai -= ai_offset;
00458     b_dim1 = *lda;
00459     b_offset = 1 + b_dim1;
00460     b -= b_offset;
00461     a_dim1 = *lda;
00462     a_offset = 1 + a_dim1;
00463     a -= a_offset;
00464     --alphar;
00465     --alphai;
00466     --beta;
00467     c_dim1 = *ldc;
00468     c_offset = 1 + c_dim1;
00469     c__ -= c_offset;
00470     --s;
00471     --work;
00472     --iwork;
00473     --bwork;
00474 
00475     /* Function Body */
00476     if (*nsize < 0) {
00477         *info = -1;
00478     } else if (*thresh < 0.f) {
00479         *info = -2;
00480     } else if (*nin <= 0) {
00481         *info = -3;
00482     } else if (*nout <= 0) {
00483         *info = -4;
00484     } else if (*lda < 1 || *lda < *nsize) {
00485         *info = -6;
00486     } else if (*ldc < 1 || *ldc < *nsize * *nsize / 2) {
00487         *info = -17;
00488     } else if (*liwork < *nsize + 6) {
00489         *info = -21;
00490     }
00491 
00492 /*     Compute workspace */
00493 /*      (Note: Comments in the code beginning "Workspace:" describe the */
00494 /*       minimal amount of workspace needed at that point in the code, */
00495 /*       as well as the preferred amount for good performance. */
00496 /*       NB refers to the optimal block size for the immediately */
00497 /*       following subroutine, as returned by ILAENV.) */
00498 
00499     minwrk = 1;
00500     if (*info == 0 && *lwork >= 1) {
00501 /*        MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 ) */
00502 /* Computing MAX */
00503         i__1 = (*nsize + 1) * 10, i__2 = *nsize * 5 * *nsize / 2;
00504         minwrk = max(i__1,i__2);
00505 
00506 /*        workspace for sggesx */
00507 
00508         maxwrk = (*nsize + 1) * 9 + *nsize * ilaenv_(&c__1, "SGEQRF", " ", 
00509                 nsize, &c__1, nsize, &c__0);
00510 /* Computing MAX */
00511         i__1 = maxwrk, i__2 = (*nsize + 1) * 9 + *nsize * ilaenv_(&c__1, 
00512                 "SORGQR", " ", nsize, &c__1, nsize, &c_n1);
00513         maxwrk = max(i__1,i__2);
00514 
00515 /*        workspace for sgesvd */
00516 
00517         bdspac = *nsize * 5 * *nsize / 2;
00518 /* Computing MAX */
00519         i__3 = *nsize * *nsize / 2;
00520         i__4 = *nsize * *nsize / 2;
00521         i__1 = maxwrk, i__2 = *nsize * 3 * *nsize / 2 + *nsize * *nsize * 
00522                 ilaenv_(&c__1, "SGEBRD", " ", &i__3, &i__4, &c_n1, &c_n1);
00523         maxwrk = max(i__1,i__2);
00524         maxwrk = max(maxwrk,bdspac);
00525 
00526         maxwrk = max(maxwrk,minwrk);
00527 
00528         work[1] = (real) maxwrk;
00529     }
00530 
00531     if (*lwork < minwrk) {
00532         *info = -19;
00533     }
00534 
00535     if (*info != 0) {
00536         i__1 = -(*info);
00537         xerbla_("SDRGSX", &i__1);
00538         return 0;
00539     }
00540 
00541 /*     Important constants */
00542 
00543     ulp = slamch_("P");
00544     ulpinv = 1.f / ulp;
00545     smlnum = slamch_("S") / ulp;
00546     bignum = 1.f / smlnum;
00547     slabad_(&smlnum, &bignum);
00548     thrsh2 = *thresh * 10.f;
00549     ntestt = 0;
00550     nerrs = 0;
00551 
00552 /*     Go to the tests for read-in matrix pairs */
00553 
00554     ifunc = 0;
00555     if (*nsize == 0) {
00556         goto L70;
00557     }
00558 
00559 /*     Test the built-in matrix pairs. */
00560 /*     Loop over different functions (IFUNC) of SGGESX, types (PRTYPE) */
00561 /*     of test matrices, different size (M+N) */
00562 
00563     prtype = 0;
00564     qba = 3;
00565     qbb = 4;
00566     weight = sqrt(ulp);
00567 
00568     for (ifunc = 0; ifunc <= 3; ++ifunc) {
00569         for (prtype = 1; prtype <= 5; ++prtype) {
00570             i__1 = *nsize - 1;
00571             for (mn_1.m = 1; mn_1.m <= i__1; ++mn_1.m) {
00572                 i__2 = *nsize - mn_1.m;
00573                 for (mn_1.n = 1; mn_1.n <= i__2; ++mn_1.n) {
00574 
00575                     weight = 1.f / weight;
00576                     mn_1.mplusn = mn_1.m + mn_1.n;
00577 
00578 /*                 Generate test matrices */
00579 
00580                     mn_1.fs = TRUE_;
00581                     mn_1.k = 0;
00582 
00583                     slaset_("Full", &mn_1.mplusn, &mn_1.mplusn, &c_b26, &
00584                             c_b26, &ai[ai_offset], lda);
00585                     slaset_("Full", &mn_1.mplusn, &mn_1.mplusn, &c_b26, &
00586                             c_b26, &bi[bi_offset], lda);
00587 
00588                     slatm5_(&prtype, &mn_1.m, &mn_1.n, &ai[ai_offset], lda, &
00589                             ai[mn_1.m + 1 + (mn_1.m + 1) * ai_dim1], lda, &ai[
00590                             (mn_1.m + 1) * ai_dim1 + 1], lda, &bi[bi_offset], 
00591                             lda, &bi[mn_1.m + 1 + (mn_1.m + 1) * bi_dim1], 
00592                             lda, &bi[(mn_1.m + 1) * bi_dim1 + 1], lda, &q[
00593                             q_offset], lda, &z__[z_offset], lda, &weight, &
00594                             qba, &qbb);
00595 
00596 /*                 Compute the Schur factorization and swapping the */
00597 /*                 m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. */
00598 /*                 Swapping is accomplished via the function SLCTSX */
00599 /*                 which is supplied below. */
00600 
00601                     if (ifunc == 0) {
00602                         *(unsigned char *)sense = 'N';
00603                     } else if (ifunc == 1) {
00604                         *(unsigned char *)sense = 'E';
00605                     } else if (ifunc == 2) {
00606                         *(unsigned char *)sense = 'V';
00607                     } else if (ifunc == 3) {
00608                         *(unsigned char *)sense = 'B';
00609                     }
00610 
00611                     slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &ai[ai_offset]
00612 , lda, &a[a_offset], lda);
00613                     slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &bi[bi_offset]
00614 , lda, &b[b_offset], lda);
00615 
00616                     sggesx_("V", "V", "S", (L_fp)slctsx_, sense, &mn_1.mplusn, 
00617                              &ai[ai_offset], lda, &bi[bi_offset], lda, &mm, &
00618                             alphar[1], &alphai[1], &beta[1], &q[q_offset], 
00619                             lda, &z__[z_offset], lda, pl, difest, &work[1], 
00620                             lwork, &iwork[1], liwork, &bwork[1], &linfo);
00621 
00622                     if (linfo != 0 && linfo != mn_1.mplusn + 2) {
00623                         result[0] = ulpinv;
00624                         io___22.ciunit = *nout;
00625                         s_wsfe(&io___22);
00626                         do_fio(&c__1, "SGGESX", (ftnlen)6);
00627                         do_fio(&c__1, (char *)&linfo, (ftnlen)sizeof(integer))
00628                                 ;
00629                         do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(
00630                                 integer));
00631                         do_fio(&c__1, (char *)&prtype, (ftnlen)sizeof(integer)
00632                                 );
00633                         e_wsfe();
00634                         *info = linfo;
00635                         goto L30;
00636                     }
00637 
00638 /*                 Compute the norm(A, B) */
00639 
00640                     slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &ai[ai_offset]
00641 , lda, &work[1], &mn_1.mplusn);
00642                     slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &bi[bi_offset]
00643 , lda, &work[mn_1.mplusn * mn_1.mplusn + 1], &
00644                             mn_1.mplusn);
00645                     i__3 = mn_1.mplusn << 1;
00646                     abnrm = slange_("Fro", &mn_1.mplusn, &i__3, &work[1], &
00647                             mn_1.mplusn, &work[1]);
00648 
00649 /*                 Do tests (1) to (4) */
00650 
00651                     sget51_(&c__1, &mn_1.mplusn, &a[a_offset], lda, &ai[
00652                             ai_offset], lda, &q[q_offset], lda, &z__[z_offset]
00653 , lda, &work[1], result);
00654                     sget51_(&c__1, &mn_1.mplusn, &b[b_offset], lda, &bi[
00655                             bi_offset], lda, &q[q_offset], lda, &z__[z_offset]
00656 , lda, &work[1], &result[1]);
00657                     sget51_(&c__3, &mn_1.mplusn, &b[b_offset], lda, &bi[
00658                             bi_offset], lda, &q[q_offset], lda, &q[q_offset], 
00659                             lda, &work[1], &result[2]);
00660                     sget51_(&c__3, &mn_1.mplusn, &b[b_offset], lda, &bi[
00661                             bi_offset], lda, &z__[z_offset], lda, &z__[
00662                             z_offset], lda, &work[1], &result[3]);
00663                     ntest = 4;
00664 
00665 /*                 Do tests (5) and (6): check Schur form of A and */
00666 /*                 compare eigenvalues with diagonals. */
00667 
00668                     temp1 = 0.f;
00669                     result[4] = 0.f;
00670                     result[5] = 0.f;
00671 
00672                     i__3 = mn_1.mplusn;
00673                     for (j = 1; j <= i__3; ++j) {
00674                         ilabad = FALSE_;
00675                         if (alphai[j] == 0.f) {
00676 /* Computing MAX */
00677                             r__7 = smlnum, r__8 = (r__2 = alphar[j], dabs(
00678                                     r__2)), r__7 = max(r__7,r__8), r__8 = (
00679                                     r__3 = ai[j + j * ai_dim1], dabs(r__3));
00680 /* Computing MAX */
00681                             r__9 = smlnum, r__10 = (r__5 = beta[j], dabs(r__5)
00682                                     ), r__9 = max(r__9,r__10), r__10 = (r__6 =
00683                                      bi[j + j * bi_dim1], dabs(r__6));
00684                             temp2 = ((r__1 = alphar[j] - ai[j + j * ai_dim1], 
00685                                     dabs(r__1)) / dmax(r__7,r__8) + (r__4 = 
00686                                     beta[j] - bi[j + j * bi_dim1], dabs(r__4))
00687                                      / dmax(r__9,r__10)) / ulp;
00688                             if (j < mn_1.mplusn) {
00689                                 if (ai[j + 1 + j * ai_dim1] != 0.f) {
00690                                     ilabad = TRUE_;
00691                                     result[4] = ulpinv;
00692                                 }
00693                             }
00694                             if (j > 1) {
00695                                 if (ai[j + (j - 1) * ai_dim1] != 0.f) {
00696                                     ilabad = TRUE_;
00697                                     result[4] = ulpinv;
00698                                 }
00699                             }
00700                         } else {
00701                             if (alphai[j] > 0.f) {
00702                                 i1 = j;
00703                             } else {
00704                                 i1 = j - 1;
00705                             }
00706                             if (i1 <= 0 || i1 >= mn_1.mplusn) {
00707                                 ilabad = TRUE_;
00708                             } else if (i1 < mn_1.mplusn - 1) {
00709                                 if (ai[i1 + 2 + (i1 + 1) * ai_dim1] != 0.f) {
00710                                     ilabad = TRUE_;
00711                                     result[4] = ulpinv;
00712                                 }
00713                             } else if (i1 > 1) {
00714                                 if (ai[i1 + (i1 - 1) * ai_dim1] != 0.f) {
00715                                     ilabad = TRUE_;
00716                                     result[4] = ulpinv;
00717                                 }
00718                             }
00719                             if (! ilabad) {
00720                                 sget53_(&ai[i1 + i1 * ai_dim1], lda, &bi[i1 + 
00721                                         i1 * bi_dim1], lda, &beta[j], &alphar[
00722                                         j], &alphai[j], &temp2, &iinfo);
00723                                 if (iinfo >= 3) {
00724                                     io___31.ciunit = *nout;
00725                                     s_wsfe(&io___31);
00726                                     do_fio(&c__1, (char *)&iinfo, (ftnlen)
00727                                             sizeof(integer));
00728                                     do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
00729                                             integer));
00730                                     do_fio(&c__1, (char *)&mn_1.mplusn, (
00731                                             ftnlen)sizeof(integer));
00732                                     do_fio(&c__1, (char *)&prtype, (ftnlen)
00733                                             sizeof(integer));
00734                                     e_wsfe();
00735                                     *info = abs(iinfo);
00736                                 }
00737                             } else {
00738                                 temp2 = ulpinv;
00739                             }
00740                         }
00741                         temp1 = dmax(temp1,temp2);
00742                         if (ilabad) {
00743                             io___32.ciunit = *nout;
00744                             s_wsfe(&io___32);
00745                             do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer))
00746                                     ;
00747                             do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)
00748                                     sizeof(integer));
00749                             do_fio(&c__1, (char *)&prtype, (ftnlen)sizeof(
00750                                     integer));
00751                             e_wsfe();
00752                         }
00753 /* L10: */
00754                     }
00755                     result[5] = temp1;
00756                     ntest += 2;
00757 
00758 /*                 Test (7) (if sorting worked) */
00759 
00760                     result[6] = 0.f;
00761                     if (linfo == mn_1.mplusn + 3) {
00762                         result[6] = ulpinv;
00763                     } else if (mm != mn_1.n) {
00764                         result[6] = ulpinv;
00765                     }
00766                     ++ntest;
00767 
00768 /*                 Test (8): compare the estimated value DIF and its */
00769 /*                 value. first, compute the exact DIF. */
00770 
00771                     result[7] = 0.f;
00772                     mn2 = mm * (mn_1.mplusn - mm) << 1;
00773                     if (ifunc >= 2 && mn2 <= *ncmax * *ncmax) {
00774 
00775 /*                    Note: for either following two causes, there are */
00776 /*                    almost same number of test cases fail the test. */
00777 
00778                         i__3 = mn_1.mplusn - mm;
00779                         slakf2_(&mm, &i__3, &ai[ai_offset], lda, &ai[mm + 1 + 
00780                                 (mm + 1) * ai_dim1], &bi[bi_offset], &bi[mm + 
00781                                 1 + (mm + 1) * bi_dim1], &c__[c_offset], ldc);
00782 
00783                         i__3 = *lwork - 2;
00784                         sgesvd_("N", "N", &mn2, &mn2, &c__[c_offset], ldc, &s[
00785                                 1], &work[1], &c__1, &work[2], &c__1, &work[3]
00786 , &i__3, info);
00787                         diftru = s[mn2];
00788 
00789                         if (difest[1] == 0.f) {
00790                             if (diftru > abnrm * ulp) {
00791                                 result[7] = ulpinv;
00792                             }
00793                         } else if (diftru == 0.f) {
00794                             if (difest[1] > abnrm * ulp) {
00795                                 result[7] = ulpinv;
00796                             }
00797                         } else if (diftru > thrsh2 * difest[1] || diftru * 
00798                                 thrsh2 < difest[1]) {
00799 /* Computing MAX */
00800                             r__1 = diftru / difest[1], r__2 = difest[1] / 
00801                                     diftru;
00802                             result[7] = dmax(r__1,r__2);
00803                         }
00804                         ++ntest;
00805                     }
00806 
00807 /*                 Test (9) */
00808 
00809                     result[8] = 0.f;
00810                     if (linfo == mn_1.mplusn + 2) {
00811                         if (diftru > abnrm * ulp) {
00812                             result[8] = ulpinv;
00813                         }
00814                         if (ifunc > 1 && difest[1] != 0.f) {
00815                             result[8] = ulpinv;
00816                         }
00817                         if (ifunc == 1 && pl[0] != 0.f) {
00818                             result[8] = ulpinv;
00819                         }
00820                         ++ntest;
00821                     }
00822 
00823                     ntestt += ntest;
00824 
00825 /*                 Print out tests which fail. */
00826 
00827                     for (j = 1; j <= 9; ++j) {
00828                         if (result[j - 1] >= *thresh) {
00829 
00830 /*                       If this is the first test to fail, */
00831 /*                       print a header to the data file. */
00832 
00833                             if (nerrs == 0) {
00834                                 io___35.ciunit = *nout;
00835                                 s_wsfe(&io___35);
00836                                 do_fio(&c__1, "SGX", (ftnlen)3);
00837                                 e_wsfe();
00838 
00839 /*                          Matrix types */
00840 
00841                                 io___36.ciunit = *nout;
00842                                 s_wsfe(&io___36);
00843                                 e_wsfe();
00844 
00845 /*                          Tests performed */
00846 
00847                                 io___37.ciunit = *nout;
00848                                 s_wsfe(&io___37);
00849                                 do_fio(&c__1, "orthogonal", (ftnlen)10);
00850                                 do_fio(&c__1, "'", (ftnlen)1);
00851                                 do_fio(&c__1, "transpose", (ftnlen)9);
00852                                 for (i__ = 1; i__ <= 4; ++i__) {
00853                                     do_fio(&c__1, "'", (ftnlen)1);
00854                                 }
00855                                 e_wsfe();
00856 
00857                             }
00858                             ++nerrs;
00859                             if (result[j - 1] < 1e4f) {
00860                                 io___39.ciunit = *nout;
00861                                 s_wsfe(&io___39);
00862                                 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)
00863                                         sizeof(integer));
00864                                 do_fio(&c__1, (char *)&prtype, (ftnlen)sizeof(
00865                                         integer));
00866                                 do_fio(&c__1, (char *)&weight, (ftnlen)sizeof(
00867                                         real));
00868                                 do_fio(&c__1, (char *)&mn_1.m, (ftnlen)sizeof(
00869                                         integer));
00870                                 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
00871                                         integer));
00872                                 do_fio(&c__1, (char *)&result[j - 1], (ftnlen)
00873                                         sizeof(real));
00874                                 e_wsfe();
00875                             } else {
00876                                 io___40.ciunit = *nout;
00877                                 s_wsfe(&io___40);
00878                                 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)
00879                                         sizeof(integer));
00880                                 do_fio(&c__1, (char *)&prtype, (ftnlen)sizeof(
00881                                         integer));
00882                                 do_fio(&c__1, (char *)&weight, (ftnlen)sizeof(
00883                                         real));
00884                                 do_fio(&c__1, (char *)&mn_1.m, (ftnlen)sizeof(
00885                                         integer));
00886                                 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(
00887                                         integer));
00888                                 do_fio(&c__1, (char *)&result[j - 1], (ftnlen)
00889                                         sizeof(real));
00890                                 e_wsfe();
00891                             }
00892                         }
00893 /* L20: */
00894                     }
00895 
00896 L30:
00897                     ;
00898                 }
00899 /* L40: */
00900             }
00901 /* L50: */
00902         }
00903 /* L60: */
00904     }
00905 
00906     goto L150;
00907 
00908 L70:
00909 
00910 /*     Read in data from file to check accuracy of condition estimation */
00911 /*     Read input data until N=0 */
00912 
00913     nptknt = 0;
00914 
00915 L80:
00916     io___42.ciunit = *nin;
00917     i__1 = s_rsle(&io___42);
00918     if (i__1 != 0) {
00919         goto L140;
00920     }
00921     i__1 = do_lio(&c__3, &c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer))
00922             ;
00923     if (i__1 != 0) {
00924         goto L140;
00925     }
00926     i__1 = e_rsle();
00927     if (i__1 != 0) {
00928         goto L140;
00929     }
00930     if (mn_1.mplusn == 0) {
00931         goto L140;
00932     }
00933     io___43.ciunit = *nin;
00934     i__1 = s_rsle(&io___43);
00935     if (i__1 != 0) {
00936         goto L140;
00937     }
00938     i__1 = do_lio(&c__3, &c__1, (char *)&mn_1.n, (ftnlen)sizeof(integer));
00939     if (i__1 != 0) {
00940         goto L140;
00941     }
00942     i__1 = e_rsle();
00943     if (i__1 != 0) {
00944         goto L140;
00945     }
00946     i__1 = mn_1.mplusn;
00947     for (i__ = 1; i__ <= i__1; ++i__) {
00948         io___44.ciunit = *nin;
00949         s_rsle(&io___44);
00950         i__2 = mn_1.mplusn;
00951         for (j = 1; j <= i__2; ++j) {
00952             do_lio(&c__4, &c__1, (char *)&ai[i__ + j * ai_dim1], (ftnlen)
00953                     sizeof(real));
00954         }
00955         e_rsle();
00956 /* L90: */
00957     }
00958     i__1 = mn_1.mplusn;
00959     for (i__ = 1; i__ <= i__1; ++i__) {
00960         io___45.ciunit = *nin;
00961         s_rsle(&io___45);
00962         i__2 = mn_1.mplusn;
00963         for (j = 1; j <= i__2; ++j) {
00964             do_lio(&c__4, &c__1, (char *)&bi[i__ + j * bi_dim1], (ftnlen)
00965                     sizeof(real));
00966         }
00967         e_rsle();
00968 /* L100: */
00969     }
00970     io___46.ciunit = *nin;
00971     s_rsle(&io___46);
00972     do_lio(&c__4, &c__1, (char *)&pltru, (ftnlen)sizeof(real));
00973     do_lio(&c__4, &c__1, (char *)&diftru, (ftnlen)sizeof(real));
00974     e_rsle();
00975 
00976     ++nptknt;
00977     mn_1.fs = TRUE_;
00978     mn_1.k = 0;
00979     mn_1.m = mn_1.mplusn - mn_1.n;
00980 
00981     slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &ai[ai_offset], lda, &a[
00982             a_offset], lda);
00983     slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &bi[bi_offset], lda, &b[
00984             b_offset], lda);
00985 
00986 /*     Compute the Schur factorization while swaping the */
00987 /*     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks. */
00988 
00989     sggesx_("V", "V", "S", (L_fp)slctsx_, "B", &mn_1.mplusn, &ai[ai_offset], 
00990             lda, &bi[bi_offset], lda, &mm, &alphar[1], &alphai[1], &beta[1], &
00991             q[q_offset], lda, &z__[z_offset], lda, pl, difest, &work[1], 
00992             lwork, &iwork[1], liwork, &bwork[1], &linfo);
00993 
00994     if (linfo != 0 && linfo != mn_1.mplusn + 2) {
00995         result[0] = ulpinv;
00996         io___48.ciunit = *nout;
00997         s_wsfe(&io___48);
00998         do_fio(&c__1, "SGGESX", (ftnlen)6);
00999         do_fio(&c__1, (char *)&linfo, (ftnlen)sizeof(integer));
01000         do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer));
01001         do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01002         e_wsfe();
01003         goto L130;
01004     }
01005 
01006 /*     Compute the norm(A, B) */
01007 /*        (should this be norm of (A,B) or (AI,BI)?) */
01008 
01009     slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &ai[ai_offset], lda, &work[1], 
01010              &mn_1.mplusn);
01011     slacpy_("Full", &mn_1.mplusn, &mn_1.mplusn, &bi[bi_offset], lda, &work[
01012             mn_1.mplusn * mn_1.mplusn + 1], &mn_1.mplusn);
01013     i__1 = mn_1.mplusn << 1;
01014     abnrm = slange_("Fro", &mn_1.mplusn, &i__1, &work[1], &mn_1.mplusn, &work[
01015             1]);
01016 
01017 /*     Do tests (1) to (4) */
01018 
01019     sget51_(&c__1, &mn_1.mplusn, &a[a_offset], lda, &ai[ai_offset], lda, &q[
01020             q_offset], lda, &z__[z_offset], lda, &work[1], result);
01021     sget51_(&c__1, &mn_1.mplusn, &b[b_offset], lda, &bi[bi_offset], lda, &q[
01022             q_offset], lda, &z__[z_offset], lda, &work[1], &result[1]);
01023     sget51_(&c__3, &mn_1.mplusn, &b[b_offset], lda, &bi[bi_offset], lda, &q[
01024             q_offset], lda, &q[q_offset], lda, &work[1], &result[2]);
01025     sget51_(&c__3, &mn_1.mplusn, &b[b_offset], lda, &bi[bi_offset], lda, &z__[
01026             z_offset], lda, &z__[z_offset], lda, &work[1], &result[3]);
01027 
01028 /*     Do tests (5) and (6): check Schur form of A and compare */
01029 /*     eigenvalues with diagonals. */
01030 
01031     ntest = 6;
01032     temp1 = 0.f;
01033     result[4] = 0.f;
01034     result[5] = 0.f;
01035 
01036     i__1 = mn_1.mplusn;
01037     for (j = 1; j <= i__1; ++j) {
01038         ilabad = FALSE_;
01039         if (alphai[j] == 0.f) {
01040 /* Computing MAX */
01041             r__7 = smlnum, r__8 = (r__2 = alphar[j], dabs(r__2)), r__7 = max(
01042                     r__7,r__8), r__8 = (r__3 = ai[j + j * ai_dim1], dabs(r__3)
01043                     );
01044 /* Computing MAX */
01045             r__9 = smlnum, r__10 = (r__5 = beta[j], dabs(r__5)), r__9 = max(
01046                     r__9,r__10), r__10 = (r__6 = bi[j + j * bi_dim1], dabs(
01047                     r__6));
01048             temp2 = ((r__1 = alphar[j] - ai[j + j * ai_dim1], dabs(r__1)) / 
01049                     dmax(r__7,r__8) + (r__4 = beta[j] - bi[j + j * bi_dim1], 
01050                     dabs(r__4)) / dmax(r__9,r__10)) / ulp;
01051             if (j < mn_1.mplusn) {
01052                 if (ai[j + 1 + j * ai_dim1] != 0.f) {
01053                     ilabad = TRUE_;
01054                     result[4] = ulpinv;
01055                 }
01056             }
01057             if (j > 1) {
01058                 if (ai[j + (j - 1) * ai_dim1] != 0.f) {
01059                     ilabad = TRUE_;
01060                     result[4] = ulpinv;
01061                 }
01062             }
01063         } else {
01064             if (alphai[j] > 0.f) {
01065                 i1 = j;
01066             } else {
01067                 i1 = j - 1;
01068             }
01069             if (i1 <= 0 || i1 >= mn_1.mplusn) {
01070                 ilabad = TRUE_;
01071             } else if (i1 < mn_1.mplusn - 1) {
01072                 if (ai[i1 + 2 + (i1 + 1) * ai_dim1] != 0.f) {
01073                     ilabad = TRUE_;
01074                     result[4] = ulpinv;
01075                 }
01076             } else if (i1 > 1) {
01077                 if (ai[i1 + (i1 - 1) * ai_dim1] != 0.f) {
01078                     ilabad = TRUE_;
01079                     result[4] = ulpinv;
01080                 }
01081             }
01082             if (! ilabad) {
01083                 sget53_(&ai[i1 + i1 * ai_dim1], lda, &bi[i1 + i1 * bi_dim1], 
01084                         lda, &beta[j], &alphar[j], &alphai[j], &temp2, &iinfo)
01085                         ;
01086                 if (iinfo >= 3) {
01087                     io___49.ciunit = *nout;
01088                     s_wsfe(&io___49);
01089                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
01090                     do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01091                     do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(
01092                             integer));
01093                     do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01094                     e_wsfe();
01095                     *info = abs(iinfo);
01096                 }
01097             } else {
01098                 temp2 = ulpinv;
01099             }
01100         }
01101         temp1 = dmax(temp1,temp2);
01102         if (ilabad) {
01103             io___50.ciunit = *nout;
01104             s_wsfe(&io___50);
01105             do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01106             do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer));
01107             do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01108             e_wsfe();
01109         }
01110 /* L110: */
01111     }
01112     result[5] = temp1;
01113 
01114 /*     Test (7) (if sorting worked)  <--------- need to be checked. */
01115 
01116     ntest = 7;
01117     result[6] = 0.f;
01118     if (linfo == mn_1.mplusn + 3) {
01119         result[6] = ulpinv;
01120     }
01121 
01122 /*     Test (8): compare the estimated value of DIF and its true value. */
01123 
01124     ntest = 8;
01125     result[7] = 0.f;
01126     if (difest[1] == 0.f) {
01127         if (diftru > abnrm * ulp) {
01128             result[7] = ulpinv;
01129         }
01130     } else if (diftru == 0.f) {
01131         if (difest[1] > abnrm * ulp) {
01132             result[7] = ulpinv;
01133         }
01134     } else if (diftru > thrsh2 * difest[1] || diftru * thrsh2 < difest[1]) {
01135 /* Computing MAX */
01136         r__1 = diftru / difest[1], r__2 = difest[1] / diftru;
01137         result[7] = dmax(r__1,r__2);
01138     }
01139 
01140 /*     Test (9) */
01141 
01142     ntest = 9;
01143     result[8] = 0.f;
01144     if (linfo == mn_1.mplusn + 2) {
01145         if (diftru > abnrm * ulp) {
01146             result[8] = ulpinv;
01147         }
01148         if (ifunc > 1 && difest[1] != 0.f) {
01149             result[8] = ulpinv;
01150         }
01151         if (ifunc == 1 && pl[0] != 0.f) {
01152             result[8] = ulpinv;
01153         }
01154     }
01155 
01156 /*     Test (10): compare the estimated value of PL and it true value. */
01157 
01158     ntest = 10;
01159     result[9] = 0.f;
01160     if (pl[0] == 0.f) {
01161         if (pltru > abnrm * ulp) {
01162             result[9] = ulpinv;
01163         }
01164     } else if (pltru == 0.f) {
01165         if (pl[0] > abnrm * ulp) {
01166             result[9] = ulpinv;
01167         }
01168     } else if (pltru > *thresh * pl[0] || pltru * *thresh < pl[0]) {
01169         result[9] = ulpinv;
01170     }
01171 
01172     ntestt += ntest;
01173 
01174 /*     Print out tests which fail. */
01175 
01176     i__1 = ntest;
01177     for (j = 1; j <= i__1; ++j) {
01178         if (result[j - 1] >= *thresh) {
01179 
01180 /*           If this is the first test to fail, */
01181 /*           print a header to the data file. */
01182 
01183             if (nerrs == 0) {
01184                 io___51.ciunit = *nout;
01185                 s_wsfe(&io___51);
01186                 do_fio(&c__1, "SGX", (ftnlen)3);
01187                 e_wsfe();
01188 
01189 /*              Matrix types */
01190 
01191                 io___52.ciunit = *nout;
01192                 s_wsfe(&io___52);
01193                 e_wsfe();
01194 
01195 /*              Tests performed */
01196 
01197                 io___53.ciunit = *nout;
01198                 s_wsfe(&io___53);
01199                 do_fio(&c__1, "orthogonal", (ftnlen)10);
01200                 do_fio(&c__1, "'", (ftnlen)1);
01201                 do_fio(&c__1, "transpose", (ftnlen)9);
01202                 for (i__ = 1; i__ <= 4; ++i__) {
01203                     do_fio(&c__1, "'", (ftnlen)1);
01204                 }
01205                 e_wsfe();
01206 
01207             }
01208             ++nerrs;
01209             if (result[j - 1] < 1e4f) {
01210                 io___54.ciunit = *nout;
01211                 s_wsfe(&io___54);
01212                 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01213                 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer));
01214                 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01215                 do_fio(&c__1, (char *)&result[j - 1], (ftnlen)sizeof(real));
01216                 e_wsfe();
01217             } else {
01218                 io___55.ciunit = *nout;
01219                 s_wsfe(&io___55);
01220                 do_fio(&c__1, (char *)&nptknt, (ftnlen)sizeof(integer));
01221                 do_fio(&c__1, (char *)&mn_1.mplusn, (ftnlen)sizeof(integer));
01222                 do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01223                 do_fio(&c__1, (char *)&result[j - 1], (ftnlen)sizeof(real));
01224                 e_wsfe();
01225             }
01226         }
01227 
01228 /* L120: */
01229     }
01230 
01231 L130:
01232     goto L80;
01233 L140:
01234 
01235 L150:
01236 
01237 /*     Summary */
01238 
01239     alasvm_("SGX", nout, &nerrs, &ntestt, &c__0);
01240 
01241     work[1] = (real) maxwrk;
01242 
01243     return 0;
01244 
01245 
01246 
01247 
01248 
01249 
01250 
01251 
01252 
01253 /*     End of SDRGSX */
01254 
01255 } /* sdrgsx_ */


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