#include "blaswrap.h" /* cget37.f -- translated by f2c (version 20061008). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib; on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */ #include "f2c.h" /* Table of constant values */ static integer c__3 = 3; static integer c__1 = 1; static integer c__6 = 6; static integer c__4 = 4; static integer c__20 = 20; static integer c__1200 = 1200; static integer c__0 = 0; /* Subroutine */ int cget37_(real *rmax, integer *lmax, integer *ninfo, integer *knt, integer *nin) { /* System generated locals */ integer i__1, i__2, i__3; real r__1, r__2; /* Builtin functions */ double sqrt(doublereal); integer s_rsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), e_rsle(void); double r_imag(complex *); /* Local variables */ static integer i__, j, m, n; static real s[20]; static complex t[400] /* was [20][20] */; static real v; static complex w[20], le[400] /* was [20][20] */, re[400] /* was [20][20] */; static real val[3], dum[1], eps, sep[20], sin__[20], tol; static complex tmp[400] /* was [20][20] */; static integer icmp; static complex cdum[1]; static integer iscl, info, lcmp[3], kmin; static real wiin[20], vmin, vmax, tnrm; static integer isrt; static real wrin[20], vmul, stmp[20]; static complex work[1200], wtmp[20]; static real wsrt[20]; extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *); static real vcmin; extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, complex *, integer *); static real sepin[20], tolin; extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, integer *); static real rwork[40]; extern /* Subroutine */ int slabad_(real *, real *); extern doublereal clange_(char *, integer *, integer *, complex *, integer *, real *); extern /* Subroutine */ int cgehrd_(integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, integer *); extern doublereal slamch_(char *); extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer *), clacpy_(char *, integer *, integer *, complex *, integer *, complex *, integer *); static logical select[20]; static real bignum; extern /* Subroutine */ int chseqr_(char *, char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, complex *, integer *, integer *), ctrevc_(char *, char *, logical *, integer *, complex *, integer *, complex *, integer *, complex *, integer *, integer *, integer *, complex *, real *, integer *), ctrsna_(char *, char *, logical *, integer *, complex *, integer *, complex *, integer *, complex *, integer *, real *, real *, integer *, integer *, complex *, integer *, real *, integer *); static real septmp[20], smlnum; /* Fortran I/O blocks */ static cilist io___5 = { 0, 0, 0, 0, 0 }; static cilist io___9 = { 0, 0, 0, 0, 0 }; static cilist io___12 = { 0, 0, 0, 0, 0 }; /* -- LAPACK test routine (version 3.1) -- Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. November 2006 Purpose ======= CGET37 tests CTRSNA, a routine for estimating condition numbers of eigenvalues and/or right eigenvectors of a matrix. The test matrices are read from a file with logical unit number NIN. Arguments ========== RMAX (output) REAL array, dimension (3) Value of the largest test ratio. RMAX(1) = largest ratio comparing different calls to CTRSNA RMAX(2) = largest error in reciprocal condition numbers taking their conditioning into account RMAX(3) = largest error in reciprocal condition numbers not taking their conditioning into account (may be larger than RMAX(2)) LMAX (output) INTEGER array, dimension (3) LMAX(i) is example number where largest test ratio RMAX(i) is achieved. Also: If CGEHRD returns INFO nonzero on example i, LMAX(1)=i If CHSEQR returns INFO nonzero on example i, LMAX(2)=i If CTRSNA returns INFO nonzero on example i, LMAX(3)=i NINFO (output) INTEGER array, dimension (3) NINFO(1) = No. of times CGEHRD returned INFO nonzero NINFO(2) = No. of times CHSEQR returned INFO nonzero NINFO(3) = No. of times CTRSNA returned INFO nonzero KNT (output) INTEGER Total number of examples tested. NIN (input) INTEGER Input logical unit number ===================================================================== Parameter adjustments */ --ninfo; --lmax; --rmax; /* Function Body */ eps = slamch_("P"); smlnum = slamch_("S") / eps; bignum = 1.f / smlnum; slabad_(&smlnum, &bignum); /* EPSIN = 2**(-24) = precision to which input data computed */ eps = dmax(eps,5.9605e-8f); rmax[1] = 0.f; rmax[2] = 0.f; rmax[3] = 0.f; lmax[1] = 0; lmax[2] = 0; lmax[3] = 0; *knt = 0; ninfo[1] = 0; ninfo[2] = 0; ninfo[3] = 0; val[0] = sqrt(smlnum); val[1] = 1.f; val[2] = sqrt(bignum); /* Read input data until N=0. Assume input eigenvalues are sorted lexicographically (increasing by real part if ISRT = 0, increasing by imaginary part if ISRT = 1) */ L10: io___5.ciunit = *nin; s_rsle(&io___5); do_lio(&c__3, &c__1, (char *)&n, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&isrt, (ftnlen)sizeof(integer)); e_rsle(); if (n == 0) { return 0; } i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { io___9.ciunit = *nin; s_rsle(&io___9); i__2 = n; for (j = 1; j <= i__2; ++j) { do_lio(&c__6, &c__1, (char *)&tmp[i__ + j * 20 - 21], (ftnlen) sizeof(complex)); } e_rsle(); /* L20: */ } i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { io___12.ciunit = *nin; s_rsle(&io___12); do_lio(&c__4, &c__1, (char *)&wrin[i__ - 1], (ftnlen)sizeof(real)); do_lio(&c__4, &c__1, (char *)&wiin[i__ - 1], (ftnlen)sizeof(real)); do_lio(&c__4, &c__1, (char *)&sin__[i__ - 1], (ftnlen)sizeof(real)); do_lio(&c__4, &c__1, (char *)&sepin[i__ - 1], (ftnlen)sizeof(real)); e_rsle(); /* L30: */ } tnrm = clange_("M", &n, &n, tmp, &c__20, rwork); for (iscl = 1; iscl <= 3; ++iscl) { /* Scale input matrix */ ++(*knt); clacpy_("F", &n, &n, tmp, &c__20, t, &c__20); vmul = val[iscl - 1]; i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { csscal_(&n, &vmul, &t[i__ * 20 - 20], &c__1); /* L40: */ } if (tnrm == 0.f) { vmul = 1.f; } /* Compute eigenvalues and eigenvectors */ i__1 = 1200 - n; cgehrd_(&n, &c__1, &n, t, &c__20, work, &work[n], &i__1, &info); if (info != 0) { lmax[1] = *knt; ++ninfo[1]; goto L260; } i__1 = n - 2; for (j = 1; j <= i__1; ++j) { i__2 = n; for (i__ = j + 2; i__ <= i__2; ++i__) { i__3 = i__ + j * 20 - 21; t[i__3].r = 0.f, t[i__3].i = 0.f; /* L50: */ } /* L60: */ } /* Compute Schur form */ chseqr_("S", "N", &n, &c__1, &n, t, &c__20, w, cdum, &c__1, work, & c__1200, &info); if (info != 0) { lmax[2] = *knt; ++ninfo[2]; goto L260; } /* Compute eigenvectors */ i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { select[i__ - 1] = TRUE_; /* L70: */ } ctrevc_("B", "A", select, &n, t, &c__20, le, &c__20, re, &c__20, &n, & m, work, rwork, &info); /* Compute condition numbers */ ctrsna_("B", "A", select, &n, t, &c__20, le, &c__20, re, &c__20, s, sep, &n, &m, work, &n, rwork, &info); if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } /* Sort eigenvalues and condition numbers lexicographically to compare with inputs */ ccopy_(&n, w, &c__1, wtmp, &c__1); if (isrt == 0) { /* Sort by increasing real part */ i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = i__ - 1; wsrt[i__ - 1] = w[i__2].r; /* L80: */ } } else { /* Sort by increasing imaginary part */ i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { wsrt[i__ - 1] = r_imag(&w[i__ - 1]); /* L90: */ } } scopy_(&n, s, &c__1, stmp, &c__1); scopy_(&n, sep, &c__1, septmp, &c__1); r__1 = 1.f / vmul; sscal_(&n, &r__1, septmp, &c__1); i__1 = n - 1; for (i__ = 1; i__ <= i__1; ++i__) { kmin = i__; vmin = wsrt[i__ - 1]; i__2 = n; for (j = i__ + 1; j <= i__2; ++j) { if (wsrt[j - 1] < vmin) { kmin = j; vmin = wsrt[j - 1]; } /* L100: */ } wsrt[kmin - 1] = wsrt[i__ - 1]; wsrt[i__ - 1] = vmin; i__2 = i__ - 1; vcmin = wtmp[i__2].r; i__2 = i__ - 1; i__3 = kmin - 1; wtmp[i__2].r = w[i__3].r, wtmp[i__2].i = w[i__3].i; i__2 = kmin - 1; wtmp[i__2].r = vcmin, wtmp[i__2].i = 0.f; vmin = stmp[kmin - 1]; stmp[kmin - 1] = stmp[i__ - 1]; stmp[i__ - 1] = vmin; vmin = septmp[kmin - 1]; septmp[kmin - 1] = septmp[i__ - 1]; septmp[i__ - 1] = vmin; /* L110: */ } /* Compare condition numbers for eigenvalues taking their condition numbers into account Computing MAX */ r__1 = (real) n * 2.f * eps * tnrm; v = dmax(r__1,smlnum); if (tnrm == 0.f) { v = 1.f; } i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (v > septmp[i__ - 1]) { tol = 1.f; } else { tol = v / septmp[i__ - 1]; } if (v > sepin[i__ - 1]) { tolin = 1.f; } else { tolin = v / sepin[i__ - 1]; } /* Computing MAX */ r__1 = tol, r__2 = smlnum / eps; tol = dmax(r__1,r__2); /* Computing MAX */ r__1 = tolin, r__2 = smlnum / eps; tolin = dmax(r__1,r__2); if (eps * (sin__[i__ - 1] - tolin) > stmp[i__ - 1] + tol) { vmax = 1.f / eps; } else if (sin__[i__ - 1] - tolin > stmp[i__ - 1] + tol) { vmax = (sin__[i__ - 1] - tolin) / (stmp[i__ - 1] + tol); } else if (sin__[i__ - 1] + tolin < eps * (stmp[i__ - 1] - tol)) { vmax = 1.f / eps; } else if (sin__[i__ - 1] + tolin < stmp[i__ - 1] - tol) { vmax = (stmp[i__ - 1] - tol) / (sin__[i__ - 1] + tolin); } else { vmax = 1.f; } if (vmax > rmax[2]) { rmax[2] = vmax; if (ninfo[2] == 0) { lmax[2] = *knt; } } /* L120: */ } /* Compare condition numbers for eigenvectors taking their condition numbers into account */ i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (v > septmp[i__ - 1] * stmp[i__ - 1]) { tol = septmp[i__ - 1]; } else { tol = v / stmp[i__ - 1]; } if (v > sepin[i__ - 1] * sin__[i__ - 1]) { tolin = sepin[i__ - 1]; } else { tolin = v / sin__[i__ - 1]; } /* Computing MAX */ r__1 = tol, r__2 = smlnum / eps; tol = dmax(r__1,r__2); /* Computing MAX */ r__1 = tolin, r__2 = smlnum / eps; tolin = dmax(r__1,r__2); if (eps * (sepin[i__ - 1] - tolin) > septmp[i__ - 1] + tol) { vmax = 1.f / eps; } else if (sepin[i__ - 1] - tolin > septmp[i__ - 1] + tol) { vmax = (sepin[i__ - 1] - tolin) / (septmp[i__ - 1] + tol); } else if (sepin[i__ - 1] + tolin < eps * (septmp[i__ - 1] - tol)) { vmax = 1.f / eps; } else if (sepin[i__ - 1] + tolin < septmp[i__ - 1] - tol) { vmax = (septmp[i__ - 1] - tol) / (sepin[i__ - 1] + tolin); } else { vmax = 1.f; } if (vmax > rmax[2]) { rmax[2] = vmax; if (ninfo[2] == 0) { lmax[2] = *knt; } } /* L130: */ } /* Compare condition numbers for eigenvalues without taking their condition numbers into account */ i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (sin__[i__ - 1] <= (real) (n << 1) * eps && stmp[i__ - 1] <= ( real) (n << 1) * eps) { vmax = 1.f; } else if (eps * sin__[i__ - 1] > stmp[i__ - 1]) { vmax = 1.f / eps; } else if (sin__[i__ - 1] > stmp[i__ - 1]) { vmax = sin__[i__ - 1] / stmp[i__ - 1]; } else if (sin__[i__ - 1] < eps * stmp[i__ - 1]) { vmax = 1.f / eps; } else if (sin__[i__ - 1] < stmp[i__ - 1]) { vmax = stmp[i__ - 1] / sin__[i__ - 1]; } else { vmax = 1.f; } if (vmax > rmax[3]) { rmax[3] = vmax; if (ninfo[3] == 0) { lmax[3] = *knt; } } /* L140: */ } /* Compare condition numbers for eigenvectors without taking their condition numbers into account */ i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (sepin[i__ - 1] <= v && septmp[i__ - 1] <= v) { vmax = 1.f; } else if (eps * sepin[i__ - 1] > septmp[i__ - 1]) { vmax = 1.f / eps; } else if (sepin[i__ - 1] > septmp[i__ - 1]) { vmax = sepin[i__ - 1] / septmp[i__ - 1]; } else if (sepin[i__ - 1] < eps * septmp[i__ - 1]) { vmax = 1.f / eps; } else if (sepin[i__ - 1] < septmp[i__ - 1]) { vmax = septmp[i__ - 1] / sepin[i__ - 1]; } else { vmax = 1.f; } if (vmax > rmax[3]) { rmax[3] = vmax; if (ninfo[3] == 0) { lmax[3] = *knt; } } /* L150: */ } /* Compute eigenvalue condition numbers only and compare */ vmax = 0.f; dum[0] = -1.f; scopy_(&n, dum, &c__0, stmp, &c__1); scopy_(&n, dum, &c__0, septmp, &c__1); ctrsna_("E", "A", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp, septmp, &n, &m, work, &n, rwork, &info) ; if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (stmp[i__ - 1] != s[i__ - 1]) { vmax = 1.f / eps; } if (septmp[i__ - 1] != dum[0]) { vmax = 1.f / eps; } /* L160: */ } /* Compute eigenvector condition numbers only and compare */ scopy_(&n, dum, &c__0, stmp, &c__1); scopy_(&n, dum, &c__0, septmp, &c__1); ctrsna_("V", "A", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp, septmp, &n, &m, work, &n, rwork, &info) ; if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (stmp[i__ - 1] != dum[0]) { vmax = 1.f / eps; } if (septmp[i__ - 1] != sep[i__ - 1]) { vmax = 1.f / eps; } /* L170: */ } /* Compute all condition numbers using SELECT and compare */ i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { select[i__ - 1] = TRUE_; /* L180: */ } scopy_(&n, dum, &c__0, stmp, &c__1); scopy_(&n, dum, &c__0, septmp, &c__1); ctrsna_("B", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp, septmp, &n, &m, work, &n, rwork, &info) ; if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (septmp[i__ - 1] != sep[i__ - 1]) { vmax = 1.f / eps; } if (stmp[i__ - 1] != s[i__ - 1]) { vmax = 1.f / eps; } /* L190: */ } /* Compute eigenvalue condition numbers using SELECT and compare */ scopy_(&n, dum, &c__0, stmp, &c__1); scopy_(&n, dum, &c__0, septmp, &c__1); ctrsna_("E", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp, septmp, &n, &m, work, &n, rwork, &info) ; if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (stmp[i__ - 1] != s[i__ - 1]) { vmax = 1.f / eps; } if (septmp[i__ - 1] != dum[0]) { vmax = 1.f / eps; } /* L200: */ } /* Compute eigenvector condition numbers using SELECT and compare */ scopy_(&n, dum, &c__0, stmp, &c__1); scopy_(&n, dum, &c__0, septmp, &c__1); ctrsna_("V", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp, septmp, &n, &m, work, &n, rwork, &info) ; if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { if (stmp[i__ - 1] != dum[0]) { vmax = 1.f / eps; } if (septmp[i__ - 1] != sep[i__ - 1]) { vmax = 1.f / eps; } /* L210: */ } if (vmax > rmax[1]) { rmax[1] = vmax; if (ninfo[1] == 0) { lmax[1] = *knt; } } /* Select second and next to last eigenvalues */ i__1 = n; for (i__ = 1; i__ <= i__1; ++i__) { select[i__ - 1] = FALSE_; /* L220: */ } icmp = 0; if (n > 1) { icmp = 1; lcmp[0] = 2; select[1] = TRUE_; ccopy_(&n, &re[20], &c__1, re, &c__1); ccopy_(&n, &le[20], &c__1, le, &c__1); } if (n > 3) { icmp = 2; lcmp[1] = n - 1; select[n - 2] = TRUE_; ccopy_(&n, &re[(n - 1) * 20 - 20], &c__1, &re[20], &c__1); ccopy_(&n, &le[(n - 1) * 20 - 20], &c__1, &le[20], &c__1); } /* Compute all selected condition numbers */ scopy_(&icmp, dum, &c__0, stmp, &c__1); scopy_(&icmp, dum, &c__0, septmp, &c__1); ctrsna_("B", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp, septmp, &n, &m, work, &n, rwork, &info) ; if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } i__1 = icmp; for (i__ = 1; i__ <= i__1; ++i__) { j = lcmp[i__ - 1]; if (septmp[i__ - 1] != sep[j - 1]) { vmax = 1.f / eps; } if (stmp[i__ - 1] != s[j - 1]) { vmax = 1.f / eps; } /* L230: */ } /* Compute selected eigenvalue condition numbers */ scopy_(&icmp, dum, &c__0, stmp, &c__1); scopy_(&icmp, dum, &c__0, septmp, &c__1); ctrsna_("E", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp, septmp, &n, &m, work, &n, rwork, &info) ; if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } i__1 = icmp; for (i__ = 1; i__ <= i__1; ++i__) { j = lcmp[i__ - 1]; if (stmp[i__ - 1] != s[j - 1]) { vmax = 1.f / eps; } if (septmp[i__ - 1] != dum[0]) { vmax = 1.f / eps; } /* L240: */ } /* Compute selected eigenvector condition numbers */ scopy_(&icmp, dum, &c__0, stmp, &c__1); scopy_(&icmp, dum, &c__0, septmp, &c__1); ctrsna_("V", "S", select, &n, t, &c__20, le, &c__20, re, &c__20, stmp, septmp, &n, &m, work, &n, rwork, &info) ; if (info != 0) { lmax[3] = *knt; ++ninfo[3]; goto L260; } i__1 = icmp; for (i__ = 1; i__ <= i__1; ++i__) { j = lcmp[i__ - 1]; if (stmp[i__ - 1] != dum[0]) { vmax = 1.f / eps; } if (septmp[i__ - 1] != sep[j - 1]) { vmax = 1.f / eps; } /* L250: */ } if (vmax > rmax[1]) { rmax[1] = vmax; if (ninfo[1] == 0) { lmax[1] = *knt; } } L260: ; } goto L10; /* End of CGET37 */ } /* cget37_ */