slasd4.c
Go to the documentation of this file.
00001 /* slasd4.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int slasd4_(integer *n, integer *i__, real *d__, real *z__, 
00017         real *delta, real *rho, real *sigma, real *work, integer *info)
00018 {
00019     /* System generated locals */
00020     integer i__1;
00021     real r__1;
00022 
00023     /* Builtin functions */
00024     double sqrt(doublereal);
00025 
00026     /* Local variables */
00027     real a, b, c__;
00028     integer j;
00029     real w, dd[3];
00030     integer ii;
00031     real dw, zz[3];
00032     integer ip1;
00033     real eta, phi, eps, tau, psi;
00034     integer iim1, iip1;
00035     real dphi, dpsi;
00036     integer iter;
00037     real temp, prew, sg2lb, sg2ub, temp1, temp2, dtiim, delsq, dtiip;
00038     integer niter;
00039     real dtisq;
00040     logical swtch;
00041     real dtnsq;
00042     extern /* Subroutine */ int slaed6_(integer *, logical *, real *, real *, 
00043             real *, real *, real *, integer *);
00044     real delsq2;
00045     extern /* Subroutine */ int slasd5_(integer *, real *, real *, real *, 
00046             real *, real *, real *);
00047     real dtnsq1;
00048     logical swtch3;
00049     extern doublereal slamch_(char *);
00050     logical orgati;
00051     real erretm, dtipsq, rhoinv;
00052 
00053 
00054 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00055 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00056 /*     November 2006 */
00057 
00058 /*     .. Scalar Arguments .. */
00059 /*     .. */
00060 /*     .. Array Arguments .. */
00061 /*     .. */
00062 
00063 /*  Purpose */
00064 /*  ======= */
00065 
00066 /*  This subroutine computes the square root of the I-th updated */
00067 /*  eigenvalue of a positive symmetric rank-one modification to */
00068 /*  a positive diagonal matrix whose entries are given as the squares */
00069 /*  of the corresponding entries in the array d, and that */
00070 
00071 /*         0 <= D(i) < D(j)  for  i < j */
00072 
00073 /*  and that RHO > 0. This is arranged by the calling routine, and is */
00074 /*  no loss in generality.  The rank-one modified system is thus */
00075 
00076 /*         diag( D ) * diag( D ) +  RHO *  Z * Z_transpose. */
00077 
00078 /*  where we assume the Euclidean norm of Z is 1. */
00079 
00080 /*  The method consists of approximating the rational functions in the */
00081 /*  secular equation by simpler interpolating rational functions. */
00082 
00083 /*  Arguments */
00084 /*  ========= */
00085 
00086 /*  N      (input) INTEGER */
00087 /*         The length of all arrays. */
00088 
00089 /*  I      (input) INTEGER */
00090 /*         The index of the eigenvalue to be computed.  1 <= I <= N. */
00091 
00092 /*  D      (input) REAL array, dimension ( N ) */
00093 /*         The original eigenvalues.  It is assumed that they are in */
00094 /*         order, 0 <= D(I) < D(J)  for I < J. */
00095 
00096 /*  Z      (input) REAL array, dimension (N) */
00097 /*         The components of the updating vector. */
00098 
00099 /*  DELTA  (output) REAL array, dimension (N) */
00100 /*         If N .ne. 1, DELTA contains (D(j) - sigma_I) in its  j-th */
00101 /*         component.  If N = 1, then DELTA(1) = 1.  The vector DELTA */
00102 /*         contains the information necessary to construct the */
00103 /*         (singular) eigenvectors. */
00104 
00105 /*  RHO    (input) REAL */
00106 /*         The scalar in the symmetric updating formula. */
00107 
00108 /*  SIGMA  (output) REAL */
00109 /*         The computed sigma_I, the I-th updated eigenvalue. */
00110 
00111 /*  WORK   (workspace) REAL array, dimension (N) */
00112 /*         If N .ne. 1, WORK contains (D(j) + sigma_I) in its  j-th */
00113 /*         component.  If N = 1, then WORK( 1 ) = 1. */
00114 
00115 /*  INFO   (output) INTEGER */
00116 /*         = 0:  successful exit */
00117 /*         > 0:  if INFO = 1, the updating process failed. */
00118 
00119 /*  Internal Parameters */
00120 /*  =================== */
00121 
00122 /*  Logical variable ORGATI (origin-at-i?) is used for distinguishing */
00123 /*  whether D(i) or D(i+1) is treated as the origin. */
00124 
00125 /*            ORGATI = .true.    origin at i */
00126 /*            ORGATI = .false.   origin at i+1 */
00127 
00128 /*  Logical variable SWTCH3 (switch-for-3-poles?) is for noting */
00129 /*  if we are working with THREE poles! */
00130 
00131 /*  MAXIT is the maximum number of iterations allowed for each */
00132 /*  eigenvalue. */
00133 
00134 /*  Further Details */
00135 /*  =============== */
00136 
00137 /*  Based on contributions by */
00138 /*     Ren-Cang Li, Computer Science Division, University of California */
00139 /*     at Berkeley, USA */
00140 
00141 /*  ===================================================================== */
00142 
00143 /*     .. Parameters .. */
00144 /*     .. */
00145 /*     .. Local Scalars .. */
00146 /*     .. */
00147 /*     .. Local Arrays .. */
00148 /*     .. */
00149 /*     .. External Subroutines .. */
00150 /*     .. */
00151 /*     .. External Functions .. */
00152 /*     .. */
00153 /*     .. Intrinsic Functions .. */
00154 /*     .. */
00155 /*     .. Executable Statements .. */
00156 
00157 /*     Since this routine is called in an inner loop, we do no argument */
00158 /*     checking. */
00159 
00160 /*     Quick return for N=1 and 2. */
00161 
00162     /* Parameter adjustments */
00163     --work;
00164     --delta;
00165     --z__;
00166     --d__;
00167 
00168     /* Function Body */
00169     *info = 0;
00170     if (*n == 1) {
00171 
00172 /*        Presumably, I=1 upon entry */
00173 
00174         *sigma = sqrt(d__[1] * d__[1] + *rho * z__[1] * z__[1]);
00175         delta[1] = 1.f;
00176         work[1] = 1.f;
00177         return 0;
00178     }
00179     if (*n == 2) {
00180         slasd5_(i__, &d__[1], &z__[1], &delta[1], rho, sigma, &work[1]);
00181         return 0;
00182     }
00183 
00184 /*     Compute machine epsilon */
00185 
00186     eps = slamch_("Epsilon");
00187     rhoinv = 1.f / *rho;
00188 
00189 /*     The case I = N */
00190 
00191     if (*i__ == *n) {
00192 
00193 /*        Initialize some basic variables */
00194 
00195         ii = *n - 1;
00196         niter = 1;
00197 
00198 /*        Calculate initial guess */
00199 
00200         temp = *rho / 2.f;
00201 
00202 /*        If ||Z||_2 is not one, then TEMP should be set to */
00203 /*        RHO * ||Z||_2^2 / TWO */
00204 
00205         temp1 = temp / (d__[*n] + sqrt(d__[*n] * d__[*n] + temp));
00206         i__1 = *n;
00207         for (j = 1; j <= i__1; ++j) {
00208             work[j] = d__[j] + d__[*n] + temp1;
00209             delta[j] = d__[j] - d__[*n] - temp1;
00210 /* L10: */
00211         }
00212 
00213         psi = 0.f;
00214         i__1 = *n - 2;
00215         for (j = 1; j <= i__1; ++j) {
00216             psi += z__[j] * z__[j] / (delta[j] * work[j]);
00217 /* L20: */
00218         }
00219 
00220         c__ = rhoinv + psi;
00221         w = c__ + z__[ii] * z__[ii] / (delta[ii] * work[ii]) + z__[*n] * z__[*
00222                 n] / (delta[*n] * work[*n]);
00223 
00224         if (w <= 0.f) {
00225             temp1 = sqrt(d__[*n] * d__[*n] + *rho);
00226             temp = z__[*n - 1] * z__[*n - 1] / ((d__[*n - 1] + temp1) * (d__[*
00227                     n] - d__[*n - 1] + *rho / (d__[*n] + temp1))) + z__[*n] * 
00228                     z__[*n] / *rho;
00229 
00230 /*           The following TAU is to approximate */
00231 /*           SIGMA_n^2 - D( N )*D( N ) */
00232 
00233             if (c__ <= temp) {
00234                 tau = *rho;
00235             } else {
00236                 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
00237                 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*
00238                         n];
00239                 b = z__[*n] * z__[*n] * delsq;
00240                 if (a < 0.f) {
00241                     tau = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
00242                 } else {
00243                     tau = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
00244                 }
00245             }
00246 
00247 /*           It can be proved that */
00248 /*               D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU <= D(N)^2+RHO */
00249 
00250         } else {
00251             delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
00252             a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
00253             b = z__[*n] * z__[*n] * delsq;
00254 
00255 /*           The following TAU is to approximate */
00256 /*           SIGMA_n^2 - D( N )*D( N ) */
00257 
00258             if (a < 0.f) {
00259                 tau = b * 2.f / (sqrt(a * a + b * 4.f * c__) - a);
00260             } else {
00261                 tau = (a + sqrt(a * a + b * 4.f * c__)) / (c__ * 2.f);
00262             }
00263 
00264 /*           It can be proved that */
00265 /*           D(N)^2 < D(N)^2+TAU < SIGMA(N)^2 < D(N)^2+RHO/2 */
00266 
00267         }
00268 
00269 /*        The following ETA is to approximate SIGMA_n - D( N ) */
00270 
00271         eta = tau / (d__[*n] + sqrt(d__[*n] * d__[*n] + tau));
00272 
00273         *sigma = d__[*n] + eta;
00274         i__1 = *n;
00275         for (j = 1; j <= i__1; ++j) {
00276             delta[j] = d__[j] - d__[*i__] - eta;
00277             work[j] = d__[j] + d__[*i__] + eta;
00278 /* L30: */
00279         }
00280 
00281 /*        Evaluate PSI and the derivative DPSI */
00282 
00283         dpsi = 0.f;
00284         psi = 0.f;
00285         erretm = 0.f;
00286         i__1 = ii;
00287         for (j = 1; j <= i__1; ++j) {
00288             temp = z__[j] / (delta[j] * work[j]);
00289             psi += z__[j] * temp;
00290             dpsi += temp * temp;
00291             erretm += psi;
00292 /* L40: */
00293         }
00294         erretm = dabs(erretm);
00295 
00296 /*        Evaluate PHI and the derivative DPHI */
00297 
00298         temp = z__[*n] / (delta[*n] * work[*n]);
00299         phi = z__[*n] * temp;
00300         dphi = temp * temp;
00301         erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
00302                 dpsi + dphi);
00303 
00304         w = rhoinv + phi + psi;
00305 
00306 /*        Test for convergence */
00307 
00308         if (dabs(w) <= eps * erretm) {
00309             goto L240;
00310         }
00311 
00312 /*        Calculate the new step */
00313 
00314         ++niter;
00315         dtnsq1 = work[*n - 1] * delta[*n - 1];
00316         dtnsq = work[*n] * delta[*n];
00317         c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
00318         a = (dtnsq + dtnsq1) * w - dtnsq * dtnsq1 * (dpsi + dphi);
00319         b = dtnsq * dtnsq1 * w;
00320         if (c__ < 0.f) {
00321             c__ = dabs(c__);
00322         }
00323         if (c__ == 0.f) {
00324             eta = *rho - *sigma * *sigma;
00325         } else if (a >= 0.f) {
00326             eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) / (
00327                     c__ * 2.f);
00328         } else {
00329             eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
00330                     r__1))));
00331         }
00332 
00333 /*        Note, eta should be positive if w is negative, and */
00334 /*        eta should be negative otherwise. However, */
00335 /*        if for some reason caused by roundoff, eta*w > 0, */
00336 /*        we simply use one Newton step instead. This way */
00337 /*        will guarantee eta*w < 0. */
00338 
00339         if (w * eta > 0.f) {
00340             eta = -w / (dpsi + dphi);
00341         }
00342         temp = eta - dtnsq;
00343         if (temp > *rho) {
00344             eta = *rho + dtnsq;
00345         }
00346 
00347         tau += eta;
00348         eta /= *sigma + sqrt(eta + *sigma * *sigma);
00349         i__1 = *n;
00350         for (j = 1; j <= i__1; ++j) {
00351             delta[j] -= eta;
00352             work[j] += eta;
00353 /* L50: */
00354         }
00355 
00356         *sigma += eta;
00357 
00358 /*        Evaluate PSI and the derivative DPSI */
00359 
00360         dpsi = 0.f;
00361         psi = 0.f;
00362         erretm = 0.f;
00363         i__1 = ii;
00364         for (j = 1; j <= i__1; ++j) {
00365             temp = z__[j] / (work[j] * delta[j]);
00366             psi += z__[j] * temp;
00367             dpsi += temp * temp;
00368             erretm += psi;
00369 /* L60: */
00370         }
00371         erretm = dabs(erretm);
00372 
00373 /*        Evaluate PHI and the derivative DPHI */
00374 
00375         temp = z__[*n] / (work[*n] * delta[*n]);
00376         phi = z__[*n] * temp;
00377         dphi = temp * temp;
00378         erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * (
00379                 dpsi + dphi);
00380 
00381         w = rhoinv + phi + psi;
00382 
00383 /*        Main loop to update the values of the array   DELTA */
00384 
00385         iter = niter + 1;
00386 
00387         for (niter = iter; niter <= 20; ++niter) {
00388 
00389 /*           Test for convergence */
00390 
00391             if (dabs(w) <= eps * erretm) {
00392                 goto L240;
00393             }
00394 
00395 /*           Calculate the new step */
00396 
00397             dtnsq1 = work[*n - 1] * delta[*n - 1];
00398             dtnsq = work[*n] * delta[*n];
00399             c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
00400             a = (dtnsq + dtnsq1) * w - dtnsq1 * dtnsq * (dpsi + dphi);
00401             b = dtnsq1 * dtnsq * w;
00402             if (a >= 0.f) {
00403                 eta = (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
00404                          (c__ * 2.f);
00405             } else {
00406                 eta = b * 2.f / (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(
00407                         r__1))));
00408             }
00409 
00410 /*           Note, eta should be positive if w is negative, and */
00411 /*           eta should be negative otherwise. However, */
00412 /*           if for some reason caused by roundoff, eta*w > 0, */
00413 /*           we simply use one Newton step instead. This way */
00414 /*           will guarantee eta*w < 0. */
00415 
00416             if (w * eta > 0.f) {
00417                 eta = -w / (dpsi + dphi);
00418             }
00419             temp = eta - dtnsq;
00420             if (temp <= 0.f) {
00421                 eta /= 2.f;
00422             }
00423 
00424             tau += eta;
00425             eta /= *sigma + sqrt(eta + *sigma * *sigma);
00426             i__1 = *n;
00427             for (j = 1; j <= i__1; ++j) {
00428                 delta[j] -= eta;
00429                 work[j] += eta;
00430 /* L70: */
00431             }
00432 
00433             *sigma += eta;
00434 
00435 /*           Evaluate PSI and the derivative DPSI */
00436 
00437             dpsi = 0.f;
00438             psi = 0.f;
00439             erretm = 0.f;
00440             i__1 = ii;
00441             for (j = 1; j <= i__1; ++j) {
00442                 temp = z__[j] / (work[j] * delta[j]);
00443                 psi += z__[j] * temp;
00444                 dpsi += temp * temp;
00445                 erretm += psi;
00446 /* L80: */
00447             }
00448             erretm = dabs(erretm);
00449 
00450 /*           Evaluate PHI and the derivative DPHI */
00451 
00452             temp = z__[*n] / (work[*n] * delta[*n]);
00453             phi = z__[*n] * temp;
00454             dphi = temp * temp;
00455             erretm = (-phi - psi) * 8.f + erretm - phi + rhoinv + dabs(tau) * 
00456                     (dpsi + dphi);
00457 
00458             w = rhoinv + phi + psi;
00459 /* L90: */
00460         }
00461 
00462 /*        Return with INFO = 1, NITER = MAXIT and not converged */
00463 
00464         *info = 1;
00465         goto L240;
00466 
00467 /*        End for the case I = N */
00468 
00469     } else {
00470 
00471 /*        The case for I < N */
00472 
00473         niter = 1;
00474         ip1 = *i__ + 1;
00475 
00476 /*        Calculate initial guess */
00477 
00478         delsq = (d__[ip1] - d__[*i__]) * (d__[ip1] + d__[*i__]);
00479         delsq2 = delsq / 2.f;
00480         temp = delsq2 / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + delsq2));
00481         i__1 = *n;
00482         for (j = 1; j <= i__1; ++j) {
00483             work[j] = d__[j] + d__[*i__] + temp;
00484             delta[j] = d__[j] - d__[*i__] - temp;
00485 /* L100: */
00486         }
00487 
00488         psi = 0.f;
00489         i__1 = *i__ - 1;
00490         for (j = 1; j <= i__1; ++j) {
00491             psi += z__[j] * z__[j] / (work[j] * delta[j]);
00492 /* L110: */
00493         }
00494 
00495         phi = 0.f;
00496         i__1 = *i__ + 2;
00497         for (j = *n; j >= i__1; --j) {
00498             phi += z__[j] * z__[j] / (work[j] * delta[j]);
00499 /* L120: */
00500         }
00501         c__ = rhoinv + psi + phi;
00502         w = c__ + z__[*i__] * z__[*i__] / (work[*i__] * delta[*i__]) + z__[
00503                 ip1] * z__[ip1] / (work[ip1] * delta[ip1]);
00504 
00505         if (w > 0.f) {
00506 
00507 /*           d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2 */
00508 
00509 /*           We choose d(i) as origin. */
00510 
00511             orgati = TRUE_;
00512             sg2lb = 0.f;
00513             sg2ub = delsq2;
00514             a = c__ * delsq + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
00515             b = z__[*i__] * z__[*i__] * delsq;
00516             if (a > 0.f) {
00517                 tau = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
00518                         r__1))));
00519             } else {
00520                 tau = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
00521                          (c__ * 2.f);
00522             }
00523 
00524 /*           TAU now is an estimation of SIGMA^2 - D( I )^2. The */
00525 /*           following, however, is the corresponding estimation of */
00526 /*           SIGMA - D( I ). */
00527 
00528             eta = tau / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + tau));
00529         } else {
00530 
00531 /*           (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2 */
00532 
00533 /*           We choose d(i+1) as origin. */
00534 
00535             orgati = FALSE_;
00536             sg2lb = -delsq2;
00537             sg2ub = 0.f;
00538             a = c__ * delsq - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
00539             b = z__[ip1] * z__[ip1] * delsq;
00540             if (a < 0.f) {
00541                 tau = b * 2.f / (a - sqrt((r__1 = a * a + b * 4.f * c__, dabs(
00542                         r__1))));
00543             } else {
00544                 tau = -(a + sqrt((r__1 = a * a + b * 4.f * c__, dabs(r__1)))) 
00545                         / (c__ * 2.f);
00546             }
00547 
00548 /*           TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The */
00549 /*           following, however, is the corresponding estimation of */
00550 /*           SIGMA - D( IP1 ). */
00551 
00552             eta = tau / (d__[ip1] + sqrt((r__1 = d__[ip1] * d__[ip1] + tau, 
00553                     dabs(r__1))));
00554         }
00555 
00556         if (orgati) {
00557             ii = *i__;
00558             *sigma = d__[*i__] + eta;
00559             i__1 = *n;
00560             for (j = 1; j <= i__1; ++j) {
00561                 work[j] = d__[j] + d__[*i__] + eta;
00562                 delta[j] = d__[j] - d__[*i__] - eta;
00563 /* L130: */
00564             }
00565         } else {
00566             ii = *i__ + 1;
00567             *sigma = d__[ip1] + eta;
00568             i__1 = *n;
00569             for (j = 1; j <= i__1; ++j) {
00570                 work[j] = d__[j] + d__[ip1] + eta;
00571                 delta[j] = d__[j] - d__[ip1] - eta;
00572 /* L140: */
00573             }
00574         }
00575         iim1 = ii - 1;
00576         iip1 = ii + 1;
00577 
00578 /*        Evaluate PSI and the derivative DPSI */
00579 
00580         dpsi = 0.f;
00581         psi = 0.f;
00582         erretm = 0.f;
00583         i__1 = iim1;
00584         for (j = 1; j <= i__1; ++j) {
00585             temp = z__[j] / (work[j] * delta[j]);
00586             psi += z__[j] * temp;
00587             dpsi += temp * temp;
00588             erretm += psi;
00589 /* L150: */
00590         }
00591         erretm = dabs(erretm);
00592 
00593 /*        Evaluate PHI and the derivative DPHI */
00594 
00595         dphi = 0.f;
00596         phi = 0.f;
00597         i__1 = iip1;
00598         for (j = *n; j >= i__1; --j) {
00599             temp = z__[j] / (work[j] * delta[j]);
00600             phi += z__[j] * temp;
00601             dphi += temp * temp;
00602             erretm += phi;
00603 /* L160: */
00604         }
00605 
00606         w = rhoinv + phi + psi;
00607 
00608 /*        W is the value of the secular function with */
00609 /*        its ii-th element removed. */
00610 
00611         swtch3 = FALSE_;
00612         if (orgati) {
00613             if (w < 0.f) {
00614                 swtch3 = TRUE_;
00615             }
00616         } else {
00617             if (w > 0.f) {
00618                 swtch3 = TRUE_;
00619             }
00620         }
00621         if (ii == 1 || ii == *n) {
00622             swtch3 = FALSE_;
00623         }
00624 
00625         temp = z__[ii] / (work[ii] * delta[ii]);
00626         dw = dpsi + dphi + temp * temp;
00627         temp = z__[ii] * temp;
00628         w += temp;
00629         erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 3.f 
00630                 + dabs(tau) * dw;
00631 
00632 /*        Test for convergence */
00633 
00634         if (dabs(w) <= eps * erretm) {
00635             goto L240;
00636         }
00637 
00638         if (w <= 0.f) {
00639             sg2lb = dmax(sg2lb,tau);
00640         } else {
00641             sg2ub = dmin(sg2ub,tau);
00642         }
00643 
00644 /*        Calculate the new step */
00645 
00646         ++niter;
00647         if (! swtch3) {
00648             dtipsq = work[ip1] * delta[ip1];
00649             dtisq = work[*i__] * delta[*i__];
00650             if (orgati) {
00651 /* Computing 2nd power */
00652                 r__1 = z__[*i__] / dtisq;
00653                 c__ = w - dtipsq * dw + delsq * (r__1 * r__1);
00654             } else {
00655 /* Computing 2nd power */
00656                 r__1 = z__[ip1] / dtipsq;
00657                 c__ = w - dtisq * dw - delsq * (r__1 * r__1);
00658             }
00659             a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
00660             b = dtipsq * dtisq * w;
00661             if (c__ == 0.f) {
00662                 if (a == 0.f) {
00663                     if (orgati) {
00664                         a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi + 
00665                                 dphi);
00666                     } else {
00667                         a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi + 
00668                                 dphi);
00669                     }
00670                 }
00671                 eta = b / a;
00672             } else if (a <= 0.f) {
00673                 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1)))) /
00674                          (c__ * 2.f);
00675             } else {
00676                 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, dabs(
00677                         r__1))));
00678             }
00679         } else {
00680 
00681 /*           Interpolation using THREE most relevant poles */
00682 
00683             dtiim = work[iim1] * delta[iim1];
00684             dtiip = work[iip1] * delta[iip1];
00685             temp = rhoinv + psi + phi;
00686             if (orgati) {
00687                 temp1 = z__[iim1] / dtiim;
00688                 temp1 *= temp1;
00689                 c__ = temp - dtiip * (dpsi + dphi) - (d__[iim1] - d__[iip1]) *
00690                          (d__[iim1] + d__[iip1]) * temp1;
00691                 zz[0] = z__[iim1] * z__[iim1];
00692                 if (dpsi < temp1) {
00693                     zz[2] = dtiip * dtiip * dphi;
00694                 } else {
00695                     zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
00696                 }
00697             } else {
00698                 temp1 = z__[iip1] / dtiip;
00699                 temp1 *= temp1;
00700                 c__ = temp - dtiim * (dpsi + dphi) - (d__[iip1] - d__[iim1]) *
00701                          (d__[iim1] + d__[iip1]) * temp1;
00702                 if (dphi < temp1) {
00703                     zz[0] = dtiim * dtiim * dpsi;
00704                 } else {
00705                     zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
00706                 }
00707                 zz[2] = z__[iip1] * z__[iip1];
00708             }
00709             zz[1] = z__[ii] * z__[ii];
00710             dd[0] = dtiim;
00711             dd[1] = delta[ii] * work[ii];
00712             dd[2] = dtiip;
00713             slaed6_(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
00714             if (*info != 0) {
00715                 goto L240;
00716             }
00717         }
00718 
00719 /*        Note, eta should be positive if w is negative, and */
00720 /*        eta should be negative otherwise. However, */
00721 /*        if for some reason caused by roundoff, eta*w > 0, */
00722 /*        we simply use one Newton step instead. This way */
00723 /*        will guarantee eta*w < 0. */
00724 
00725         if (w * eta >= 0.f) {
00726             eta = -w / dw;
00727         }
00728         if (orgati) {
00729             temp1 = work[*i__] * delta[*i__];
00730             temp = eta - temp1;
00731         } else {
00732             temp1 = work[ip1] * delta[ip1];
00733             temp = eta - temp1;
00734         }
00735         if (temp > sg2ub || temp < sg2lb) {
00736             if (w < 0.f) {
00737                 eta = (sg2ub - tau) / 2.f;
00738             } else {
00739                 eta = (sg2lb - tau) / 2.f;
00740             }
00741         }
00742 
00743         tau += eta;
00744         eta /= *sigma + sqrt(*sigma * *sigma + eta);
00745 
00746         prew = w;
00747 
00748         *sigma += eta;
00749         i__1 = *n;
00750         for (j = 1; j <= i__1; ++j) {
00751             work[j] += eta;
00752             delta[j] -= eta;
00753 /* L170: */
00754         }
00755 
00756 /*        Evaluate PSI and the derivative DPSI */
00757 
00758         dpsi = 0.f;
00759         psi = 0.f;
00760         erretm = 0.f;
00761         i__1 = iim1;
00762         for (j = 1; j <= i__1; ++j) {
00763             temp = z__[j] / (work[j] * delta[j]);
00764             psi += z__[j] * temp;
00765             dpsi += temp * temp;
00766             erretm += psi;
00767 /* L180: */
00768         }
00769         erretm = dabs(erretm);
00770 
00771 /*        Evaluate PHI and the derivative DPHI */
00772 
00773         dphi = 0.f;
00774         phi = 0.f;
00775         i__1 = iip1;
00776         for (j = *n; j >= i__1; --j) {
00777             temp = z__[j] / (work[j] * delta[j]);
00778             phi += z__[j] * temp;
00779             dphi += temp * temp;
00780             erretm += phi;
00781 /* L190: */
00782         }
00783 
00784         temp = z__[ii] / (work[ii] * delta[ii]);
00785         dw = dpsi + dphi + temp * temp;
00786         temp = z__[ii] * temp;
00787         w = rhoinv + phi + psi + temp;
00788         erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 3.f 
00789                 + dabs(tau) * dw;
00790 
00791         if (w <= 0.f) {
00792             sg2lb = dmax(sg2lb,tau);
00793         } else {
00794             sg2ub = dmin(sg2ub,tau);
00795         }
00796 
00797         swtch = FALSE_;
00798         if (orgati) {
00799             if (-w > dabs(prew) / 10.f) {
00800                 swtch = TRUE_;
00801             }
00802         } else {
00803             if (w > dabs(prew) / 10.f) {
00804                 swtch = TRUE_;
00805             }
00806         }
00807 
00808 /*        Main loop to update the values of the array   DELTA and WORK */
00809 
00810         iter = niter + 1;
00811 
00812         for (niter = iter; niter <= 20; ++niter) {
00813 
00814 /*           Test for convergence */
00815 
00816             if (dabs(w) <= eps * erretm) {
00817                 goto L240;
00818             }
00819 
00820 /*           Calculate the new step */
00821 
00822             if (! swtch3) {
00823                 dtipsq = work[ip1] * delta[ip1];
00824                 dtisq = work[*i__] * delta[*i__];
00825                 if (! swtch) {
00826                     if (orgati) {
00827 /* Computing 2nd power */
00828                         r__1 = z__[*i__] / dtisq;
00829                         c__ = w - dtipsq * dw + delsq * (r__1 * r__1);
00830                     } else {
00831 /* Computing 2nd power */
00832                         r__1 = z__[ip1] / dtipsq;
00833                         c__ = w - dtisq * dw - delsq * (r__1 * r__1);
00834                     }
00835                 } else {
00836                     temp = z__[ii] / (work[ii] * delta[ii]);
00837                     if (orgati) {
00838                         dpsi += temp * temp;
00839                     } else {
00840                         dphi += temp * temp;
00841                     }
00842                     c__ = w - dtisq * dpsi - dtipsq * dphi;
00843                 }
00844                 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
00845                 b = dtipsq * dtisq * w;
00846                 if (c__ == 0.f) {
00847                     if (a == 0.f) {
00848                         if (! swtch) {
00849                             if (orgati) {
00850                                 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * 
00851                                         (dpsi + dphi);
00852                             } else {
00853                                 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (
00854                                         dpsi + dphi);
00855                             }
00856                         } else {
00857                             a = dtisq * dtisq * dpsi + dtipsq * dtipsq * dphi;
00858                         }
00859                     }
00860                     eta = b / a;
00861                 } else if (a <= 0.f) {
00862                     eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, dabs(r__1))
00863                             )) / (c__ * 2.f);
00864                 } else {
00865                     eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, 
00866                             dabs(r__1))));
00867                 }
00868             } else {
00869 
00870 /*              Interpolation using THREE most relevant poles */
00871 
00872                 dtiim = work[iim1] * delta[iim1];
00873                 dtiip = work[iip1] * delta[iip1];
00874                 temp = rhoinv + psi + phi;
00875                 if (swtch) {
00876                     c__ = temp - dtiim * dpsi - dtiip * dphi;
00877                     zz[0] = dtiim * dtiim * dpsi;
00878                     zz[2] = dtiip * dtiip * dphi;
00879                 } else {
00880                     if (orgati) {
00881                         temp1 = z__[iim1] / dtiim;
00882                         temp1 *= temp1;
00883                         temp2 = (d__[iim1] - d__[iip1]) * (d__[iim1] + d__[
00884                                 iip1]) * temp1;
00885                         c__ = temp - dtiip * (dpsi + dphi) - temp2;
00886                         zz[0] = z__[iim1] * z__[iim1];
00887                         if (dpsi < temp1) {
00888                             zz[2] = dtiip * dtiip * dphi;
00889                         } else {
00890                             zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
00891                         }
00892                     } else {
00893                         temp1 = z__[iip1] / dtiip;
00894                         temp1 *= temp1;
00895                         temp2 = (d__[iip1] - d__[iim1]) * (d__[iim1] + d__[
00896                                 iip1]) * temp1;
00897                         c__ = temp - dtiim * (dpsi + dphi) - temp2;
00898                         if (dphi < temp1) {
00899                             zz[0] = dtiim * dtiim * dpsi;
00900                         } else {
00901                             zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
00902                         }
00903                         zz[2] = z__[iip1] * z__[iip1];
00904                     }
00905                 }
00906                 dd[0] = dtiim;
00907                 dd[1] = delta[ii] * work[ii];
00908                 dd[2] = dtiip;
00909                 slaed6_(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
00910                 if (*info != 0) {
00911                     goto L240;
00912                 }
00913             }
00914 
00915 /*           Note, eta should be positive if w is negative, and */
00916 /*           eta should be negative otherwise. However, */
00917 /*           if for some reason caused by roundoff, eta*w > 0, */
00918 /*           we simply use one Newton step instead. This way */
00919 /*           will guarantee eta*w < 0. */
00920 
00921             if (w * eta >= 0.f) {
00922                 eta = -w / dw;
00923             }
00924             if (orgati) {
00925                 temp1 = work[*i__] * delta[*i__];
00926                 temp = eta - temp1;
00927             } else {
00928                 temp1 = work[ip1] * delta[ip1];
00929                 temp = eta - temp1;
00930             }
00931             if (temp > sg2ub || temp < sg2lb) {
00932                 if (w < 0.f) {
00933                     eta = (sg2ub - tau) / 2.f;
00934                 } else {
00935                     eta = (sg2lb - tau) / 2.f;
00936                 }
00937             }
00938 
00939             tau += eta;
00940             eta /= *sigma + sqrt(*sigma * *sigma + eta);
00941 
00942             *sigma += eta;
00943             i__1 = *n;
00944             for (j = 1; j <= i__1; ++j) {
00945                 work[j] += eta;
00946                 delta[j] -= eta;
00947 /* L200: */
00948             }
00949 
00950             prew = w;
00951 
00952 /*           Evaluate PSI and the derivative DPSI */
00953 
00954             dpsi = 0.f;
00955             psi = 0.f;
00956             erretm = 0.f;
00957             i__1 = iim1;
00958             for (j = 1; j <= i__1; ++j) {
00959                 temp = z__[j] / (work[j] * delta[j]);
00960                 psi += z__[j] * temp;
00961                 dpsi += temp * temp;
00962                 erretm += psi;
00963 /* L210: */
00964             }
00965             erretm = dabs(erretm);
00966 
00967 /*           Evaluate PHI and the derivative DPHI */
00968 
00969             dphi = 0.f;
00970             phi = 0.f;
00971             i__1 = iip1;
00972             for (j = *n; j >= i__1; --j) {
00973                 temp = z__[j] / (work[j] * delta[j]);
00974                 phi += z__[j] * temp;
00975                 dphi += temp * temp;
00976                 erretm += phi;
00977 /* L220: */
00978             }
00979 
00980             temp = z__[ii] / (work[ii] * delta[ii]);
00981             dw = dpsi + dphi + temp * temp;
00982             temp = z__[ii] * temp;
00983             w = rhoinv + phi + psi + temp;
00984             erretm = (phi - psi) * 8.f + erretm + rhoinv * 2.f + dabs(temp) * 
00985                     3.f + dabs(tau) * dw;
00986             if (w * prew > 0.f && dabs(w) > dabs(prew) / 10.f) {
00987                 swtch = ! swtch;
00988             }
00989 
00990             if (w <= 0.f) {
00991                 sg2lb = dmax(sg2lb,tau);
00992             } else {
00993                 sg2ub = dmin(sg2ub,tau);
00994             }
00995 
00996 /* L230: */
00997         }
00998 
00999 /*        Return with INFO = 1, NITER = MAXIT and not converged */
01000 
01001         *info = 1;
01002 
01003     }
01004 
01005 L240:
01006     return 0;
01007 
01008 /*     End of SLASD4 */
01009 
01010 } /* slasd4_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:11