src/gromacs/gmxlib/splitter.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <assert.h>
  41 #include <stdio.h>
  42 #include <string.h>
  43
  44 #include "sysstuff.h"
  45 #include "macros.h"
  46 #include "smalloc.h"
  47 #include "typedefs.h"
  48 #include "mshift.h"
  49 #include "invblock.h"
  50 #include "txtdump.h"
  51 #include <math.h>
  52 #include "gmx_fatal.h"
  53 #include "splitter.h"
  54
  55 typedef struct {
  56   int     nr;
  57   t_iatom *ia;
  58 } t_sf;
  59
  60 static t_sf *init_sf(int nr)
  61 {
  62   t_sf *sf;
  63   int  i;
  64
  65   snew(sf,nr);
  66   for(i=0; (i<nr); i++) {
  67     sf[i].nr=0;
  68     sf[i].ia=NULL;
  69   }
  70
  71   return sf;
  72 }
  73
  74 static void done_sf(int nr,t_sf *sf)
  75 {
  76   int i;
  77
  78   for(i=0; (i<nr); i++) {
  79     sf[i].nr=0;
  80     sfree(sf[i].ia);
  81     sf[i].ia=NULL;
  82   }
  83   sfree(sf);
  84 }
  85
  86 static void push_sf(t_sf *sf,int nr,t_iatom ia[])
  87 {
  88   int i;
  89
  90   srenew(sf->ia,sf->nr+nr);
  91   for(i=0; (i<nr); i++)
  92     sf->ia[sf->nr+i]=ia[i];
  93   sf->nr+=nr;
  94 }
  95
  96 static int min_nodeid(int nr,atom_id list[],int hid[])
  97 {
  98   int i,nodeid,minnodeid;
  99
 100   if (nr <= 0)
 101     gmx_incons("Invalid node number");
 102   minnodeid=hid[list[0]];
 103   for (i=1; (i<nr); i++)
 104     if ((nodeid=hid[list[i]]) < minnodeid)
 105       minnodeid=nodeid;
 106
 107   return minnodeid;
 108 }
 109
 110
 111 static void split_force2(t_inputrec *ir, int nnodes,int hid[],int ftype,t_ilist *ilist,
 112                                                  int *multinr,
 113                                                  int *constr_min_nodeid,int * constr_max_nodeid,
 114                                                  int *left_range, int *right_range)
 115 {
 116   int     i,j,k,type,nodeid,nratoms,tnr;
 117   int     nvsite_constr;
 118   t_iatom ai,aj;
 119   int     node_low_ai,node_low_aj,node_high_ai,node_high_aj;
 120   int     node_low,node_high;
 121   int     nodei,nodej;
 122   t_sf    *sf;
 123   int     nextra;
 124
 125   sf = init_sf(nnodes);
 126
 127   node_high = node_low = 0;
 128   nextra    = 0;
 129
 130   /* Walk along all the bonded forces, find the appropriate node
 131    * to calc it on, and add it to that nodes list.
 132    */
 133   for (i=0; i<ilist->nr; i+=(1+nratoms))
 134   {
 135     type    = ilist->iatoms[i];
 136     nratoms = interaction_function[ftype].nratoms;
 137
 138     if (ftype == F_CONSTR)
 139         {
 140                 ai = ilist->iatoms[i+1];
 141                 aj = ilist->iatoms[i+2];
 142
 143                 nodei = hid[ai];
 144                 nodej = hid[aj];
 145                 nodeid = nodei;
 146
 147                 if(ir->eConstrAlg == econtLINCS)
 148                 {
 149                         node_low_ai   = constr_min_nodeid[ai];
 150                         node_low_aj   = constr_min_nodeid[aj];
 151                         node_high_ai  = constr_max_nodeid[ai];
 152                         node_high_aj  = constr_max_nodeid[aj];
 153
 154                         node_low      = min(node_low_ai,node_low_aj);
 155                         node_high     = max(node_high_ai,node_high_aj);
 156
 157             if (node_high-nodei > 1 || nodei-node_low > 1 ||
 158                 node_high-nodej > 1 || nodej-node_low > 1 )
 159                         {
 160                                 gmx_fatal(FARGS,"Constraint dependencies further away than next-neighbor\n"
 161                                                   "in particle decomposition. Constraint between atoms %d--%d evaluated\n"
 162                                                   "on node %d and %d, but atom %d has connections within %d bonds (lincs_order)\n"
 163                                                   "of node %d, and atom %d has connections within %d bonds of node %d.\n"
 164                                                   "Reduce the # nodes, lincs_order, or\n"
 165                                                   "try domain decomposition.",ai,aj,nodei,nodej,ai,ir->nProjOrder,node_low,aj,ir->nProjOrder,node_high);
 166                         }
 167
 168                         if (node_low < nodei || node_low < nodej)
 169                         {
 170                                 right_range[node_low] = max(right_range[node_low],aj);
 171                         }
 172                         if (node_high > nodei || node_high > nodej)
 173                         {
 174                                 left_range[node_high] = min(left_range[node_high],ai);
 175                         }
 176                 }
 177                 else
 178                 {
 179                         /* Shake */
 180                         if (hid[ilist->iatoms[i+2]] != nodei)
 181                                 gmx_fatal(FARGS,"Shake block crossing node boundaries\n"
 182                                                   "constraint between atoms (%d,%d) (try LINCS instead!)",
 183                                                   ilist->iatoms[i+1]+1,ilist->iatoms[i+2]+1);
 184                 }
 185     }
 186     else if (ftype == F_SETTLE)
 187         {
 188       /* Only the first particle is stored for settles ... */
 189       ai=ilist->iatoms[i+1];
 190       nodeid=hid[ai];
 191       if (nodeid != hid[ilist->iatoms[i+2]] ||
 192           nodeid != hid[ilist->iatoms[i+3]])
 193         gmx_fatal(FARGS,"Settle block crossing node boundaries\n"
 194                   "constraint between atoms %d, %d, %d)",
 195               ai,ilist->iatoms[i+2],ilist->iatoms[i+3]);
 196     }
 197     else if(interaction_function[ftype].flags & IF_VSITE)
 198         {
 199       /* Virtual sites are special, since we need to pre-communicate
 200            * their coordinates to construct vsites before then main
 201            * coordinate communication.
 202            * Vsites can have constructing atoms both larger and smaller than themselves.
 203            * To minimize communication and book-keeping, each vsite is constructed on
 204            * the home node of the atomnr of the vsite.
 205            * Since the vsite coordinates too have to be communicated to the next node,
 206            * we need to
 207            *
 208            * 1. Pre-communicate coordinates of constructing atoms
 209            * 2. Construct the vsite
 210            * 3. Perform main coordinate communication
 211            *
 212            * Note that this has change from gromacs 4.0 and earlier, where the vsite
 213            * was constructed on the home node of the lowest index of any of the constructing
 214            * atoms and the vsite itself.
 215            */
 216
 217                 if (ftype==F_VSITE2)
 218                         nvsite_constr=2;
 219                 else if(ftype==F_VSITE4FD || ftype==F_VSITE4FDN)
 220                         nvsite_constr=4;
 221                 else
 222                         nvsite_constr=3;
 223
 224                 /* Vsites are constructed on the home node of the actual site to save communication
 225                  * and simplify the book-keeping.
 226                  */
 227                 nodeid=hid[ilist->iatoms[i+1]];
 228
 229                 for(k=2;k<nvsite_constr+2;k++)
 230                 {
 231                         if(hid[ilist->iatoms[i+k]]<(nodeid-1) ||
 232                            hid[ilist->iatoms[i+k]]>(nodeid+1))
 233                                 gmx_fatal(FARGS,"Virtual site %d and its constructing"
 234                                                   " atoms are not on the same or adjacent\n"
 235                                                   " nodes. This is necessary to avoid a lot\n"
 236                                                   " of extra communication. The easiest way"
 237                                                   " to ensure this is to place virtual sites\n"
 238                                                   " close to the constructing atoms.\n"
 239                                                   " Sorry, but you will have to rework your topology!\n",
 240                                                   ilist->iatoms[i+1]);
 241                 }
 242     }
 243         else
 244     {
 245                 nodeid=min_nodeid(nratoms,&ilist->iatoms[i+1],hid);
 246         }
 247
 248         if (ftype == F_CONSTR && ir->eConstrAlg == econtLINCS)
 249         {
 250                 push_sf(&(sf[nodeid]),nratoms+1,&(ilist->iatoms[i]));
 251
 252                 if(node_low<nodeid)
 253                 {
 254                         push_sf(&(sf[node_low]),nratoms+1,&(ilist->iatoms[i]));
 255                         nextra+=nratoms+1;
 256                 }
 257                 if(node_high>nodeid)
 258                 {
 259                         push_sf(&(sf[node_high]),nratoms+1,&(ilist->iatoms[i]));
 260                         nextra+=nratoms+1;
 261                 }
 262         }
 263         else
 264         {
 265                 push_sf(&(sf[nodeid]),nratoms+1,&(ilist->iatoms[i]));
 266         }
 267   }
 268
 269   if(nextra>0)
 270   {
 271           ilist->nr += nextra;
 272           srenew(ilist->iatoms,ilist->nr);
 273   }
 274
 275   tnr=0;
 276   for(nodeid=0; (nodeid<nnodes); nodeid++) {
 277     for (i=0; (i<sf[nodeid].nr); i++)
 278       ilist->iatoms[tnr++]=sf[nodeid].ia[i];
 279
 280     multinr[nodeid]=(nodeid==0) ? 0 : multinr[nodeid-1];
 281     multinr[nodeid]+=sf[nodeid].nr;
 282   }
 283
 284         if (tnr != ilist->nr)
 285                 gmx_incons("Splitting forces over processors");
 286
 287   done_sf(nnodes,sf);
 288 }
 289
 290 static int *home_index(int nnodes,t_block *cgs,int *multinr)
 291 {
 292   /* This routine determines the node id for each particle */
 293   int *hid;
 294   int nodeid,j0,j1,j,k;
 295
 296   snew(hid,cgs->index[cgs->nr]);
 297   /* Initiate to -1 to make it possible to check afterwards,
 298    * all hid's should be set in the loop below
 299    */
 300   for(k=0; (k<cgs->index[cgs->nr]); k++)
 301     hid[k]=-1;
 302
 303   /* loop over nodes */
 304   for(nodeid=0; (nodeid<nnodes); nodeid++) {
 305     j0 = (nodeid==0) ? 0 : multinr[nodeid-1];
 306     j1 = multinr[nodeid];
 307
 308     /* j0 and j1 are the boundariesin the index array */
 309     for(j=j0; (j<j1); j++) {
 310       for(k=cgs->index[j]; (k<cgs->index[j+1]); k++) {
 311         hid[k]=nodeid;
 312       }
 313     }
 314   }
 315   /* Now verify that all hid's are not -1 */
 316   for(k=0; (k<cgs->index[cgs->nr]); k++)
 317     if (hid[k] == -1)
 318       gmx_fatal(FARGS,"hid[%d] = -1, cgs->nr = %d, natoms = %d",
 319                   k,cgs->nr,cgs->index[cgs->nr]);
 320
 321   return hid;
 322 }
 323
 324 typedef struct {
 325   int atom,ic,is;
 326 } t_border;
 327
 328 void set_bor(t_border *b,int atom,int ic,int is)
 329 {
 330   if (debug)
 331     fprintf(debug,"border @ atom %5d [ ic = %5d,  is = %5d ]\n",atom,ic,is);
 332   b->atom = atom;
 333   b->ic   = ic;
 334   b->is   = is;
 335 }
 336
 337 static gmx_bool is_bor(atom_id ai[],int i)
 338 {
 339   return ((ai[i] != ai[i-1]) || ((ai[i] == NO_ATID) && (ai[i-1] == NO_ATID)));
 340 }
 341
 342 static t_border *mk_border(FILE *fp,int natom,atom_id *invcgs,
 343                            atom_id *invshk,int *nb)
 344 {
 345   t_border *bor;
 346   atom_id  *sbor,*cbor;
 347   int      i,j,is,ic,ns,nc,nbor;
 348
 349   if (debug) {
 350     for(i=0; (i<natom); i++) {
 351       fprintf(debug,"atom: %6d  cgindex: %6d  shkindex: %6d\n",
 352                           i, invcgs[i], invshk[i]);
 353     }
 354   }
 355
 356   snew(sbor,natom+1);
 357   snew(cbor,natom+1);
 358   ns = nc = 1;
 359   for(i=1; (i<natom); i++) {
 360     if (is_bor(invcgs,i))
 361       cbor[nc++] = i;
 362     if (is_bor(invshk,i))
 363       sbor[ns++] = i;
 364   }
 365   sbor[ns] = 0;
 366   cbor[nc] = 0;
 367   if (fp)
 368   {
 369           fprintf(fp,"There are %d charge group borders",nc);
 370           if(invshk!=NULL)
 371           {
 372                   fprintf(fp," and %d shake borders",ns);
 373           }
 374           fprintf(fp,".\n");
 375   }
 376   snew(bor,max(nc,ns));
 377   ic = is = nbor = 0;
 378   while ((ic < nc) || (is < ns)) {
 379     if (sbor[is] == cbor[ic]) {
 380       set_bor(&(bor[nbor]),cbor[ic],ic,is);
 381       nbor++;
 382       if (ic < nc) ic++;
 383       if (is < ns) is++;
 384     }
 385     else if (cbor[ic] > sbor[is]) {
 386       if (is == ns) {
 387         set_bor(&(bor[nbor]),cbor[ic],ic,is);
 388         nbor++;
 389         if (ic < nc) ic++;
 390       }
 391       else if (is < ns)
 392         is++;
 393     }
 394     else if (ic < nc)
 395       ic++;
 396     else
 397       is++;/*gmx_fatal(FARGS,"Can't happen is=%d, ic=%d (%s, %d)",
 398              is,ic,__FILE__,__LINE__);*/
 399   }
 400   if (fp)
 401     fprintf(fp,"There are %d total borders\n",nbor);
 402
 403   if (debug) {
 404     fprintf(debug,"There are %d actual bor entries\n",nbor);
 405     for(i=0; (i<nbor); i++)
 406       fprintf(debug,"bor[%5d] = atom: %d  ic: %d  is: %d\n",i,
 407               bor[i].atom,bor[i].ic,bor[i].is);
 408   }
 409
 410   *nb = nbor;
 411
 412   return bor;
 413 }
 414
 415 static void split_blocks(FILE *fp,t_inputrec *ir, int nnodes,
 416                          t_block *cgs,t_blocka *sblock,real capacity[],
 417                          int *multinr_cgs)
 418 {
 419   int      natoms,*maxatom;
 420   int      i,ii,ai,b0,b1;
 421   int      nodeid,last_shk,nbor;
 422   t_border *border;
 423   double   tload,tcap;
 424
 425   gmx_bool    bSHK;
 426   atom_id *shknum,*cgsnum;
 427
 428   natoms = cgs->index[cgs->nr];
 429
 430   if (NULL != debug) {
 431     pr_block(debug,0,"cgs",cgs,TRUE);
 432     pr_blocka(debug,0,"sblock",sblock,TRUE);
 433     fflush(debug);
 434   }
 435
 436   cgsnum = make_invblock(cgs,natoms+1);
 437   shknum = make_invblocka(sblock,natoms+1);
 438   border = mk_border(fp,natoms,cgsnum,shknum,&nbor);
 439
 440   snew(maxatom,nnodes);
 441   tload  = capacity[0]*natoms;
 442   tcap   = 1.0;
 443   nodeid = 0;
 444   /* Start at bor is 1, to force the first block on the first processor */
 445   for(i=0; (i<nbor) && (tload < natoms); i++) {
 446     if(i<(nbor-1))
 447       b1=border[i+1].atom;
 448     else
 449       b1=natoms;
 450
 451     b0 = border[i].atom;
 452
 453     if ((fabs(b0-tload)<fabs(b1-tload))) {
 454       /* New nodeid time */
 455       multinr_cgs[nodeid] = border[i].ic;
 456       maxatom[nodeid]     = b0;
 457       tcap -= capacity[nodeid];
 458       nodeid++;
 459
 460       /* Recompute target load */
 461       tload = b0 + (natoms-b0)*capacity[nodeid]/tcap;
 462
 463     if(debug)
 464           printf("tload: %g tcap: %g  nodeid: %d\n",tload,tcap,nodeid);
 465     }
 466   }
 467   /* Now the last one... */
 468   while (nodeid < nnodes) {
 469     multinr_cgs[nodeid] = cgs->nr;
 470     /* Store atom number, see above */
 471     maxatom[nodeid]     = natoms;
 472     nodeid++;
 473   }
 474   if (nodeid != nnodes) {
 475     gmx_fatal(FARGS,"nodeid = %d, nnodes = %d, file %s, line %d",
 476                 nodeid,nnodes,__FILE__,__LINE__);
 477   }
 478
 479   for(i=nnodes-1; (i>0); i--)
 480     maxatom[i]-=maxatom[i-1];
 481
 482   if (fp) {
 483     fprintf(fp,"Division over nodes in atoms:\n");
 484     for(i=0; (i<nnodes); i++)
 485       fprintf(fp," %7d",maxatom[i]);
 486     fprintf(fp,"\n");
 487   }
 488
 489   sfree(maxatom);
 490   sfree(shknum);
 491   sfree(cgsnum);
 492   sfree(border);
 493 }
 494
 495 typedef struct {
 496   int atom,sid;
 497 } t_sid;
 498
 499 static int sid_comp(const void *a,const void *b)
 500 {
 501   t_sid *sa,*sb;
 502   int   dd;
 503
 504   sa=(t_sid *)a;
 505   sb=(t_sid *)b;
 506
 507   dd=sa->sid-sb->sid;
 508   if (dd == 0)
 509     return (sa->atom-sb->atom);
 510   else
 511     return dd;
 512 }
 513
 514 static int mk_grey(int nnodes,egCol egc[],t_graph *g,int *AtomI,
 515                    int maxsid,t_sid sid[])
 516 {
 517   int  j,ng,ai,aj,g0;
 518
 519   ng=0;
 520   ai=*AtomI;
 521
 522   g0=g->at_start;
 523   /* Loop over all the bonds */
 524   for(j=0; (j<g->nedge[ai]); j++) {
 525     aj=g->edge[ai][j]-g0;
 526     /* If there is a white one, make it gray and set pbc */
 527     if (egc[aj] == egcolWhite) {
 528       if (aj < *AtomI)
 529         *AtomI = aj;
 530       egc[aj] = egcolGrey;
 531
 532       /* Check whether this one has been set before... */
 533       range_check(aj+g0,0,maxsid);
 534       range_check(ai+g0,0,maxsid);
 535       if (sid[aj+g0].sid != -1)
 536         gmx_fatal(FARGS,"sid[%d]=%d, sid[%d]=%d, file %s, line %d",
 537                     ai,sid[ai+g0].sid,aj,sid[aj+g0].sid,__FILE__,__LINE__);
 538       else {
 539         sid[aj+g0].sid  = sid[ai+g0].sid;
 540         sid[aj+g0].atom = aj+g0;
 541       }
 542       ng++;
 543     }
 544   }
 545   return ng;
 546 }
 547
 548 static int first_colour(int fC,egCol Col,t_graph *g,egCol egc[])
 549 /* Return the first node with colour Col starting at fC.
 550  * return -1 if none found.
 551  */
 552 {
 553   int i;
 554
 555   for(i=fC; (i<g->nnodes); i++)
 556     if ((g->nedge[i] > 0) && (egc[i]==Col))
 557       return i;
 558
 559   return -1;
 560 }
 561
 562 static int mk_sblocks(FILE *fp,t_graph *g,int maxsid,t_sid sid[])
 563 {
 564   int    ng,nnodes;
 565   int    nW,nG,nB;              /* Number of Grey, Black, White */
 566   int    fW,fG;                 /* First of each category       */
 567   egCol  *egc=NULL;             /* The colour of each node      */
 568   int    g0,nblock;
 569
 570   if (!g->nbound)
 571     return 0;
 572
 573   nblock=0;
 574
 575   nnodes=g->nnodes;
 576   snew(egc,nnodes);
 577
 578   g0=g->at_start;
 579   nW=g->nbound;
 580   nG=0;
 581   nB=0;
 582
 583   fW=0;
 584
 585   /* We even have a loop invariant:
 586    * nW+nG+nB == g->nbound
 587    */
 588
 589   if (fp)
 590     fprintf(fp,"Walking down the molecule graph to make constraint-blocks\n");
 591
 592   while (nW > 0) {
 593     /* Find the first white, this will allways be a larger
 594      * number than before, because no nodes are made white
 595      * in the loop
 596      */
 597     if ((fW=first_colour(fW,egcolWhite,g,egc)) == -1)
 598       gmx_fatal(FARGS,"No WHITE nodes found while nW=%d\n",nW);
 599
 600     /* Make the first white node grey, and set the block number */
 601     egc[fW]        = egcolGrey;
 602     range_check(fW+g0,0,maxsid);
 603     sid[fW+g0].sid = nblock++;
 604     nG++;
 605     nW--;
 606
 607     /* Initial value for the first grey */
 608     fG=fW;
 609
 610     if (debug)
 611       fprintf(debug,"Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
 612               nW,nG,nB,nW+nG+nB);
 613
 614     while (nG > 0) {
 615       if ((fG=first_colour(fG,egcolGrey,g,egc)) == -1)
 616         gmx_fatal(FARGS,"No GREY nodes found while nG=%d\n",nG);
 617
 618       /* Make the first grey node black */
 619       egc[fG]=egcolBlack;
 620       nB++;
 621       nG--;
 622
 623       /* Make all the neighbours of this black node grey
 624        * and set their block number
 625        */
 626       ng=mk_grey(nnodes,egc,g,&fG,maxsid,sid);
 627       /* ng is the number of white nodes made grey */
 628       nG+=ng;
 629       nW-=ng;
 630     }
 631   }
 632   sfree(egc);
 633
 634   if (debug)
 635     fprintf(debug,"Found %d shake blocks\n",nblock);
 636
 637   return nblock;
 638 }
 639
 640
 641 typedef struct {
 642   int first,last,sid;
 643 } t_merge_sid;
 644
 645 static int ms_comp(const void *a, const void *b)
 646 {
 647   t_merge_sid *ma = (t_merge_sid *)a;
 648   t_merge_sid *mb = (t_merge_sid *)b;
 649   int d;
 650
 651   d = ma->first-mb->first;
 652   if (d == 0)
 653     return ma->last-mb->last;
 654   else
 655     return d;
 656 }
 657
 658 static int merge_sid(int i0,int at_start,int at_end,int nsid,t_sid sid[],
 659                      t_blocka *sblock)
 660 {
 661   int  i,j,k,n,isid,ndel;
 662   t_merge_sid *ms;
 663   int  nChanged;
 664
 665   /* We try to remdy the following problem:
 666    * Atom: 1  2  3  4  5 6 7 8 9 10
 667    * Sid:  0 -1  1  0 -1 1 1 1 1 1
 668    */
 669
 670   /* Determine first and last atom in each shake ID */
 671   snew(ms,nsid);
 672
 673   for(k=0; (k<nsid); k++) {
 674     ms[k].first = at_end+1;
 675     ms[k].last  = -1;
 676     ms[k].sid   = k;
 677   }
 678   for(i=at_start; (i<at_end); i++) {
 679     isid = sid[i].sid;
 680     range_check(isid,-1,nsid);
 681     if (isid >= 0) {
 682       ms[isid].first = min(ms[isid].first,sid[i].atom);
 683       ms[isid].last  = max(ms[isid].last,sid[i].atom);
 684     }
 685   }
 686   qsort(ms,nsid,sizeof(ms[0]),ms_comp);
 687
 688   /* Now merge the overlapping ones */
 689   ndel = 0;
 690   for(k=0; (k<nsid); ) {
 691     for(j=k+1; (j<nsid); ) {
 692       if (ms[j].first <= ms[k].last) {
 693         ms[k].last  = max(ms[k].last,ms[j].last);
 694         ms[k].first = min(ms[k].first,ms[j].first);
 695         ms[j].sid = -1;
 696         ndel++;
 697         j++;
 698       }
 699       else {
 700         k = j;
 701         j = k+1;
 702       }
 703     }
 704     if (j == nsid)
 705       k++;
 706   }
 707   for(k=0; (k<nsid); k++) {
 708     while ((k < nsid-1) && (ms[k].sid == -1)) {
 709       for(j=k+1; (j<nsid); j++) {
 710         memcpy(&(ms[j-1]),&(ms[j]),sizeof(ms[0]));
 711       }
 712       nsid--;
 713     }
 714   }
 715
 716   for(k=at_start; (k<at_end); k++) {
 717     sid[k].atom = k;
 718     sid[k].sid = -1;
 719   }
 720   sblock->nr  = nsid;
 721   sblock->nalloc_index = sblock->nr+1;
 722   snew(sblock->index,sblock->nalloc_index);
 723   sblock->nra = at_end - at_start;
 724   sblock->nalloc_a = sblock->nra;
 725   snew(sblock->a,sblock->nalloc_a);
 726   sblock->index[0] = 0;
 727   for(k=n=0; (k<nsid); k++) {
 728     sblock->index[k+1] = sblock->index[k] + ms[k].last - ms[k].first+1;
 729     for(j=ms[k].first; (j<=ms[k].last); j++) {
 730       range_check(n,0,sblock->nra);
 731       sblock->a[n++] = j;
 732       range_check(j,0,at_end);
 733       if (sid[j].sid == -1)
 734         sid[j].sid = k;
 735       else
 736         fprintf(stderr,"Double sids (%d, %d) for atom %d\n",sid[j].sid,k,j);
 737     }
 738   }
 739   assert(k == nsid);
 740   /* Removed 2007-09-04
 741   sblock->index[k+1] = natoms;
 742   for(k=0; (k<natoms); k++)
 743     if (sid[k].sid == -1)
 744       sblock->a[n++] = k;
 745   assert(n == natoms);
 746   */
 747   sblock->nra = n;
 748   assert(sblock->index[k] == sblock->nra);
 749   sfree(ms);
 750
 751   return nsid;
 752 }
 753
 754 void gen_sblocks(FILE *fp,int at_start,int at_end,
 755                                  t_idef *idef,t_blocka *sblock,
 756                                  gmx_bool bSettle)
 757 {
 758   t_graph *g;
 759   int     i,i0,j,k,istart,n;
 760   t_sid   *sid;
 761   int     isid,nsid;
 762
 763   g=mk_graph(NULL,idef,at_start,at_end,TRUE,bSettle);
 764   if (debug)
 765     p_graph(debug,"Graaf Dracula",g);
 766   snew(sid,at_end);
 767   for(i=at_start; (i<at_end); i++) {
 768     sid[i].atom =  i;
 769     sid[i].sid  = -1;
 770   }
 771   nsid=mk_sblocks(fp,g,at_end,sid);
 772
 773   if (!nsid)
 774     return;
 775
 776   /* Now sort the shake blocks... */
 777   qsort(sid+at_start,at_end-at_start,(size_t)sizeof(sid[0]),sid_comp);
 778
 779   if (debug) {
 780     fprintf(debug,"Sorted shake block\n");
 781     for(i=at_start; (i<at_end); i++)
 782       fprintf(debug,"sid[%5d] = atom:%5d sid:%5d\n",i,sid[i].atom,sid[i].sid);
 783   }
 784   /* Now check how many are NOT -1, i.e. how many have to be shaken */
 785   for(i0=at_start; (i0<at_end); i0++)
 786     if (sid[i0].sid > -1)
 787       break;
 788
 789   /* Now we have the sids that have to be shaken. We'll check the min and
 790    * max atom numbers and this determines the shake block. DvdS 2007-07-19.
 791    * For the purpose of making boundaries all atoms in between need to be
 792    * part of the shake block too. There may be cases where blocks overlap
 793    * and they will have to be merged.
 794    */
 795   nsid = merge_sid(i0,at_start,at_end,nsid,sid,sblock);
 796   /* Now sort the shake blocks again... */
 797   /*qsort(sid,natoms,(size_t)sizeof(sid[0]),sid_comp);*/
 798
 799   /* Fill the sblock struct */
 800   /*  sblock->nr  = nsid;
 801   sblock->nra = natoms;
 802   srenew(sblock->a,sblock->nra);
 803   srenew(sblock->index,sblock->nr+1);
 804
 805   i    = i0;
 806   isid = sid[i].sid;
 807   n    = k = 0;
 808   sblock->index[n++]=k;
 809   while (i < natoms) {
 810     istart = sid[i].atom;
 811     while ((i<natoms-1) && (sid[i+1].sid == isid))
 812     i++;*/
 813     /* After while: we found a new block, or are thru with the atoms */
 814   /*    for(j=istart; (j<=sid[i].atom); j++,k++)
 815       sblock->a[k]=j;
 816     sblock->index[n] = k;
 817     if (i < natoms-1)
 818       n++;
 819     if (n > nsid)
 820       gmx_fatal(FARGS,"Death Horror: nsid = %d, n= %d",nsid,n);
 821     i++;
 822     isid = sid[i].sid;
 823   }
 824   */
 825   sfree(sid);
 826   /* Due to unknown reason this free generates a problem sometimes */
 827   done_graph(g);
 828   sfree(g);
 829   if (debug)
 830     fprintf(debug,"Done gen_sblocks\n");
 831 }
 832
 833 static t_blocka block2blocka(t_block *block)
 834 {
 835   t_blocka blocka;
 836   int i;
 837
 838   blocka.nr = block->nr;
 839   blocka.nalloc_index = blocka.nr + 1;
 840   snew(blocka.index,blocka.nalloc_index);
 841   for(i=0; i<=block->nr; i++)
 842     blocka.index[i] = block->index[i];
 843   blocka.nra = block->index[block->nr];
 844   blocka.nalloc_a = blocka.nra;
 845   snew(blocka.a,blocka.nalloc_a);
 846   for(i=0; i<blocka.nra; i++)
 847     blocka.a[i] = i;
 848
 849   return blocka;
 850 }
 851
 852 typedef struct
 853 {
 854         int   nconstr;
 855         int   index[10];
 856 } pd_constraintlist_t;
 857
 858
 859 static void
 860 find_constraint_range_recursive(pd_constraintlist_t *  constraintlist,
 861                                                                 int                    thisatom,
 862                                                                 int                    depth,
 863                                                                 int *                  min_atomid,
 864                                                                 int *                  max_atomid)
 865 {
 866         int i,j;
 867         int nconstr;
 868
 869         for(i=0;i<constraintlist[thisatom].nconstr;i++)
 870         {
 871                 j = constraintlist[thisatom].index[i];
 872
 873                 *min_atomid = (j<*min_atomid) ? j : *min_atomid;
 874                 *max_atomid = (j>*max_atomid) ? j : *max_atomid;
 875
 876                 if(depth>0)
 877                 {
 878                         find_constraint_range_recursive(constraintlist,j,depth-1,min_atomid,max_atomid);
 879                 }
 880         }
 881 }
 882
 883 static void
 884 pd_determine_constraints_range(t_inputrec *      ir,
 885                                                            int               natoms,
 886                                                            t_ilist *         ilist,
 887                                                            int               hid[],
 888                                                            int *             min_nodeid,
 889                                                            int *             max_nodeid)
 890 {
 891         int i,j,k;
 892         int nratoms;
 893         int depth;
 894         int ai,aj;
 895         int min_atomid,max_atomid;
 896         pd_constraintlist_t *constraintlist;
 897
 898         nratoms = interaction_function[F_CONSTR].nratoms;
 899         depth   = ir->nProjOrder;
 900
 901         snew(constraintlist,natoms);
 902
 903         /* Make a list of all the connections */
 904         for(i=0;i<ilist->nr;i+=nratoms+1)
 905         {
 906                 ai=ilist->iatoms[i+1];
 907                 aj=ilist->iatoms[i+2];
 908                 constraintlist[ai].index[constraintlist[ai].nconstr++]=aj;
 909                 constraintlist[aj].index[constraintlist[aj].nconstr++]=ai;
 910         }
 911
 912         for(i=0;i<natoms;i++)
 913         {
 914                 min_atomid = i;
 915                 max_atomid = i;
 916
 917                 find_constraint_range_recursive(constraintlist,i,depth,&min_atomid,&max_atomid);
 918
 919                 min_nodeid[i] = hid[min_atomid];
 920                 max_nodeid[i] = hid[max_atomid];
 921         }
 922         sfree(constraintlist);
 923 }
 924
 925
 926 void split_top(FILE *fp,int nnodes,gmx_localtop_t *top,t_inputrec *ir,t_block *mols,
 927                real *capacity,int *multinr_cgs,int **multinr_nre, int *left_range, int * right_range)
 928 {
 929   int     natoms,i,j,k,mj,atom,maxatom,sstart,send,bstart,nodeid;
 930   t_blocka sblock;
 931   int     *homeind;
 932   int ftype,nvsite_constr,nra,nrd;
 933   t_iatom   *ia;
 934   int minhome,ihome,minidx;
 935         int *constr_min_nodeid;
 936         int *constr_max_nodeid;
 937
 938   if (nnodes <= 1)
 939     return;
 940
 941   natoms = mols->index[mols->nr];
 942
 943   if (fp)
 944     fprintf(fp,"splitting topology...\n");
 945
 946 /*#define MOL_BORDER*/
 947 /*Removed the above to allow splitting molecules with h-bond constraints
 948   over processors. The results in DP are the same. */
 949   init_blocka(&sblock);
 950   if(ir->eConstrAlg != econtLINCS)
 951   {
 952 #ifndef MOL_BORDER
 953   /* Make a special shake block that includes settles */
 954           gen_sblocks(fp,0,natoms,&top->idef,&sblock,TRUE);
 955 #else
 956           sblock = block2blocka(mols);
 957 #endif
 958   }
 959
 960   split_blocks(fp,ir,nnodes,&top->cgs,&sblock,capacity,multinr_cgs);
 961
 962   homeind=home_index(nnodes,&top->cgs,multinr_cgs);
 963
 964   snew(constr_min_nodeid,natoms);
 965   snew(constr_max_nodeid,natoms);
 966
 967   if(top->idef.il[F_CONSTR].nr>0)
 968   {
 969           pd_determine_constraints_range(ir,natoms,&top->idef.il[F_CONSTR],homeind,constr_min_nodeid,constr_max_nodeid);
 970   }
 971   else
 972   {
 973           /* Not 100% necessary, but it is a bad habit to have uninitialized arrays around... */
 974           for(i=0;i<natoms;i++)
 975           {
 976                   constr_min_nodeid[i] = constr_max_nodeid[i] = homeind[i];
 977           }
 978   }
 979
 980   /* Default limits (no communication) for PD constraints */
 981   left_range[0] = 0;
 982   for(i=1;i<nnodes;i++)
 983   {
 984         left_range[i]    = top->cgs.index[multinr_cgs[i-1]];
 985         right_range[i-1] = left_range[i]-1;
 986   }
 987   right_range[nnodes-1] = top->cgs.index[multinr_cgs[nnodes-1]]-1;
 988
 989   for(j=0; (j<F_NRE); j++)
 990     split_force2(ir, nnodes,homeind,j,&top->idef.il[j],multinr_nre[j],constr_min_nodeid,constr_max_nodeid,
 991                                  left_range,right_range);
 992
 993   sfree(constr_min_nodeid);
 994   sfree(constr_max_nodeid);
 995
 996   sfree(homeind);
 997   done_blocka(&sblock);
 998 }
 999