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


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