1 /* ./src_f77/slaed6.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
slaed6_(integer * kniter,logical * orgati,real * rho,real * d__,real * z__,real * finit,real * tau,integer * info)8 /* Subroutine */ int slaed6_(integer *kniter, logical *orgati, real *rho,
9 real *d__, real *z__, real *finit, real *tau, integer *info)
10 {
11 /* Initialized data */
12
13 static logical first = TRUE_;
14
15 /* System generated locals */
16 integer i__1;
17 real r__1, r__2, r__3, r__4;
18
19 /* Builtin functions */
20 double sqrt(doublereal), log(doublereal), pow_ri(real *, integer *);
21
22 /* Local variables */
23 static real a, b, c__, f;
24 static integer i__;
25 static real fc, df, ddf, eta, eps, base;
26 static integer iter;
27 static real temp, temp1, temp2, temp3, temp4;
28 static logical scale;
29 static integer niter;
30 static real small1, small2, sminv1, sminv2, dscale[3], sclfac;
31 extern doublereal slamch_(char *, ftnlen);
32 static real zscale[3], erretm, sclinv;
33
34
35 /* -- LAPACK routine (version 3.0) -- */
36 /* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, */
37 /* Courant Institute, NAG Ltd., and Rice University */
38 /* June 30, 1999 */
39
40 /* .. Scalar Arguments .. */
41 /* .. */
42 /* .. Array Arguments .. */
43 /* .. */
44
45 /* Purpose */
46 /* ======= */
47
48 /* SLAED6 computes the positive or negative root (closest to the origin) */
49 /* of */
50 /* z(1) z(2) z(3) */
51 /* f(x) = rho + --------- + ---------- + --------- */
52 /* d(1)-x d(2)-x d(3)-x */
53
54 /* It is assumed that */
55
56 /* if ORGATI = .true. the root is between d(2) and d(3); */
57 /* otherwise it is between d(1) and d(2) */
58
59 /* This routine will be called by SLAED4 when necessary. In most cases, */
60 /* the root sought is the smallest in magnitude, though it might not be */
61 /* in some extremely rare situations. */
62
63 /* Arguments */
64 /* ========= */
65
66 /* KNITER (input) INTEGER */
67 /* Refer to SLAED4 for its significance. */
68
69 /* ORGATI (input) LOGICAL */
70 /* If ORGATI is true, the needed root is between d(2) and */
71 /* d(3); otherwise it is between d(1) and d(2). See */
72 /* SLAED4 for further details. */
73
74 /* RHO (input) REAL */
75 /* Refer to the equation f(x) above. */
76
77 /* D (input) REAL array, dimension (3) */
78 /* D satisfies d(1) < d(2) < d(3). */
79
80 /* Z (input) REAL array, dimension (3) */
81 /* Each of the elements in z must be positive. */
82
83 /* FINIT (input) REAL */
84 /* The value of f at 0. It is more accurate than the one */
85 /* evaluated inside this routine (if someone wants to do */
86 /* so). */
87
88 /* TAU (output) REAL */
89 /* The root of the equation f(x). */
90
91 /* INFO (output) INTEGER */
92 /* = 0: successful exit */
93 /* > 0: if INFO = 1, failure to converge */
94
95 /* Further Details */
96 /* =============== */
97
98 /* Based on contributions by */
99 /* Ren-Cang Li, Computer Science Division, University of California */
100 /* at Berkeley, USA */
101
102 /* ===================================================================== */
103
104 /* .. Parameters .. */
105 /* .. */
106 /* .. External Functions .. */
107 /* .. */
108 /* .. Local Arrays .. */
109 /* .. */
110 /* .. Local Scalars .. */
111 /* .. */
112 /* .. Save statement .. */
113 /* .. */
114 /* .. Intrinsic Functions .. */
115 /* .. */
116 /* .. Data statements .. */
117 /* Parameter adjustments */
118 --z__;
119 --d__;
120
121 /* Function Body */
122 /* .. */
123 /* .. Executable Statements .. */
124
125 *info = 0;
126
127 niter = 1;
128 *tau = 0.f;
129 if (*kniter == 2) {
130 if (*orgati) {
131 temp = (d__[3] - d__[2]) / 2.f;
132 c__ = *rho + z__[1] / (d__[1] - d__[2] - temp);
133 a = c__ * (d__[2] + d__[3]) + z__[2] + z__[3];
134 b = c__ * d__[2] * d__[3] + z__[2] * d__[3] + z__[3] * d__[2];
135 } else {
136 temp = (d__[1] - d__[2]) / 2.f;
137 c__ = *rho + z__[3] / (d__[3] - d__[2] - temp);
138 a = c__ * (d__[1] + d__[2]) + z__[1] + z__[2];
139 b = c__ * d__[1] * d__[2] + z__[1] * d__[2] + z__[2] * d__[1];
140 }
141 /* Computing MAX */
142 r__1 = dabs(a), r__2 = dabs(b), r__1 = max(r__1,r__2), r__2 = dabs(
143 c__);
144 temp = dmax(r__1,r__2);
145 a /= temp;
146 b /= temp;
147 c__ /= temp;
148 if (c__ == 0.f) {
149 *tau = b / a;
150 } else if (a <= 0.f) {
151 *tau = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) / (
152 c__ * 2.f);
153 } else {
154 *tau = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
155 r__1))));
156 }
157 temp = *rho + z__[1] / (d__[1] - *tau) + z__[2] / (d__[2] - *tau) +
158 z__[3] / (d__[3] - *tau);
159 if (dabs(*finit) <= dabs(temp)) {
160 *tau = 0.f;
161 }
162 }
163
164 /* On first call to routine, get machine parameters for */
165 /* possible scaling to avoid overflow */
166
167 if (first) {
168 eps = slamch_("Epsilon", (ftnlen)7);
169 base = slamch_("Base", (ftnlen)4);
170 i__1 = (integer) (log(slamch_("SafMin", (ftnlen)6)) / log(base) / 3.f)
171 ;
172 small1 = pow_ri(&base, &i__1);
173 sminv1 = 1.f / small1;
174 small2 = small1 * small1;
175 sminv2 = sminv1 * sminv1;
176 first = FALSE_;
177 }
178
179 /* Determine if scaling of inputs necessary to avoid overflow */
180 /* when computing 1/TEMP**3 */
181
182 if (*orgati) {
183 /* Computing MIN */
184 r__3 = (r__1 = d__[2] - *tau, dabs(r__1)), r__4 = (r__2 = d__[3] - *
185 tau, dabs(r__2));
186 temp = dmin(r__3,r__4);
187 } else {
188 /* Computing MIN */
189 r__3 = (r__1 = d__[1] - *tau, dabs(r__1)), r__4 = (r__2 = d__[2] - *
190 tau, dabs(r__2));
191 temp = dmin(r__3,r__4);
192 }
193 scale = FALSE_;
194 if (temp <= small1) {
195 scale = TRUE_;
196 if (temp <= small2) {
197
198 /* Scale up by power of radix nearest 1/SAFMIN**(2/3) */
199
200 sclfac = sminv2;
201 sclinv = small2;
202 } else {
203
204 /* Scale up by power of radix nearest 1/SAFMIN**(1/3) */
205
206 sclfac = sminv1;
207 sclinv = small1;
208 }
209
210 /* Scaling up safe because D, Z, TAU scaled elsewhere to be O(1) */
211
212 for (i__ = 1; i__ <= 3; ++i__) {
213 dscale[i__ - 1] = d__[i__] * sclfac;
214 zscale[i__ - 1] = z__[i__] * sclfac;
215 /* L10: */
216 }
217 *tau *= sclfac;
218 } else {
219
220 /* Copy D and Z to DSCALE and ZSCALE */
221
222 for (i__ = 1; i__ <= 3; ++i__) {
223 dscale[i__ - 1] = d__[i__];
224 zscale[i__ - 1] = z__[i__];
225 /* L20: */
226 }
227 }
228
229 fc = 0.f;
230 df = 0.f;
231 ddf = 0.f;
232 for (i__ = 1; i__ <= 3; ++i__) {
233 temp = 1.f / (dscale[i__ - 1] - *tau);
234 temp1 = zscale[i__ - 1] * temp;
235 temp2 = temp1 * temp;
236 temp3 = temp2 * temp;
237 fc += temp1 / dscale[i__ - 1];
238 df += temp2;
239 ddf += temp3;
240 /* L30: */
241 }
242 f = *finit + *tau * fc;
243
244 if (dabs(f) <= 0.f) {
245 goto L60;
246 }
247
248 /* Iteration begins */
249
250 /* It is not hard to see that */
251
252 /* 1) Iterations will go up monotonically */
253 /* if FINIT < 0; */
254
255 /* 2) Iterations will go down monotonically */
256 /* if FINIT > 0. */
257
258 iter = niter + 1;
259
260 for (niter = iter; niter <= 20; ++niter) {
261
262 if (*orgati) {
263 temp1 = dscale[1] - *tau;
264 temp2 = dscale[2] - *tau;
265 } else {
266 temp1 = dscale[0] - *tau;
267 temp2 = dscale[1] - *tau;
268 }
269 a = (temp1 + temp2) * f - temp1 * temp2 * df;
270 b = temp1 * temp2 * f;
271 c__ = f - (temp1 + temp2) * df + temp1 * temp2 * ddf;
272 /* Computing MAX */
273 r__1 = dabs(a), r__2 = dabs(b), r__1 = max(r__1,r__2), r__2 = dabs(
274 c__);
275 temp = dmax(r__1,r__2);
276 a /= temp;
277 b /= temp;
278 c__ /= temp;
279 if (c__ == 0.f) {
280 eta = b / a;
281 } else if (a <= 0.f) {
282 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) / (
283 c__ * 2.f);
284 } else {
285 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
286 r__1))));
287 }
288 if (f * eta >= 0.f) {
289 eta = -f / df;
290 }
291
292 temp = eta + *tau;
293 if (*orgati) {
294 if (eta > 0.f && temp >= dscale[2]) {
295 eta = (dscale[2] - *tau) / 2.f;
296 }
297 if (eta < 0.f && temp <= dscale[1]) {
298 eta = (dscale[1] - *tau) / 2.f;
299 }
300 } else {
301 if (eta > 0.f && temp >= dscale[1]) {
302 eta = (dscale[1] - *tau) / 2.f;
303 }
304 if (eta < 0.f && temp <= dscale[0]) {
305 eta = (dscale[0] - *tau) / 2.f;
306 }
307 }
308 *tau += eta;
309
310 fc = 0.f;
311 erretm = 0.f;
312 df = 0.f;
313 ddf = 0.f;
314 for (i__ = 1; i__ <= 3; ++i__) {
315 temp = 1.f / (dscale[i__ - 1] - *tau);
316 temp1 = zscale[i__ - 1] * temp;
317 temp2 = temp1 * temp;
318 temp3 = temp2 * temp;
319 temp4 = temp1 / dscale[i__ - 1];
320 fc += temp4;
321 erretm += dabs(temp4);
322 df += temp2;
323 ddf += temp3;
324 /* L40: */
325 }
326 f = *finit + *tau * fc;
327 erretm = (dabs(*finit) + dabs(*tau) * erretm) * 8.f + dabs(*tau) * df;
328 if (dabs(f) <= eps * erretm) {
329 goto L60;
330 }
331 /* L50: */
332 }
333 *info = 1;
334 L60:
335
336 /* Undo scaling */
337
338 if (scale) {
339 *tau *= sclinv;
340 }
341 return 0;
342
343 /* End of SLAED6 */
344
345 } /* slaed6_ */
346
347