cget23.c
Go to the documentation of this file.
00001 /* cget23.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__4 = 4;
00020 
00021 /* Subroutine */ int cget23_(logical *comp, integer *isrt, char *balanc, 
00022         integer *jtype, real *thresh, integer *iseed, integer *nounit, 
00023         integer *n, complex *a, integer *lda, complex *h__, complex *w, 
00024         complex *w1, complex *vl, integer *ldvl, complex *vr, integer *ldvr, 
00025         complex *lre, integer *ldlre, real *rcondv, real *rcndv1, real *
00026         rcdvin, real *rconde, real *rcnde1, real *rcdein, real *scale, real *
00027         scale1, real *result, complex *work, integer *lwork, real *rwork, 
00028         integer *info)
00029 {
00030     /* Initialized data */
00031 
00032     static char sens[1*2] = "N" "V";
00033 
00034     /* Format strings */
00035     static char fmt_9998[] = "(\002 CGET23: \002,a,\002 returned INFO=\002,i"
00036             "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, BALANC = "
00037             "\002,a,\002, ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
00038     static char fmt_9999[] = "(\002 CGET23: \002,a,\002 returned INFO=\002,i"
00039             "6,\002.\002,/9x,\002N=\002,i6,\002, INPUT EXAMPLE NUMBER = \002,"
00040             "i4)";
00041 
00042     /* System generated locals */
00043     integer a_dim1, a_offset, h_dim1, h_offset, lre_dim1, lre_offset, vl_dim1,
00044              vl_offset, vr_dim1, vr_offset, i__1, i__2, i__3, i__4, i__5;
00045     real r__1, r__2, r__3, r__4, r__5;
00046 
00047     /* Builtin functions */
00048     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00049     double c_abs(complex *), r_imag(complex *);
00050 
00051     /* Local variables */
00052     integer i__, j;
00053     real v;
00054     integer jj, ihi, ilo;
00055     real eps, res[2], tol, ulp, vmx;
00056     integer ihi1, ilo1;
00057     complex cdum[1];
00058     integer kmin;
00059     complex ctmp;
00060     real vmax, tnrm, vrmx, vtst;
00061     extern /* Subroutine */ int cget22_(char *, char *, char *, integer *, 
00062             complex *, integer *, complex *, integer *, complex *, complex *, 
00063             real *, real *);
00064     logical balok, nobal;
00065     real abnrm;
00066     extern logical lsame_(char *, char *);
00067     integer iinfo;
00068     char sense[1];
00069     integer isens;
00070     real tolin, abnrm1;
00071     extern doublereal scnrm2_(integer *, complex *, integer *), slamch_(char *
00072 );
00073     extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 
00074             *, integer *, complex *, integer *), xerbla_(char *, 
00075             integer *), cgeevx_(char *, char *, char *, char *, 
00076             integer *, complex *, integer *, complex *, complex *, integer *, 
00077             complex *, integer *, integer *, integer *, real *, real *, real *
00078 , real *, complex *, integer *, real *, integer *);
00079     integer isensm;
00080     real vricmp, vrimin, smlnum, ulpinv;
00081 
00082     /* Fortran I/O blocks */
00083     static cilist io___14 = { 0, 0, 0, fmt_9998, 0 };
00084     static cilist io___15 = { 0, 0, 0, fmt_9999, 0 };
00085     static cilist io___28 = { 0, 0, 0, fmt_9998, 0 };
00086     static cilist io___29 = { 0, 0, 0, fmt_9999, 0 };
00087     static cilist io___30 = { 0, 0, 0, fmt_9998, 0 };
00088     static cilist io___31 = { 0, 0, 0, fmt_9999, 0 };
00089     static cilist io___32 = { 0, 0, 0, fmt_9998, 0 };
00090     static cilist io___33 = { 0, 0, 0, fmt_9999, 0 };
00091     static cilist io___34 = { 0, 0, 0, fmt_9999, 0 };
00092 
00093 
00094 
00095 /*  -- LAPACK test routine (version 3.1) -- */
00096 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00097 /*     November 2006 */
00098 
00099 /*     .. Scalar Arguments .. */
00100 /*     .. */
00101 /*     .. Array Arguments .. */
00102 /*     .. */
00103 
00104 /*  Purpose */
00105 /*  ======= */
00106 
00107 /*     CGET23  checks the nonsymmetric eigenvalue problem driver CGEEVX. */
00108 /*     If COMP = .FALSE., the first 8 of the following tests will be */
00109 /*     performed on the input matrix A, and also test 9 if LWORK is */
00110 /*     sufficiently large. */
00111 /*     if COMP is .TRUE. all 11 tests will be performed. */
00112 
00113 /*     (1)     | A * VR - VR * W | / ( n |A| ulp ) */
00114 
00115 /*       Here VR is the matrix of unit right eigenvectors. */
00116 /*       W is a diagonal matrix with diagonal entries W(j). */
00117 
00118 /*     (2)     | A**H * VL - VL * W**H | / ( n |A| ulp ) */
00119 
00120 /*       Here VL is the matrix of unit left eigenvectors, A**H is the */
00121 /*       conjugate transpose of A, and W is as above. */
00122 
00123 /*     (3)     | |VR(i)| - 1 | / ulp and largest component real */
00124 
00125 /*       VR(i) denotes the i-th column of VR. */
00126 
00127 /*     (4)     | |VL(i)| - 1 | / ulp and largest component real */
00128 
00129 /*       VL(i) denotes the i-th column of VL. */
00130 
00131 /*     (5)     0 if W(full) = W(partial), 1/ulp otherwise */
00132 
00133 /*       W(full) denotes the eigenvalues computed when VR, VL, RCONDV */
00134 /*       and RCONDE are also computed, and W(partial) denotes the */
00135 /*       eigenvalues computed when only some of VR, VL, RCONDV, and */
00136 /*       RCONDE are computed. */
00137 
00138 /*     (6)     0 if VR(full) = VR(partial), 1/ulp otherwise */
00139 
00140 /*       VR(full) denotes the right eigenvectors computed when VL, RCONDV */
00141 /*       and RCONDE are computed, and VR(partial) denotes the result */
00142 /*       when only some of VL and RCONDV are computed. */
00143 
00144 /*     (7)     0 if VL(full) = VL(partial), 1/ulp otherwise */
00145 
00146 /*       VL(full) denotes the left eigenvectors computed when VR, RCONDV */
00147 /*       and RCONDE are computed, and VL(partial) denotes the result */
00148 /*       when only some of VR and RCONDV are computed. */
00149 
00150 /*     (8)     0 if SCALE, ILO, IHI, ABNRM (full) = */
00151 /*                  SCALE, ILO, IHI, ABNRM (partial) */
00152 /*             1/ulp otherwise */
00153 
00154 /*       SCALE, ILO, IHI and ABNRM describe how the matrix is balanced. */
00155 /*       (full) is when VR, VL, RCONDE and RCONDV are also computed, and */
00156 /*       (partial) is when some are not computed. */
00157 
00158 /*     (9)     0 if RCONDV(full) = RCONDV(partial), 1/ulp otherwise */
00159 
00160 /*       RCONDV(full) denotes the reciprocal condition numbers of the */
00161 /*       right eigenvectors computed when VR, VL and RCONDE are also */
00162 /*       computed. RCONDV(partial) denotes the reciprocal condition */
00163 /*       numbers when only some of VR, VL and RCONDE are computed. */
00164 
00165 /*    (10)     |RCONDV - RCDVIN| / cond(RCONDV) */
00166 
00167 /*       RCONDV is the reciprocal right eigenvector condition number */
00168 /*       computed by CGEEVX and RCDVIN (the precomputed true value) */
00169 /*       is supplied as input. cond(RCONDV) is the condition number of */
00170 /*       RCONDV, and takes errors in computing RCONDV into account, so */
00171 /*       that the resulting quantity should be O(ULP). cond(RCONDV) is */
00172 /*       essentially given by norm(A)/RCONDE. */
00173 
00174 /*    (11)     |RCONDE - RCDEIN| / cond(RCONDE) */
00175 
00176 /*       RCONDE is the reciprocal eigenvalue condition number */
00177 /*       computed by CGEEVX and RCDEIN (the precomputed true value) */
00178 /*       is supplied as input.  cond(RCONDE) is the condition number */
00179 /*       of RCONDE, and takes errors in computing RCONDE into account, */
00180 /*       so that the resulting quantity should be O(ULP). cond(RCONDE) */
00181 /*       is essentially given by norm(A)/RCONDV. */
00182 
00183 /*  Arguments */
00184 /*  ========= */
00185 
00186 /*  COMP    (input) LOGICAL */
00187 /*          COMP describes which input tests to perform: */
00188 /*            = .FALSE. if the computed condition numbers are not to */
00189 /*                      be tested against RCDVIN and RCDEIN */
00190 /*            = .TRUE.  if they are to be compared */
00191 
00192 /*  ISRT    (input) INTEGER */
00193 /*          If COMP = .TRUE., ISRT indicates in how the eigenvalues */
00194 /*          corresponding to values in RCDVIN and RCDEIN are ordered: */
00195 /*            = 0 means the eigenvalues are sorted by */
00196 /*                increasing real part */
00197 /*            = 1 means the eigenvalues are sorted by */
00198 /*                increasing imaginary part */
00199 /*          If COMP = .FALSE., ISRT is not referenced. */
00200 
00201 /*  BALANC  (input) CHARACTER */
00202 /*          Describes the balancing option to be tested. */
00203 /*            = 'N' for no permuting or diagonal scaling */
00204 /*            = 'P' for permuting but no diagonal scaling */
00205 /*            = 'S' for no permuting but diagonal scaling */
00206 /*            = 'B' for permuting and diagonal scaling */
00207 
00208 /*  JTYPE   (input) INTEGER */
00209 /*          Type of input matrix. Used to label output if error occurs. */
00210 
00211 /*  THRESH  (input) REAL */
00212 /*          A test will count as "failed" if the "error", computed as */
00213 /*          described above, exceeds THRESH.  Note that the error */
00214 /*          is scaled to be O(1), so THRESH should be a reasonably */
00215 /*          small multiple of 1, e.g., 10 or 100.  In particular, */
00216 /*          it should not depend on the precision (single vs. double) */
00217 /*          or the size of the matrix.  It must be at least zero. */
00218 
00219 /*  ISEED   (input) INTEGER array, dimension (4) */
00220 /*          If COMP = .FALSE., the random number generator seed */
00221 /*          used to produce matrix. */
00222 /*          If COMP = .TRUE., ISEED(1) = the number of the example. */
00223 /*          Used to label output if error occurs. */
00224 
00225 /*  NOUNIT  (input) INTEGER */
00226 /*          The FORTRAN unit number for printing out error messages */
00227 /*          (e.g., if a routine returns INFO not equal to 0.) */
00228 
00229 /*  N       (input) INTEGER */
00230 /*          The dimension of A. N must be at least 0. */
00231 
00232 /*  A       (input/output) COMPLEX array, dimension (LDA,N) */
00233 /*          Used to hold the matrix whose eigenvalues are to be */
00234 /*          computed. */
00235 
00236 /*  LDA     (input) INTEGER */
00237 /*          The leading dimension of A, and H. LDA must be at */
00238 /*          least 1 and at least N. */
00239 
00240 /*  H       (workspace) COMPLEX array, dimension (LDA,N) */
00241 /*          Another copy of the test matrix A, modified by CGEEVX. */
00242 
00243 /*  W       (workspace) COMPLEX array, dimension (N) */
00244 /*          Contains the eigenvalues of A. */
00245 
00246 /*  W1      (workspace) COMPLEX array, dimension (N) */
00247 /*          Like W, this array contains the eigenvalues of A, */
00248 /*          but those computed when CGEEVX only computes a partial */
00249 /*          eigendecomposition, i.e. not the eigenvalues and left */
00250 /*          and right eigenvectors. */
00251 
00252 /*  VL      (workspace) COMPLEX array, dimension (LDVL,N) */
00253 /*          VL holds the computed left eigenvectors. */
00254 
00255 /*  LDVL    (input) INTEGER */
00256 /*          Leading dimension of VL. Must be at least max(1,N). */
00257 
00258 /*  VR      (workspace) COMPLEX array, dimension (LDVR,N) */
00259 /*          VR holds the computed right eigenvectors. */
00260 
00261 /*  LDVR    (input) INTEGER */
00262 /*          Leading dimension of VR. Must be at least max(1,N). */
00263 
00264 /*  LRE     (workspace) COMPLEX array, dimension (LDLRE,N) */
00265 /*          LRE holds the computed right or left eigenvectors. */
00266 
00267 /*  LDLRE   (input) INTEGER */
00268 /*          Leading dimension of LRE. Must be at least max(1,N). */
00269 
00270 /*  RCONDV  (workspace) REAL array, dimension (N) */
00271 /*          RCONDV holds the computed reciprocal condition numbers */
00272 /*          for eigenvectors. */
00273 
00274 /*  RCNDV1  (workspace) REAL array, dimension (N) */
00275 /*          RCNDV1 holds more computed reciprocal condition numbers */
00276 /*          for eigenvectors. */
00277 
00278 /*  RCDVIN  (input) REAL array, dimension (N) */
00279 /*          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal */
00280 /*          condition numbers for eigenvectors to be compared with */
00281 /*          RCONDV. */
00282 
00283 /*  RCONDE  (workspace) REAL array, dimension (N) */
00284 /*          RCONDE holds the computed reciprocal condition numbers */
00285 /*          for eigenvalues. */
00286 
00287 /*  RCNDE1  (workspace) REAL array, dimension (N) */
00288 /*          RCNDE1 holds more computed reciprocal condition numbers */
00289 /*          for eigenvalues. */
00290 
00291 /*  RCDEIN  (input) REAL array, dimension (N) */
00292 /*          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal */
00293 /*          condition numbers for eigenvalues to be compared with */
00294 /*          RCONDE. */
00295 
00296 /*  SCALE   (workspace) REAL array, dimension (N) */
00297 /*          Holds information describing balancing of matrix. */
00298 
00299 /*  SCALE1  (workspace) REAL array, dimension (N) */
00300 /*          Holds information describing balancing of matrix. */
00301 
00302 /*  RESULT  (output) REAL array, dimension (11) */
00303 /*          The values computed by the 11 tests described above. */
00304 /*          The values are currently limited to 1/ulp, to avoid */
00305 /*          overflow. */
00306 
00307 /*  WORK    (workspace) COMPLEX array, dimension (LWORK) */
00308 
00309 /*  LWORK   (input) INTEGER */
00310 /*          The number of entries in WORK.  This must be at least */
00311 /*          2*N, and 2*N+N**2 if tests 9, 10 or 11 are to be performed. */
00312 
00313 /*  RWORK   (workspace) REAL array, dimension (2*N) */
00314 
00315 /*  INFO    (output) INTEGER */
00316 /*          If 0,  successful exit. */
00317 /*          If <0, input parameter -INFO had an incorrect value. */
00318 /*          If >0, CGEEVX returned an error code, the absolute */
00319 /*                 value of which is returned. */
00320 
00321 /*  ===================================================================== */
00322 
00323 /*     .. Parameters .. */
00324 /*     .. */
00325 /*     .. Local Scalars .. */
00326 /*     .. */
00327 /*     .. Local Arrays .. */
00328 /*     .. */
00329 /*     .. External Functions .. */
00330 /*     .. */
00331 /*     .. External Subroutines .. */
00332 /*     .. */
00333 /*     .. Intrinsic Functions .. */
00334 /*     .. */
00335 /*     .. Data statements .. */
00336     /* Parameter adjustments */
00337     --iseed;
00338     h_dim1 = *lda;
00339     h_offset = 1 + h_dim1;
00340     h__ -= h_offset;
00341     a_dim1 = *lda;
00342     a_offset = 1 + a_dim1;
00343     a -= a_offset;
00344     --w;
00345     --w1;
00346     vl_dim1 = *ldvl;
00347     vl_offset = 1 + vl_dim1;
00348     vl -= vl_offset;
00349     vr_dim1 = *ldvr;
00350     vr_offset = 1 + vr_dim1;
00351     vr -= vr_offset;
00352     lre_dim1 = *ldlre;
00353     lre_offset = 1 + lre_dim1;
00354     lre -= lre_offset;
00355     --rcondv;
00356     --rcndv1;
00357     --rcdvin;
00358     --rconde;
00359     --rcnde1;
00360     --rcdein;
00361     --scale;
00362     --scale1;
00363     --result;
00364     --work;
00365     --rwork;
00366 
00367     /* Function Body */
00368 /*     .. */
00369 /*     .. Executable Statements .. */
00370 
00371 /*     Check for errors */
00372 
00373     nobal = lsame_(balanc, "N");
00374     balok = nobal || lsame_(balanc, "P") || lsame_(
00375             balanc, "S") || lsame_(balanc, "B");
00376     *info = 0;
00377     if (*isrt != 0 && *isrt != 1) {
00378         *info = -2;
00379     } else if (! balok) {
00380         *info = -3;
00381     } else if (*thresh < 0.f) {
00382         *info = -5;
00383     } else if (*nounit <= 0) {
00384         *info = -7;
00385     } else if (*n < 0) {
00386         *info = -8;
00387     } else if (*lda < 1 || *lda < *n) {
00388         *info = -10;
00389     } else if (*ldvl < 1 || *ldvl < *n) {
00390         *info = -15;
00391     } else if (*ldvr < 1 || *ldvr < *n) {
00392         *info = -17;
00393     } else if (*ldlre < 1 || *ldlre < *n) {
00394         *info = -19;
00395     } else if (*lwork < *n << 1 || *comp && *lwork < (*n << 1) + *n * *n) {
00396         *info = -30;
00397     }
00398 
00399     if (*info != 0) {
00400         i__1 = -(*info);
00401         xerbla_("CGET23", &i__1);
00402         return 0;
00403     }
00404 
00405 /*     Quick return if nothing to do */
00406 
00407     for (i__ = 1; i__ <= 11; ++i__) {
00408         result[i__] = -1.f;
00409 /* L10: */
00410     }
00411 
00412     if (*n == 0) {
00413         return 0;
00414     }
00415 
00416 /*     More Important constants */
00417 
00418     ulp = slamch_("Precision");
00419     smlnum = slamch_("S");
00420     ulpinv = 1.f / ulp;
00421 
00422 /*     Compute eigenvalues and eigenvectors, and test them */
00423 
00424     if (*lwork >= (*n << 1) + *n * *n) {
00425         *(unsigned char *)sense = 'B';
00426         isensm = 2;
00427     } else {
00428         *(unsigned char *)sense = 'E';
00429         isensm = 1;
00430     }
00431     clacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00432     cgeevx_(balanc, "V", "V", sense, n, &h__[h_offset], lda, &w[1], &vl[
00433             vl_offset], ldvl, &vr[vr_offset], ldvr, &ilo, &ihi, &scale[1], &
00434             abnrm, &rconde[1], &rcondv[1], &work[1], lwork, &rwork[1], &iinfo);
00435     if (iinfo != 0) {
00436         result[1] = ulpinv;
00437         if (*jtype != 22) {
00438             io___14.ciunit = *nounit;
00439             s_wsfe(&io___14);
00440             do_fio(&c__1, "CGEEVX1", (ftnlen)7);
00441             do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00442             do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00443             do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00444             do_fio(&c__1, balanc, (ftnlen)1);
00445             do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00446             e_wsfe();
00447         } else {
00448             io___15.ciunit = *nounit;
00449             s_wsfe(&io___15);
00450             do_fio(&c__1, "CGEEVX1", (ftnlen)7);
00451             do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00452             do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00453             do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00454             e_wsfe();
00455         }
00456         *info = abs(iinfo);
00457         return 0;
00458     }
00459 
00460 /*     Do Test (1) */
00461 
00462     cget22_("N", "N", "N", n, &a[a_offset], lda, &vr[vr_offset], ldvr, &w[1], 
00463             &work[1], &rwork[1], res);
00464     result[1] = res[0];
00465 
00466 /*     Do Test (2) */
00467 
00468     cget22_("C", "N", "C", n, &a[a_offset], lda, &vl[vl_offset], ldvl, &w[1], 
00469             &work[1], &rwork[1], res);
00470     result[2] = res[0];
00471 
00472 /*     Do Test (3) */
00473 
00474     i__1 = *n;
00475     for (j = 1; j <= i__1; ++j) {
00476         tnrm = scnrm2_(n, &vr[j * vr_dim1 + 1], &c__1);
00477 /* Computing MAX */
00478 /* Computing MIN */
00479         r__4 = ulpinv, r__5 = (r__1 = tnrm - 1.f, dabs(r__1)) / ulp;
00480         r__2 = result[3], r__3 = dmin(r__4,r__5);
00481         result[3] = dmax(r__2,r__3);
00482         vmx = 0.f;
00483         vrmx = 0.f;
00484         i__2 = *n;
00485         for (jj = 1; jj <= i__2; ++jj) {
00486             vtst = c_abs(&vr[jj + j * vr_dim1]);
00487             if (vtst > vmx) {
00488                 vmx = vtst;
00489             }
00490             i__3 = jj + j * vr_dim1;
00491             if (r_imag(&vr[jj + j * vr_dim1]) == 0.f && (r__1 = vr[i__3].r, 
00492                     dabs(r__1)) > vrmx) {
00493                 i__4 = jj + j * vr_dim1;
00494                 vrmx = (r__2 = vr[i__4].r, dabs(r__2));
00495             }
00496 /* L20: */
00497         }
00498         if (vrmx / vmx < 1.f - ulp * 2.f) {
00499             result[3] = ulpinv;
00500         }
00501 /* L30: */
00502     }
00503 
00504 /*     Do Test (4) */
00505 
00506     i__1 = *n;
00507     for (j = 1; j <= i__1; ++j) {
00508         tnrm = scnrm2_(n, &vl[j * vl_dim1 + 1], &c__1);
00509 /* Computing MAX */
00510 /* Computing MIN */
00511         r__4 = ulpinv, r__5 = (r__1 = tnrm - 1.f, dabs(r__1)) / ulp;
00512         r__2 = result[4], r__3 = dmin(r__4,r__5);
00513         result[4] = dmax(r__2,r__3);
00514         vmx = 0.f;
00515         vrmx = 0.f;
00516         i__2 = *n;
00517         for (jj = 1; jj <= i__2; ++jj) {
00518             vtst = c_abs(&vl[jj + j * vl_dim1]);
00519             if (vtst > vmx) {
00520                 vmx = vtst;
00521             }
00522             i__3 = jj + j * vl_dim1;
00523             if (r_imag(&vl[jj + j * vl_dim1]) == 0.f && (r__1 = vl[i__3].r, 
00524                     dabs(r__1)) > vrmx) {
00525                 i__4 = jj + j * vl_dim1;
00526                 vrmx = (r__2 = vl[i__4].r, dabs(r__2));
00527             }
00528 /* L40: */
00529         }
00530         if (vrmx / vmx < 1.f - ulp * 2.f) {
00531             result[4] = ulpinv;
00532         }
00533 /* L50: */
00534     }
00535 
00536 /*     Test for all options of computing condition numbers */
00537 
00538     i__1 = isensm;
00539     for (isens = 1; isens <= i__1; ++isens) {
00540 
00541         *(unsigned char *)sense = *(unsigned char *)&sens[isens - 1];
00542 
00543 /*        Compute eigenvalues only, and test them */
00544 
00545         clacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00546         cgeevx_(balanc, "N", "N", sense, n, &h__[h_offset], lda, &w1[1], cdum, 
00547                  &c__1, cdum, &c__1, &ilo1, &ihi1, &scale1[1], &abnrm1, &
00548                 rcnde1[1], &rcndv1[1], &work[1], lwork, &rwork[1], &iinfo);
00549         if (iinfo != 0) {
00550             result[1] = ulpinv;
00551             if (*jtype != 22) {
00552                 io___28.ciunit = *nounit;
00553                 s_wsfe(&io___28);
00554                 do_fio(&c__1, "CGEEVX2", (ftnlen)7);
00555                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00556                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00557                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00558                 do_fio(&c__1, balanc, (ftnlen)1);
00559                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00560                 e_wsfe();
00561             } else {
00562                 io___29.ciunit = *nounit;
00563                 s_wsfe(&io___29);
00564                 do_fio(&c__1, "CGEEVX2", (ftnlen)7);
00565                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00566                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00567                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00568                 e_wsfe();
00569             }
00570             *info = abs(iinfo);
00571             goto L190;
00572         }
00573 
00574 /*        Do Test (5) */
00575 
00576         i__2 = *n;
00577         for (j = 1; j <= i__2; ++j) {
00578             i__3 = j;
00579             i__4 = j;
00580             if (w[i__3].r != w1[i__4].r || w[i__3].i != w1[i__4].i) {
00581                 result[5] = ulpinv;
00582             }
00583 /* L60: */
00584         }
00585 
00586 /*        Do Test (8) */
00587 
00588         if (! nobal) {
00589             i__2 = *n;
00590             for (j = 1; j <= i__2; ++j) {
00591                 if (scale[j] != scale1[j]) {
00592                     result[8] = ulpinv;
00593                 }
00594 /* L70: */
00595             }
00596             if (ilo != ilo1) {
00597                 result[8] = ulpinv;
00598             }
00599             if (ihi != ihi1) {
00600                 result[8] = ulpinv;
00601             }
00602             if (abnrm != abnrm1) {
00603                 result[8] = ulpinv;
00604             }
00605         }
00606 
00607 /*        Do Test (9) */
00608 
00609         if (isens == 2 && *n > 1) {
00610             i__2 = *n;
00611             for (j = 1; j <= i__2; ++j) {
00612                 if (rcondv[j] != rcndv1[j]) {
00613                     result[9] = ulpinv;
00614                 }
00615 /* L80: */
00616             }
00617         }
00618 
00619 /*        Compute eigenvalues and right eigenvectors, and test them */
00620 
00621         clacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00622         cgeevx_(balanc, "N", "V", sense, n, &h__[h_offset], lda, &w1[1], cdum, 
00623                  &c__1, &lre[lre_offset], ldlre, &ilo1, &ihi1, &scale1[1], &
00624                 abnrm1, &rcnde1[1], &rcndv1[1], &work[1], lwork, &rwork[1], &
00625                 iinfo);
00626         if (iinfo != 0) {
00627             result[1] = ulpinv;
00628             if (*jtype != 22) {
00629                 io___30.ciunit = *nounit;
00630                 s_wsfe(&io___30);
00631                 do_fio(&c__1, "CGEEVX3", (ftnlen)7);
00632                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00633                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00634                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00635                 do_fio(&c__1, balanc, (ftnlen)1);
00636                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00637                 e_wsfe();
00638             } else {
00639                 io___31.ciunit = *nounit;
00640                 s_wsfe(&io___31);
00641                 do_fio(&c__1, "CGEEVX3", (ftnlen)7);
00642                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00643                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00644                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00645                 e_wsfe();
00646             }
00647             *info = abs(iinfo);
00648             goto L190;
00649         }
00650 
00651 /*        Do Test (5) again */
00652 
00653         i__2 = *n;
00654         for (j = 1; j <= i__2; ++j) {
00655             i__3 = j;
00656             i__4 = j;
00657             if (w[i__3].r != w1[i__4].r || w[i__3].i != w1[i__4].i) {
00658                 result[5] = ulpinv;
00659             }
00660 /* L90: */
00661         }
00662 
00663 /*        Do Test (6) */
00664 
00665         i__2 = *n;
00666         for (j = 1; j <= i__2; ++j) {
00667             i__3 = *n;
00668             for (jj = 1; jj <= i__3; ++jj) {
00669                 i__4 = j + jj * vr_dim1;
00670                 i__5 = j + jj * lre_dim1;
00671                 if (vr[i__4].r != lre[i__5].r || vr[i__4].i != lre[i__5].i) {
00672                     result[6] = ulpinv;
00673                 }
00674 /* L100: */
00675             }
00676 /* L110: */
00677         }
00678 
00679 /*        Do Test (8) again */
00680 
00681         if (! nobal) {
00682             i__2 = *n;
00683             for (j = 1; j <= i__2; ++j) {
00684                 if (scale[j] != scale1[j]) {
00685                     result[8] = ulpinv;
00686                 }
00687 /* L120: */
00688             }
00689             if (ilo != ilo1) {
00690                 result[8] = ulpinv;
00691             }
00692             if (ihi != ihi1) {
00693                 result[8] = ulpinv;
00694             }
00695             if (abnrm != abnrm1) {
00696                 result[8] = ulpinv;
00697             }
00698         }
00699 
00700 /*        Do Test (9) again */
00701 
00702         if (isens == 2 && *n > 1) {
00703             i__2 = *n;
00704             for (j = 1; j <= i__2; ++j) {
00705                 if (rcondv[j] != rcndv1[j]) {
00706                     result[9] = ulpinv;
00707                 }
00708 /* L130: */
00709             }
00710         }
00711 
00712 /*        Compute eigenvalues and left eigenvectors, and test them */
00713 
00714         clacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00715         cgeevx_(balanc, "V", "N", sense, n, &h__[h_offset], lda, &w1[1], &lre[
00716                 lre_offset], ldlre, cdum, &c__1, &ilo1, &ihi1, &scale1[1], &
00717                 abnrm1, &rcnde1[1], &rcndv1[1], &work[1], lwork, &rwork[1], &
00718                 iinfo);
00719         if (iinfo != 0) {
00720             result[1] = ulpinv;
00721             if (*jtype != 22) {
00722                 io___32.ciunit = *nounit;
00723                 s_wsfe(&io___32);
00724                 do_fio(&c__1, "CGEEVX4", (ftnlen)7);
00725                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00726                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00727                 do_fio(&c__1, (char *)&(*jtype), (ftnlen)sizeof(integer));
00728                 do_fio(&c__1, balanc, (ftnlen)1);
00729                 do_fio(&c__4, (char *)&iseed[1], (ftnlen)sizeof(integer));
00730                 e_wsfe();
00731             } else {
00732                 io___33.ciunit = *nounit;
00733                 s_wsfe(&io___33);
00734                 do_fio(&c__1, "CGEEVX4", (ftnlen)7);
00735                 do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00736                 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00737                 do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00738                 e_wsfe();
00739             }
00740             *info = abs(iinfo);
00741             goto L190;
00742         }
00743 
00744 /*        Do Test (5) again */
00745 
00746         i__2 = *n;
00747         for (j = 1; j <= i__2; ++j) {
00748             i__3 = j;
00749             i__4 = j;
00750             if (w[i__3].r != w1[i__4].r || w[i__3].i != w1[i__4].i) {
00751                 result[5] = ulpinv;
00752             }
00753 /* L140: */
00754         }
00755 
00756 /*        Do Test (7) */
00757 
00758         i__2 = *n;
00759         for (j = 1; j <= i__2; ++j) {
00760             i__3 = *n;
00761             for (jj = 1; jj <= i__3; ++jj) {
00762                 i__4 = j + jj * vl_dim1;
00763                 i__5 = j + jj * lre_dim1;
00764                 if (vl[i__4].r != lre[i__5].r || vl[i__4].i != lre[i__5].i) {
00765                     result[7] = ulpinv;
00766                 }
00767 /* L150: */
00768             }
00769 /* L160: */
00770         }
00771 
00772 /*        Do Test (8) again */
00773 
00774         if (! nobal) {
00775             i__2 = *n;
00776             for (j = 1; j <= i__2; ++j) {
00777                 if (scale[j] != scale1[j]) {
00778                     result[8] = ulpinv;
00779                 }
00780 /* L170: */
00781             }
00782             if (ilo != ilo1) {
00783                 result[8] = ulpinv;
00784             }
00785             if (ihi != ihi1) {
00786                 result[8] = ulpinv;
00787             }
00788             if (abnrm != abnrm1) {
00789                 result[8] = ulpinv;
00790             }
00791         }
00792 
00793 /*        Do Test (9) again */
00794 
00795         if (isens == 2 && *n > 1) {
00796             i__2 = *n;
00797             for (j = 1; j <= i__2; ++j) {
00798                 if (rcondv[j] != rcndv1[j]) {
00799                     result[9] = ulpinv;
00800                 }
00801 /* L180: */
00802             }
00803         }
00804 
00805 L190:
00806 
00807 /* L200: */
00808         ;
00809     }
00810 
00811 /*     If COMP, compare condition numbers to precomputed ones */
00812 
00813     if (*comp) {
00814         clacpy_("F", n, n, &a[a_offset], lda, &h__[h_offset], lda);
00815         cgeevx_("N", "V", "V", "B", n, &h__[h_offset], lda, &w[1], &vl[
00816                 vl_offset], ldvl, &vr[vr_offset], ldvr, &ilo, &ihi, &scale[1], 
00817                  &abnrm, &rconde[1], &rcondv[1], &work[1], lwork, &rwork[1], &
00818                 iinfo);
00819         if (iinfo != 0) {
00820             result[1] = ulpinv;
00821             io___34.ciunit = *nounit;
00822             s_wsfe(&io___34);
00823             do_fio(&c__1, "CGEEVX5", (ftnlen)7);
00824             do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00825             do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00826             do_fio(&c__1, (char *)&iseed[1], (ftnlen)sizeof(integer));
00827             e_wsfe();
00828             *info = abs(iinfo);
00829             goto L250;
00830         }
00831 
00832 /*        Sort eigenvalues and condition numbers lexicographically */
00833 /*        to compare with inputs */
00834 
00835         i__1 = *n - 1;
00836         for (i__ = 1; i__ <= i__1; ++i__) {
00837             kmin = i__;
00838             if (*isrt == 0) {
00839                 i__2 = i__;
00840                 vrimin = w[i__2].r;
00841             } else {
00842                 vrimin = r_imag(&w[i__]);
00843             }
00844             i__2 = *n;
00845             for (j = i__ + 1; j <= i__2; ++j) {
00846                 if (*isrt == 0) {
00847                     i__3 = j;
00848                     vricmp = w[i__3].r;
00849                 } else {
00850                     vricmp = r_imag(&w[j]);
00851                 }
00852                 if (vricmp < vrimin) {
00853                     kmin = j;
00854                     vrimin = vricmp;
00855                 }
00856 /* L210: */
00857             }
00858             i__2 = kmin;
00859             ctmp.r = w[i__2].r, ctmp.i = w[i__2].i;
00860             i__2 = kmin;
00861             i__3 = i__;
00862             w[i__2].r = w[i__3].r, w[i__2].i = w[i__3].i;
00863             i__2 = i__;
00864             w[i__2].r = ctmp.r, w[i__2].i = ctmp.i;
00865             vrimin = rconde[kmin];
00866             rconde[kmin] = rconde[i__];
00867             rconde[i__] = vrimin;
00868             vrimin = rcondv[kmin];
00869             rcondv[kmin] = rcondv[i__];
00870             rcondv[i__] = vrimin;
00871 /* L220: */
00872         }
00873 
00874 /*        Compare condition numbers for eigenvectors */
00875 /*        taking their condition numbers into account */
00876 
00877         result[10] = 0.f;
00878         eps = dmax(5.9605e-8f,ulp);
00879 /* Computing MAX */
00880         r__1 = (real) (*n) * eps * abnrm;
00881         v = dmax(r__1,smlnum);
00882         if (abnrm == 0.f) {
00883             v = 1.f;
00884         }
00885         i__1 = *n;
00886         for (i__ = 1; i__ <= i__1; ++i__) {
00887             if (v > rcondv[i__] * rconde[i__]) {
00888                 tol = rcondv[i__];
00889             } else {
00890                 tol = v / rconde[i__];
00891             }
00892             if (v > rcdvin[i__] * rcdein[i__]) {
00893                 tolin = rcdvin[i__];
00894             } else {
00895                 tolin = v / rcdein[i__];
00896             }
00897 /* Computing MAX */
00898             r__1 = tol, r__2 = smlnum / eps;
00899             tol = dmax(r__1,r__2);
00900 /* Computing MAX */
00901             r__1 = tolin, r__2 = smlnum / eps;
00902             tolin = dmax(r__1,r__2);
00903             if (eps * (rcdvin[i__] - tolin) > rcondv[i__] + tol) {
00904                 vmax = 1.f / eps;
00905             } else if (rcdvin[i__] - tolin > rcondv[i__] + tol) {
00906                 vmax = (rcdvin[i__] - tolin) / (rcondv[i__] + tol);
00907             } else if (rcdvin[i__] + tolin < eps * (rcondv[i__] - tol)) {
00908                 vmax = 1.f / eps;
00909             } else if (rcdvin[i__] + tolin < rcondv[i__] - tol) {
00910                 vmax = (rcondv[i__] - tol) / (rcdvin[i__] + tolin);
00911             } else {
00912                 vmax = 1.f;
00913             }
00914             result[10] = dmax(result[10],vmax);
00915 /* L230: */
00916         }
00917 
00918 /*        Compare condition numbers for eigenvalues */
00919 /*        taking their condition numbers into account */
00920 
00921         result[11] = 0.f;
00922         i__1 = *n;
00923         for (i__ = 1; i__ <= i__1; ++i__) {
00924             if (v > rcondv[i__]) {
00925                 tol = 1.f;
00926             } else {
00927                 tol = v / rcondv[i__];
00928             }
00929             if (v > rcdvin[i__]) {
00930                 tolin = 1.f;
00931             } else {
00932                 tolin = v / rcdvin[i__];
00933             }
00934 /* Computing MAX */
00935             r__1 = tol, r__2 = smlnum / eps;
00936             tol = dmax(r__1,r__2);
00937 /* Computing MAX */
00938             r__1 = tolin, r__2 = smlnum / eps;
00939             tolin = dmax(r__1,r__2);
00940             if (eps * (rcdein[i__] - tolin) > rconde[i__] + tol) {
00941                 vmax = 1.f / eps;
00942             } else if (rcdein[i__] - tolin > rconde[i__] + tol) {
00943                 vmax = (rcdein[i__] - tolin) / (rconde[i__] + tol);
00944             } else if (rcdein[i__] + tolin < eps * (rconde[i__] - tol)) {
00945                 vmax = 1.f / eps;
00946             } else if (rcdein[i__] + tolin < rconde[i__] - tol) {
00947                 vmax = (rconde[i__] - tol) / (rcdein[i__] + tolin);
00948             } else {
00949                 vmax = 1.f;
00950             }
00951             result[11] = dmax(result[11],vmax);
00952 /* L240: */
00953         }
00954 L250:
00955 
00956         ;
00957     }
00958 
00959 
00960     return 0;
00961 
00962 /*     End of CGET23 */
00963 
00964 } /* cget23_ */


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