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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlasq6_(integer *i0, integer *n0, doublereal *z__,
   4         integer *pp, doublereal *dmin__, doublereal *dmin1, doublereal *dmin2,
   5          doublereal *dn, doublereal *dnm1, doublereal *dnm2)
   6 {
   7     /* System generated locals */
   8     integer i__1;
   9     doublereal d__1, d__2;
  10
  11     /* Local variables */
  12     doublereal d__;
  13     integer j4, j4p2;
  14     doublereal emin, temp;
  15     extern doublereal dlamch_(char *);
  16     doublereal safmin;
  17
  18
  19 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  20 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  21 /*     November 2006 */
  22
  23 /*     .. Scalar Arguments .. */
  24 /*     .. */
  25 /*     .. Array Arguments .. */
  26 /*     .. */
  27
  28 /*  Purpose */
  29 /*  ======= */
  30
  31 /*  DLASQ6 computes one dqd (shift equal to zero) transform in */
  32 /*  ping-pong form, with protection against underflow and overflow. */
  33
  34 /*  Arguments */
  35 /*  ========= */
  36
  37 /*  I0    (input) INTEGER */
  38 /*        First index. */
  39
  40 /*  N0    (input) INTEGER */
  41 /*        Last index. */
  42
  43 /*  Z     (input) DOUBLE PRECISION array, dimension ( 4*N ) */
  44 /*        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid */
  45 /*        an extra argument. */
  46
  47 /*  PP    (input) INTEGER */
  48 /*        PP=0 for ping, PP=1 for pong. */
  49
  50 /*  DMIN  (output) DOUBLE PRECISION */
  51 /*        Minimum value of d. */
  52
  53 /*  DMIN1 (output) DOUBLE PRECISION */
  54 /*        Minimum value of d, excluding D( N0 ). */
  55
  56 /*  DMIN2 (output) DOUBLE PRECISION */
  57 /*        Minimum value of d, excluding D( N0 ) and D( N0-1 ). */
  58
  59 /*  DN    (output) DOUBLE PRECISION */
  60 /*        d(N0), the last value of d. */
  61
  62 /*  DNM1  (output) DOUBLE PRECISION */
  63 /*        d(N0-1). */
  64
  65 /*  DNM2  (output) DOUBLE PRECISION */
  66 /*        d(N0-2). */
  67
  68 /*  ===================================================================== */
  69
  70 /*     .. Parameter .. */
  71 /*     .. */
  72 /*     .. Local Scalars .. */
  73 /*     .. */
  74 /*     .. External Function .. */
  75 /*     .. */
  76 /*     .. Intrinsic Functions .. */
  77 /*     .. */
  78 /*     .. Executable Statements .. */
  79
  80     /* Parameter adjustments */
  81     --z__;
  82
  83     /* Function Body */
  84     if (*n0 - *i0 - 1 <= 0) {
  85         return 0;
  86     }
  87
  88     safmin = dlamch_("Safe minimum");
  89     j4 = (*i0 << 2) + *pp - 3;
  90     emin = z__[j4 + 4];
  91     d__ = z__[j4];
  92     *dmin__ = d__;
  93
  94     if (*pp == 0) {
  95         i__1 = *n0 - 3 << 2;
  96         for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
  97             z__[j4 - 2] = d__ + z__[j4 - 1];
  98             if (z__[j4 - 2] == 0.) {
  99                 z__[j4] = 0.;
 100                 d__ = z__[j4 + 1];
 101                 *dmin__ = d__;
 102                 emin = 0.;
 103             } else if (safmin * z__[j4 + 1] < z__[j4 - 2] && safmin * z__[j4
 104                     - 2] < z__[j4 + 1]) {
 105                 temp = z__[j4 + 1] / z__[j4 - 2];
 106                 z__[j4] = z__[j4 - 1] * temp;
 107                 d__ *= temp;
 108             } else {
 109                 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
 110                 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]);
 111             }
 112             *dmin__ = min(*dmin__,d__);
 113 /* Computing MIN */
 114             d__1 = emin, d__2 = z__[j4];
 115             emin = min(d__1,d__2);
 116 /* L10: */
 117         }
 118     } else {
 119         i__1 = *n0 - 3 << 2;
 120         for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
 121             z__[j4 - 3] = d__ + z__[j4];
 122             if (z__[j4 - 3] == 0.) {
 123                 z__[j4 - 1] = 0.;
 124                 d__ = z__[j4 + 2];
 125                 *dmin__ = d__;
 126                 emin = 0.;
 127             } else if (safmin * z__[j4 + 2] < z__[j4 - 3] && safmin * z__[j4
 128                     - 3] < z__[j4 + 2]) {
 129                 temp = z__[j4 + 2] / z__[j4 - 3];
 130                 z__[j4 - 1] = z__[j4] * temp;
 131                 d__ *= temp;
 132             } else {
 133                 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
 134                 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]);
 135             }
 136             *dmin__ = min(*dmin__,d__);
 137 /* Computing MIN */
 138             d__1 = emin, d__2 = z__[j4 - 1];
 139             emin = min(d__1,d__2);
 140 /* L20: */
 141         }
 142     }
 143
 144 /*     Unroll last two steps. */
 145
 146     *dnm2 = d__;
 147     *dmin2 = *dmin__;
 148     j4 = (*n0 - 2 << 2) - *pp;
 149     j4p2 = j4 + (*pp << 1) - 1;
 150     z__[j4 - 2] = *dnm2 + z__[j4p2];
 151     if (z__[j4 - 2] == 0.) {
 152         z__[j4] = 0.;
 153         *dnm1 = z__[j4p2 + 2];
 154         *dmin__ = *dnm1;
 155         emin = 0.;
 156     } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] <
 157             z__[j4p2 + 2]) {
 158         temp = z__[j4p2 + 2] / z__[j4 - 2];
 159         z__[j4] = z__[j4p2] * temp;
 160         *dnm1 = *dnm2 * temp;
 161     } else {
 162         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
 163         *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]);
 164     }
 165     *dmin__ = min(*dmin__,*dnm1);
 166
 167     *dmin1 = *dmin__;
 168     j4 += 4;
 169     j4p2 = j4 + (*pp << 1) - 1;
 170     z__[j4 - 2] = *dnm1 + z__[j4p2];
 171     if (z__[j4 - 2] == 0.) {
 172         z__[j4] = 0.;
 173         *dn = z__[j4p2 + 2];
 174         *dmin__ = *dn;
 175         emin = 0.;
 176     } else if (safmin * z__[j4p2 + 2] < z__[j4 - 2] && safmin * z__[j4 - 2] <
 177             z__[j4p2 + 2]) {
 178         temp = z__[j4p2 + 2] / z__[j4 - 2];
 179         z__[j4] = z__[j4p2] * temp;
 180         *dn = *dnm1 * temp;
 181     } else {
 182         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
 183         *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]);
 184     }
 185     *dmin__ = min(*dmin__,*dn);
 186
 187     z__[j4 + 2] = *dn;
 188     z__[(*n0 << 2) - *pp] = emin;
 189     return 0;
 190
 191 /*     End of DLASQ6 */
 192
 193 } /* dlasq6_ */