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


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