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

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