src/gmxlib/gmx_sort.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * This file is part of Gromacs        Copyright (c) 1991-2010
   5  * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
   6  * Copyright (c) 2012, by the GROMACS development team, led by
   7  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   8  * others, as listed in the AUTHORS file in the top-level source
   9  * directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdlib.h>
  43 #include "gmx_sort.h"
  44
  45
  46 static void
  47 qsort_swapfunc(void *        a,
  48                void *        b,
  49                size_t        n,
  50                int           swaptype)
  51 {
  52     int *  ia;
  53     int *  ib;
  54     int    itmp;
  55
  56     char *  ca;
  57     char *  cb;
  58     char    ctmp;
  59
  60     if (swaptype <= 1)
  61     {
  62         ia = (int *)a;
  63         ib = (int *)b;
  64         for( ; n > 0; ia += 1, ib += 1, n -= sizeof(int))
  65         {
  66             itmp      = *ia;
  67             *ia       = *ib;
  68             *ib       = itmp;
  69         }
  70     }
  71     else
  72     {
  73         ca = (char *)a;
  74         cb = (char *)b;
  75         for( ; n > 0; ca += 1, cb += 1, n -= 1)
  76         {
  77             ctmp       = *ca;
  78             *ca        = *cb;
  79             *cb        = ctmp;
  80         }
  81     }
  82 }
  83
  84
  85 static void *
  86 qsort_med3(void *        a,
  87            void *        b,
  88            void *        c,
  89            int          (*compar) (const void *a, const void *b))
  90 {
  91     if(compar(a,b) < 0)
  92     {
  93         if(compar(b,c) < 0)
  94             return b;
  95         else if(compar(a,c) < 0)
  96             return c;
  97         else
  98             return a;
  99     }
 100     else
 101     {
 102         if(compar(b,c) > 0)
 103             return b;
 104         else if(compar(a,c) > 0)
 105             return c;
 106         else
 107             return a;
 108     }
 109 }
 110
 111
 112 void
 113 gmx_qsort(void *           base,
 114           size_t           nmemb,
 115           size_t           size,
 116           int            (*compar)(const void *, const void *))
 117 {
 118 #define QSORT_EXCH(a, b, t) (t = a, a = b, b = t)
 119 #define QSORT_SWAP(a, b) swaptype != 0 ? qsort_swapfunc(a, b, size, swaptype) : \
 120                             (void)QSORT_EXCH(*(int *)(a), *(int *)(b), t)
 121
 122     char *pa, *pb, *pc, *pd, *pl, *pm, *pn, *pv, *cbase;
 123     int r, swaptype;
 124     int t, v;
 125     size_t s, st;
 126
 127     cbase = (char *)base;
 128
 129     swaptype = (size_t)(cbase - (char *)0) % sizeof(int) || size % sizeof(int) ? 2 : size == sizeof(int)? 0 : 1;
 130
 131     if (nmemb < 7)
 132     {
 133         /* Insertion sort on smallest arrays */
 134         for (pm = cbase + size; pm < cbase + nmemb*size; pm += size)
 135         {
 136             for (pl = pm; (pl > cbase) && compar((void *)(pl-size),(void *) pl) > 0; pl -= size)
 137             {
 138                 QSORT_SWAP(pl, pl-size);
 139             }
 140         }
 141         return;
 142     }
 143
 144     /* Small arrays, middle element */
 145     pm = cbase + (nmemb/2)*size;
 146
 147     if (nmemb > 7)
 148     {
 149         pl = cbase;
 150         pn = cbase + (nmemb-1)*size;
 151         if (nmemb > 40)
 152         {
 153             /* Big arrays, pseudomedian of 9 */
 154             s = (nmemb/8)*size;
 155             pl = (char *)qsort_med3((void *)pl, (void *)((size_t)pl+s), (void *)((size_t)pl+2*s), compar);
 156             pm = (char *)qsort_med3((void *)((size_t)pm-s), (void *)pm, (void *)((size_t)pm+s), compar);
 157             pn = (char *)qsort_med3((void *)((size_t)pn-2*s), (void *)((size_t)pn-s), (void *)pn, compar);
 158         }
 159         /* Mid-size, med of 3 */
 160         pm = (char *)qsort_med3((void *)pl, (void *)pm, (void *)pn, compar);
 161     }
 162
 163     /* pv points to partition value */
 164     if (swaptype != 0)
 165     {
 166         pv = cbase;
 167         QSORT_SWAP(pv, pm);
 168     }
 169     else
 170     {
 171         pv = (char*)(void*)&v;
 172         v = *(int *)pm;
 173     }
 174
 175     pa = pb = cbase;
 176     pc = pd = cbase + (nmemb-1)*size;
 177
 178     for (;;)
 179     {
 180         while (pb <= pc && (r = compar((void *)pb,(void *) pv)) <= 0)
 181         {
 182             if (r == 0)
 183             {
 184                 QSORT_SWAP(pa, pb);
 185                 pa += size;
 186             }
 187             pb += size;
 188         }
 189         while (pc >= pb && (r = compar((void *)pc,(void *) pv)) >= 0)
 190         {
 191             if (r == 0)
 192             {
 193                 QSORT_SWAP(pc, pd);
 194                 pd -= size;
 195             }
 196             pc -= size;
 197         }
 198         if (pb > pc)
 199         {
 200             break;
 201         }
 202         QSORT_SWAP(pb, pc);
 203         pb += size;
 204         pc -= size;
 205     }
 206     pn = cbase + nmemb*size;
 207
 208     s  = pa-cbase;
 209     st = pb-pa;
 210     if(st<s)
 211     {
 212         s = st;
 213     }
 214
 215     if(s>0)
 216     {
 217         qsort_swapfunc(cbase, pb-s, s, swaptype);
 218     }
 219
 220     s = pd-pc;
 221     st = pn-pd-size;
 222     if(st<s)
 223     {
 224         s = st;
 225     }
 226
 227     if(s>0)
 228     {
 229         qsort_swapfunc(pb, pn-s, s, swaptype);
 230     }
 231
 232     if ((s = pb-pa) > size)
 233     {
 234         gmx_qsort(cbase, s/size, size, compar);
 235     }
 236
 237     if ((s = pd-pc) > size)
 238     {
 239         gmx_qsort(pn-s, s/size, size, compar);
 240     }
 241
 242 #undef QSORT_EXCH
 243 #undef QSORT_SWAP
 244
 245     return;
 246 }
 247
 248
 249