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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlasrt_(char *id, integer *n, doublereal *d__, integer *
   4         info)
   5 {
   6     /* System generated locals */
   7     integer i__1, i__2;
   8
   9     /* Local variables */
  10     integer i__, j;
  11     doublereal d1, d2, d3;
  12     integer dir;
  13     doublereal tmp;
  14     integer endd;
  15     extern logical lsame_(char *, char *);
  16     integer stack[64]   /* was [2][32] */;
  17     doublereal dmnmx;
  18     integer start;
  19     extern /* Subroutine */ int xerbla_(char *, integer *);
  20     integer stkpnt;
  21
  22
  23 /*  -- LAPACK routine (version 3.1) -- */
  24 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  25 /*     November 2006 */
  26
  27 /*     .. Scalar Arguments .. */
  28 /*     .. */
  29 /*     .. Array Arguments .. */
  30 /*     .. */
  31
  32 /*  Purpose */
  33 /*  ======= */
  34
  35 /*  Sort the numbers in D in increasing order (if ID = 'I') or */
  36 /*  in decreasing order (if ID = 'D' ). */
  37
  38 /*  Use Quick Sort, reverting to Insertion sort on arrays of */
  39 /*  size <= 20. Dimension of STACK limits N to about 2**32. */
  40
  41 /*  Arguments */
  42 /*  ========= */
  43
  44 /*  ID      (input) CHARACTER*1 */
  45 /*          = 'I': sort D in increasing order; */
  46 /*          = 'D': sort D in decreasing order. */
  47
  48 /*  N       (input) INTEGER */
  49 /*          The length of the array D. */
  50
  51 /*  D       (input/output) DOUBLE PRECISION array, dimension (N) */
  52 /*          On entry, the array to be sorted. */
  53 /*          On exit, D has been sorted into increasing order */
  54 /*          (D(1) <= ... <= D(N) ) or into decreasing order */
  55 /*          (D(1) >= ... >= D(N) ), depending on ID. */
  56
  57 /*  INFO    (output) INTEGER */
  58 /*          = 0:  successful exit */
  59 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
  60
  61 /*  ===================================================================== */
  62
  63 /*     .. Parameters .. */
  64 /*     .. */
  65 /*     .. Local Scalars .. */
  66 /*     .. */
  67 /*     .. Local Arrays .. */
  68 /*     .. */
  69 /*     .. External Functions .. */
  70 /*     .. */
  71 /*     .. External Subroutines .. */
  72 /*     .. */
  73 /*     .. Executable Statements .. */
  74
  75 /*     Test the input paramters. */
  76
  77     /* Parameter adjustments */
  78     --d__;
  79
  80     /* Function Body */
  81     *info = 0;
  82     dir = -1;
  83     if (lsame_(id, "D")) {
  84         dir = 0;
  85     } else if (lsame_(id, "I")) {
  86         dir = 1;
  87     }
  88     if (dir == -1) {
  89         *info = -1;
  90     } else if (*n < 0) {
  91         *info = -2;
  92     }
  93     if (*info != 0) {
  94         i__1 = -(*info);
  95         xerbla_("DLASRT", &i__1);
  96         return 0;
  97     }
  98
  99 /*     Quick return if possible */
 100
 101     if (*n <= 1) {
 102         return 0;
 103     }
 104
 105     stkpnt = 1;
 106     stack[0] = 1;
 107     stack[1] = *n;
 108 L10:
 109     start = stack[(stkpnt << 1) - 2];
 110     endd = stack[(stkpnt << 1) - 1];
 111     --stkpnt;
 112     if (endd - start <= 20 && endd - start > 0) {
 113
 114 /*        Do Insertion sort on D( START:ENDD ) */
 115
 116         if (dir == 0) {
 117
 118 /*           Sort into decreasing order */
 119
 120             i__1 = endd;
 121             for (i__ = start + 1; i__ <= i__1; ++i__) {
 122                 i__2 = start + 1;
 123                 for (j = i__; j >= i__2; --j) {
 124                     if (d__[j] > d__[j - 1]) {
 125                         dmnmx = d__[j];
 126                         d__[j] = d__[j - 1];
 127                         d__[j - 1] = dmnmx;
 128                     } else {
 129                         goto L30;
 130                     }
 131 /* L20: */
 132                 }
 133 L30:
 134                 ;
 135             }
 136
 137         } else {
 138
 139 /*           Sort into increasing order */
 140
 141             i__1 = endd;
 142             for (i__ = start + 1; i__ <= i__1; ++i__) {
 143                 i__2 = start + 1;
 144                 for (j = i__; j >= i__2; --j) {
 145                     if (d__[j] < d__[j - 1]) {
 146                         dmnmx = d__[j];
 147                         d__[j] = d__[j - 1];
 148                         d__[j - 1] = dmnmx;
 149                     } else {
 150                         goto L50;
 151                     }
 152 /* L40: */
 153                 }
 154 L50:
 155                 ;
 156             }
 157
 158         }
 159
 160     } else if (endd - start > 20) {
 161
 162 /*        Partition D( START:ENDD ) and stack parts, largest one first */
 163
 164 /*        Choose partition entry as median of 3 */
 165
 166         d1 = d__[start];
 167         d2 = d__[endd];
 168         i__ = (start + endd) / 2;
 169         d3 = d__[i__];
 170         if (d1 < d2) {
 171             if (d3 < d1) {
 172                 dmnmx = d1;
 173             } else if (d3 < d2) {
 174                 dmnmx = d3;
 175             } else {
 176                 dmnmx = d2;
 177             }
 178         } else {
 179             if (d3 < d2) {
 180                 dmnmx = d2;
 181             } else if (d3 < d1) {
 182                 dmnmx = d3;
 183             } else {
 184                 dmnmx = d1;
 185             }
 186         }
 187
 188         if (dir == 0) {
 189
 190 /*           Sort into decreasing order */
 191
 192             i__ = start - 1;
 193             j = endd + 1;
 194 L60:
 195 L70:
 196             --j;
 197             if (d__[j] < dmnmx) {
 198                 goto L70;
 199             }
 200 L80:
 201             ++i__;
 202             if (d__[i__] > dmnmx) {
 203                 goto L80;
 204             }
 205             if (i__ < j) {
 206                 tmp = d__[i__];
 207                 d__[i__] = d__[j];
 208                 d__[j] = tmp;
 209                 goto L60;
 210             }
 211             if (j - start > endd - j - 1) {
 212                 ++stkpnt;
 213                 stack[(stkpnt << 1) - 2] = start;
 214                 stack[(stkpnt << 1) - 1] = j;
 215                 ++stkpnt;
 216                 stack[(stkpnt << 1) - 2] = j + 1;
 217                 stack[(stkpnt << 1) - 1] = endd;
 218             } else {
 219                 ++stkpnt;
 220                 stack[(stkpnt << 1) - 2] = j + 1;
 221                 stack[(stkpnt << 1) - 1] = endd;
 222                 ++stkpnt;
 223                 stack[(stkpnt << 1) - 2] = start;
 224                 stack[(stkpnt << 1) - 1] = j;
 225             }
 226         } else {
 227
 228 /*           Sort into increasing order */
 229
 230             i__ = start - 1;
 231             j = endd + 1;
 232 L90:
 233 L100:
 234             --j;
 235             if (d__[j] > dmnmx) {
 236                 goto L100;
 237             }
 238 L110:
 239             ++i__;
 240             if (d__[i__] < dmnmx) {
 241                 goto L110;
 242             }
 243             if (i__ < j) {
 244                 tmp = d__[i__];
 245                 d__[i__] = d__[j];
 246                 d__[j] = tmp;
 247                 goto L90;
 248             }
 249             if (j - start > endd - j - 1) {
 250                 ++stkpnt;
 251                 stack[(stkpnt << 1) - 2] = start;
 252                 stack[(stkpnt << 1) - 1] = j;
 253                 ++stkpnt;
 254                 stack[(stkpnt << 1) - 2] = j + 1;
 255                 stack[(stkpnt << 1) - 1] = endd;
 256             } else {
 257                 ++stkpnt;
 258                 stack[(stkpnt << 1) - 2] = j + 1;
 259                 stack[(stkpnt << 1) - 1] = endd;
 260                 ++stkpnt;
 261                 stack[(stkpnt << 1) - 2] = start;
 262                 stack[(stkpnt << 1) - 1] = j;
 263             }
 264         }
 265     }
 266     if (stkpnt > 0) {
 267         goto L10;
 268     }
 269     return 0;
 270
 271 /*     End of DLASRT */
 272
 273 } /* dlasrt_ */