2 #include "../gmx_lapack.h"
3 #include "lapack_limits.h"
5 #include "types/simple.h"
8 F77_FUNC(dlasd4,DLASD4)(int *n,
27 double eta, phi, eps, tau, psi;
31 double temp, prew, sg2lb, sg2ub, temp1, temp2, dtiim, delsq,
37 double delsq2, dtnsq1;
40 double erretm, dtipsq, rhoinv;
50 *sigma = sqrt(d__[1] * d__[1] + *rho * z__[1] * z__[1]);
56 F77_FUNC(dlasd5,DLASD5)(i__, &d__[1], &z__[1], &delta[1], rho, sigma, &work[1]);
70 temp1 = temp / (d__[*n] + sqrt(d__[*n] * d__[*n] + temp));
72 for (j = 1; j <= i__1; ++j) {
73 work[j] = d__[j] + d__[*n] + temp1;
74 delta[j] = d__[j] - d__[*n] - temp1;
79 for (j = 1; j <= i__1; ++j) {
80 psi += z__[j] * z__[j] / (delta[j] * work[j]);
84 w = c__ + z__[ii] * z__[ii] / (delta[ii] * work[ii]) + z__[*n] * z__[*
85 n] / (delta[*n] * work[*n]);
88 temp1 = sqrt(d__[*n] * d__[*n] + *rho);
89 temp = z__[*n - 1] * z__[*n - 1] / ((d__[*n - 1] + temp1) * (d__[*
90 n] - d__[*n - 1] + *rho / (d__[*n] + temp1))) + z__[*n] *
96 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
97 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*
99 b = z__[*n] * z__[*n] * delsq;
101 tau = b * 2. / (sqrt(a * a + b * 4. * c__) - a);
103 tau = (a + sqrt(a * a + b * 4. * c__)) / (c__ * 2.);
108 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
109 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
110 b = z__[*n] * z__[*n] * delsq;
113 tau = b * 2. / (sqrt(a * a + b * 4. * c__) - a);
115 tau = (a + sqrt(a * a + b * 4. * c__)) / (c__ * 2.);
120 eta = tau / (d__[*n] + sqrt(d__[*n] * d__[*n] + tau));
122 *sigma = d__[*n] + eta;
124 for (j = 1; j <= i__1; ++j) {
125 delta[j] = d__[j] - d__[*i__] - eta;
126 work[j] = d__[j] + d__[*i__] + eta;
133 for (j = 1; j <= i__1; ++j) {
134 temp = z__[j] / (delta[j] * work[j]);
135 psi += z__[j] * temp;
139 erretm = fabs(erretm);
141 temp = z__[*n] / (delta[*n] * work[*n]);
142 phi = z__[*n] * temp;
144 erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (dpsi
147 w = rhoinv + phi + psi;
149 if (fabs(w) <= eps * erretm) {
154 dtnsq1 = work[*n - 1] * delta[*n - 1];
155 dtnsq = work[*n] * delta[*n];
156 c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
157 a = (dtnsq + dtnsq1) * w - dtnsq * dtnsq1 * (dpsi + dphi);
158 b = dtnsq * dtnsq1 * w;
162 if ( fabs(c__)<GMX_DOUBLE_MIN) {
163 eta = *rho - *sigma * *sigma;
164 } else if (a >= 0.) {
165 eta = (a + sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
167 eta = b * 2. / (a - sqrt(fabs(a * a - b * 4. * c__)));
171 eta = -w / (dpsi + dphi);
179 eta /= *sigma + sqrt(eta + *sigma * *sigma);
181 for (j = 1; j <= i__1; ++j) {
192 for (j = 1; j <= i__1; ++j) {
193 temp = z__[j] / (work[j] * delta[j]);
194 psi += z__[j] * temp;
198 erretm = fabs(erretm);
200 temp = z__[*n] / (work[*n] * delta[*n]);
201 phi = z__[*n] * temp;
203 erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (dpsi
206 w = rhoinv + phi + psi;
210 for (niter = iter; niter <= 20; ++niter) {
212 if (fabs(w) <= eps * erretm) {
215 dtnsq1 = work[*n - 1] * delta[*n - 1];
216 dtnsq = work[*n] * delta[*n];
217 c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
218 a = (dtnsq + dtnsq1) * w - dtnsq1 * dtnsq * (dpsi + dphi);
219 b = dtnsq1 * dtnsq * w;
221 eta = (a + sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
223 eta = b * 2. / (a - sqrt(fabs(a * a - b * 4. * c__)));
227 eta = -w / (dpsi + dphi);
235 eta /= *sigma + sqrt(eta + *sigma * *sigma);
237 for (j = 1; j <= i__1; ++j) {
248 for (j = 1; j <= i__1; ++j) {
249 temp = z__[j] / (work[j] * delta[j]);
250 psi += z__[j] * temp;
254 erretm = fabs(erretm);
256 temp = z__[*n] / (work[*n] * delta[*n]);
257 phi = z__[*n] * temp;
259 erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (
262 w = rhoinv + phi + psi;
273 delsq = (d__[ip1] - d__[*i__]) * (d__[ip1] + d__[*i__]);
275 temp = delsq2 / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + delsq2));
277 for (j = 1; j <= i__1; ++j) {
278 work[j] = d__[j] + d__[*i__] + temp;
279 delta[j] = d__[j] - d__[*i__] - temp;
284 for (j = 1; j <= i__1; ++j) {
285 psi += z__[j] * z__[j] / (work[j] * delta[j]);
290 for (j = *n; j >= i__1; --j) {
291 phi += z__[j] * z__[j] / (work[j] * delta[j]);
293 c__ = rhoinv + psi + phi;
294 w = c__ + z__[*i__] * z__[*i__] / (work[*i__] * delta[*i__]) + z__[
295 ip1] * z__[ip1] / (work[ip1] * delta[ip1]);
302 a = c__ * delsq + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
303 b = z__[*i__] * z__[*i__] * delsq;
305 tau = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
307 tau = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
309 eta = tau / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + tau));
315 a = c__ * delsq - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
316 b = z__[ip1] * z__[ip1] * delsq;
318 tau = b * 2. / (a - sqrt(fabs(a * a + b * 4. * c__)));
320 tau = -(a + sqrt(fabs(a * a + b * 4. * c__))) / (c__ * 2.);
322 eta = tau / (d__[ip1] + sqrt(fabs(d__[ip1] * d__[ip1] + tau)));
327 *sigma = d__[*i__] + eta;
329 for (j = 1; j <= i__1; ++j) {
330 work[j] = d__[j] + d__[*i__] + eta;
331 delta[j] = d__[j] - d__[*i__] - eta;
335 *sigma = d__[ip1] + eta;
337 for (j = 1; j <= i__1; ++j) {
338 work[j] = d__[j] + d__[ip1] + eta;
339 delta[j] = d__[j] - d__[ip1] - eta;
349 for (j = 1; j <= i__1; ++j) {
350 temp = z__[j] / (work[j] * delta[j]);
351 psi += z__[j] * temp;
355 erretm = fabs(erretm);
360 for (j = *n; j >= i__1; --j) {
361 temp = z__[j] / (work[j] * delta[j]);
362 phi += z__[j] * temp;
367 w = rhoinv + phi + psi;
379 if (ii == 1 || ii == *n) {
383 temp = z__[ii] / (work[ii] * delta[ii]);
384 dw = dpsi + dphi + temp * temp;
385 temp = z__[ii] * temp;
387 erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3. +
390 if (fabs(w) <= eps * erretm) {
395 sg2lb = (sg2lb > tau) ? sg2lb : tau;
397 sg2ub = (sg2ub < tau) ? sg2ub : tau;
402 dtipsq = work[ip1] * delta[ip1];
403 dtisq = work[*i__] * delta[*i__];
405 d__1 = z__[*i__] / dtisq;
406 c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
408 d__1 = z__[ip1] / dtipsq;
409 c__ = w - dtisq * dw - delsq * (d__1 * d__1);
411 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
412 b = dtipsq * dtisq * w;
413 if ( fabs(c__)<GMX_DOUBLE_MIN) {
414 if ( fabs(a)<GMX_DOUBLE_MIN) {
416 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi +
419 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi +
424 } else if (a <= 0.) {
425 eta = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
427 eta = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
431 dtiim = work[iim1] * delta[iim1];
432 dtiip = work[iip1] * delta[iip1];
433 temp = rhoinv + psi + phi;
435 temp1 = z__[iim1] / dtiim;
437 c__ = temp - dtiip * (dpsi + dphi) - (d__[iim1] - d__[iip1]) *
438 (d__[iim1] + d__[iip1]) * temp1;
439 zz[0] = z__[iim1] * z__[iim1];
441 zz[2] = dtiip * dtiip * dphi;
443 zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
446 temp1 = z__[iip1] / dtiip;
448 c__ = temp - dtiim * (dpsi + dphi) - (d__[iip1] - d__[iim1]) *
449 (d__[iim1] + d__[iip1]) * temp1;
451 zz[0] = dtiim * dtiim * dpsi;
453 zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
455 zz[2] = z__[iip1] * z__[iip1];
457 zz[1] = z__[ii] * z__[ii];
459 dd[1] = delta[ii] * work[ii];
461 F77_FUNC(dlaed6,DLAED6)(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
471 temp1 = work[*i__] * delta[*i__];
474 temp1 = work[ip1] * delta[ip1];
477 if (temp > sg2ub || temp < sg2lb) {
479 eta = (sg2ub - tau) / 2.;
481 eta = (sg2lb - tau) / 2.;
486 eta /= *sigma + sqrt(*sigma * *sigma + eta);
492 for (j = 1; j <= i__1; ++j) {
501 for (j = 1; j <= i__1; ++j) {
502 temp = z__[j] / (work[j] * delta[j]);
503 psi += z__[j] * temp;
507 erretm = fabs(erretm);
512 for (j = *n; j >= i__1; --j) {
513 temp = z__[j] / (work[j] * delta[j]);
514 phi += z__[j] * temp;
519 temp = z__[ii] / (work[ii] * delta[ii]);
520 dw = dpsi + dphi + temp * temp;
521 temp = z__[ii] * temp;
522 w = rhoinv + phi + psi + temp;
523 erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3. +
527 sg2lb = (sg2lb > tau) ? sg2lb : tau;
529 sg2ub = (sg2ub < tau) ? sg2ub : tau;
534 if (-w > fabs(prew) / 10.) {
538 if (w > fabs(prew) / 10.) {
545 for (niter = iter; niter <= 20; ++niter) {
547 if (fabs(w) <= eps * erretm) {
552 dtipsq = work[ip1] * delta[ip1];
553 dtisq = work[*i__] * delta[*i__];
556 d__1 = z__[*i__] / dtisq;
557 c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
559 d__1 = z__[ip1] / dtipsq;
560 c__ = w - dtisq * dw - delsq * (d__1 * d__1);
563 temp = z__[ii] / (work[ii] * delta[ii]);
569 c__ = w - dtisq * dpsi - dtipsq * dphi;
571 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
572 b = dtipsq * dtisq * w;
573 if (fabs(c__)<GMX_DOUBLE_MIN) {
574 if (fabs(a)<GMX_DOUBLE_MIN) {
577 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq *
580 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (
584 a = dtisq * dtisq * dpsi + dtipsq * dtipsq * dphi;
588 } else if (a <= 0.) {
589 eta = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
591 eta = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
595 dtiim = work[iim1] * delta[iim1];
596 dtiip = work[iip1] * delta[iip1];
597 temp = rhoinv + psi + phi;
599 c__ = temp - dtiim * dpsi - dtiip * dphi;
600 zz[0] = dtiim * dtiim * dpsi;
601 zz[2] = dtiip * dtiip * dphi;
604 temp1 = z__[iim1] / dtiim;
606 temp2 = (d__[iim1] - d__[iip1]) * (d__[iim1] + d__[
608 c__ = temp - dtiip * (dpsi + dphi) - temp2;
609 zz[0] = z__[iim1] * z__[iim1];
611 zz[2] = dtiip * dtiip * dphi;
613 zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
616 temp1 = z__[iip1] / dtiip;
618 temp2 = (d__[iip1] - d__[iim1]) * (d__[iim1] + d__[
620 c__ = temp - dtiim * (dpsi + dphi) - temp2;
622 zz[0] = dtiim * dtiim * dpsi;
624 zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
626 zz[2] = z__[iip1] * z__[iip1];
630 dd[1] = delta[ii] * work[ii];
632 F77_FUNC(dlaed6,DLAED6)(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
642 temp1 = work[*i__] * delta[*i__];
645 temp1 = work[ip1] * delta[ip1];
648 if (temp > sg2ub || temp < sg2lb) {
650 eta = (sg2ub - tau) / 2.;
652 eta = (sg2lb - tau) / 2.;
657 eta /= *sigma + sqrt(*sigma * *sigma + eta);
661 for (j = 1; j <= i__1; ++j) {
672 for (j = 1; j <= i__1; ++j) {
673 temp = z__[j] / (work[j] * delta[j]);
674 psi += z__[j] * temp;
678 erretm = fabs(erretm);
683 for (j = *n; j >= i__1; --j) {
684 temp = z__[j] / (work[j] * delta[j]);
685 phi += z__[j] * temp;
690 temp = z__[ii] / (work[ii] * delta[ii]);
691 dw = dpsi + dphi + temp * temp;
692 temp = z__[ii] * temp;
693 w = rhoinv + phi + psi + temp;
694 erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3.
696 if (w * prew > 0. && fabs(w) > fabs(prew) / 10.) {
701 sg2lb = (sg2lb > tau) ? sg2lb : tau;
703 sg2ub = (sg2ub < tau) ? sg2ub : tau;