src/kernel/sorting.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <limits.h>
  40 #include "sysstuff.h"
  41 #include "smalloc.h"
  42 #include "sorting.h"
  43
  44 /*****************************************************************************
  45  *                                                                           *
  46  *                   Block sorting on coordinates                            *
  47  *                                                                           *
  48  *****************************************************************************/
  49
  50 static rvec *make_xblock(t_block *block,rvec x[])
  51 {
  52   int i,j,k,nr,n;
  53   rvec *xblock;
  54
  55   nr=block->nr;
  56   snew(xblock,nr);
  57   for (i=0; i<nr; i++)
  58     {
  59       for (j=0; j<DIM; j++) xblock[i][j]=0.0;
  60       for (j=block->index[i]; j<(int)(block->index[i+1]); j++)
  61         for (k=0; k<DIM; k++) xblock[i][k]+=x[j][k];
  62       n=block->index[i+1]-block->index[i];
  63       for (k=0; k<DIM; k++) xblock[i][k]/=n;
  64     }
  65   return xblock;
  66 }
  67
  68 static rvec *xblock; /* just global to bcomp1, used in qsort */
  69
  70 static int bomp1(const void *p1,const void *p2)
  71 {
  72   int i,i1,i2;
  73
  74   i1=*(int *)p1;
  75   i2=*(int *)p2;
  76   for (i=0; i<DIM; i++)
  77     if (xblock[i1][i]<xblock[i2][i]) return -1;
  78     else if (xblock[i1][i]>xblock[i2][i]) return 1;
  79   return 0;
  80 }
  81
  82 void sort_xblock(t_block *block,rvec x[],int renum[])
  83 {
  84   int i,nr,*invnum;
  85
  86   nr=block->nr;
  87   snew(invnum,nr);
  88   xblock=make_xblock(block,x);
  89   for (i=0; i<nr; i++) invnum[i]=i;
  90   qsort((void *)invnum,nr,(size_t)sizeof(invnum[0]),bomp1);
  91   for (i=0; i<nr; i++) renum[invnum[i]]=i;
  92   sfree(xblock);
  93   sfree(invnum);
  94 }
  95
  96 /*****************************************************************************
  97  *                                                                           *
  98  *                        Bond list sorting                                  *
  99  *                                                                           *
 100  *****************************************************************************/
 101
 102 static int bcomp2(const void *p1,const void *p2)
 103 {
 104   int done;
 105
 106   if ((((atom_id *)p1)[0])!=(((atom_id *)p2)[0]))
 107     done=((((atom_id *)p1)[0])-(((atom_id *)p2)[0]));
 108   else
 109     done=((((atom_id *)p1)[1])-(((atom_id *)p2)[1]));
 110 #ifdef DEBUG
 111   printf("bcomp2: [%d,%d] with [%d,%d] result %d\n",
 112           ((atom_id *)p1)[0],((atom_id *)p1)[1],
 113           ((atom_id *)p2)[0],((atom_id *)p2)[1],done);
 114 #endif
 115   return done;
 116 }
 117
 118 void sort_bond_list(t_bond bonds[],int nr)
 119 {
 120   qsort((void *)bonds,nr,(size_t)sizeof(bonds[0]),bcomp2);
 121 }