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

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