1 /* ../netlib/slasd4.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 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 */
3 #include "FLA_f2c.h" /* > \brief \b SLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modification to a positive diagonal matrix. Used by sbdsdc. */
4 /* =========== DOCUMENTATION =========== */
5 /* Online html documentation available at */
6 /* http://www.netlib.org/lapack/explore-html/ */
7 /* > \htmlonly */
8 /* > Download SLASD4 + dependencies */
9 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd4. f"> */
10 /* > [TGZ]</a> */
11 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd4. f"> */
12 /* > [ZIP]</a> */
13 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd4. f"> */
14 /* > [TXT]</a> */
15 /* > \endhtmlonly */
16 /* Definition: */
17 /* =========== */
18 /* SUBROUTINE SLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO ) */
19 /* .. Scalar Arguments .. */
20 /* INTEGER I, INFO, N */
21 /* REAL RHO, SIGMA */
22 /* .. */
23 /* .. Array Arguments .. */
24 /* REAL D( * ), DELTA( * ), WORK( * ), Z( * ) */
25 /* .. */
26 /* > \par Purpose: */
27 /* ============= */
28 /* > */
29 /* > \verbatim */
30 /* > */
31 /* > This subroutine computes the square root of the I-th updated */
32 /* > eigenvalue of a positive symmetric rank-one modification to */
33 /* > a positive diagonal matrix whose entries are given as the squares */
34 /* > of the corresponding entries in the array d, and that */
35 /* > */
36 /* > 0 <= D(i) < D(j) for i < j */
37 /* > */
38 /* > and that RHO > 0. This is arranged by the calling routine, and is */
39 /* > no loss in generality. The rank-one modified system is thus */
40 /* > */
41 /* > diag( D ) * diag( D ) + RHO * Z * Z_transpose. */
42 /* > */
43 /* > where we assume the Euclidean norm of Z is 1. */
44 /* > */
45 /* > The method consists of approximating the rational functions in the */
46 /* > secular equation by simpler interpolating rational functions. */
47 /* > \endverbatim */
48 /* Arguments: */
49 /* ========== */
50 /* > \param[in] N */
51 /* > \verbatim */
52 /* > N is INTEGER */
53 /* > The length of all arrays. */
54 /* > \endverbatim */
55 /* > */
56 /* > \param[in] I */
57 /* > \verbatim */
58 /* > I is INTEGER */
59 /* > The index of the eigenvalue to be computed. 1 <= I <= N. */
60 /* > \endverbatim */
61 /* > */
62 /* > \param[in] D */
63 /* > \verbatim */
64 /* > D is REAL array, dimension ( N ) */
65 /* > The original eigenvalues. It is assumed that they are in */
66 /* > order, 0 <= D(I) < D(J) for I < J. */
67 /* > \endverbatim */
68 /* > */
69 /* > \param[in] Z */
70 /* > \verbatim */
71 /* > Z is REAL array, dimension ( N ) */
72 /* > The components of the updating vector. */
73 /* > \endverbatim */
74 /* > */
75 /* > \param[out] DELTA */
76 /* > \verbatim */
77 /* > DELTA is REAL array, dimension ( N ) */
78 /* > If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th */
79 /* > component. If N = 1, then DELTA(1) = 1. The vector DELTA */
80 /* > contains the information necessary to construct the */
81 /* > (singular) eigenvectors. */
82 /* > \endverbatim */
83 /* > */
84 /* > \param[in] RHO */
85 /* > \verbatim */
86 /* > RHO is REAL */
87 /* > The scalar in the symmetric updating formula. */
88 /* > \endverbatim */
89 /* > */
90 /* > \param[out] SIGMA */
91 /* > \verbatim */
92 /* > SIGMA is REAL */
93 /* > The computed sigma_I, the I-th updated eigenvalue. */
94 /* > \endverbatim */
95 /* > */
96 /* > \param[out] WORK */
97 /* > \verbatim */
98 /* > WORK is REAL array, dimension ( N ) */
99 /* > If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th */
100 /* > component. If N = 1, then WORK( 1 ) = 1. */
101 /* > \endverbatim */
102 /* > */
103 /* > \param[out] INFO */
104 /* > \verbatim */
105 /* > INFO is INTEGER */
106 /* > = 0: successful exit */
107 /* > > 0: if INFO = 1, the updating process failed. */
108 /* > \endverbatim */
109 /* > \par Internal Parameters: */
110 /* ========================= */
111 /* > */
112 /* > \verbatim */
113 /* > Logical variable ORGATI (origin-at-i?) is used for distinguishing */
114 /* > whether D(i) or D(i+1) is treated as the origin. */
115 /* > */
116 /* > ORGATI = .true. origin at i */
117 /* > ORGATI = .false. origin at i+1 */
118 /* > */
119 /* > Logical variable SWTCH3 (switch-for-3-poles?) is for noting */
120 /* > if we are working with THREE poles! */
121 /* > */
122 /* > MAXIT is the maximum number of iterations allowed for each */
123 /* > eigenvalue. */
124 /* > \endverbatim */
125 /* Authors: */
126 /* ======== */
127 /* > \author Univ. of Tennessee */
128 /* > \author Univ. of California Berkeley */
129 /* > \author Univ. of Colorado Denver */
130 /* > \author NAG Ltd. */
131 /* > \date November 2013 */
132 /* > \ingroup auxOTHERauxiliary */
133 /* > \par Contributors: */
134 /* ================== */
135 /* > */
136 /* > Ren-Cang Li, Computer Science Division, University of California */
137 /* > at Berkeley, USA */
138 /* > */
139 /* ===================================================================== */
140 /* Subroutine */
slasd4_(integer * n,integer * i__,real * d__,real * z__,real * delta,real * rho,real * sigma,real * work,integer * info)141 int slasd4_(integer *n, integer *i__, real *d__, real *z__, real *delta, real *rho, real *sigma, real *work, integer *info)
142 {
143 /* System generated locals */
144 integer i__1;
145 real r__1;
146 /* Builtin functions */
147 double sqrt(doublereal);
148 /* Local variables */
149 real a, b, c__;
150 integer j;
151 real w, dd[3];
152 integer ii;
153 real dw, zz[3];
154 integer ip1;
155 real sq2, eta, phi, eps, tau, psi;
156 integer iim1, iip1;
157 real tau2, dphi, sglb, dpsi, sgub;
158 integer iter;
159 real temp, prew, temp1, temp2, dtiim, delsq, dtiip;
160 integer niter;
161 real dtisq;
162 logical swtch;
163 real dtnsq;
164 extern /* Subroutine */
165 int slaed6_(integer *, logical *, real *, real *, real *, real *, real *, integer *);
166 real delsq2;
167 extern /* Subroutine */
168 int slasd5_(integer *, real *, real *, real *, real *, real *, real *);
169 real dtnsq1;
170 logical swtch3;
171 extern real slamch_(char *);
172 logical orgati;
173 real erretm, dtipsq, rhoinv;
174 logical geomavg;
175 /* -- LAPACK auxiliary routine (version 3.5.0) -- */
176 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
177 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
178 /* November 2013 */
179 /* .. Scalar Arguments .. */
180 /* .. */
181 /* .. Array Arguments .. */
182 /* .. */
183 /* ===================================================================== */
184 /* .. Parameters .. */
185 /* .. */
186 /* .. Local Scalars .. */
187 /* .. */
188 /* .. Local Arrays .. */
189 /* .. */
190 /* .. External Subroutines .. */
191 /* .. */
192 /* .. External Functions .. */
193 /* .. */
194 /* .. Intrinsic Functions .. */
195 /* .. */
196 /* .. Executable Statements .. */
197 /* Since this routine is called in an inner loop, we do no argument */
198 /* checking. */
199 /* Quick return for N=1 and 2. */
200 /* Parameter adjustments */
201 --work;
202 --delta;
203 --z__;
204 --d__;
205 /* Function Body */
206 *info = 0;
207 if (*n == 1)
208 {
209 /* Presumably, I=1 upon entry */
210 *sigma = sqrt(d__[1] * d__[1] + *rho * z__[1] * z__[1]);
211 delta[1] = 1.f;
212 work[1] = 1.f;
213 return 0;
214 }
215 if (*n == 2)
216 {
217 slasd5_(i__, &d__[1], &z__[1], &delta[1], rho, sigma, &work[1]);
218 return 0;
219 }
220 /* Compute machine epsilon */
221 eps = slamch_("Epsilon");
222 rhoinv = 1.f / *rho;
223 tau2 = 0.f;
224 /* The case I = N */
225 if (*i__ == *n)
226 {
227 /* Initialize some basic variables */
228 ii = *n - 1;
229 niter = 1;
230 /* Calculate initial guess */
231 temp = *rho / 2.f;
232 /* If ||Z||_2 is not one, then TEMP should be set to */
233 /* RHO * ||Z||_2^2 / TWO */
234 temp1 = temp / (d__[*n] + sqrt(d__[*n] * d__[*n] + temp));
235 i__1 = *n;
236 for (j = 1;
237 j <= i__1;
238 ++j)
239 {
240 work[j] = d__[j] + d__[*n] + temp1;
241 delta[j] = d__[j] - d__[*n] - temp1;
242 /* L10: */
243 }
244 psi = 0.f;
245 i__1 = *n - 2;
246 for (j = 1;
247 j <= i__1;
248 ++j)
249 {
250 psi += z__[j] * z__[j] / (delta[j] * work[j]);
251 /* L20: */
252 }
253 c__ = rhoinv + psi;
254 w = c__ + z__[ii] * z__[ii] / (delta[ii] * work[ii]) + z__[*n] * z__[* n] / (delta[*n] * work[*n]);
255 if (w <= 0.f)
256 {
257 temp1 = sqrt(d__[*n] * d__[*n] + *rho);
258 temp = z__[*n - 1] * z__[*n - 1] / ((d__[*n - 1] + temp1) * (d__[* n] - d__[*n - 1] + *rho / (d__[*n] + temp1))) + z__[*n] * z__[*n] / *rho;
259 /* The following TAU2 is to approximate */
260 /* SIGMA_n^2 - D( N )*D( N ) */
261 if (c__ <= temp)
262 {
263 tau = *rho;
264 }
265 else
266 {
267 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
268 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[* n];
269 b = z__[*n] * z__[*n] * delsq;
270 if (a < 0.f)
271 {
272 tau2 = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
273 }
274 else
275 {
276 tau2 = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
277 }
278 tau = tau2 / (d__[*n] + sqrt(d__[*n] * d__[*n] + tau2));
279 }
280 /* It can be proved that */
281 /* D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU2 <= D(N)^2+RHO */
282 }
283 else
284 {
285 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
286 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
287 b = z__[*n] * z__[*n] * delsq;
288 /* The following TAU2 is to approximate */
289 /* SIGMA_n^2 - D( N )*D( N ) */
290 if (a < 0.f)
291 {
292 tau2 = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
293 }
294 else
295 {
296 tau2 = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
297 }
298 tau = tau2 / (d__[*n] + sqrt(d__[*n] * d__[*n] + tau2));
299 /* It can be proved that */
300 /* D(N)^2 < D(N)^2+TAU2 < SIGMA(N)^2 < D(N)^2+RHO/2 */
301 }
302 /* The following TAU is to approximate SIGMA_n - D( N ) */
303 /* TAU = TAU2 / ( D( N )+SQRT( D( N )*D( N )+TAU2 ) ) */
304 *sigma = d__[*n] + tau;
305 i__1 = *n;
306 for (j = 1;
307 j <= i__1;
308 ++j)
309 {
310 delta[j] = d__[j] - d__[*n] - tau;
311 work[j] = d__[j] + d__[*n] + tau;
312 /* L30: */
313 }
314 /* Evaluate PSI and the derivative DPSI */
315 dpsi = 0.f;
316 psi = 0.f;
317 erretm = 0.f;
318 i__1 = ii;
319 for (j = 1;
320 j <= i__1;
321 ++j)
322 {
323 temp = z__[j] / (delta[j] * work[j]);
324 psi += z__[j] * temp;
325 dpsi += temp * temp;
326 erretm += psi;
327 /* L40: */
328 }
329 erretm = f2c_abs(erretm);
330 /* Evaluate PHI and the derivative DPHI */
331 temp = z__[*n] / (delta[*n] * work[*n]);
332 phi = z__[*n] * temp;
333 dphi = temp * temp;
334 erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv;
335 /* $ + ABS( TAU2 )*( DPSI+DPHI ) */
336 w = rhoinv + phi + psi;
337 /* Test for convergence */
338 if (f2c_abs(w) <= eps * erretm)
339 {
340 goto L240;
341 }
342 /* Calculate the new step */
343 ++niter;
344 dtnsq1 = work[*n - 1] * delta[*n - 1];
345 dtnsq = work[*n] * delta[*n];
346 c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
347 a = (dtnsq + dtnsq1) * w - dtnsq * dtnsq1 * (dpsi + dphi);
348 b = dtnsq * dtnsq1 * w;
349 if (c__ < 0.f)
350 {
351 c__ = f2c_abs(c__);
352 }
353 if (c__ == 0.f)
354 {
355 eta = *rho - *sigma * *sigma;
356 }
357 else if (a >= 0.f)
358 {
359 eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1)))) / ( c__ * 2.f);
360 }
361 else
362 {
363 eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1) )));
364 }
365 /* Note, eta should be positive if w is negative, and */
366 /* eta should be negative otherwise. However, */
367 /* if for some reason caused by roundoff, eta*w > 0, */
368 /* we simply use one Newton step instead. This way */
369 /* will guarantee eta*w < 0. */
370 if (w * eta > 0.f)
371 {
372 eta = -w / (dpsi + dphi);
373 }
374 temp = eta - dtnsq;
375 if (temp > *rho)
376 {
377 eta = *rho + dtnsq;
378 }
379 eta /= *sigma + sqrt(eta + *sigma * *sigma);
380 tau += eta;
381 *sigma += eta;
382 i__1 = *n;
383 for (j = 1;
384 j <= i__1;
385 ++j)
386 {
387 delta[j] -= eta;
388 work[j] += eta;
389 /* L50: */
390 }
391 /* Evaluate PSI and the derivative DPSI */
392 dpsi = 0.f;
393 psi = 0.f;
394 erretm = 0.f;
395 i__1 = ii;
396 for (j = 1;
397 j <= i__1;
398 ++j)
399 {
400 temp = z__[j] / (work[j] * delta[j]);
401 psi += z__[j] * temp;
402 dpsi += temp * temp;
403 erretm += psi;
404 /* L60: */
405 }
406 erretm = f2c_abs(erretm);
407 /* Evaluate PHI and the derivative DPHI */
408 tau2 = work[*n] * delta[*n];
409 temp = z__[*n] / tau2;
410 phi = z__[*n] * temp;
411 dphi = temp * temp;
412 erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv;
413 /* $ + ABS( TAU2 )*( DPSI+DPHI ) */
414 w = rhoinv + phi + psi;
415 /* Main loop to update the values of the array DELTA */
416 iter = niter + 1;
417 for (niter = iter;
418 niter <= 400;
419 ++niter)
420 {
421 /* Test for convergence */
422 if (f2c_abs(w) <= eps * erretm)
423 {
424 goto L240;
425 }
426 /* Calculate the new step */
427 dtnsq1 = work[*n - 1] * delta[*n - 1];
428 dtnsq = work[*n] * delta[*n];
429 c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
430 a = (dtnsq + dtnsq1) * w - dtnsq1 * dtnsq * (dpsi + dphi);
431 b = dtnsq1 * dtnsq * w;
432 if (a >= 0.f)
433 {
434 eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1)))) / (c__ * 2.f);
435 }
436 else
437 {
438 eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs( r__1))));
439 }
440 /* Note, eta should be positive if w is negative, and */
441 /* eta should be negative otherwise. However, */
442 /* if for some reason caused by roundoff, eta*w > 0, */
443 /* we simply use one Newton step instead. This way */
444 /* will guarantee eta*w < 0. */
445 if (w * eta > 0.f)
446 {
447 eta = -w / (dpsi + dphi);
448 }
449 temp = eta - dtnsq;
450 if (temp <= 0.f)
451 {
452 eta /= 2.f;
453 }
454 eta /= *sigma + sqrt(eta + *sigma * *sigma);
455 tau += eta;
456 *sigma += eta;
457 i__1 = *n;
458 for (j = 1;
459 j <= i__1;
460 ++j)
461 {
462 delta[j] -= eta;
463 work[j] += eta;
464 /* L70: */
465 }
466 /* Evaluate PSI and the derivative DPSI */
467 dpsi = 0.f;
468 psi = 0.f;
469 erretm = 0.f;
470 i__1 = ii;
471 for (j = 1;
472 j <= i__1;
473 ++j)
474 {
475 temp = z__[j] / (work[j] * delta[j]);
476 psi += z__[j] * temp;
477 dpsi += temp * temp;
478 erretm += psi;
479 /* L80: */
480 }
481 erretm = f2c_abs(erretm);
482 /* Evaluate PHI and the derivative DPHI */
483 tau2 = work[*n] * delta[*n];
484 temp = z__[*n] / tau2;
485 phi = z__[*n] * temp;
486 dphi = temp * temp;
487 erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv;
488 /* $ + ABS( TAU2 )*( DPSI+DPHI ) */
489 w = rhoinv + phi + psi;
490 /* L90: */
491 }
492 /* Return with INFO = 1, NITER = MAXIT and not converged */
493 *info = 1;
494 goto L240;
495 /* End for the case I = N */
496 }
497 else
498 {
499 /* The case for I < N */
500 niter = 1;
501 ip1 = *i__ + 1;
502 /* Calculate initial guess */
503 delsq = (d__[ip1] - d__[*i__]) * (d__[ip1] + d__[*i__]);
504 delsq2 = delsq / 2.f;
505 sq2 = sqrt((d__[*i__] * d__[*i__] + d__[ip1] * d__[ip1]) / 2.f);
506 temp = delsq2 / (d__[*i__] + sq2);
507 i__1 = *n;
508 for (j = 1;
509 j <= i__1;
510 ++j)
511 {
512 work[j] = d__[j] + d__[*i__] + temp;
513 delta[j] = d__[j] - d__[*i__] - temp;
514 /* L100: */
515 }
516 psi = 0.f;
517 i__1 = *i__ - 1;
518 for (j = 1;
519 j <= i__1;
520 ++j)
521 {
522 psi += z__[j] * z__[j] / (work[j] * delta[j]);
523 /* L110: */
524 }
525 phi = 0.f;
526 i__1 = *i__ + 2;
527 for (j = *n;
528 j >= i__1;
529 --j)
530 {
531 phi += z__[j] * z__[j] / (work[j] * delta[j]);
532 /* L120: */
533 }
534 c__ = rhoinv + psi + phi;
535 w = c__ + z__[*i__] * z__[*i__] / (work[*i__] * delta[*i__]) + z__[ ip1] * z__[ip1] / (work[ip1] * delta[ip1]);
536 geomavg = FALSE_;
537 if (w > 0.f)
538 {
539 /* d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2 */
540 /* We choose d(i) as origin. */
541 orgati = TRUE_;
542 ii = *i__;
543 sglb = 0.f;
544 sgub = delsq2 / (d__[*i__] + sq2);
545 a = c__ * delsq + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
546 b = z__[*i__] * z__[*i__] * delsq;
547 if (a > 0.f)
548 {
549 tau2 = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs( r__1))));
550 }
551 else
552 {
553 tau2 = (a - sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1)))) / (c__ * 2.f);
554 }
555 /* TAU2 now is an estimation of SIGMA^2 - D( I )^2. The */
556 /* following, however, is the corresponding estimation of */
557 /* SIGMA - D( I ). */
558 tau = tau2 / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + tau2));
559 temp = sqrt(eps);
560 if (d__[*i__] <= temp * d__[ip1] && (r__1 = z__[*i__], f2c_abs(r__1)) <= temp && d__[*i__] > 0.f)
561 {
562 /* Computing MIN */
563 r__1 = d__[*i__] * 10.f;
564 tau = min(r__1,sgub);
565 geomavg = TRUE_;
566 }
567 }
568 else
569 {
570 /* (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2 */
571 /* We choose d(i+1) as origin. */
572 orgati = FALSE_;
573 ii = ip1;
574 sglb = -delsq2 / (d__[ii] + sq2);
575 sgub = 0.f;
576 a = c__ * delsq - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
577 b = z__[ip1] * z__[ip1] * delsq;
578 if (a < 0.f)
579 {
580 tau2 = b * 2.f / (a - sqrt((r__1 = a * a + b * 4.f * c__, f2c_abs( r__1))));
581 }
582 else
583 {
584 tau2 = -(a + sqrt((r__1 = a * a + b * 4.f * c__, f2c_abs(r__1)))) / (c__ * 2.f);
585 }
586 /* TAU2 now is an estimation of SIGMA^2 - D( IP1 )^2. The */
587 /* following, however, is the corresponding estimation of */
588 /* SIGMA - D( IP1 ). */
589 tau = tau2 / (d__[ip1] + sqrt((r__1 = d__[ip1] * d__[ip1] + tau2, f2c_abs(r__1))));
590 }
591 *sigma = d__[ii] + tau;
592 i__1 = *n;
593 for (j = 1;
594 j <= i__1;
595 ++j)
596 {
597 work[j] = d__[j] + d__[ii] + tau;
598 delta[j] = d__[j] - d__[ii] - tau;
599 /* L130: */
600 }
601 iim1 = ii - 1;
602 iip1 = ii + 1;
603 /* Evaluate PSI and the derivative DPSI */
604 dpsi = 0.f;
605 psi = 0.f;
606 erretm = 0.f;
607 i__1 = iim1;
608 for (j = 1;
609 j <= i__1;
610 ++j)
611 {
612 temp = z__[j] / (work[j] * delta[j]);
613 psi += z__[j] * temp;
614 dpsi += temp * temp;
615 erretm += psi;
616 /* L150: */
617 }
618 erretm = f2c_abs(erretm);
619 /* Evaluate PHI and the derivative DPHI */
620 dphi = 0.f;
621 phi = 0.f;
622 i__1 = iip1;
623 for (j = *n;
624 j >= i__1;
625 --j)
626 {
627 temp = z__[j] / (work[j] * delta[j]);
628 phi += z__[j] * temp;
629 dphi += temp * temp;
630 erretm += phi;
631 /* L160: */
632 }
633 w = rhoinv + phi + psi;
634 /* W is the value of the secular function with */
635 /* its ii-th element removed. */
636 swtch3 = FALSE_;
637 if (orgati)
638 {
639 if (w < 0.f)
640 {
641 swtch3 = TRUE_;
642 }
643 }
644 else
645 {
646 if (w > 0.f)
647 {
648 swtch3 = TRUE_;
649 }
650 }
651 if (ii == 1 || ii == *n)
652 {
653 swtch3 = FALSE_;
654 }
655 temp = z__[ii] / (work[ii] * delta[ii]);
656 dw = dpsi + dphi + temp * temp;
657 temp = z__[ii] * temp;
658 w += temp;
659 erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + f2c_abs(temp) * 3.f;
660 /* $ + ABS( TAU2 )*DW */
661 /* Test for convergence */
662 if (f2c_abs(w) <= eps * erretm)
663 {
664 goto L240;
665 }
666 if (w <= 0.f)
667 {
668 sglb = max(sglb,tau);
669 }
670 else
671 {
672 sgub = min(sgub,tau);
673 }
674 /* Calculate the new step */
675 ++niter;
676 if (! swtch3)
677 {
678 dtipsq = work[ip1] * delta[ip1];
679 dtisq = work[*i__] * delta[*i__];
680 if (orgati)
681 {
682 /* Computing 2nd power */
683 r__1 = z__[*i__] / dtisq;
684 c__ = w - dtipsq * dw + delsq * (r__1 * r__1);
685 }
686 else
687 {
688 /* Computing 2nd power */
689 r__1 = z__[ip1] / dtipsq;
690 c__ = w - dtisq * dw - delsq * (r__1 * r__1);
691 }
692 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
693 b = dtipsq * dtisq * w;
694 if (c__ == 0.f)
695 {
696 if (a == 0.f)
697 {
698 if (orgati)
699 {
700 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi + dphi);
701 }
702 else
703 {
704 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi + dphi);
705 }
706 }
707 eta = b / a;
708 }
709 else if (a <= 0.f)
710 {
711 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1)))) / (c__ * 2.f);
712 }
713 else
714 {
715 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs( r__1))));
716 }
717 }
718 else
719 {
720 /* Interpolation using THREE most relevant poles */
721 dtiim = work[iim1] * delta[iim1];
722 dtiip = work[iip1] * delta[iip1];
723 temp = rhoinv + psi + phi;
724 if (orgati)
725 {
726 temp1 = z__[iim1] / dtiim;
727 temp1 *= temp1;
728 c__ = temp - dtiip * (dpsi + dphi) - (d__[iim1] - d__[iip1]) * (d__[iim1] + d__[iip1]) * temp1;
729 zz[0] = z__[iim1] * z__[iim1];
730 if (dpsi < temp1)
731 {
732 zz[2] = dtiip * dtiip * dphi;
733 }
734 else
735 {
736 zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
737 }
738 }
739 else
740 {
741 temp1 = z__[iip1] / dtiip;
742 temp1 *= temp1;
743 c__ = temp - dtiim * (dpsi + dphi) - (d__[iip1] - d__[iim1]) * (d__[iim1] + d__[iip1]) * temp1;
744 if (dphi < temp1)
745 {
746 zz[0] = dtiim * dtiim * dpsi;
747 }
748 else
749 {
750 zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
751 }
752 zz[2] = z__[iip1] * z__[iip1];
753 }
754 zz[1] = z__[ii] * z__[ii];
755 dd[0] = dtiim;
756 dd[1] = delta[ii] * work[ii];
757 dd[2] = dtiip;
758 slaed6_(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
759 if (*info != 0)
760 {
761 /* If INFO is not 0, i.e., SLAED6 failed, switch back */
762 /* to 2 pole interpolation. */
763 swtch3 = FALSE_;
764 *info = 0;
765 dtipsq = work[ip1] * delta[ip1];
766 dtisq = work[*i__] * delta[*i__];
767 if (orgati)
768 {
769 /* Computing 2nd power */
770 r__1 = z__[*i__] / dtisq;
771 c__ = w - dtipsq * dw + delsq * (r__1 * r__1);
772 }
773 else
774 {
775 /* Computing 2nd power */
776 r__1 = z__[ip1] / dtipsq;
777 c__ = w - dtisq * dw - delsq * (r__1 * r__1);
778 }
779 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
780 b = dtipsq * dtisq * w;
781 if (c__ == 0.f)
782 {
783 if (a == 0.f)
784 {
785 if (orgati)
786 {
787 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * ( dpsi + dphi);
788 }
789 else
790 {
791 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi + dphi);
792 }
793 }
794 eta = b / a;
795 }
796 else if (a <= 0.f)
797 {
798 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1))) ) / (c__ * 2.f);
799 }
800 else
801 {
802 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1))));
803 }
804 }
805 }
806 /* Note, eta should be positive if w is negative, and */
807 /* eta should be negative otherwise. However, */
808 /* if for some reason caused by roundoff, eta*w > 0, */
809 /* we simply use one Newton step instead. This way */
810 /* will guarantee eta*w < 0. */
811 if (w * eta >= 0.f)
812 {
813 eta = -w / dw;
814 }
815 eta /= *sigma + sqrt(*sigma * *sigma + eta);
816 temp = tau + eta;
817 if (temp > sgub || temp < sglb)
818 {
819 if (w < 0.f)
820 {
821 eta = (sgub - tau) / 2.f;
822 }
823 else
824 {
825 eta = (sglb - tau) / 2.f;
826 }
827 if (geomavg)
828 {
829 if (w < 0.f)
830 {
831 if (tau > 0.f)
832 {
833 eta = sqrt(sgub * tau) - tau;
834 }
835 }
836 else
837 {
838 if (sglb > 0.f)
839 {
840 eta = sqrt(sglb * tau) - tau;
841 }
842 }
843 }
844 }
845 prew = w;
846 tau += eta;
847 *sigma += eta;
848 i__1 = *n;
849 for (j = 1;
850 j <= i__1;
851 ++j)
852 {
853 work[j] += eta;
854 delta[j] -= eta;
855 /* L170: */
856 }
857 /* Evaluate PSI and the derivative DPSI */
858 dpsi = 0.f;
859 psi = 0.f;
860 erretm = 0.f;
861 i__1 = iim1;
862 for (j = 1;
863 j <= i__1;
864 ++j)
865 {
866 temp = z__[j] / (work[j] * delta[j]);
867 psi += z__[j] * temp;
868 dpsi += temp * temp;
869 erretm += psi;
870 /* L180: */
871 }
872 erretm = f2c_abs(erretm);
873 /* Evaluate PHI and the derivative DPHI */
874 dphi = 0.f;
875 phi = 0.f;
876 i__1 = iip1;
877 for (j = *n;
878 j >= i__1;
879 --j)
880 {
881 temp = z__[j] / (work[j] * delta[j]);
882 phi += z__[j] * temp;
883 dphi += temp * temp;
884 erretm += phi;
885 /* L190: */
886 }
887 tau2 = work[ii] * delta[ii];
888 temp = z__[ii] / tau2;
889 dw = dpsi + dphi + temp * temp;
890 temp = z__[ii] * temp;
891 w = rhoinv + phi + psi + temp;
892 erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + f2c_abs(temp) * 3.f;
893 /* $ + ABS( TAU2 )*DW */
894 swtch = FALSE_;
895 if (orgati)
896 {
897 if (-w > f2c_abs(prew) / 10.f)
898 {
899 swtch = TRUE_;
900 }
901 }
902 else
903 {
904 if (w > f2c_abs(prew) / 10.f)
905 {
906 swtch = TRUE_;
907 }
908 }
909 /* Main loop to update the values of the array DELTA and WORK */
910 iter = niter + 1;
911 for (niter = iter;
912 niter <= 400;
913 ++niter)
914 {
915 /* Test for convergence */
916 if (f2c_abs(w) <= eps * erretm)
917 {
918 /* $ .OR. (SGUB-SGLB).LE.EIGHT*ABS(SGUB+SGLB) ) THEN */
919 goto L240;
920 }
921 if (w <= 0.f)
922 {
923 sglb = max(sglb,tau);
924 }
925 else
926 {
927 sgub = min(sgub,tau);
928 }
929 /* Calculate the new step */
930 if (! swtch3)
931 {
932 dtipsq = work[ip1] * delta[ip1];
933 dtisq = work[*i__] * delta[*i__];
934 if (! swtch)
935 {
936 if (orgati)
937 {
938 /* Computing 2nd power */
939 r__1 = z__[*i__] / dtisq;
940 c__ = w - dtipsq * dw + delsq * (r__1 * r__1);
941 }
942 else
943 {
944 /* Computing 2nd power */
945 r__1 = z__[ip1] / dtipsq;
946 c__ = w - dtisq * dw - delsq * (r__1 * r__1);
947 }
948 }
949 else
950 {
951 temp = z__[ii] / (work[ii] * delta[ii]);
952 if (orgati)
953 {
954 dpsi += temp * temp;
955 }
956 else
957 {
958 dphi += temp * temp;
959 }
960 c__ = w - dtisq * dpsi - dtipsq * dphi;
961 }
962 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
963 b = dtipsq * dtisq * w;
964 if (c__ == 0.f)
965 {
966 if (a == 0.f)
967 {
968 if (! swtch)
969 {
970 if (orgati)
971 {
972 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi + dphi);
973 }
974 else
975 {
976 a = z__[ip1] * z__[ip1] + dtisq * dtisq * ( dpsi + dphi);
977 }
978 }
979 else
980 {
981 a = dtisq * dtisq * dpsi + dtipsq * dtipsq * dphi;
982 }
983 }
984 eta = b / a;
985 }
986 else if (a <= 0.f)
987 {
988 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1))) ) / (c__ * 2.f);
989 }
990 else
991 {
992 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1))));
993 }
994 }
995 else
996 {
997 /* Interpolation using THREE most relevant poles */
998 dtiim = work[iim1] * delta[iim1];
999 dtiip = work[iip1] * delta[iip1];
1000 temp = rhoinv + psi + phi;
1001 if (swtch)
1002 {
1003 c__ = temp - dtiim * dpsi - dtiip * dphi;
1004 zz[0] = dtiim * dtiim * dpsi;
1005 zz[2] = dtiip * dtiip * dphi;
1006 }
1007 else
1008 {
1009 if (orgati)
1010 {
1011 temp1 = z__[iim1] / dtiim;
1012 temp1 *= temp1;
1013 temp2 = (d__[iim1] - d__[iip1]) * (d__[iim1] + d__[ iip1]) * temp1;
1014 c__ = temp - dtiip * (dpsi + dphi) - temp2;
1015 zz[0] = z__[iim1] * z__[iim1];
1016 if (dpsi < temp1)
1017 {
1018 zz[2] = dtiip * dtiip * dphi;
1019 }
1020 else
1021 {
1022 zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
1023 }
1024 }
1025 else
1026 {
1027 temp1 = z__[iip1] / dtiip;
1028 temp1 *= temp1;
1029 temp2 = (d__[iip1] - d__[iim1]) * (d__[iim1] + d__[ iip1]) * temp1;
1030 c__ = temp - dtiim * (dpsi + dphi) - temp2;
1031 if (dphi < temp1)
1032 {
1033 zz[0] = dtiim * dtiim * dpsi;
1034 }
1035 else
1036 {
1037 zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
1038 }
1039 zz[2] = z__[iip1] * z__[iip1];
1040 }
1041 }
1042 dd[0] = dtiim;
1043 dd[1] = delta[ii] * work[ii];
1044 dd[2] = dtiip;
1045 slaed6_(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
1046 if (*info != 0)
1047 {
1048 /* If INFO is not 0, i.e., SLAED6 failed, switch */
1049 /* back to two pole interpolation */
1050 swtch3 = FALSE_;
1051 *info = 0;
1052 dtipsq = work[ip1] * delta[ip1];
1053 dtisq = work[*i__] * delta[*i__];
1054 if (! swtch)
1055 {
1056 if (orgati)
1057 {
1058 /* Computing 2nd power */
1059 r__1 = z__[*i__] / dtisq;
1060 c__ = w - dtipsq * dw + delsq * (r__1 * r__1);
1061 }
1062 else
1063 {
1064 /* Computing 2nd power */
1065 r__1 = z__[ip1] / dtipsq;
1066 c__ = w - dtisq * dw - delsq * (r__1 * r__1);
1067 }
1068 }
1069 else
1070 {
1071 temp = z__[ii] / (work[ii] * delta[ii]);
1072 if (orgati)
1073 {
1074 dpsi += temp * temp;
1075 }
1076 else
1077 {
1078 dphi += temp * temp;
1079 }
1080 c__ = w - dtisq * dpsi - dtipsq * dphi;
1081 }
1082 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
1083 b = dtipsq * dtisq * w;
1084 if (c__ == 0.f)
1085 {
1086 if (a == 0.f)
1087 {
1088 if (! swtch)
1089 {
1090 if (orgati)
1091 {
1092 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi + dphi);
1093 }
1094 else
1095 {
1096 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi + dphi);
1097 }
1098 }
1099 else
1100 {
1101 a = dtisq * dtisq * dpsi + dtipsq * dtipsq * dphi;
1102 }
1103 }
1104 eta = b / a;
1105 }
1106 else if (a <= 0.f)
1107 {
1108 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs( r__1)))) / (c__ * 2.f);
1109 }
1110 else
1111 {
1112 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, f2c_abs(r__1))));
1113 }
1114 }
1115 }
1116 /* Note, eta should be positive if w is negative, and */
1117 /* eta should be negative otherwise. However, */
1118 /* if for some reason caused by roundoff, eta*w > 0, */
1119 /* we simply use one Newton step instead. This way */
1120 /* will guarantee eta*w < 0. */
1121 if (w * eta >= 0.f)
1122 {
1123 eta = -w / dw;
1124 }
1125 eta /= *sigma + sqrt(*sigma * *sigma + eta);
1126 temp = tau + eta;
1127 if (temp > sgub || temp < sglb)
1128 {
1129 if (w < 0.f)
1130 {
1131 eta = (sgub - tau) / 2.f;
1132 }
1133 else
1134 {
1135 eta = (sglb - tau) / 2.f;
1136 }
1137 if (geomavg)
1138 {
1139 if (w < 0.f)
1140 {
1141 if (tau > 0.f)
1142 {
1143 eta = sqrt(sgub * tau) - tau;
1144 }
1145 }
1146 else
1147 {
1148 if (sglb > 0.f)
1149 {
1150 eta = sqrt(sglb * tau) - tau;
1151 }
1152 }
1153 }
1154 }
1155 prew = w;
1156 tau += eta;
1157 *sigma += eta;
1158 i__1 = *n;
1159 for (j = 1;
1160 j <= i__1;
1161 ++j)
1162 {
1163 work[j] += eta;
1164 delta[j] -= eta;
1165 /* L200: */
1166 }
1167 /* Evaluate PSI and the derivative DPSI */
1168 dpsi = 0.f;
1169 psi = 0.f;
1170 erretm = 0.f;
1171 i__1 = iim1;
1172 for (j = 1;
1173 j <= i__1;
1174 ++j)
1175 {
1176 temp = z__[j] / (work[j] * delta[j]);
1177 psi += z__[j] * temp;
1178 dpsi += temp * temp;
1179 erretm += psi;
1180 /* L210: */
1181 }
1182 erretm = f2c_abs(erretm);
1183 /* Evaluate PHI and the derivative DPHI */
1184 dphi = 0.f;
1185 phi = 0.f;
1186 i__1 = iip1;
1187 for (j = *n;
1188 j >= i__1;
1189 --j)
1190 {
1191 temp = z__[j] / (work[j] * delta[j]);
1192 phi += z__[j] * temp;
1193 dphi += temp * temp;
1194 erretm += phi;
1195 /* L220: */
1196 }
1197 tau2 = work[ii] * delta[ii];
1198 temp = z__[ii] / tau2;
1199 dw = dpsi + dphi + temp * temp;
1200 temp = z__[ii] * temp;
1201 w = rhoinv + phi + psi + temp;
1202 erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + f2c_abs(temp) * 3.f;
1203 /* $ + ABS( TAU2 )*DW */
1204 if (w * prew > 0.f && f2c_abs(w) > f2c_abs(prew) / 10.f)
1205 {
1206 swtch = ! swtch;
1207 }
1208 /* L230: */
1209 }
1210 /* Return with INFO = 1, NITER = MAXIT and not converged */
1211 *info = 1;
1212 }
1213 L240:
1214 return 0;
1215 /* End of SLASD4 */
1216 }
1217 /* slasd4_ */
1218