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