git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlasq4.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlasq4_(integer *i0, integer *n0, doublereal *z__,
   4         integer *pp, integer *n0in, doublereal *dmin__, doublereal *dmin1,
   5         doublereal *dmin2, doublereal *dn, doublereal *dn1, doublereal *dn2,
   6         doublereal *tau, integer *ttype)
   7 {
   8     /* Initialized data */
   9
  10     static doublereal g = 0.;
  11
  12     /* System generated locals */
  13     integer i__1;
  14     doublereal d__1, d__2;
  15
  16     /* Builtin functions */
  17     double sqrt(doublereal);
  18
  19     /* Local variables */
  20     doublereal s, a2, b1, b2;
  21     integer i4, nn, np;
  22     doublereal gam, gap1, gap2;
  23
  24
  25 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  26 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  27 /*     November 2006 */
  28
  29 /*     .. Scalar Arguments .. */
  30 /*     .. */
  31 /*     .. Array Arguments .. */
  32 /*     .. */
  33
  34 /*  Purpose */
  35 /*  ======= */
  36
  37 /*  DLASQ4 computes an approximation TAU to the smallest eigenvalue */
  38 /*  using values of d from the previous transform. */
  39
  40 /*  I0    (input) INTEGER */
  41 /*        First index. */
  42
  43 /*  N0    (input) INTEGER */
  44 /*        Last index. */
  45
  46 /*  Z     (input) DOUBLE PRECISION array, dimension ( 4*N ) */
  47 /*        Z holds the qd array. */
  48
  49 /*  PP    (input) INTEGER */
  50 /*        PP=0 for ping, PP=1 for pong. */
  51
  52 /*  N0IN  (input) INTEGER */
  53 /*        The value of N0 at start of EIGTEST. */
  54
  55 /*  DMIN  (input) DOUBLE PRECISION */
  56 /*        Minimum value of d. */
  57
  58 /*  DMIN1 (input) DOUBLE PRECISION */
  59 /*        Minimum value of d, excluding D( N0 ). */
  60
  61 /*  DMIN2 (input) DOUBLE PRECISION */
  62 /*        Minimum value of d, excluding D( N0 ) and D( N0-1 ). */
  63
  64 /*  DN    (input) DOUBLE PRECISION */
  65 /*        d(N) */
  66
  67 /*  DN1   (input) DOUBLE PRECISION */
  68 /*        d(N-1) */
  69
  70 /*  DN2   (input) DOUBLE PRECISION */
  71 /*        d(N-2) */
  72
  73 /*  TAU   (output) DOUBLE PRECISION */
  74 /*        This is the shift. */
  75
  76 /*  TTYPE (output) INTEGER */
  77 /*        Shift type. */
  78
  79 /*  Further Details */
  80 /*  =============== */
  81 /*  CNST1 = 9/16 */
  82
  83 /*  ===================================================================== */
  84
  85 /*     .. Parameters .. */
  86 /*     .. */
  87 /*     .. Local Scalars .. */
  88 /*     .. */
  89 /*     .. Intrinsic Functions .. */
  90 /*     .. */
  91 /*     .. Save statement .. */
  92 /*     .. */
  93 /*     .. Data statement .. */
  94     /* Parameter adjustments */
  95     --z__;
  96
  97     /* Function Body */
  98 /*     .. */
  99 /*     .. Executable Statements .. */
 100
 101 /*     A negative DMIN forces the shift to take that absolute value */
 102 /*     TTYPE records the type of shift. */
 103
 104     if (*dmin__ <= 0.) {
 105         *tau = -(*dmin__);
 106         *ttype = -1;
 107         return 0;
 108     }
 109
 110     nn = (*n0 << 2) + *pp;
 111     if (*n0in == *n0) {
 112
 113 /*        No eigenvalues deflated. */
 114
 115         if (*dmin__ == *dn || *dmin__ == *dn1) {
 116
 117             b1 = sqrt(z__[nn - 3]) * sqrt(z__[nn - 5]);
 118             b2 = sqrt(z__[nn - 7]) * sqrt(z__[nn - 9]);
 119             a2 = z__[nn - 7] + z__[nn - 5];
 120
 121 /*           Cases 2 and 3. */
 122
 123             if (*dmin__ == *dn && *dmin1 == *dn1) {
 124                 gap2 = *dmin2 - a2 - *dmin2 * .25;
 125                 if (gap2 > 0. && gap2 > b2) {
 126                     gap1 = a2 - *dn - b2 / gap2 * b2;
 127                 } else {
 128                     gap1 = a2 - *dn - (b1 + b2);
 129                 }
 130                 if (gap1 > 0. && gap1 > b1) {
 131 /* Computing MAX */
 132                     d__1 = *dn - b1 / gap1 * b1, d__2 = *dmin__ * .5;
 133                     s = max(d__1,d__2);
 134                     *ttype = -2;
 135                 } else {
 136                     s = 0.;
 137                     if (*dn > b1) {
 138                         s = *dn - b1;
 139                     }
 140                     if (a2 > b1 + b2) {
 141 /* Computing MIN */
 142                         d__1 = s, d__2 = a2 - (b1 + b2);
 143                         s = min(d__1,d__2);
 144                     }
 145 /* Computing MAX */
 146                     d__1 = s, d__2 = *dmin__ * .333;
 147                     s = max(d__1,d__2);
 148                     *ttype = -3;
 149                 }
 150             } else {
 151
 152 /*              Case 4. */
 153
 154                 *ttype = -4;
 155                 s = *dmin__ * .25;
 156                 if (*dmin__ == *dn) {
 157                     gam = *dn;
 158                     a2 = 0.;
 159                     if (z__[nn - 5] > z__[nn - 7]) {
 160                         return 0;
 161                     }
 162                     b2 = z__[nn - 5] / z__[nn - 7];
 163                     np = nn - 9;
 164                 } else {
 165                     np = nn - (*pp << 1);
 166                     b2 = z__[np - 2];
 167                     gam = *dn1;
 168                     if (z__[np - 4] > z__[np - 2]) {
 169                         return 0;
 170                     }
 171                     a2 = z__[np - 4] / z__[np - 2];
 172                     if (z__[nn - 9] > z__[nn - 11]) {
 173                         return 0;
 174                     }
 175                     b2 = z__[nn - 9] / z__[nn - 11];
 176                     np = nn - 13;
 177                 }
 178
 179 /*              Approximate contribution to norm squared from I < NN-1. */
 180
 181                 a2 += b2;
 182                 i__1 = (*i0 << 2) - 1 + *pp;
 183                 for (i4 = np; i4 >= i__1; i4 += -4) {
 184                     if (b2 == 0.) {
 185                         goto L20;
 186                     }
 187                     b1 = b2;
 188                     if (z__[i4] > z__[i4 - 2]) {
 189                         return 0;
 190                     }
 191                     b2 *= z__[i4] / z__[i4 - 2];
 192                     a2 += b2;
 193                     if (max(b2,b1) * 100. < a2 || .563 < a2) {
 194                         goto L20;
 195                     }
 196 /* L10: */
 197                 }
 198 L20:
 199                 a2 *= 1.05;
 200
 201 /*              Rayleigh quotient residual bound. */
 202
 203                 if (a2 < .563) {
 204                     s = gam * (1. - sqrt(a2)) / (a2 + 1.);
 205                 }
 206             }
 207         } else if (*dmin__ == *dn2) {
 208
 209 /*           Case 5. */
 210
 211             *ttype = -5;
 212             s = *dmin__ * .25;
 213
 214 /*           Compute contribution to norm squared from I > NN-2. */
 215
 216             np = nn - (*pp << 1);
 217             b1 = z__[np - 2];
 218             b2 = z__[np - 6];
 219             gam = *dn2;
 220             if (z__[np - 8] > b2 || z__[np - 4] > b1) {
 221                 return 0;
 222             }
 223             a2 = z__[np - 8] / b2 * (z__[np - 4] / b1 + 1.);
 224
 225 /*           Approximate contribution to norm squared from I < NN-2. */
 226
 227             if (*n0 - *i0 > 2) {
 228                 b2 = z__[nn - 13] / z__[nn - 15];
 229                 a2 += b2;
 230                 i__1 = (*i0 << 2) - 1 + *pp;
 231                 for (i4 = nn - 17; i4 >= i__1; i4 += -4) {
 232                     if (b2 == 0.) {
 233                         goto L40;
 234                     }
 235                     b1 = b2;
 236                     if (z__[i4] > z__[i4 - 2]) {
 237                         return 0;
 238                     }
 239                     b2 *= z__[i4] / z__[i4 - 2];
 240                     a2 += b2;
 241                     if (max(b2,b1) * 100. < a2 || .563 < a2) {
 242                         goto L40;
 243                     }
 244 /* L30: */
 245                 }
 246 L40:
 247                 a2 *= 1.05;
 248             }
 249
 250             if (a2 < .563) {
 251                 s = gam * (1. - sqrt(a2)) / (a2 + 1.);
 252             }
 253         } else {
 254
 255 /*           Case 6, no information to guide us. */
 256
 257             if (*ttype == -6) {
 258                 g += (1. - g) * .333;
 259             } else if (*ttype == -18) {
 260                 g = .083250000000000005;
 261             } else {
 262                 g = .25;
 263             }
 264             s = g * *dmin__;
 265             *ttype = -6;
 266         }
 267
 268     } else if (*n0in == *n0 + 1) {
 269
 270 /*        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. */
 271
 272         if (*dmin1 == *dn1 && *dmin2 == *dn2) {
 273
 274 /*           Cases 7 and 8. */
 275
 276             *ttype = -7;
 277             s = *dmin1 * .333;
 278             if (z__[nn - 5] > z__[nn - 7]) {
 279                 return 0;
 280             }
 281             b1 = z__[nn - 5] / z__[nn - 7];
 282             b2 = b1;
 283             if (b2 == 0.) {
 284                 goto L60;
 285             }
 286             i__1 = (*i0 << 2) - 1 + *pp;
 287             for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
 288                 a2 = b1;
 289                 if (z__[i4] > z__[i4 - 2]) {
 290                     return 0;
 291                 }
 292                 b1 *= z__[i4] / z__[i4 - 2];
 293                 b2 += b1;
 294                 if (max(b1,a2) * 100. < b2) {
 295                     goto L60;
 296                 }
 297 /* L50: */
 298             }
 299 L60:
 300             b2 = sqrt(b2 * 1.05);
 301 /* Computing 2nd power */
 302             d__1 = b2;
 303             a2 = *dmin1 / (d__1 * d__1 + 1.);
 304             gap2 = *dmin2 * .5 - a2;
 305             if (gap2 > 0. && gap2 > b2 * a2) {
 306 /* Computing MAX */
 307                 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
 308                 s = max(d__1,d__2);
 309             } else {
 310 /* Computing MAX */
 311                 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
 312                 s = max(d__1,d__2);
 313                 *ttype = -8;
 314             }
 315         } else {
 316
 317 /*           Case 9. */
 318
 319             s = *dmin1 * .25;
 320             if (*dmin1 == *dn1) {
 321                 s = *dmin1 * .5;
 322             }
 323             *ttype = -9;
 324         }
 325
 326     } else if (*n0in == *n0 + 2) {
 327
 328 /*        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN. */
 329
 330 /*        Cases 10 and 11. */
 331
 332         if (*dmin2 == *dn2 && z__[nn - 5] * 2. < z__[nn - 7]) {
 333             *ttype = -10;
 334             s = *dmin2 * .333;
 335             if (z__[nn - 5] > z__[nn - 7]) {
 336                 return 0;
 337             }
 338             b1 = z__[nn - 5] / z__[nn - 7];
 339             b2 = b1;
 340             if (b2 == 0.) {
 341                 goto L80;
 342             }
 343             i__1 = (*i0 << 2) - 1 + *pp;
 344             for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
 345                 if (z__[i4] > z__[i4 - 2]) {
 346                     return 0;
 347                 }
 348                 b1 *= z__[i4] / z__[i4 - 2];
 349                 b2 += b1;
 350                 if (b1 * 100. < b2) {
 351                     goto L80;
 352                 }
 353 /* L70: */
 354             }
 355 L80:
 356             b2 = sqrt(b2 * 1.05);
 357 /* Computing 2nd power */
 358             d__1 = b2;
 359             a2 = *dmin2 / (d__1 * d__1 + 1.);
 360             gap2 = z__[nn - 7] + z__[nn - 9] - sqrt(z__[nn - 11]) * sqrt(z__[
 361                     nn - 9]) - a2;
 362             if (gap2 > 0. && gap2 > b2 * a2) {
 363 /* Computing MAX */
 364                 d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
 365                 s = max(d__1,d__2);
 366             } else {
 367 /* Computing MAX */
 368                 d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
 369                 s = max(d__1,d__2);
 370             }
 371         } else {
 372             s = *dmin2 * .25;
 373             *ttype = -11;
 374         }
 375     } else if (*n0in > *n0 + 2) {
 376
 377 /*        Case 12, more than two eigenvalues deflated. No information. */
 378
 379         s = 0.;
 380         *ttype = -12;
 381     }
 382
 383     *tau = s;
 384     return 0;
 385
 386 /*     End of DLASQ4 */
 387
 388 } /* dlasq4_ */