dlasq4.c
Go to the documentation of this file.
00001 /* dlasq4.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 dlasq4_(integer *i0, integer *n0, doublereal *z__, 
00017         integer *pp, integer *n0in, doublereal *dmin__, doublereal *dmin1, 
00018         doublereal *dmin2, doublereal *dn, doublereal *dn1, doublereal *dn2, 
00019         doublereal *tau, integer *ttype, doublereal *g)
00020 {
00021     /* System generated locals */
00022     integer i__1;
00023     doublereal d__1, d__2;
00024 
00025     /* Builtin functions */
00026     double sqrt(doublereal);
00027 
00028     /* Local variables */
00029     doublereal s, a2, b1, b2;
00030     integer i4, nn, np;
00031     doublereal gam, gap1, gap2;
00032 
00033 
00034 /*  -- LAPACK routine (version 3.2)                                    -- */
00035 
00036 /*  -- Contributed by Osni Marques of the Lawrence Berkeley National   -- */
00037 /*  -- Laboratory and Beresford Parlett of the Univ. of California at  -- */
00038 /*  -- Berkeley                                                        -- */
00039 /*  -- November 2008                                                   -- */
00040 
00041 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00042 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00043 
00044 /*     .. Scalar Arguments .. */
00045 /*     .. */
00046 /*     .. Array Arguments .. */
00047 /*     .. */
00048 
00049 /*  Purpose */
00050 /*  ======= */
00051 
00052 /*  DLASQ4 computes an approximation TAU to the smallest eigenvalue */
00053 /*  using values of d from the previous transform. */
00054 
00055 /*  I0    (input) INTEGER */
00056 /*        First index. */
00057 
00058 /*  N0    (input) INTEGER */
00059 /*        Last index. */
00060 
00061 /*  Z     (input) DOUBLE PRECISION array, dimension ( 4*N ) */
00062 /*        Z holds the qd array. */
00063 
00064 /*  PP    (input) INTEGER */
00065 /*        PP=0 for ping, PP=1 for pong. */
00066 
00067 /*  NOIN  (input) INTEGER */
00068 /*        The value of N0 at start of EIGTEST. */
00069 
00070 /*  DMIN  (input) DOUBLE PRECISION */
00071 /*        Minimum value of d. */
00072 
00073 /*  DMIN1 (input) DOUBLE PRECISION */
00074 /*        Minimum value of d, excluding D( N0 ). */
00075 
00076 /*  DMIN2 (input) DOUBLE PRECISION */
00077 /*        Minimum value of d, excluding D( N0 ) and D( N0-1 ). */
00078 
00079 /*  DN    (input) DOUBLE PRECISION */
00080 /*        d(N) */
00081 
00082 /*  DN1   (input) DOUBLE PRECISION */
00083 /*        d(N-1) */
00084 
00085 /*  DN2   (input) DOUBLE PRECISION */
00086 /*        d(N-2) */
00087 
00088 /*  TAU   (output) DOUBLE PRECISION */
00089 /*        This is the shift. */
00090 
00091 /*  TTYPE (output) INTEGER */
00092 /*        Shift type. */
00093 
00094 /*  G     (input/output) REAL */
00095 /*        G is passed as an argument in order to save its value between */
00096 /*        calls to DLASQ4. */
00097 
00098 /*  Further Details */
00099 /*  =============== */
00100 /*  CNST1 = 9/16 */
00101 
00102 /*  ===================================================================== */
00103 
00104 /*     .. Parameters .. */
00105 /*     .. */
00106 /*     .. Local Scalars .. */
00107 /*     .. */
00108 /*     .. Intrinsic Functions .. */
00109 /*     .. */
00110 /*     .. Executable Statements .. */
00111 
00112 /*     A negative DMIN forces the shift to take that absolute value */
00113 /*     TTYPE records the type of shift. */
00114 
00115     /* Parameter adjustments */
00116     --z__;
00117 
00118     /* Function Body */
00119     if (*dmin__ <= 0.) {
00120         *tau = -(*dmin__);
00121         *ttype = -1;
00122         return 0;
00123     }
00124 
00125     nn = (*n0 << 2) + *pp;
00126     if (*n0in == *n0) {
00127 
00128 /*        No eigenvalues deflated. */
00129 
00130         if (*dmin__ == *dn || *dmin__ == *dn1) {
00131 
00132             b1 = sqrt(z__[nn - 3]) * sqrt(z__[nn - 5]);
00133             b2 = sqrt(z__[nn - 7]) * sqrt(z__[nn - 9]);
00134             a2 = z__[nn - 7] + z__[nn - 5];
00135 
00136 /*           Cases 2 and 3. */
00137 
00138             if (*dmin__ == *dn && *dmin1 == *dn1) {
00139                 gap2 = *dmin2 - a2 - *dmin2 * .25;
00140                 if (gap2 > 0. && gap2 > b2) {
00141                     gap1 = a2 - *dn - b2 / gap2 * b2;
00142                 } else {
00143                     gap1 = a2 - *dn - (b1 + b2);
00144                 }
00145                 if (gap1 > 0. && gap1 > b1) {
00146 /* Computing MAX */
00147                     d__1 = *dn - b1 / gap1 * b1, d__2 = *dmin__ * .5;
00148                     s = max(d__1,d__2);
00149                     *ttype = -2;
00150                 } else {
00151                     s = 0.;
00152                     if (*dn > b1) {
00153                         s = *dn - b1;
00154                     }
00155                     if (a2 > b1 + b2) {
00156 /* Computing MIN */
00157                         d__1 = s, d__2 = a2 - (b1 + b2);
00158                         s = min(d__1,d__2);
00159                     }
00160 /* Computing MAX */
00161                     d__1 = s, d__2 = *dmin__ * .333;
00162                     s = max(d__1,d__2);
00163                     *ttype = -3;
00164                 }
00165             } else {
00166 
00167 /*              Case 4. */
00168 
00169                 *ttype = -4;
00170                 s = *dmin__ * .25;
00171                 if (*dmin__ == *dn) {
00172                     gam = *dn;
00173                     a2 = 0.;
00174                     if (z__[nn - 5] > z__[nn - 7]) {
00175                         return 0;
00176                     }
00177                     b2 = z__[nn - 5] / z__[nn - 7];
00178                     np = nn - 9;
00179                 } else {
00180                     np = nn - (*pp << 1);
00181                     b2 = z__[np - 2];
00182                     gam = *dn1;
00183                     if (z__[np - 4] > z__[np - 2]) {
00184                         return 0;
00185                     }
00186                     a2 = z__[np - 4] / z__[np - 2];
00187                     if (z__[nn - 9] > z__[nn - 11]) {
00188                         return 0;
00189                     }
00190                     b2 = z__[nn - 9] / z__[nn - 11];
00191                     np = nn - 13;
00192                 }
00193 
00194 /*              Approximate contribution to norm squared from I < NN-1. */
00195 
00196                 a2 += b2;
00197                 i__1 = (*i0 << 2) - 1 + *pp;
00198                 for (i4 = np; i4 >= i__1; i4 += -4) {
00199                     if (b2 == 0.) {
00200                         goto L20;
00201                     }
00202                     b1 = b2;
00203                     if (z__[i4] > z__[i4 - 2]) {
00204                         return 0;
00205                     }
00206                     b2 *= z__[i4] / z__[i4 - 2];
00207                     a2 += b2;
00208                     if (max(b2,b1) * 100. < a2 || .563 < a2) {
00209                         goto L20;
00210                     }
00211 /* L10: */
00212                 }
00213 L20:
00214                 a2 *= 1.05;
00215 
00216 /*              Rayleigh quotient residual bound. */
00217 
00218                 if (a2 < .563) {
00219                     s = gam * (1. - sqrt(a2)) / (a2 + 1.);
00220                 }
00221             }
00222         } else if (*dmin__ == *dn2) {
00223 
00224 /*           Case 5. */
00225 
00226             *ttype = -5;
00227             s = *dmin__ * .25;
00228 
00229 /*           Compute contribution to norm squared from I > NN-2. */
00230 
00231             np = nn - (*pp << 1);
00232             b1 = z__[np - 2];
00233             b2 = z__[np - 6];
00234             gam = *dn2;
00235             if (z__[np - 8] > b2 || z__[np - 4] > b1) {
00236                 return 0;
00237             }
00238             a2 = z__[np - 8] / b2 * (z__[np - 4] / b1 + 1.);
00239 
00240 /*           Approximate contribution to norm squared from I < NN-2. */
00241 
00242             if (*n0 - *i0 > 2) {
00243                 b2 = z__[nn - 13] / z__[nn - 15];
00244                 a2 += b2;
00245                 i__1 = (*i0 << 2) - 1 + *pp;
00246                 for (i4 = nn - 17; i4 >= i__1; i4 += -4) {
00247                     if (b2 == 0.) {
00248                         goto L40;
00249                     }
00250                     b1 = b2;
00251                     if (z__[i4] > z__[i4 - 2]) {
00252                         return 0;
00253                     }
00254                     b2 *= z__[i4] / z__[i4 - 2];
00255                     a2 += b2;
00256                     if (max(b2,b1) * 100. < a2 || .563 < a2) {
00257                         goto L40;
00258                     }
00259 /* L30: */
00260                 }
00261 L40:
00262                 a2 *= 1.05;
00263             }
00264 
00265             if (a2 < .563) {
00266                 s = gam * (1. - sqrt(a2)) / (a2 + 1.);
00267             }
00268         } else {
00269 
00270 /*           Case 6, no information to guide us. */
00271 
00272             if (*ttype == -6) {
00273                 *g += (1. - *g) * .333;
00274             } else if (*ttype == -18) {
00275                 *g = .083250000000000005;
00276             } else {
00277                 *g = .25;
00278             }
00279             s = *g * *dmin__;
00280             *ttype = -6;
00281         }
00282 
00283     } else if (*n0in == *n0 + 1) {
00284 
00285 /*        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. */
00286 
00287         if (*dmin1 == *dn1 && *dmin2 == *dn2) {
00288 
00289 /*           Cases 7 and 8. */
00290 
00291             *ttype = -7;
00292             s = *dmin1 * .333;
00293             if (z__[nn - 5] > z__[nn - 7]) {
00294                 return 0;
00295             }
00296             b1 = z__[nn - 5] / z__[nn - 7];
00297             b2 = b1;
00298             if (b2 == 0.) {
00299                 goto L60;
00300             }
00301             i__1 = (*i0 << 2) - 1 + *pp;
00302             for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
00303                 a2 = b1;
00304                 if (z__[i4] > z__[i4 - 2]) {
00305                     return 0;
00306                 }
00307                 b1 *= z__[i4] / z__[i4 - 2];
00308                 b2 += b1;
00309                 if (max(b1,a2) * 100. < b2) {
00310                     goto L60;
00311                 }
00312 /* L50: */
00313             }
00314 L60:
00315             b2 = sqrt(b2 * 1.05);
00316 /* Computing 2nd power */
00317             d__1 = b2;
00318             a2 = *dmin1 / (d__1 * d__1 + 1.);
00319             gap2 = *dmin2 * .5 - a2;
00320             if (gap2 > 0. && gap2 > b2 * a2) {
00321 /* Computing MAX */
00322                 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
00323                 s = max(d__1,d__2);
00324             } else {
00325 /* Computing MAX */
00326                 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
00327                 s = max(d__1,d__2);
00328                 *ttype = -8;
00329             }
00330         } else {
00331 
00332 /*           Case 9. */
00333 
00334             s = *dmin1 * .25;
00335             if (*dmin1 == *dn1) {
00336                 s = *dmin1 * .5;
00337             }
00338             *ttype = -9;
00339         }
00340 
00341     } else if (*n0in == *n0 + 2) {
00342 
00343 /*        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN. */
00344 
00345 /*        Cases 10 and 11. */
00346 
00347         if (*dmin2 == *dn2 && z__[nn - 5] * 2. < z__[nn - 7]) {
00348             *ttype = -10;
00349             s = *dmin2 * .333;
00350             if (z__[nn - 5] > z__[nn - 7]) {
00351                 return 0;
00352             }
00353             b1 = z__[nn - 5] / z__[nn - 7];
00354             b2 = b1;
00355             if (b2 == 0.) {
00356                 goto L80;
00357             }
00358             i__1 = (*i0 << 2) - 1 + *pp;
00359             for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
00360                 if (z__[i4] > z__[i4 - 2]) {
00361                     return 0;
00362                 }
00363                 b1 *= z__[i4] / z__[i4 - 2];
00364                 b2 += b1;
00365                 if (b1 * 100. < b2) {
00366                     goto L80;
00367                 }
00368 /* L70: */
00369             }
00370 L80:
00371             b2 = sqrt(b2 * 1.05);
00372 /* Computing 2nd power */
00373             d__1 = b2;
00374             a2 = *dmin2 / (d__1 * d__1 + 1.);
00375             gap2 = z__[nn - 7] + z__[nn - 9] - sqrt(z__[nn - 11]) * sqrt(z__[
00376                     nn - 9]) - a2;
00377             if (gap2 > 0. && gap2 > b2 * a2) {
00378 /* Computing MAX */
00379                 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
00380                 s = max(d__1,d__2);
00381             } else {
00382 /* Computing MAX */
00383                 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
00384                 s = max(d__1,d__2);
00385             }
00386         } else {
00387             s = *dmin2 * .25;
00388             *ttype = -11;
00389         }
00390     } else if (*n0in > *n0 + 2) {
00391 
00392 /*        Case 12, more than two eigenvalues deflated. No information. */
00393 
00394         s = 0.;
00395         *ttype = -12;
00396     }
00397 
00398     *tau = s;
00399     return 0;
00400 
00401 /*     End of DLASQ4 */
00402 
00403 } /* dlasq4_ */

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:46