zdrges.c
Go to the documentation of this file.
00001 /* zdrges.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 doublecomplex c_b1 = {0.,0.};
00019 static integer c__1 = 1;
00020 static integer c_n1 = -1;
00021 static integer c__2 = 2;
00022 static doublereal c_b29 = 1.;
00023 static integer c__3 = 3;
00024 static integer c__4 = 4;
00025 static integer c__0 = 0;
00026 
00027 /* Subroutine */ int zdrges_(integer *nsizes, integer *nn, integer *ntypes, 
00028         logical *dotype, integer *iseed, doublereal *thresh, integer *nounit, 
00029         doublecomplex *a, integer *lda, doublecomplex *b, doublecomplex *s, 
00030         doublecomplex *t, doublecomplex *q, integer *ldq, doublecomplex *z__, 
00031         doublecomplex *alpha, doublecomplex *beta, doublecomplex *work, 
00032         integer *lwork, doublereal *rwork, doublereal *result, logical *bwork, 
00033          integer *info)
00034 {
00035     /* Initialized data */
00036 
00037     static integer kclass[26] = { 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,
00038             2,2,2,3 };
00039     static integer kbmagn[26] = { 1,1,1,1,1,1,1,1,3,2,3,2,2,3,1,1,1,1,1,1,1,3,
00040             2,3,2,1 };
00041     static integer ktrian[26] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,
00042             1,1,1,1 };
00043     static logical lasign[26] = { FALSE_,FALSE_,FALSE_,FALSE_,FALSE_,FALSE_,
00044             TRUE_,FALSE_,TRUE_,TRUE_,FALSE_,FALSE_,TRUE_,TRUE_,TRUE_,FALSE_,
00045             TRUE_,FALSE_,FALSE_,FALSE_,TRUE_,TRUE_,TRUE_,TRUE_,TRUE_,FALSE_ };
00046     static logical lbsign[26] = { FALSE_,FALSE_,FALSE_,FALSE_,FALSE_,FALSE_,
00047             FALSE_,TRUE_,FALSE_,FALSE_,TRUE_,TRUE_,FALSE_,FALSE_,TRUE_,FALSE_,
00048             TRUE_,FALSE_,FALSE_,FALSE_,FALSE_,FALSE_,FALSE_,FALSE_,FALSE_,
00049             FALSE_ };
00050     static integer kz1[6] = { 0,1,2,1,3,3 };
00051     static integer kz2[6] = { 0,0,1,2,1,1 };
00052     static integer kadd[6] = { 0,0,0,0,3,2 };
00053     static integer katype[26] = { 0,1,0,1,2,3,4,1,4,4,1,1,4,4,4,2,4,5,8,7,9,4,
00054             4,4,4,0 };
00055     static integer kbtype[26] = { 0,0,1,1,2,-3,1,4,1,1,4,4,1,1,-4,2,-4,8,8,8,
00056             8,8,8,8,8,0 };
00057     static integer kazero[26] = { 1,1,1,1,1,1,2,1,2,2,1,1,2,2,3,1,3,5,5,5,5,3,
00058             3,3,3,1 };
00059     static integer kbzero[26] = { 1,1,1,1,1,1,1,2,1,1,2,2,1,1,4,1,4,6,6,6,6,4,
00060             4,4,4,1 };
00061     static integer kamagn[26] = { 1,1,1,1,1,1,1,1,2,3,2,3,2,3,1,1,1,1,1,1,1,2,
00062             3,3,2,1 };
00063 
00064     /* Format strings */
00065     static char fmt_9999[] = "(\002 ZDRGES: \002,a,\002 returned INFO=\002,i"
00066             "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, ISEED="
00067             "(\002,4(i4,\002,\002),i5,\002)\002)";
00068     static char fmt_9998[] = "(\002 ZDRGES: S not in Schur form at eigenvalu"
00069             "e \002,i6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, "
00070             "ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
00071     static char fmt_9997[] = "(/1x,a3,\002 -- Complex Generalized Schur from"
00072             " problem \002,\002driver\002)";
00073     static char fmt_9996[] = "(\002 Matrix types (see ZDRGES for details):"
00074             " \002)";
00075     static char fmt_9995[] = "(\002 Special Matrices:\002,23x,\002(J'=transp"
00076             "osed Jordan block)\002,/\002   1=(0,0)  2=(I,0)  3=(0,I)  4=(I,I"
00077             ")  5=(J',J')  \002,\0026=(diag(J',I), diag(I,J'))\002,/\002 Diag"
00078             "onal Matrices:  ( \002,\002D=diag(0,1,2,...) )\002,/\002   7=(D,"
00079             "I)   9=(large*D, small*I\002,\002)  11=(large*I, small*D)  13=(l"
00080             "arge*D, large*I)\002,/\002   8=(I,D)  10=(small*D, large*I)  12="
00081             "(small*I, large*D) \002,\002 14=(small*D, small*I)\002,/\002  15"
00082             "=(D, reversed D)\002)";
00083     static char fmt_9994[] = "(\002 Matrices Rotated by Random \002,a,\002 M"
00084             "atrices U, V:\002,/\002  16=Transposed Jordan Blocks            "
00085             " 19=geometric \002,\002alpha, beta=0,1\002,/\002  17=arithm. alp"
00086             "ha&beta             \002,\002      20=arithmetic alpha, beta=0,"
00087             "1\002,/\002  18=clustered \002,\002alpha, beta=0,1            21"
00088             "=random alpha, beta=0,1\002,/\002 Large & Small Matrices:\002,"
00089             "/\002  22=(large, small)   \002,\00223=(small,large)    24=(smal"
00090             "l,small)    25=(large,large)\002,/\002  26=random O(1) matrices"
00091             ".\002)";
00092     static char fmt_9993[] = "(/\002 Tests performed:  (S is Schur, T is tri"
00093             "angular, \002,\002Q and Z are \002,a,\002,\002,/19x,\002l and r "
00094             "are the appropriate left and right\002,/19x,\002eigenvectors, re"
00095             "sp., a is alpha, b is beta, and\002,/19x,a,\002 means \002,a,"
00096             "\002.)\002,/\002 Without ordering: \002,/\002  1 = | A - Q S "
00097             "Z\002,a,\002 | / ( |A| n ulp )      2 = | B - Q T Z\002,a,\002 |"
00098             " / ( |B| n ulp )\002,/\002  3 = | I - QQ\002,a,\002 | / ( n ulp "
00099             ")             4 = | I - ZZ\002,a,\002 | / ( n ulp )\002,/\002  5"
00100             " = A is in Schur form S\002,/\002  6 = difference between (alpha"
00101             ",beta)\002,\002 and diagonals of (S,T)\002,/\002 With ordering:"
00102             " \002,/\002  7 = | (A,B) - Q (S,T) Z\002,a,\002 | / ( |(A,B)| n "
00103             "ulp )\002,/\002  8 = | I - QQ\002,a,\002 | / ( n ulp )          "
00104             "   9 = | I - ZZ\002,a,\002 | / ( n ulp )\002,/\002 10 = A is in "
00105             "Schur form S\002,/\002 11 = difference between (alpha,beta) and "
00106             "diagonals\002,\002 of (S,T)\002,/\002 12 = SDIM is the correct n"
00107             "umber of \002,\002selected eigenvalues\002,/)";
00108     static char fmt_9992[] = "(\002 Matrix order=\002,i5,\002, type=\002,i2"
00109             ",\002, seed=\002,4(i4,\002,\002),\002 result \002,i2,\002 is\002"
00110             ",0p,f8.2)";
00111     static char fmt_9991[] = "(\002 Matrix order=\002,i5,\002, type=\002,i2"
00112             ",\002, seed=\002,4(i4,\002,\002),\002 result \002,i2,\002 is\002"
00113             ",1p,d10.3)";
00114 
00115     /* System generated locals */
00116     integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, s_dim1, 
00117             s_offset, t_dim1, t_offset, z_dim1, z_offset, i__1, i__2, i__3, 
00118             i__4, i__5, i__6, i__7, i__8, i__9, i__10, i__11;
00119     doublereal d__1, d__2, d__3, d__4, d__5, d__6, d__7, d__8, d__9, d__10, 
00120             d__11, d__12, d__13, d__14, d__15, d__16;
00121     doublecomplex z__1, z__2, z__3, z__4;
00122 
00123     /* Builtin functions */
00124     double d_sign(doublereal *, doublereal *), z_abs(doublecomplex *);
00125     void d_cnjg(doublecomplex *, doublecomplex *);
00126     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00127     double d_imag(doublecomplex *);
00128 
00129     /* Local variables */
00130     integer i__, j, n, n1, jc, nb, in, jr;
00131     doublereal ulp;
00132     integer iadd, sdim, nmax, rsub;
00133     char sort[1];
00134     doublereal temp1, temp2;
00135     logical badnn;
00136     integer iinfo;
00137     doublereal rmagn[4];
00138     doublecomplex ctemp;
00139     extern /* Subroutine */ int zget51_(integer *, integer *, doublecomplex *, 
00140              integer *, doublecomplex *, integer *, doublecomplex *, integer *
00141 , doublecomplex *, integer *, doublecomplex *, doublereal *, 
00142             doublereal *), zgges_(char *, char *, char *, L_fp, integer *, 
00143             doublecomplex *, integer *, doublecomplex *, integer *, integer *, 
00144              doublecomplex *, doublecomplex *, doublecomplex *, integer *, 
00145             doublecomplex *, integer *, doublecomplex *, integer *, 
00146             doublereal *, logical *, integer *);
00147     integer nmats, jsize;
00148     extern /* Subroutine */ int zget54_(integer *, doublecomplex *, integer *, 
00149              doublecomplex *, integer *, doublecomplex *, integer *, 
00150             doublecomplex *, integer *, doublecomplex *, integer *, 
00151             doublecomplex *, integer *, doublecomplex *, doublereal *);
00152     integer nerrs, jtype, ntest, isort;
00153     extern /* Subroutine */ int dlabad_(doublereal *, doublereal *), zlatm4_(
00154             integer *, integer *, integer *, integer *, logical *, doublereal 
00155             *, doublereal *, doublereal *, integer *, integer *, 
00156             doublecomplex *, integer *);
00157     logical ilabad;
00158     extern doublereal dlamch_(char *);
00159     extern /* Subroutine */ int zunm2r_(char *, char *, integer *, integer *, 
00160             integer *, doublecomplex *, integer *, doublecomplex *, 
00161             doublecomplex *, integer *, doublecomplex *, integer *);
00162     doublereal safmin, safmax;
00163     integer knteig, ioldsd[4];
00164     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00165             integer *, integer *);
00166     extern /* Subroutine */ int alasvm_(char *, integer *, integer *, integer 
00167             *, integer *), xerbla_(char *, integer *), 
00168             zlarfg_(integer *, doublecomplex *, doublecomplex *, integer *, 
00169             doublecomplex *);
00170     extern /* Double Complex */ VOID zlarnd_(doublecomplex *, integer *, 
00171             integer *);
00172     extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 
00173             doublecomplex *, integer *, doublecomplex *, integer *), 
00174             zlaset_(char *, integer *, integer *, doublecomplex *, 
00175             doublecomplex *, doublecomplex *, integer *);
00176     extern logical zlctes_(doublecomplex *, doublecomplex *);
00177     integer minwrk, maxwrk;
00178     doublereal ulpinv;
00179     integer mtypes, ntestt;
00180 
00181     /* Fortran I/O blocks */
00182     static cilist io___41 = { 0, 0, 0, fmt_9999, 0 };
00183     static cilist io___47 = { 0, 0, 0, fmt_9999, 0 };
00184     static cilist io___51 = { 0, 0, 0, fmt_9998, 0 };
00185     static cilist io___53 = { 0, 0, 0, fmt_9997, 0 };
00186     static cilist io___54 = { 0, 0, 0, fmt_9996, 0 };
00187     static cilist io___55 = { 0, 0, 0, fmt_9995, 0 };
00188     static cilist io___56 = { 0, 0, 0, fmt_9994, 0 };
00189     static cilist io___57 = { 0, 0, 0, fmt_9993, 0 };
00190     static cilist io___58 = { 0, 0, 0, fmt_9992, 0 };
00191     static cilist io___59 = { 0, 0, 0, fmt_9991, 0 };
00192 
00193 
00194 
00195 /*  -- LAPACK test routine (version 3.1.1) -- */
00196 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00197 /*     February 2007 */
00198 
00199 /*     .. Scalar Arguments .. */
00200 /*     .. */
00201 /*     .. Array Arguments .. */
00202 /*     .. */
00203 
00204 /*  Purpose */
00205 /*  ======= */
00206 
00207 /*  ZDRGES checks the nonsymmetric generalized eigenvalue (Schur form) */
00208 /*  problem driver ZGGES. */
00209 
00210 /*  ZGGES factors A and B as Q*S*Z'  and Q*T*Z' , where ' means conjugate */
00211 /*  transpose, S and T are  upper triangular (i.e., in generalized Schur */
00212 /*  form), and Q and Z are unitary. It also computes the generalized */
00213 /*  eigenvalues (alpha(j),beta(j)), j=1,...,n.  Thus, */
00214 /*  w(j) = alpha(j)/beta(j) is a root of the characteristic equation */
00215 
00216 /*                  det( A - w(j) B ) = 0 */
00217 
00218 /*  Optionally it also reorder the eigenvalues so that a selected */
00219 /*  cluster of eigenvalues appears in the leading diagonal block of the */
00220 /*  Schur forms. */
00221 
00222 /*  When ZDRGES is called, a number of matrix "sizes" ("N's") and a */
00223 /*  number of matrix "TYPES" are specified.  For each size ("N") */
00224 /*  and each TYPE of matrix, a pair of matrices (A, B) will be generated */
00225 /*  and used for testing. For each matrix pair, the following 13 tests */
00226 /*  will be performed and compared with the threshhold THRESH except */
00227 /*  the tests (5), (11) and (13). */
00228 
00229 
00230 /*  (1)   | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues) */
00231 
00232 
00233 /*  (2)   | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues) */
00234 
00235 
00236 /*  (3)   | I - QQ' | / ( n ulp ) (no sorting of eigenvalues) */
00237 
00238 
00239 /*  (4)   | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues) */
00240 
00241 /*  (5)   if A is in Schur form (i.e. triangular form) (no sorting of */
00242 /*        eigenvalues) */
00243 
00244 /*  (6)   if eigenvalues = diagonal elements of the Schur form (S, T), */
00245 /*        i.e., test the maximum over j of D(j)  where: */
00246 
00247 /*                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)| */
00248 /*            D(j) = ------------------------ + ----------------------- */
00249 /*                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|) */
00250 
00251 /*        (no sorting of eigenvalues) */
00252 
00253 /*  (7)   | (A,B) - Q (S,T) Z' | / ( |(A,B)| n ulp ) */
00254 /*        (with sorting of eigenvalues). */
00255 
00256 /*  (8)   | I - QQ' | / ( n ulp ) (with sorting of eigenvalues). */
00257 
00258 /*  (9)   | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues). */
00259 
00260 /*  (10)  if A is in Schur form (i.e. quasi-triangular form) */
00261 /*        (with sorting of eigenvalues). */
00262 
00263 /*  (11)  if eigenvalues = diagonal elements of the Schur form (S, T), */
00264 /*        i.e. test the maximum over j of D(j)  where: */
00265 
00266 /*                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)| */
00267 /*            D(j) = ------------------------ + ----------------------- */
00268 /*                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|) */
00269 
00270 /*        (with sorting of eigenvalues). */
00271 
00272 /*  (12)  if sorting worked and SDIM is the number of eigenvalues */
00273 /*        which were CELECTed. */
00274 
00275 /*  Test Matrices */
00276 /*  ============= */
00277 
00278 /*  The sizes of the test matrices are specified by an array */
00279 /*  NN(1:NSIZES); the value of each element NN(j) specifies one size. */
00280 /*  The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if */
00281 /*  DOTYPE(j) is .TRUE., then matrix type "j" will be generated. */
00282 /*  Currently, the list of possible types is: */
00283 
00284 /*  (1)  ( 0, 0 )         (a pair of zero matrices) */
00285 
00286 /*  (2)  ( I, 0 )         (an identity and a zero matrix) */
00287 
00288 /*  (3)  ( 0, I )         (an identity and a zero matrix) */
00289 
00290 /*  (4)  ( I, I )         (a pair of identity matrices) */
00291 
00292 /*          t   t */
00293 /*  (5)  ( J , J  )       (a pair of transposed Jordan blocks) */
00294 
00295 /*                                      t                ( I   0  ) */
00296 /*  (6)  ( X, Y )         where  X = ( J   0  )  and Y = (      t ) */
00297 /*                                   ( 0   I  )          ( 0   J  ) */
00298 /*                        and I is a k x k identity and J a (k+1)x(k+1) */
00299 /*                        Jordan block; k=(N-1)/2 */
00300 
00301 /*  (7)  ( D, I )         where D is diag( 0, 1,..., N-1 ) (a diagonal */
00302 /*                        matrix with those diagonal entries.) */
00303 /*  (8)  ( I, D ) */
00304 
00305 /*  (9)  ( big*D, small*I ) where "big" is near overflow and small=1/big */
00306 
00307 /*  (10) ( small*D, big*I ) */
00308 
00309 /*  (11) ( big*I, small*D ) */
00310 
00311 /*  (12) ( small*I, big*D ) */
00312 
00313 /*  (13) ( big*D, big*I ) */
00314 
00315 /*  (14) ( small*D, small*I ) */
00316 
00317 /*  (15) ( D1, D2 )        where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and */
00318 /*                         D2 is diag( 0, N-3, N-4,..., 1, 0, 0 ) */
00319 /*            t   t */
00320 /*  (16) Q ( J , J ) Z     where Q and Z are random orthogonal matrices. */
00321 
00322 /*  (17) Q ( T1, T2 ) Z    where T1 and T2 are upper triangular matrices */
00323 /*                         with random O(1) entries above the diagonal */
00324 /*                         and diagonal entries diag(T1) = */
00325 /*                         ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) = */
00326 /*                         ( 0, N-3, N-4,..., 1, 0, 0 ) */
00327 
00328 /*  (18) Q ( T1, T2 ) Z    diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 ) */
00329 /*                         diag(T2) = ( 0, 1, 0, 1,..., 1, 0 ) */
00330 /*                         s = machine precision. */
00331 
00332 /*  (19) Q ( T1, T2 ) Z    diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 ) */
00333 /*                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 ) */
00334 
00335 /*                                                         N-5 */
00336 /*  (20) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, 1, a, ..., a   =s, 0 ) */
00337 /*                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) */
00338 
00339 /*  (21) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 ) */
00340 /*                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 ) */
00341 /*                         where r1,..., r(N-4) are random. */
00342 
00343 /*  (22) Q ( big*T1, small*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) */
00344 /*                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 ) */
00345 
00346 /*  (23) Q ( small*T1, big*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) */
00347 /*                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 ) */
00348 
00349 /*  (24) Q ( small*T1, small*T2 ) Z  diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) */
00350 /*                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 ) */
00351 
00352 /*  (25) Q ( big*T1, big*T2 ) Z      diag(T1) = ( 0, 0, 1, ..., N-3, 0 ) */
00353 /*                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 ) */
00354 
00355 /*  (26) Q ( T1, T2 ) Z     where T1 and T2 are random upper-triangular */
00356 /*                          matrices. */
00357 
00358 
00359 /*  Arguments */
00360 /*  ========= */
00361 
00362 /*  NSIZES  (input) INTEGER */
00363 /*          The number of sizes of matrices to use.  If it is zero, */
00364 /*          DDRGES does nothing.  NSIZES >= 0. */
00365 
00366 /*  NN      (input) INTEGER array, dimension (NSIZES) */
00367 /*          An array containing the sizes to be used for the matrices. */
00368 /*          Zero values will be skipped.  NN >= 0. */
00369 
00370 /*  NTYPES  (input) INTEGER */
00371 /*          The number of elements in DOTYPE.   If it is zero, DDRGES */
00372 /*          does nothing.  It must be at least zero.  If it is MAXTYP+1 */
00373 /*          and NSIZES is 1, then an additional type, MAXTYP+1 is */
00374 /*          defined, which is to use whatever matrix is in A on input. */
00375 /*          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and */
00376 /*          DOTYPE(MAXTYP+1) is .TRUE. . */
00377 
00378 /*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) */
00379 /*          If DOTYPE(j) is .TRUE., then for each size in NN a */
00380 /*          matrix of that size and of type j will be generated. */
00381 /*          If NTYPES is smaller than the maximum number of types */
00382 /*          defined (PARAMETER MAXTYP), then types NTYPES+1 through */
00383 /*          MAXTYP will not be generated. If NTYPES is larger */
00384 /*          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) */
00385 /*          will be ignored. */
00386 
00387 /*  ISEED   (input/output) INTEGER array, dimension (4) */
00388 /*          On entry ISEED specifies the seed of the random number */
00389 /*          generator. The array elements should be between 0 and 4095; */
00390 /*          if not they will be reduced mod 4096. Also, ISEED(4) must */
00391 /*          be odd.  The random number generator uses a linear */
00392 /*          congruential sequence limited to small integers, and so */
00393 /*          should produce machine independent random numbers. The */
00394 /*          values of ISEED are changed on exit, and can be used in the */
00395 /*          next call to DDRGES to continue the same random number */
00396 /*          sequence. */
00397 
00398 /*  THRESH  (input) DOUBLE PRECISION */
00399 /*          A test will count as "failed" if the "error", computed as */
00400 /*          described above, exceeds THRESH.  Note that the error is */
00401 /*          scaled to be O(1), so THRESH should be a reasonably small */
00402 /*          multiple of 1, e.g., 10 or 100.  In particular, it should */
00403 /*          not depend on the precision (single vs. double) or the size */
00404 /*          of the matrix.  THRESH >= 0. */
00405 
00406 /*  NOUNIT  (input) INTEGER */
00407 /*          The FORTRAN unit number for printing out error messages */
00408 /*          (e.g., if a routine returns IINFO not equal to 0.) */
00409 
00410 /*  A       (input/workspace) COMPLEX*16 array, dimension(LDA, max(NN)) */
00411 /*          Used to hold the original A matrix.  Used as input only */
00412 /*          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and */
00413 /*          DOTYPE(MAXTYP+1)=.TRUE. */
00414 
00415 /*  LDA     (input) INTEGER */
00416 /*          The leading dimension of A, B, S, and T. */
00417 /*          It must be at least 1 and at least max( NN ). */
00418 
00419 /*  B       (input/workspace) COMPLEX*16 array, dimension(LDA, max(NN)) */
00420 /*          Used to hold the original B matrix.  Used as input only */
00421 /*          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and */
00422 /*          DOTYPE(MAXTYP+1)=.TRUE. */
00423 
00424 /*  S       (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) */
00425 /*          The Schur form matrix computed from A by ZGGES.  On exit, S */
00426 /*          contains the Schur form matrix corresponding to the matrix */
00427 /*          in A. */
00428 
00429 /*  T       (workspace) COMPLEX*16 array, dimension (LDA, max(NN)) */
00430 /*          The upper triangular matrix computed from B by ZGGES. */
00431 
00432 /*  Q       (workspace) COMPLEX*16 array, dimension (LDQ, max(NN)) */
00433 /*          The (left) orthogonal matrix computed by ZGGES. */
00434 
00435 /*  LDQ     (input) INTEGER */
00436 /*          The leading dimension of Q and Z. It must */
00437 /*          be at least 1 and at least max( NN ). */
00438 
00439 /*  Z       (workspace) COMPLEX*16 array, dimension( LDQ, max(NN) ) */
00440 /*          The (right) orthogonal matrix computed by ZGGES. */
00441 
00442 /*  ALPHA   (workspace) COMPLEX*16 array, dimension (max(NN)) */
00443 /*  BETA    (workspace) COMPLEX*16 array, dimension (max(NN)) */
00444 /*          The generalized eigenvalues of (A,B) computed by ZGGES. */
00445 /*          ALPHA(k) / BETA(k) is the k-th generalized eigenvalue of A */
00446 /*          and B. */
00447 
00448 /*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK) */
00449 
00450 /*  LWORK   (input) INTEGER */
00451 /*          The dimension of the array WORK.  LWORK >= 3*N*N. */
00452 
00453 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension ( 8*N ) */
00454 /*          Real workspace. */
00455 
00456 /*  RESULT  (output) DOUBLE PRECISION array, dimension (15) */
00457 /*          The values computed by the tests described above. */
00458 /*          The values are currently limited to 1/ulp, to avoid overflow. */
00459 
00460 /*  BWORK   (workspace) LOGICAL array, dimension (N) */
00461 
00462 /*  INFO    (output) INTEGER */
00463 /*          = 0:  successful exit */
00464 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00465 /*          > 0:  A routine returned an error code.  INFO is the */
00466 /*                absolute value of the INFO value returned. */
00467 
00468 /*  ===================================================================== */
00469 
00470 /*     .. Parameters .. */
00471 /*     .. */
00472 /*     .. Local Scalars .. */
00473 /*     .. */
00474 /*     .. Local Arrays .. */
00475 /*     .. */
00476 /*     .. External Functions .. */
00477 /*     .. */
00478 /*     .. External Subroutines .. */
00479 /*     .. */
00480 /*     .. Intrinsic Functions .. */
00481 /*     .. */
00482 /*     .. Statement Functions .. */
00483 /*     .. */
00484 /*     .. Statement Function definitions .. */
00485 /*     .. */
00486 /*     .. Data statements .. */
00487     /* Parameter adjustments */
00488     --nn;
00489     --dotype;
00490     --iseed;
00491     t_dim1 = *lda;
00492     t_offset = 1 + t_dim1;
00493     t -= t_offset;
00494     s_dim1 = *lda;
00495     s_offset = 1 + s_dim1;
00496     s -= s_offset;
00497     b_dim1 = *lda;
00498     b_offset = 1 + b_dim1;
00499     b -= b_offset;
00500     a_dim1 = *lda;
00501     a_offset = 1 + a_dim1;
00502     a -= a_offset;
00503     z_dim1 = *ldq;
00504     z_offset = 1 + z_dim1;
00505     z__ -= z_offset;
00506     q_dim1 = *ldq;
00507     q_offset = 1 + q_dim1;
00508     q -= q_offset;
00509     --alpha;
00510     --beta;
00511     --work;
00512     --rwork;
00513     --result;
00514     --bwork;
00515 
00516     /* Function Body */
00517 /*     .. */
00518 /*     .. Executable Statements .. */
00519 
00520 /*     Check for errors */
00521 
00522     *info = 0;
00523 
00524     badnn = FALSE_;
00525     nmax = 1;
00526     i__1 = *nsizes;
00527     for (j = 1; j <= i__1; ++j) {
00528 /* Computing MAX */
00529         i__2 = nmax, i__3 = nn[j];
00530         nmax = max(i__2,i__3);
00531         if (nn[j] < 0) {
00532             badnn = TRUE_;
00533         }
00534 /* L10: */
00535     }
00536 
00537     if (*nsizes < 0) {
00538         *info = -1;
00539     } else if (badnn) {
00540         *info = -2;
00541     } else if (*ntypes < 0) {
00542         *info = -3;
00543     } else if (*thresh < 0.) {
00544         *info = -6;
00545     } else if (*lda <= 1 || *lda < nmax) {
00546         *info = -9;
00547     } else if (*ldq <= 1 || *ldq < nmax) {
00548         *info = -14;
00549     }
00550 
00551 /*     Compute workspace */
00552 /*      (Note: Comments in the code beginning "Workspace:" describe the */
00553 /*       minimal amount of workspace needed at that point in the code, */
00554 /*       as well as the preferred amount for good performance. */
00555 /*       NB refers to the optimal block size for the immediately */
00556 /*       following subroutine, as returned by ILAENV. */
00557 
00558     minwrk = 1;
00559     if (*info == 0 && *lwork >= 1) {
00560         minwrk = nmax * 3 * nmax;
00561 /* Computing MAX */
00562         i__1 = 1, i__2 = ilaenv_(&c__1, "ZGEQRF", " ", &nmax, &nmax, &c_n1, &
00563                 c_n1), i__1 = max(i__1,i__2), i__2 = 
00564                 ilaenv_(&c__1, "ZUNMQR", "LC", &nmax, &nmax, &nmax, &c_n1), i__1 = max(i__1,i__2), i__2 = ilaenv_(&
00565                 c__1, "ZUNGQR", " ", &nmax, &nmax, &nmax, &c_n1);
00566         nb = max(i__1,i__2);
00567 /* Computing MAX */
00568         i__1 = nmax + nmax * nb, i__2 = nmax * 3 * nmax;
00569         maxwrk = max(i__1,i__2);
00570         work[1].r = (doublereal) maxwrk, work[1].i = 0.;
00571     }
00572 
00573     if (*lwork < minwrk) {
00574         *info = -19;
00575     }
00576 
00577     if (*info != 0) {
00578         i__1 = -(*info);
00579         xerbla_("ZDRGES", &i__1);
00580         return 0;
00581     }
00582 
00583 /*     Quick return if possible */
00584 
00585     if (*nsizes == 0 || *ntypes == 0) {
00586         return 0;
00587     }
00588 
00589     ulp = dlamch_("Precision");
00590     safmin = dlamch_("Safe minimum");
00591     safmin /= ulp;
00592     safmax = 1. / safmin;
00593     dlabad_(&safmin, &safmax);
00594     ulpinv = 1. / ulp;
00595 
00596 /*     The values RMAGN(2:3) depend on N, see below. */
00597 
00598     rmagn[0] = 0.;
00599     rmagn[1] = 1.;
00600 
00601 /*     Loop over matrix sizes */
00602 
00603     ntestt = 0;
00604     nerrs = 0;
00605     nmats = 0;
00606 
00607     i__1 = *nsizes;
00608     for (jsize = 1; jsize <= i__1; ++jsize) {
00609         n = nn[jsize];
00610         n1 = max(1,n);
00611         rmagn[2] = safmax * ulp / (doublereal) n1;
00612         rmagn[3] = safmin * ulpinv * (doublereal) n1;
00613 
00614         if (*nsizes != 1) {
00615             mtypes = min(26,*ntypes);
00616         } else {
00617             mtypes = min(27,*ntypes);
00618         }
00619 
00620 /*        Loop over matrix types */
00621 
00622         i__2 = mtypes;
00623         for (jtype = 1; jtype <= i__2; ++jtype) {
00624             if (! dotype[jtype]) {
00625                 goto L180;
00626             }
00627             ++nmats;
00628             ntest = 0;
00629 
00630 /*           Save ISEED in case of an error. */
00631 
00632             for (j = 1; j <= 4; ++j) {
00633                 ioldsd[j - 1] = iseed[j];
00634 /* L20: */
00635             }
00636 
00637 /*           Initialize RESULT */
00638 
00639             for (j = 1; j <= 13; ++j) {
00640                 result[j] = 0.;
00641 /* L30: */
00642             }
00643 
00644 /*           Generate test matrices A and B */
00645 
00646 /*           Description of control parameters: */
00647 
00648 /*           KZLASS: =1 means w/o rotation, =2 means w/ rotation, */
00649 /*                   =3 means random. */
00650 /*           KATYPE: the "type" to be passed to ZLATM4 for computing A. */
00651 /*           KAZERO: the pattern of zeros on the diagonal for A: */
00652 /*                   =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ), */
00653 /*                   =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ), */
00654 /*                   =6: ( 0, 1, 0, xxx, 0 ).  (xxx means a string of */
00655 /*                   non-zero entries.) */
00656 /*           KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1), */
00657 /*                   =2: large, =3: small. */
00658 /*           LASIGN: .TRUE. if the diagonal elements of A are to be */
00659 /*                   multiplied by a random magnitude 1 number. */
00660 /*           KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B. */
00661 /*           KTRIAN: =0: don't fill in the upper triangle, =1: do. */
00662 /*           KZ1, KZ2, KADD: used to implement KAZERO and KBZERO. */
00663 /*           RMAGN: used to implement KAMAGN and KBMAGN. */
00664 
00665             if (mtypes > 26) {
00666                 goto L110;
00667             }
00668             iinfo = 0;
00669             if (kclass[jtype - 1] < 3) {
00670 
00671 /*              Generate A (w/o rotation) */
00672 
00673                 if ((i__3 = katype[jtype - 1], abs(i__3)) == 3) {
00674                     in = ((n - 1) / 2 << 1) + 1;
00675                     if (in != n) {
00676                         zlaset_("Full", &n, &n, &c_b1, &c_b1, &a[a_offset], 
00677                                 lda);
00678                     }
00679                 } else {
00680                     in = n;
00681                 }
00682                 zlatm4_(&katype[jtype - 1], &in, &kz1[kazero[jtype - 1] - 1], 
00683                         &kz2[kazero[jtype - 1] - 1], &lasign[jtype - 1], &
00684                         rmagn[kamagn[jtype - 1]], &ulp, &rmagn[ktrian[jtype - 
00685                         1] * kamagn[jtype - 1]], &c__2, &iseed[1], &a[
00686                         a_offset], lda);
00687                 iadd = kadd[kazero[jtype - 1] - 1];
00688                 if (iadd > 0 && iadd <= n) {
00689                     i__3 = iadd + iadd * a_dim1;
00690                     i__4 = kamagn[jtype - 1];
00691                     a[i__3].r = rmagn[i__4], a[i__3].i = 0.;
00692                 }
00693 
00694 /*              Generate B (w/o rotation) */
00695 
00696                 if ((i__3 = kbtype[jtype - 1], abs(i__3)) == 3) {
00697                     in = ((n - 1) / 2 << 1) + 1;
00698                     if (in != n) {
00699                         zlaset_("Full", &n, &n, &c_b1, &c_b1, &b[b_offset], 
00700                                 lda);
00701                     }
00702                 } else {
00703                     in = n;
00704                 }
00705                 zlatm4_(&kbtype[jtype - 1], &in, &kz1[kbzero[jtype - 1] - 1], 
00706                         &kz2[kbzero[jtype - 1] - 1], &lbsign[jtype - 1], &
00707                         rmagn[kbmagn[jtype - 1]], &c_b29, &rmagn[ktrian[jtype 
00708                         - 1] * kbmagn[jtype - 1]], &c__2, &iseed[1], &b[
00709                         b_offset], lda);
00710                 iadd = kadd[kbzero[jtype - 1] - 1];
00711                 if (iadd != 0 && iadd <= n) {
00712                     i__3 = iadd + iadd * b_dim1;
00713                     i__4 = kbmagn[jtype - 1];
00714                     b[i__3].r = rmagn[i__4], b[i__3].i = 0.;
00715                 }
00716 
00717                 if (kclass[jtype - 1] == 2 && n > 0) {
00718 
00719 /*                 Include rotations */
00720 
00721 /*                 Generate Q, Z as Householder transformations times */
00722 /*                 a diagonal matrix. */
00723 
00724                     i__3 = n - 1;
00725                     for (jc = 1; jc <= i__3; ++jc) {
00726                         i__4 = n;
00727                         for (jr = jc; jr <= i__4; ++jr) {
00728                             i__5 = jr + jc * q_dim1;
00729                             zlarnd_(&z__1, &c__3, &iseed[1]);
00730                             q[i__5].r = z__1.r, q[i__5].i = z__1.i;
00731                             i__5 = jr + jc * z_dim1;
00732                             zlarnd_(&z__1, &c__3, &iseed[1]);
00733                             z__[i__5].r = z__1.r, z__[i__5].i = z__1.i;
00734 /* L40: */
00735                         }
00736                         i__4 = n + 1 - jc;
00737                         zlarfg_(&i__4, &q[jc + jc * q_dim1], &q[jc + 1 + jc * 
00738                                 q_dim1], &c__1, &work[jc]);
00739                         i__4 = (n << 1) + jc;
00740                         i__5 = jc + jc * q_dim1;
00741                         d__2 = q[i__5].r;
00742                         d__1 = d_sign(&c_b29, &d__2);
00743                         work[i__4].r = d__1, work[i__4].i = 0.;
00744                         i__4 = jc + jc * q_dim1;
00745                         q[i__4].r = 1., q[i__4].i = 0.;
00746                         i__4 = n + 1 - jc;
00747                         zlarfg_(&i__4, &z__[jc + jc * z_dim1], &z__[jc + 1 + 
00748                                 jc * z_dim1], &c__1, &work[n + jc]);
00749                         i__4 = n * 3 + jc;
00750                         i__5 = jc + jc * z_dim1;
00751                         d__2 = z__[i__5].r;
00752                         d__1 = d_sign(&c_b29, &d__2);
00753                         work[i__4].r = d__1, work[i__4].i = 0.;
00754                         i__4 = jc + jc * z_dim1;
00755                         z__[i__4].r = 1., z__[i__4].i = 0.;
00756 /* L50: */
00757                     }
00758                     zlarnd_(&z__1, &c__3, &iseed[1]);
00759                     ctemp.r = z__1.r, ctemp.i = z__1.i;
00760                     i__3 = n + n * q_dim1;
00761                     q[i__3].r = 1., q[i__3].i = 0.;
00762                     i__3 = n;
00763                     work[i__3].r = 0., work[i__3].i = 0.;
00764                     i__3 = n * 3;
00765                     d__1 = z_abs(&ctemp);
00766                     z__1.r = ctemp.r / d__1, z__1.i = ctemp.i / d__1;
00767                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00768                     zlarnd_(&z__1, &c__3, &iseed[1]);
00769                     ctemp.r = z__1.r, ctemp.i = z__1.i;
00770                     i__3 = n + n * z_dim1;
00771                     z__[i__3].r = 1., z__[i__3].i = 0.;
00772                     i__3 = n << 1;
00773                     work[i__3].r = 0., work[i__3].i = 0.;
00774                     i__3 = n << 2;
00775                     d__1 = z_abs(&ctemp);
00776                     z__1.r = ctemp.r / d__1, z__1.i = ctemp.i / d__1;
00777                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00778 
00779 /*                 Apply the diagonal matrices */
00780 
00781                     i__3 = n;
00782                     for (jc = 1; jc <= i__3; ++jc) {
00783                         i__4 = n;
00784                         for (jr = 1; jr <= i__4; ++jr) {
00785                             i__5 = jr + jc * a_dim1;
00786                             i__6 = (n << 1) + jr;
00787                             d_cnjg(&z__3, &work[n * 3 + jc]);
00788                             z__2.r = work[i__6].r * z__3.r - work[i__6].i * 
00789                                     z__3.i, z__2.i = work[i__6].r * z__3.i + 
00790                                     work[i__6].i * z__3.r;
00791                             i__7 = jr + jc * a_dim1;
00792                             z__1.r = z__2.r * a[i__7].r - z__2.i * a[i__7].i, 
00793                                     z__1.i = z__2.r * a[i__7].i + z__2.i * a[
00794                                     i__7].r;
00795                             a[i__5].r = z__1.r, a[i__5].i = z__1.i;
00796                             i__5 = jr + jc * b_dim1;
00797                             i__6 = (n << 1) + jr;
00798                             d_cnjg(&z__3, &work[n * 3 + jc]);
00799                             z__2.r = work[i__6].r * z__3.r - work[i__6].i * 
00800                                     z__3.i, z__2.i = work[i__6].r * z__3.i + 
00801                                     work[i__6].i * z__3.r;
00802                             i__7 = jr + jc * b_dim1;
00803                             z__1.r = z__2.r * b[i__7].r - z__2.i * b[i__7].i, 
00804                                     z__1.i = z__2.r * b[i__7].i + z__2.i * b[
00805                                     i__7].r;
00806                             b[i__5].r = z__1.r, b[i__5].i = z__1.i;
00807 /* L60: */
00808                         }
00809 /* L70: */
00810                     }
00811                     i__3 = n - 1;
00812                     zunm2r_("L", "N", &n, &n, &i__3, &q[q_offset], ldq, &work[
00813                             1], &a[a_offset], lda, &work[(n << 1) + 1], &
00814                             iinfo);
00815                     if (iinfo != 0) {
00816                         goto L100;
00817                     }
00818                     i__3 = n - 1;
00819                     zunm2r_("R", "C", &n, &n, &i__3, &z__[z_offset], ldq, &
00820                             work[n + 1], &a[a_offset], lda, &work[(n << 1) + 
00821                             1], &iinfo);
00822                     if (iinfo != 0) {
00823                         goto L100;
00824                     }
00825                     i__3 = n - 1;
00826                     zunm2r_("L", "N", &n, &n, &i__3, &q[q_offset], ldq, &work[
00827                             1], &b[b_offset], lda, &work[(n << 1) + 1], &
00828                             iinfo);
00829                     if (iinfo != 0) {
00830                         goto L100;
00831                     }
00832                     i__3 = n - 1;
00833                     zunm2r_("R", "C", &n, &n, &i__3, &z__[z_offset], ldq, &
00834                             work[n + 1], &b[b_offset], lda, &work[(n << 1) + 
00835                             1], &iinfo);
00836                     if (iinfo != 0) {
00837                         goto L100;
00838                     }
00839                 }
00840             } else {
00841 
00842 /*              Random matrices */
00843 
00844                 i__3 = n;
00845                 for (jc = 1; jc <= i__3; ++jc) {
00846                     i__4 = n;
00847                     for (jr = 1; jr <= i__4; ++jr) {
00848                         i__5 = jr + jc * a_dim1;
00849                         i__6 = kamagn[jtype - 1];
00850                         zlarnd_(&z__2, &c__4, &iseed[1]);
00851                         z__1.r = rmagn[i__6] * z__2.r, z__1.i = rmagn[i__6] * 
00852                                 z__2.i;
00853                         a[i__5].r = z__1.r, a[i__5].i = z__1.i;
00854                         i__5 = jr + jc * b_dim1;
00855                         i__6 = kbmagn[jtype - 1];
00856                         zlarnd_(&z__2, &c__4, &iseed[1]);
00857                         z__1.r = rmagn[i__6] * z__2.r, z__1.i = rmagn[i__6] * 
00858                                 z__2.i;
00859                         b[i__5].r = z__1.r, b[i__5].i = z__1.i;
00860 /* L80: */
00861                     }
00862 /* L90: */
00863                 }
00864             }
00865 
00866 L100:
00867 
00868             if (iinfo != 0) {
00869                 io___41.ciunit = *nounit;
00870                 s_wsfe(&io___41);
00871                 do_fio(&c__1, "Generator", (ftnlen)9);
00872                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00873                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00874                 do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00875                 do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer));
00876                 e_wsfe();
00877                 *info = abs(iinfo);
00878                 return 0;
00879             }
00880 
00881 L110:
00882 
00883             for (i__ = 1; i__ <= 13; ++i__) {
00884                 result[i__] = -1.;
00885 /* L120: */
00886             }
00887 
00888 /*           Test with and without sorting of eigenvalues */
00889 
00890             for (isort = 0; isort <= 1; ++isort) {
00891                 if (isort == 0) {
00892                     *(unsigned char *)sort = 'N';
00893                     rsub = 0;
00894                 } else {
00895                     *(unsigned char *)sort = 'S';
00896                     rsub = 5;
00897                 }
00898 
00899 /*              Call ZGGES to compute H, T, Q, Z, alpha, and beta. */
00900 
00901                 zlacpy_("Full", &n, &n, &a[a_offset], lda, &s[s_offset], lda);
00902                 zlacpy_("Full", &n, &n, &b[b_offset], lda, &t[t_offset], lda);
00903                 ntest = rsub + 1 + isort;
00904                 result[rsub + 1 + isort] = ulpinv;
00905                 zgges_("V", "V", sort, (L_fp)zlctes_, &n, &s[s_offset], lda, &
00906                         t[t_offset], lda, &sdim, &alpha[1], &beta[1], &q[
00907                         q_offset], ldq, &z__[z_offset], ldq, &work[1], lwork, 
00908                         &rwork[1], &bwork[1], &iinfo);
00909                 if (iinfo != 0 && iinfo != n + 2) {
00910                     result[rsub + 1 + isort] = ulpinv;
00911                     io___47.ciunit = *nounit;
00912                     s_wsfe(&io___47);
00913                     do_fio(&c__1, "ZGGES", (ftnlen)5);
00914                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00915                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00916                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00917                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
00918                             ;
00919                     e_wsfe();
00920                     *info = abs(iinfo);
00921                     goto L160;
00922                 }
00923 
00924                 ntest = rsub + 4;
00925 
00926 /*              Do tests 1--4 (or tests 7--9 when reordering ) */
00927 
00928                 if (isort == 0) {
00929                     zget51_(&c__1, &n, &a[a_offset], lda, &s[s_offset], lda, &
00930                             q[q_offset], ldq, &z__[z_offset], ldq, &work[1], &
00931                             rwork[1], &result[1]);
00932                     zget51_(&c__1, &n, &b[b_offset], lda, &t[t_offset], lda, &
00933                             q[q_offset], ldq, &z__[z_offset], ldq, &work[1], &
00934                             rwork[1], &result[2]);
00935                 } else {
00936                     zget54_(&n, &a[a_offset], lda, &b[b_offset], lda, &s[
00937                             s_offset], lda, &t[t_offset], lda, &q[q_offset], 
00938                             ldq, &z__[z_offset], ldq, &work[1], &result[rsub 
00939                             + 2]);
00940                 }
00941 
00942                 zget51_(&c__3, &n, &b[b_offset], lda, &t[t_offset], lda, &q[
00943                         q_offset], ldq, &q[q_offset], ldq, &work[1], &rwork[1]
00944 , &result[rsub + 3]);
00945                 zget51_(&c__3, &n, &b[b_offset], lda, &t[t_offset], lda, &z__[
00946                         z_offset], ldq, &z__[z_offset], ldq, &work[1], &rwork[
00947                         1], &result[rsub + 4]);
00948 
00949 /*              Do test 5 and 6 (or Tests 10 and 11 when reordering): */
00950 /*              check Schur form of A and compare eigenvalues with */
00951 /*              diagonals. */
00952 
00953                 ntest = rsub + 6;
00954                 temp1 = 0.;
00955 
00956                 i__3 = n;
00957                 for (j = 1; j <= i__3; ++j) {
00958                     ilabad = FALSE_;
00959                     i__4 = j;
00960                     i__5 = j + j * s_dim1;
00961                     z__2.r = alpha[i__4].r - s[i__5].r, z__2.i = alpha[i__4]
00962                             .i - s[i__5].i;
00963                     z__1.r = z__2.r, z__1.i = z__2.i;
00964                     i__6 = j;
00965                     i__7 = j + j * t_dim1;
00966                     z__4.r = beta[i__6].r - t[i__7].r, z__4.i = beta[i__6].i 
00967                             - t[i__7].i;
00968                     z__3.r = z__4.r, z__3.i = z__4.i;
00969 /* Computing MAX */
00970                     i__8 = j;
00971                     i__9 = j + j * s_dim1;
00972                     d__13 = safmin, d__14 = (d__1 = alpha[i__8].r, abs(d__1)) 
00973                             + (d__2 = d_imag(&alpha[j]), abs(d__2)), d__13 = 
00974                             max(d__13,d__14), d__14 = (d__3 = s[i__9].r, abs(
00975                             d__3)) + (d__4 = d_imag(&s[j + j * s_dim1]), abs(
00976                             d__4));
00977 /* Computing MAX */
00978                     i__10 = j;
00979                     i__11 = j + j * t_dim1;
00980                     d__15 = safmin, d__16 = (d__5 = beta[i__10].r, abs(d__5)) 
00981                             + (d__6 = d_imag(&beta[j]), abs(d__6)), d__15 = 
00982                             max(d__15,d__16), d__16 = (d__7 = t[i__11].r, abs(
00983                             d__7)) + (d__8 = d_imag(&t[j + j * t_dim1]), abs(
00984                             d__8));
00985                     temp2 = (((d__9 = z__1.r, abs(d__9)) + (d__10 = d_imag(&
00986                             z__1), abs(d__10))) / max(d__13,d__14) + ((d__11 =
00987                              z__3.r, abs(d__11)) + (d__12 = d_imag(&z__3), 
00988                             abs(d__12))) / max(d__15,d__16)) / ulp;
00989 
00990                     if (j < n) {
00991                         i__4 = j + 1 + j * s_dim1;
00992                         if (s[i__4].r != 0. || s[i__4].i != 0.) {
00993                             ilabad = TRUE_;
00994                             result[rsub + 5] = ulpinv;
00995                         }
00996                     }
00997                     if (j > 1) {
00998                         i__4 = j + (j - 1) * s_dim1;
00999                         if (s[i__4].r != 0. || s[i__4].i != 0.) {
01000                             ilabad = TRUE_;
01001                             result[rsub + 5] = ulpinv;
01002                         }
01003                     }
01004                     temp1 = max(temp1,temp2);
01005                     if (ilabad) {
01006                         io___51.ciunit = *nounit;
01007                         s_wsfe(&io___51);
01008                         do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
01009                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01010                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
01011                                 ;
01012                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
01013                                 integer));
01014                         e_wsfe();
01015                     }
01016 /* L130: */
01017                 }
01018                 result[rsub + 6] = temp1;
01019 
01020                 if (isort >= 1) {
01021 
01022 /*                 Do test 12 */
01023 
01024                     ntest = 12;
01025                     result[12] = 0.;
01026                     knteig = 0;
01027                     i__3 = n;
01028                     for (i__ = 1; i__ <= i__3; ++i__) {
01029                         if (zlctes_(&alpha[i__], &beta[i__])) {
01030                             ++knteig;
01031                         }
01032 /* L140: */
01033                     }
01034                     if (sdim != knteig) {
01035                         result[13] = ulpinv;
01036                     }
01037                 }
01038 
01039 /* L150: */
01040             }
01041 
01042 /*           End of Loop -- Check for RESULT(j) > THRESH */
01043 
01044 L160:
01045 
01046             ntestt += ntest;
01047 
01048 /*           Print out tests which fail. */
01049 
01050             i__3 = ntest;
01051             for (jr = 1; jr <= i__3; ++jr) {
01052                 if (result[jr] >= *thresh) {
01053 
01054 /*                 If this is the first test to fail, */
01055 /*                 print a header to the data file. */
01056 
01057                     if (nerrs == 0) {
01058                         io___53.ciunit = *nounit;
01059                         s_wsfe(&io___53);
01060                         do_fio(&c__1, "ZGS", (ftnlen)3);
01061                         e_wsfe();
01062 
01063 /*                    Matrix types */
01064 
01065                         io___54.ciunit = *nounit;
01066                         s_wsfe(&io___54);
01067                         e_wsfe();
01068                         io___55.ciunit = *nounit;
01069                         s_wsfe(&io___55);
01070                         e_wsfe();
01071                         io___56.ciunit = *nounit;
01072                         s_wsfe(&io___56);
01073                         do_fio(&c__1, "Unitary", (ftnlen)7);
01074                         e_wsfe();
01075 
01076 /*                    Tests performed */
01077 
01078                         io___57.ciunit = *nounit;
01079                         s_wsfe(&io___57);
01080                         do_fio(&c__1, "unitary", (ftnlen)7);
01081                         do_fio(&c__1, "'", (ftnlen)1);
01082                         do_fio(&c__1, "transpose", (ftnlen)9);
01083                         for (j = 1; j <= 8; ++j) {
01084                             do_fio(&c__1, "'", (ftnlen)1);
01085                         }
01086                         e_wsfe();
01087 
01088                     }
01089                     ++nerrs;
01090                     if (result[jr] < 1e4) {
01091                         io___58.ciunit = *nounit;
01092                         s_wsfe(&io___58);
01093                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01094                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
01095                                 ;
01096                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
01097                                 integer));
01098                         do_fio(&c__1, (char *)&jr, (ftnlen)sizeof(integer));
01099                         do_fio(&c__1, (char *)&result[jr], (ftnlen)sizeof(
01100                                 doublereal));
01101                         e_wsfe();
01102                     } else {
01103                         io___59.ciunit = *nounit;
01104                         s_wsfe(&io___59);
01105                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
01106                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
01107                                 ;
01108                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
01109                                 integer));
01110                         do_fio(&c__1, (char *)&jr, (ftnlen)sizeof(integer));
01111                         do_fio(&c__1, (char *)&result[jr], (ftnlen)sizeof(
01112                                 doublereal));
01113                         e_wsfe();
01114                     }
01115                 }
01116 /* L170: */
01117             }
01118 
01119 L180:
01120             ;
01121         }
01122 /* L190: */
01123     }
01124 
01125 /*     Summary */
01126 
01127     alasvm_("ZGS", nounit, &nerrs, &ntestt, &c__0);
01128 
01129     work[1].r = (doublereal) maxwrk, work[1].i = 0.;
01130 
01131     return 0;
01132 
01133 
01134 
01135 
01136 
01137 
01138 
01139 /*     End of ZDRGES */
01140 
01141 } /* zdrges_ */


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