src/gmxlib/gmx_lapack/dlasrt2.c

   1 #include "gmx_lapack.h"
   2
   3 void F77_FUNC(dlasrt2,DLASRT2)(const char *id,
   4               int *n,
   5               double *d__,
   6               int * key,
   7               int *info)
   8 {
   9     int i__1, i__2;
  10
  11     int i__, j;
  12     double d1, d2, d3;
  13     int dir;
  14     double tmp;
  15     int endd;
  16     int stack[64];
  17     double dmnmx;
  18     int start;
  19     int tmpkey, stkpnt;
  20
  21     --key;
  22     --d__;
  23
  24     *info = 0;
  25     dir = -1;
  26     if (*id=='D' || *id=='d')
  27         dir = 0;
  28     else if (*id=='I' || *id=='i')
  29         dir = 1;
  30
  31     if (dir == -1) {
  32         *info = -1;
  33     } else if (*n < 0) {
  34         *info = -2;
  35     }
  36     if (*info != 0) {
  37         i__1 = -(*info);
  38         return;
  39     }
  40
  41     if (*n <= 1) {
  42         return;
  43     }
  44
  45     stkpnt = 1;
  46     stack[0] = 1;
  47     stack[1] = *n;
  48 L10:
  49     start = stack[(stkpnt << 1) - 2];
  50     endd = stack[(stkpnt << 1) - 1];
  51     --stkpnt;
  52     if (endd - start > 0) {
  53
  54         if (dir == 0) {
  55
  56             i__1 = endd;
  57             for (i__ = start + 1; i__ <= i__1; ++i__) {
  58                 i__2 = start + 1;
  59                 for (j = i__; j >= i__2; --j) {
  60                     if (d__[j] > d__[j - 1]) {
  61                         dmnmx = d__[j];
  62                         d__[j] = d__[j - 1];
  63                         d__[j - 1] = dmnmx;
  64                         tmpkey = key[j];
  65                         key[j] = key[j - 1];
  66                         key[j - 1] = tmpkey;
  67                     } else {
  68                         break;
  69                     }
  70                 }
  71             }
  72
  73         } else {
  74
  75             i__1 = endd;
  76             for (i__ = start + 1; i__ <= i__1; ++i__) {
  77                 i__2 = start + 1;
  78                 for (j = i__; j >= i__2; --j) {
  79                     if (d__[j] < d__[j - 1]) {
  80                         dmnmx = d__[j];
  81                         d__[j] = d__[j - 1];
  82                         d__[j - 1] = dmnmx;
  83                         tmpkey = key[j];
  84                         key[j] = key[j - 1];
  85                         key[j - 1] = tmpkey;
  86                     } else {
  87                         break;
  88                     }
  89                 }
  90             }
  91
  92         }
  93
  94     } else if (endd - start > 20) {
  95
  96         d1 = d__[start];
  97         d2 = d__[endd];
  98         i__ = (start + endd) / 2;
  99         d3 = d__[i__];
 100         if (d1 < d2) {
 101             if (d3 < d1) {
 102                 dmnmx = d1;
 103             } else if (d3 < d2) {
 104                 dmnmx = d3;
 105             } else {
 106                 dmnmx = d2;
 107             }
 108         } else {
 109             if (d3 < d2) {
 110                 dmnmx = d2;
 111             } else if (d3 < d1) {
 112                 dmnmx = d3;
 113             } else {
 114                 dmnmx = d1;
 115             }
 116         }
 117
 118         if (dir == 0) {
 119
 120             i__ = start - 1;
 121             j = endd + 1;
 122 L60:
 123 L70:
 124             --j;
 125             if (d__[j] < dmnmx) {
 126                 goto L70;
 127             }
 128 L80:
 129             ++i__;
 130             if (d__[i__] > dmnmx) {
 131                 goto L80;
 132             }
 133             if (i__ < j) {
 134                 tmp = d__[i__];
 135                 d__[i__] = d__[j];
 136                 d__[j] = tmp;
 137                 tmpkey = key[j];
 138                 key[j] = key[i__];
 139                 key[i__] = tmpkey;
 140                 goto L60;
 141             }
 142             if (j - start > endd - j - 1) {
 143                 ++stkpnt;
 144                 stack[(stkpnt << 1) - 2] = start;
 145                 stack[(stkpnt << 1) - 1] = j;
 146                 ++stkpnt;
 147                 stack[(stkpnt << 1) - 2] = j + 1;
 148                 stack[(stkpnt << 1) - 1] = endd;
 149             } else {
 150                 ++stkpnt;
 151                 stack[(stkpnt << 1) - 2] = j + 1;
 152                 stack[(stkpnt << 1) - 1] = endd;
 153                 ++stkpnt;
 154                 stack[(stkpnt << 1) - 2] = start;
 155                 stack[(stkpnt << 1) - 1] = j;
 156             }
 157         } else {
 158
 159             i__ = start - 1;
 160             j = endd + 1;
 161 L90:
 162 L100:
 163             --j;
 164             if (d__[j] > dmnmx) {
 165                 goto L100;
 166             }
 167 L110:
 168             ++i__;
 169             if (d__[i__] < dmnmx) {
 170                 goto L110;
 171             }
 172             if (i__ < j) {
 173                 tmp = d__[i__];
 174                 d__[i__] = d__[j];
 175                 d__[j] = tmp;
 176                 tmpkey = key[j];
 177                 key[j] = key[i__];
 178                 key[i__] = tmpkey;
 179                 goto L90;
 180             }
 181             if (j - start > endd - j - 1) {
 182                 ++stkpnt;
 183                 stack[(stkpnt << 1) - 2] = start;
 184                 stack[(stkpnt << 1) - 1] = j;
 185                 ++stkpnt;
 186                 stack[(stkpnt << 1) - 2] = j + 1;
 187                 stack[(stkpnt << 1) - 1] = endd;
 188             } else {
 189                 ++stkpnt;
 190                 stack[(stkpnt << 1) - 2] = j + 1;
 191                 stack[(stkpnt << 1) - 1] = endd;
 192                 ++stkpnt;
 193                 stack[(stkpnt << 1) - 2] = start;
 194                 stack[(stkpnt << 1) - 1] = j;
 195             }
 196         }
 197     }
 198     if (stkpnt > 0) {
 199         goto L10;
 200     }
 201
 202     return;
 203 }