src/gromacs/gmxlib/splitter.c

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