Project

General

Profile

splitter.c

spliiter.c source filee - David van der Spoel, 01/04/2008 09:48 AM

 
1
/*
2
 * $Id: splitter.c,v 1.35.2.5 2007/10/20 07:52:23 spoel Exp $
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.3.2
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-2007, 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
 * Groningen Machine for Chemical Simulation
35
 */
36
#ifdef HAVE_CONFIG_H
37
#include <config.h>
38
#endif
39

    
40
#include <stdio.h>
41
#include <string.h>
42

    
43
#include "sysstuff.h"
44
#include "macros.h"
45
#include "smalloc.h"
46
#include "assert.h"
47
#include "typedefs.h"
48
#include "mshift.h"
49
#include "invblock.h"
50
#include "txtdump.h"
51
#include "math.h"
52
#include "splitter.h"
53

    
54
typedef struct {
55
  int     nr;
56
  t_iatom *ia;
57
} t_sf;
58

    
59
static void init_sf(int nr,t_sf sf[])
60
{
61
  int i;
62

    
63
  for(i=0; (i<nr); i++) {
64
    sf[i].nr=0;
65
    sf[i].ia=NULL;
66
  }
67
}
68

    
69
static void done_sf(int nr,t_sf sf[])
70
{
71
  int i;
72

    
73
  for(i=0; (i<nr); i++) {
74
    sf[i].nr=0;
75
    sfree(sf[i].ia);
76
    sf[i].ia=NULL;
77
  }
78
}
79

    
80
static void push_sf(t_sf *sf,int nr,t_iatom ia[])
81
{
82
  int i;
83

    
84
  srenew(sf->ia,sf->nr+nr);
85
  for(i=0; (i<nr); i++)
86
    sf->ia[sf->nr+i]=ia[i];
87
  sf->nr+=nr;
88
}
89

    
90
static int min_nodeid(int nr,atom_id list[],int hid[])
91
{
92
  int i,nodeid,minnodeid;
93

    
94
  if (nr <= 0)
95
    gmx_incons("Invalid node number");
96
  minnodeid=hid[list[0]];
97
  for (i=1; (i<nr); i++) 
98
    if ((nodeid=hid[list[i]]) < minnodeid) 
99
      minnodeid=nodeid;
100

    
101
  return minnodeid;
102
}
103

    
104
static void split_force2(int nnodes,int hid[],t_idef *idef,t_ilist *ilist)
105
{
106
  int     i,j,k,type,ftype,nodeid,nratoms,tnr;
107
  int     nvsite_constr,tmpid;
108
  t_iatom ai;
109
  t_sf    sf[MAXNODES];
110
  
111
  init_sf(MAXNODES,sf);
112

    
113
  /* Walk along all the bonded forces, find the appropriate node
114
   * to calc it on, and add it to that nodes list.
115
   */
116
  for (i=0; i<ilist->nr; i+=(1+nratoms)) {
117
    type    = ilist->iatoms[i];
118
    ftype   = idef->functype[type];
119
    nratoms = interaction_function[ftype].nratoms;
120

    
121
    if (ftype == F_SHAKE) {
122
      /* SPECIAL CASE: All Atoms must have the same home node! */
123
      nodeid=hid[ilist->iatoms[i+1]];
124
      if (hid[ilist->iatoms[i+2]] != nodeid) 
125
        gmx_fatal(FARGS,"Shake block crossing node boundaries\n"
126
                  "constraint between atoms (%d,%d)",
127
                  ilist->iatoms[i+1]+1,ilist->iatoms[i+2]+1);
128
    }
129
    else if (ftype == F_SETTLE) {
130
      /* Only the first particle is stored for settles ... */
131
      ai=ilist->iatoms[i+1];
132
      nodeid=hid[ai];
133
      if ((nodeid != hid[ai+1]) ||
134
          (nodeid != hid[ai+2]))
135
        gmx_fatal(FARGS,"Settle block crossing node boundaries\n"
136
                  "constraint between atoms (%d-%d)",ai+1,ai+2+1);
137
    }
138
    else if(interaction_function[ftype].flags & IF_VSITE) {
139
      /* Virtual sites should be constructed on the home node */
140
  
141
      if (ftype==F_VSITE2)
142
        nvsite_constr=2;
143
      else if(ftype==F_VSITE4FD)
144
        nvsite_constr=4;
145
      else
146
        nvsite_constr=3;
147
      
148
      tmpid=hid[ilist->iatoms[i+1]];
149
      
150
      for(k=2;k<nvsite_constr+2;k++) {
151
        if(hid[ilist->iatoms[i+k]]<(tmpid-1) ||
152
           hid[ilist->iatoms[i+k]]>(tmpid+1))
153
          gmx_fatal(FARGS,"Virtual site %d and its constructing"
154
                      " atoms are not on the same or adjacent\n" 
155
                      " nodes. This is necessary to avoid a lot\n"
156
                      " of extra communication. The easiest way"
157
                      " to ensure this is to place virtual sites\n"
158
                      " close to the constructing atoms.\n"
159
                      " Sorry, but you will have to rework your topology!\n",
160
                      ilist->iatoms[i+1]);
161
      }
162
      nodeid=min_nodeid(nratoms,&ilist->iatoms[i+1],hid);
163
    } else
164
      nodeid=min_nodeid(nratoms,&ilist->iatoms[i+1],hid);
165

    
166
    /* Add it to the list */
167
    push_sf(&(sf[nodeid]),nratoms+1,&(ilist->iatoms[i]));
168
    }
169
  tnr=0;
170
  for(nodeid=0; (nodeid<MAXNODES); nodeid++) {
171
    for (i=0; (i<sf[nodeid].nr); i++) 
172
      ilist->iatoms[tnr++]=sf[nodeid].ia[i];
173

    
174
    ilist->multinr[nodeid]=(nodeid==0) ? 0 : ilist->multinr[nodeid-1];
175
    ilist->multinr[nodeid]+=sf[nodeid].nr;
176
  }
177
  if (tnr != ilist->nr)
178
    gmx_incons("Splitting forces over processors");
179
  done_sf(MAXNODES,sf);
180
}
181

    
182
static int *home_index(int nnodes,t_block *cgs)
183
{
184
  /* This routine determines the node id for each particle */
185
  int *hid;
186
  int nodeid,j0,j1,j,k,ak;
187
  
188
  snew(hid,cgs->nra);
189
  /* Initiate to -1 to make it possible to check afterwards,
190
   * all hid's should be set in the loop below
191
   */
192
  for(k=0; (k<cgs->nra); k++)
193
    hid[k]=-1;
194
    
195
  /* loop over nodes */
196
  for(nodeid=0; (nodeid<nnodes); nodeid++) {
197
    j0 = (nodeid==0) ? 0 : cgs->multinr[nodeid-1];
198
    j1 = cgs->multinr[nodeid];
199
    
200
    /* j0 and j1 are the boundariesin the index array */
201
    for(j=j0; (j<j1); j++) {
202
      for(k=cgs->index[j]; (k<cgs->index[j+1]); k++) {
203
        ak=cgs->a[k];
204
        hid[ak]=nodeid;
205
      }
206
    }
207
  }
208
  /* Now verify that all hid's are not -1 */
209
  for(k=0; (k<cgs->nra); k++)
210
    if (hid[k] == -1)
211
      gmx_fatal(FARGS,"hid[%d] = -1, cgs->nr = %d, cgs->nra = %d",
212
                  k,cgs->nr,cgs->nra);
213
  
214
  return hid;
215
}
216

    
217
static int max_index(int start,t_block *b)
218
{
219
  int k,k0,k1;
220
  int mi=-1;
221

    
222
  if (start < b->nr) {  
223
    k0=b->index[start];
224
    k1=b->index[start+1];
225
    if (k1 > k0) {
226
      mi=b->a[k0];
227
      for(k=k0+1; (k<k1); k++)
228
        mi=max(mi,b->a[k]);
229
    }
230
  }
231
  return mi;
232
}
233

    
234
static int min_index(int start,t_block *b)
235
{
236
  int k,k0,k1;
237
  int mi=INT_MAX;
238

    
239
  if (start < b->nr) {  
240
    k0=b->index[start];
241
    k1=b->index[start+1];
242
    if (k1 > k0) {
243
      mi=b->a[k0];
244
      for(k=k0+1; (k<k1); k++)
245
        mi=min(mi,b->a[k]);
246
    }
247
  }
248
  return mi;
249
}
250

    
251
typedef struct {
252
  int atom,ic,is;
253
} t_border;
254

    
255
void set_bor(t_border *b,int atom,int ic,int is)
256
{
257
  if (debug)
258
    fprintf(debug,"border @ atom %5d [ ic = %5d,  is = %5d ]\n",atom,ic,is);
259
  b->atom = atom;
260
  b->ic   = ic;
261
  b->is   = is;
262
}
263

    
264
static bool is_bor(atom_id ai[],int i)
265
{
266
  return ((ai[i] != ai[i-1]) || ((ai[i] == NO_ATID) && (ai[i-1] == NO_ATID)));
267
}
268

    
269
t_border *mk_border(bool bVerbose,int natom,atom_id *invcgs,
270
                    atom_id *invshk,int *nb)
271
{
272
  t_border *bor;
273
  atom_id  *sbor,*cbor;
274
  int      i,j,is,ic,ns,nc,nbor;
275

    
276
  if (debug) {
277
    for(i=0; (i<natom); i++) {
278
      fprintf(debug,"atom: %6d  cgindex: %6d  shkindex: %6d\n",
279
              i,invcgs[i],invshk[i]);
280
    }
281
  }
282
    
283
  snew(sbor,natom+1);
284
  snew(cbor,natom+1);
285
  ns = nc = 1;
286
  for(i=1; (i<natom); i++) {
287
    if (is_bor(invcgs,i)) 
288
      cbor[nc++] = i;
289
    if (is_bor(invshk,i)) 
290
      sbor[ns++] = i;
291
  }
292
  sbor[ns] = 0;
293
  cbor[nc] = 0;
294
  fprintf(stderr,"There are %d charge group borders and %d shake borders\n",
295
          nc,ns);
296
  snew(bor,max(nc,ns));
297
  ic = is = nbor = 0;
298
  while ((ic < nc) || (is < ns)) {
299
    if (sbor[is] == cbor[ic]) {
300
      set_bor(&(bor[nbor]),cbor[ic],ic,is);
301
      nbor++;
302
      if (ic < nc) ic++;
303
      if (is < ns) is++;
304
    }
305
    else if (cbor[ic] > sbor[is]) {
306
      if (is == ns) {
307
        set_bor(&(bor[nbor]),cbor[ic],ic,is);
308
        nbor++;
309
        if (ic < nc) ic++;
310
      }
311
      else if (is < ns) 
312
        is++;
313
    }
314
    else {
315
      if (ic == nc) {
316
        set_bor(&(bor[nbor]),cbor[ic],ic,is);
317
        nbor++;
318
        if (is < ns) is++;
319
      }
320
      else if (ic < nc) 
321
        ic++;
322
    }
323
    /*if (ic < nc)
324
      ic++;
325
      else
326
      is++;*//*gmx_fatal(FARGS,"Can't happen is=%d, ic=%d (%s, %d)",
327
               is,ic,__FILE__,__LINE__);*/
328
  }
329
  fprintf(stderr,"There are %d total borders\n",nbor);
330

    
331
  if (debug) {
332
    fprintf(debug,"There are %d actual bor entries\n",nbor);
333
    for(i=0; (i<nbor); i++) 
334
      fprintf(debug,"bor[%5d] = atom: %d  ic: %d  is: %d\n",i,
335
              bor[i].atom,bor[i].ic,bor[i].is);
336
  }
337
  
338
  *nb = nbor;
339
  
340
  return bor;
341
}
342

    
343
static void split_blocks(bool bVerbose,int nnodes,
344
                         t_block *cgs,t_block *sblock,real capacity[])
345
{
346
  int      maxatom[MAXNODES];
347
  int      i,ii,ai,b0,b1;
348
  int      nodeid,last_shk,nbor;
349
  t_border *border;
350
  double   tload,tcap;
351
  
352
  bool    bSHK;
353
  atom_id *shknum,*cgsnum;
354
  
355
  if (debug) {
356
    pr_block(debug,0,"cgs",cgs,TRUE);
357
    pr_block(debug,0,"sblock",sblock,TRUE);
358
  }
359

    
360
  shknum = make_invblock(sblock,cgs->nra+1);
361
  cgsnum = make_invblock(cgs,cgs->nra+1);
362
  border = mk_border(bVerbose,cgs->nra,cgsnum,shknum,&nbor);
363

    
364
  tload  = capacity[0]*cgs->nra;
365
  tcap   = 1.0;
366
  nodeid = 0;
367
  /* Start at bor is 1, to force the first block on the first processor */
368
  for(i=0; (i<nbor) && (tload < cgs->nra); i++) {
369
    if(i<(nbor-1)) 
370
      b1=border[i+1].atom;
371
    else
372
      b1=cgs->nra;
373

    
374
    b0 = border[i].atom;
375
    
376
    if ((fabs(b0-tload)<fabs(b1-tload))) {
377
      /* New nodeid time */
378
      cgs->multinr[nodeid]    = border[i].ic;
379
      /* Store the atom number here, has to be processed later */
380
      sblock->multinr[nodeid] = border[i].atom;
381
      maxatom[nodeid]         = b0;
382
      tcap -= capacity[nodeid];
383
      nodeid++;
384
      
385
      /* Recompute target load */
386
      tload = b0 +
387
        (cgs->nra-b0)*capacity[nodeid]/tcap;
388

    
389
      if (debug)
390
        fprintf(debug,"tload: %g tcap: %g  nodeid: %d\n",tload,tcap,nodeid);
391
    } 
392
  }
393
  /* Now the last one... */
394
  while (nodeid < nnodes) {
395
    cgs->multinr[nodeid]    = cgs->nr;
396
    /* Store atom number, see above */
397
    sblock->multinr[nodeid] = cgs->nra;
398
    maxatom[nodeid]         = cgs->nra;
399
    nodeid++;
400
  }
401
  if (nodeid != nnodes) {
402
    gmx_fatal(FARGS,"nodeid = %d, nnodes = %d, file %s, line %d",
403
                nodeid,nnodes,__FILE__,__LINE__);
404
  }
405

    
406
  if (bVerbose) {
407
    for(i=nnodes-1; (i>0); i--)
408
      maxatom[i]-=maxatom[i-1];
409
    fprintf(stderr,"Division over nodes in atoms:\n");
410
    for(i=0; (i<nnodes); i++)
411
      fprintf(stderr," %7d",maxatom[i]);
412
    fprintf(stderr,"\n");
413
  }
414
  sfree(shknum);
415
  sfree(cgsnum);
416
  sfree(border);
417
}
418

    
419
static void def_mnr(int nr,int mnr[])
420
{
421
  int i;
422

    
423
  for (i=0; (i<MAXNODES); i++) 
424
    mnr[i]=0;
425
  mnr[0]=nr;
426
}
427

    
428
typedef struct {
429
  int atom,sid;
430
} t_sid;
431

    
432
static int sid_comp(const void *a,const void *b)
433
{
434
  t_sid *sa,*sb;
435
  int   dd;
436
  
437
  sa=(t_sid *)a;
438
  sb=(t_sid *)b;
439
  
440
  dd=sa->sid-sb->sid;
441
  if (dd == 0)
442
    return (sa->atom-sb->atom);
443
  else
444
    return dd;
445
}
446

    
447
static int mk_grey(int nnodes,egCol egc[],t_graph *g,int *AtomI,
448
                   t_sid sid[])
449
{
450
  int  j,ng,ai,aj,g0;
451

    
452
  ng=0;
453
  ai=*AtomI;
454
  
455
  g0=g->start;
456
  /* Loop over all the bonds */
457
  for(j=0; (j<g->nedge[ai]); j++) {
458
    aj=g->edge[ai][j]-g0;
459
    /* If there is a white one, make it gray and set pbc */
460
    if (egc[aj] == egcolWhite) {
461
      if (aj < *AtomI)
462
        *AtomI = aj;
463
      egc[aj] = egcolGrey;
464

    
465
      /* Check whether this one has been set before... */
466
      if (sid[aj+g0].sid != -1) 
467
        gmx_fatal(FARGS,"sid[%d]=%d, sid[%d]=%d, file %s, line %d",
468
                    ai,sid[ai+g0].sid,aj,sid[aj+g0].sid,__FILE__,__LINE__);
469
      else {
470
        sid[aj+g0].sid  = sid[ai+g0].sid;
471
        sid[aj+g0].atom = aj+g0;
472
      }
473
      ng++;
474
    }
475
  }
476
  return ng;
477
}
478

    
479
static int first_colour(int fC,egCol Col,t_graph *g,egCol egc[])
480
/* Return the first node with colour Col starting at fC.
481
 * return -1 if none found.
482
 */
483
{
484
  int i;
485
  
486
  for(i=fC; (i<g->nnodes); i++)
487
    if ((g->nedge[i] > 0) && (egc[i]==Col))
488
      return i;
489
  
490
  return -1;
491
}
492

    
493
static int mk_sblocks(bool bVerbose,t_graph *g,t_sid sid[])
494
{
495
  int    ng,nnodes;
496
  int    nW,nG,nB;                /* Number of Grey, Black, White        */
497
  int    fW,fG;                        /* First of each category        */
498
  egCol  *egc=NULL;                /* The colour of each node        */
499
  int    g0,nblock;
500
  
501
  if (!g->nbound) 
502
    return 0;
503

    
504
  nblock=0;
505
  
506
  nnodes=g->nnodes;
507
  snew(egc,nnodes);
508
  
509
  g0=g->start;
510
  nW=g->nbound;
511
  nG=0;
512
  nB=0;
513

    
514
  fW=0;
515

    
516
  /* We even have a loop invariant:
517
   * nW+nG+nB == g->nbound
518
   */
519
  
520
  if (bVerbose)
521
    fprintf(stderr,"Walking down the molecule graph to make shake-blocks\n");
522

    
523
  while (nW > 0) {
524
    /* Find the first white, this will allways be a larger
525
     * number than before, because no nodes are made white
526
     * in the loop
527
     */
528
    if ((fW=first_colour(fW,egcolWhite,g,egc)) == -1) 
529
      gmx_fatal(FARGS,"No WHITE nodes found while nW=%d\n",nW);
530
    
531
    /* Make the first white node grey, and set the block number */
532
    egc[fW]        = egcolGrey;
533
    sid[fW+g0].sid = nblock++;
534
    nG++;
535
    nW--;
536

    
537
    /* Initial value for the first grey */
538
    fG=fW;
539

    
540
    if (debug)
541
      fprintf(debug,"Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
542
              nW,nG,nB,nW+nG+nB);
543
              
544
    while (nG > 0) {
545
      if ((fG=first_colour(fG,egcolGrey,g,egc)) == -1)
546
        gmx_fatal(FARGS,"No GREY nodes found while nG=%d\n",nG);
547
      
548
      /* Make the first grey node black */
549
      egc[fG]=egcolBlack;
550
      nB++;
551
      nG--;
552

    
553
      /* Make all the neighbours of this black node grey
554
       * and set their block number
555
       */
556
      ng=mk_grey(nnodes,egc,g,&fG,sid);
557
      /* ng is the number of white nodes made grey */
558
      nG+=ng;
559
      nW-=ng;
560
    }
561
  }
562
  sfree(egc);
563

    
564
  if (debug)
565
    fprintf(debug,"Found %d shake blocks\n",nblock);
566
    
567
  return nblock;
568
}
569

    
570
typedef struct {
571
  int first,last,sid;
572
} t_merge_sid;
573

    
574
static int ms_comp(const void *a, const void *b)
575
{
576
  t_merge_sid *ma = (t_merge_sid *)a;
577
  t_merge_sid *mb = (t_merge_sid *)b;
578
  int d;
579
  
580
  d = ma->first-mb->first;
581
  if (d == 0)
582
    return ma->last-mb->last;
583
  else
584
    return d;
585
}
586

    
587
static int merge_sid(int i0,int natoms,int nsid,t_sid sid[],
588
                     t_block *sblock)
589
{
590
  int  i,j,k,n,isid,ndel;
591
  t_merge_sid *ms;
592
  int  nChanged;
593

    
594
  /* We try to remdy the following problem:
595
   * Atom: 1  2  3  4  5 6 7 8 9 10
596
   * Sid:  0 -1  1  0 -1 1 1 1 1 1
597
   */
598
   
599
  /* Determine first and last atom in each shake ID */
600
  snew(ms,nsid);
601

    
602
  for(k=0; (k<nsid); k++) {
603
    ms[k].first = natoms+1;
604
    ms[k].last  = -1;
605
    ms[k].sid   = k;
606
  }
607
  for(i=0; (i<natoms); i++) {
608
    isid = sid[i].sid;
609
    if (isid >= 0) {
610
      ms[isid].first = min(ms[isid].first,sid[i].atom);
611
      ms[isid].last  = max(ms[isid].last,sid[i].atom);
612
    }
613
  }
614
  qsort(ms,nsid,sizeof(ms[0]),ms_comp);
615
  
616
  /* Now merge the overlapping ones */
617
  ndel = 0;
618
  for(k=0; (k<nsid); ) {
619
    for(j=k+1; (j<nsid); ) {
620
      if (ms[j].first <= ms[k].last) {
621
        ms[k].last  = max(ms[k].last,ms[j].last);
622
        ms[k].first = min(ms[k].first,ms[j].first);
623
        ms[j].sid = -1;
624
        ndel++;
625
        j++;
626
      }
627
      else {
628
        k = j;
629
        j = k+1;
630
      }
631
    }
632
    if (j == nsid)
633
      k++;
634
  }
635
  for(k=0; (k<nsid); k++) {
636
    while ((k < nsid-1) && (ms[k].sid == -1)) {
637
      for(j=k+1; (j<nsid); j++) {
638
        memcpy(&(ms[j-1]),&(ms[j]),sizeof(ms[0]));
639
      }
640
      nsid--;
641
    }
642
  }
643
  while ((nsid > 0) && (ms[nsid-1].sid == -1)) 
644
    nsid--;
645
    
646
  for(k=0; (k<natoms); k++) {
647
    sid[k].atom = k;
648
    sid[k].sid = -1;
649
  }
650
  sblock->nr  = nsid;
651
  sblock->nra = natoms;
652
  srenew(sblock->a,sblock->nra);
653
  srenew(sblock->index,sblock->nr+2);
654
  sblock->index[0] = 0;
655
  for(k=n=0; (k<nsid); k++) {
656
    sblock->index[k+1] = sblock->index[k] + ms[k].last - ms[k].first+1;
657
    for(j=ms[k].first; (j<=ms[k].last); j++) {
658
      sblock->a[n++] = j;
659
      if (sid[j].sid == -1)
660
        sid[j].sid = k;
661
      else 
662
        fprintf(stderr,"Double sids (%d, %d) for atom %d\n",sid[j].sid,k,j);
663
    }
664
  }
665
  assert(k == nsid);
666
  sblock->index[k+1] = natoms;
667
  for(k=0; (k<natoms); k++)
668
    if (sid[k].sid == -1)
669
      sblock->a[n++] = k;
670
  assert(n == natoms);
671
  
672
  sfree(ms);
673
  
674
  return nsid;
675
}
676

    
677
void gen_sblocks(bool bVerbose,int natoms,t_idef *idef,t_block *sblock,
678
                 bool bSettle)
679
{
680
  t_graph *g;
681
  int     i,i0,j,k,istart,n;
682
  t_sid   *sid;
683
  int     isid,nsid;
684
  
685
  g=mk_graph(idef,natoms,TRUE,bSettle);
686
  if (bVerbose && debug)
687
    p_graph(debug,"Graaf Dracula",g);
688
  snew(sid,natoms);
689
  for(i=0; (i<natoms); i++) {
690
    sid[i].atom =  i;
691
    sid[i].sid  = -1;
692
  }
693
  nsid=mk_sblocks(bVerbose,g,sid);
694

    
695
  if (!nsid)
696
    return;
697
    
698
  /* Now sort the shake blocks... */
699
  qsort(sid,natoms,(size_t)sizeof(sid[0]),sid_comp);
700
  
701
  if (debug) {
702
    fprintf(debug,"Sorted shake block\n");
703
    for(i=0; (i<natoms); i++) 
704
      fprintf(debug,"sid[%5d] = atom:%5d sid:%5d\n",i,sid[i].atom,sid[i].sid);
705
  }
706
  /* Now check how many are NOT -1, i.e. how many have to be shaken */
707
  for(i0=0; (i0<natoms); i0++)
708
    if (sid[i0].sid > -1)
709
      break;
710
  
711
  /* Now we have the sids that have to be shaken. We'll check the min and
712
   * max atom numbers and this determines the shake block. DvdS 2007-07-19.
713
   * For the purpose of making boundaries all atoms in between need to be
714
   * part of the shake block too. There may be cases where blocks overlap
715
   * and they will have to be merged.
716
   */
717
  nsid = merge_sid(i0,natoms,nsid,sid,sblock);
718
  /* Now sort the shake blocks again... */
719
  /*qsort(sid,natoms,(size_t)sizeof(sid[0]),sid_comp);*/
720
  
721
  /* Fill the sblock struct */    
722
  /*  sblock->nr  = nsid;
723
  sblock->nra = natoms;
724
  srenew(sblock->a,sblock->nra);
725
  srenew(sblock->index,sblock->nr+1);
726
  
727
  i    = i0;
728
  isid = sid[i].sid;
729
  n    = k = 0;
730
  sblock->index[n++]=k;
731
  while (i < natoms) {
732
    istart = sid[i].atom;
733
    while ((i<natoms-1) && (sid[i+1].sid == isid))
734
    i++;*/
735
    /* After while: we found a new block, or are thru with the atoms */
736
  /*    for(j=istart; (j<=sid[i].atom); j++,k++)
737
      sblock->a[k]=j;
738
    sblock->index[n] = k;
739
    if (i < natoms-1)
740
      n++;
741
    if (n > nsid)
742
      gmx_fatal(FARGS,"Death Horror: nsid = %d, n= %d",nsid,n);
743
    i++;
744
    isid = sid[i].sid;
745
  }
746
  */
747
  sfree(sid);
748
  /* Due to unknown reason this free generates a problem sometimes */
749
  done_graph(g);
750
  sfree(g); 
751
  if (debug)
752
    fprintf(debug,"Done gen_sblocks\n");
753
}
754

    
755
void split_top(bool bVerbose,int nnodes,t_topology *top,real
756
               *capacity)
757
{
758
  int     i,j,k,mj,atom,maxatom;
759
  t_block sblock;
760
  int     *homeind;
761
  atom_id *sblinv;
762
  int ftype,nvsite_constr,nra,nrd;
763
  t_iatom   *ia;
764
  int minhome,ihome,minidx;
765
  
766
  if ((bVerbose) && (nnodes>1))
767
    fprintf(stderr,"splitting topology...\n");
768
  
769
  for(j=0; (j<F_NRE); j++)  
770
    def_mnr(top->idef.il[j].nr,top->idef.il[j].multinr);
771
  def_mnr(top->atoms.excl.nr,top->atoms.excl.multinr);
772

    
773
  for (j=0; j<ebNR; j++) 
774
    def_mnr(top->blocks[j].nr,top->blocks[j].multinr);
775

    
776
  if (nnodes > 1) {
777
    /* Make a special shake block that includes settles */
778
    init_block(&sblock);
779
    gen_sblocks(bVerbose,top->atoms.nr,&top->idef,&sblock,TRUE);
780

    
781
    split_blocks(bVerbose,nnodes,&(top->blocks[ebCGS]),&sblock,capacity);
782
    
783
    /* Now transform atom numbers to real inverted shake blocks */
784
    sblinv = make_invblock(&(top->blocks[ebSBLOCKS]),top->atoms.nr+1);
785
    for(j=0; (j<MAXNODES); j++) {
786
      atom = sblock.multinr[j];
787
      mj   = NO_ATID;
788
      for(k=(j == 0) ? 0 : sblock.multinr[j-1]; (k<atom); k++)
789
        if (sblinv[k] != NO_ATID)
790
          mj = max(mj,(int)sblinv[k]);
791
      if (mj == NO_ATID) 
792
        mj = (j == 0) ? -1 : top->blocks[ebSBLOCKS].multinr[j-1]-1;
793
      
794
      top->blocks[ebSBLOCKS].multinr[j] = mj+1;
795
    }
796
    sfree(sblinv);
797
    
798
    homeind=home_index(nnodes,&(top->blocks[ebCGS]));
799
    
800
    for(j=0; (j<F_NRE); j++)
801
      split_force2(nnodes,homeind,&top->idef,&top->idef.il[j]);
802
    sfree(homeind);
803
    done_block(&sblock);
804
  }
805
}
806