1 /* slasq3.f -- translated by f2c (version 20061008).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
10 http://www.netlib.org/f2c/libf2c.zip
16 /* Subroutine */ int slasq3_(integer *i0, integer *n0, real *z__, integer *pp,
17 real *dmin__, real *sigma, real *desig, real *qmax, integer *nfail,
18 integer *iter, integer *ndiv, logical *ieee, integer *ttype, real *
19 dmin1, real *dmin2, real *dn, real *dn1, real *dn2, real *g, real *
22 /* System generated locals */
26 /* Builtin functions */
27 double sqrt(doublereal);
35 extern /* Subroutine */ int slasq4_(integer *, integer *, real *, integer
36 *, integer *, real *, real *, real *, real *, real *, real *,
37 real *, integer *, real *), slasq5_(integer *, integer *, real *,
38 integer *, real *, real *, real *, real *, real *, real *, real *,
39 logical *), slasq6_(integer *, integer *, real *, integer *,
40 real *, real *, real *, real *, real *, real *);
41 extern doublereal slamch_(char *);
42 extern logical sisnan_(real *);
45 /* -- LAPACK routine (version 3.2) -- */
47 /* -- Contributed by Osni Marques of the Lawrence Berkeley National -- */
48 /* -- Laboratory and Beresford Parlett of the Univ. of California at -- */
50 /* -- November 2008 -- */
52 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
53 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
55 /* .. Scalar Arguments .. */
57 /* .. Array Arguments .. */
63 /* SLASQ3 checks for deflation, computes a shift (TAU) and calls dqds. */
64 /* In case of failure it changes shifts, and tries again until output */
70 /* I0 (input) INTEGER */
73 /* N0 (input) INTEGER */
76 /* Z (input) REAL array, dimension ( 4*N ) */
77 /* Z holds the qd array. */
79 /* PP (input/output) INTEGER */
80 /* PP=0 for ping, PP=1 for pong. */
81 /* PP=2 indicates that flipping was applied to the Z array */
82 /* and that the initial tests for deflation should not be */
85 /* DMIN (output) REAL */
86 /* Minimum value of d. */
88 /* SIGMA (output) REAL */
89 /* Sum of shifts used in current segment. */
91 /* DESIG (input/output) REAL */
92 /* Lower order part of SIGMA */
94 /* QMAX (input) REAL */
95 /* Maximum value of q. */
97 /* NFAIL (output) INTEGER */
98 /* Number of times shift was too big. */
100 /* ITER (output) INTEGER */
101 /* Number of iterations. */
103 /* NDIV (output) INTEGER */
104 /* Number of divisions. */
106 /* IEEE (input) LOGICAL */
107 /* Flag for IEEE or non IEEE arithmetic (passed to SLASQ5). */
109 /* TTYPE (input/output) INTEGER */
112 /* DMIN1, DMIN2, DN, DN1, DN2, G, TAU (input/output) REAL */
113 /* These are passed as arguments in order to save their values */
114 /* between calls to SLASQ3. */
116 /* ===================================================================== */
118 /* .. Parameters .. */
120 /* .. Local Scalars .. */
122 /* .. External Subroutines .. */
124 /* .. External Function .. */
126 /* .. Intrinsic Functions .. */
128 /* .. Executable Statements .. */
130 /* Parameter adjustments */
135 eps = slamch_("Precision");
137 /* Computing 2nd power */
141 /* Check for deflation. */
151 nn = (*n0 << 2) + *pp;
152 if (*n0 == *i0 + 1) {
156 /* Check whether E(N0-1) is negligible, 1 eigenvalue. */
158 if (z__[nn - 5] > tol2 * (*sigma + z__[nn - 3]) && z__[nn - (*pp << 1) -
159 4] > tol2 * z__[nn - 7]) {
165 z__[(*n0 << 2) - 3] = z__[(*n0 << 2) + *pp - 3] + *sigma;
169 /* Check whether E(N0-2) is negligible, 2 eigenvalues. */
173 if (z__[nn - 9] > tol2 * *sigma && z__[nn - (*pp << 1) - 8] > tol2 * z__[
180 if (z__[nn - 3] > z__[nn - 7]) {
182 z__[nn - 3] = z__[nn - 7];
185 if (z__[nn - 5] > z__[nn - 3] * tol2) {
186 t = (z__[nn - 7] - z__[nn - 3] + z__[nn - 5]) * .5f;
187 s = z__[nn - 3] * (z__[nn - 5] / t);
189 s = z__[nn - 3] * (z__[nn - 5] / (t * (sqrt(s / t + 1.f) + 1.f)));
191 s = z__[nn - 3] * (z__[nn - 5] / (t + sqrt(t) * sqrt(t + s)));
193 t = z__[nn - 7] + (s + z__[nn - 5]);
194 z__[nn - 3] *= z__[nn - 7] / t;
197 z__[(*n0 << 2) - 7] = z__[nn - 7] + *sigma;
198 z__[(*n0 << 2) - 3] = z__[nn - 3] + *sigma;
207 /* Reverse the qd-array, if warranted. */
209 if (*dmin__ <= 0.f || *n0 < n0in) {
210 if (z__[(*i0 << 2) + *pp - 3] * 1.5f < z__[(*n0 << 2) + *pp - 3]) {
211 ipn4 = *i0 + *n0 << 2;
212 i__1 = *i0 + *n0 - 1 << 1;
213 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
215 z__[j4 - 3] = z__[ipn4 - j4 - 3];
216 z__[ipn4 - j4 - 3] = temp;
218 z__[j4 - 2] = z__[ipn4 - j4 - 2];
219 z__[ipn4 - j4 - 2] = temp;
221 z__[j4 - 1] = z__[ipn4 - j4 - 5];
222 z__[ipn4 - j4 - 5] = temp;
224 z__[j4] = z__[ipn4 - j4 - 4];
225 z__[ipn4 - j4 - 4] = temp;
228 if (*n0 - *i0 <= 4) {
229 z__[(*n0 << 2) + *pp - 1] = z__[(*i0 << 2) + *pp - 1];
230 z__[(*n0 << 2) - *pp] = z__[(*i0 << 2) - *pp];
233 r__1 = *dmin2, r__2 = z__[(*n0 << 2) + *pp - 1];
234 *dmin2 = dmin(r__1,r__2);
236 r__1 = z__[(*n0 << 2) + *pp - 1], r__2 = z__[(*i0 << 2) + *pp - 1]
237 , r__1 = min(r__1,r__2), r__2 = z__[(*i0 << 2) + *pp + 3];
238 z__[(*n0 << 2) + *pp - 1] = dmin(r__1,r__2);
240 r__1 = z__[(*n0 << 2) - *pp], r__2 = z__[(*i0 << 2) - *pp], r__1 =
241 min(r__1,r__2), r__2 = z__[(*i0 << 2) - *pp + 4];
242 z__[(*n0 << 2) - *pp] = dmin(r__1,r__2);
244 r__1 = *qmax, r__2 = z__[(*i0 << 2) + *pp - 3], r__1 = max(r__1,
245 r__2), r__2 = z__[(*i0 << 2) + *pp + 1];
246 *qmax = dmax(r__1,r__2);
251 /* Choose a shift. */
253 slasq4_(i0, n0, &z__[1], pp, &n0in, dmin__, dmin1, dmin2, dn, dn1, dn2,
256 /* Call dqds until DMIN > 0. */
260 slasq5_(i0, n0, &z__[1], pp, tau, dmin__, dmin1, dmin2, dn, dn1, dn2,
263 *ndiv += *n0 - *i0 + 2;
268 if (*dmin__ >= 0.f && *dmin1 > 0.f) {
274 } else if (*dmin__ < 0.f && *dmin1 > 0.f && z__[(*n0 - 1 << 2) - *pp] <
275 tol * (*sigma + *dn1) && dabs(*dn) < tol * *sigma) {
277 /* Convergence hidden by negative DN. */
279 z__[(*n0 - 1 << 2) - *pp + 2] = 0.f;
282 } else if (*dmin__ < 0.f) {
284 /* TAU too big. Select new TAU and try again. */
289 /* Failed twice. Play it safe. */
292 } else if (*dmin1 > 0.f) {
294 /* Late failure. Gives excellent shift. */
296 *tau = (*tau + *dmin__) * (1.f - eps * 2.f);
300 /* Early failure. Divide by 4. */
306 } else if (sisnan_(dmin__)) {
318 /* Possible underflow. Play it safe. */
323 /* Risk of underflow. */
326 slasq6_(i0, n0, &z__[1], pp, dmin__, dmin1, dmin2, dn, dn1, dn2);
327 *ndiv += *n0 - *i0 + 2;
335 *desig -= t - *sigma;
338 *desig = *sigma - (t - *tau) + *desig;