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

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