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


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