1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This file is part of Gromacs Copyright (c) 1991-2008
5 * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
16 * Gnomes, ROck Monsters And Chili Sauce
33 /* Margin for setting up the DD grid */
34 #define DD_GRID_MARGIN_PRES_SCALE 1.05
36 static int factorize(int n,int **fac,int **mfac)
40 /* Decompose n in factors */
49 if (ndiv == 0 || (*fac)[ndiv-1] != d)
63 static bool fits_pme_ratio(int nnodes,int npme,float ratio)
65 return ((double)npme/(double)nnodes > 0.95*ratio);
68 static bool fits_pp_pme_perf(FILE *fplog,
69 t_inputrec *ir,matrix box,gmx_mtop_t *mtop,
70 int nnodes,int npme,float ratio)
72 int ndiv,*div,*mdiv,ldiv;
74 ndiv = factorize(nnodes-npme,&div,&mdiv);
78 /* The check below gives a reasonable division:
79 * factor 5 allowed at 5 or more PP nodes,
80 * factor 7 allowed at 49 or more PP nodes.
82 if (ldiv > 3 + (int)(pow(nnodes-npme,1.0/3.0) + 0.5))
87 /* Does this division gives a reasonable PME load? */
88 return fits_pme_ratio(nnodes,npme,ratio);
91 static int guess_npme(FILE *fplog,gmx_mtop_t *mtop,t_inputrec *ir,matrix box,
98 ratio = pme_load_estimate(mtop,ir,box);
102 fprintf(fplog,"Guess for relative PME load: %.2f\n",ratio);
105 /* We assume the optimal node ratio is close to the load ratio.
106 * The communication load is neglected,
107 * but (hopefully) this will balance out between PP and PME.
110 if (!fits_pme_ratio(nnodes,nnodes/2,ratio))
112 /* We would need more than nnodes/2 PME only nodes,
113 * which is not possible. Since the PME load is very high,
114 * we will not loose much performance when all nodes do PME.
120 /* First try to find npme as a factor of nnodes up to nnodes/3.
121 * We start with a minimum PME node fraction of 1/16
122 * and avoid ratios which lead to large prime factors in nnodes-npme.
124 npme = (nnodes + 15)/16;
125 while (npme <= nnodes/3) {
126 if (nnodes % npme == 0)
128 /* Note that fits_perf might change the PME grid,
129 * in the current implementation it does not.
131 if (fits_pp_pme_perf(fplog,ir,box,mtop,nnodes,npme,ratio))
140 /* Try any possible number for npme */
142 while (npme <= nnodes/2)
144 /* Note that fits_perf may change the PME grid */
145 if (fits_pp_pme_perf(fplog,ir,box,mtop,nnodes,npme,ratio))
154 gmx_fatal(FARGS,"Could not find an appropriate number of separate PME nodes. i.e. >= %5f*#nodes (%d) and <= #nodes/2 (%d) and reasonable performance wise (grid_x=%d, grid_y=%d).\n"
155 "Use the -npme option of mdrun or change the number of processors or the PME grid dimensions, see the manual for details.",
156 ratio,(int)(0.95*ratio*nnodes+0.5),nnodes/2,ir->nkx,ir->nky);
157 /* Keep the compiler happy */
165 "Will use %d particle-particle and %d PME only nodes\n"
166 "This is a guess, check the performance at the end of the log file\n",
170 "Will use %d particle-particle and %d PME only nodes\n"
171 "This is a guess, check the performance at the end of the log file\n",
178 static int lcd(int n1,int n2)
183 for(i=2; (i<=n1 && i<=n2); i++)
185 if (n1 % i == 0 && n2 % i == 0)
194 static int div_up(int n,int f)
196 return (n + f - 1)/f;
199 real comm_box_frac(ivec dd_nc,real cutoff,gmx_ddbox_t *ddbox)
207 bt[i] = ddbox->box_size[i]*ddbox->skew_fac[i];
208 nw[i] = dd_nc[i]*cutoff/bt[i];
219 for(j=i+1; j<DIM; j++)
223 comm_vol += nw[i]*nw[j]*M_PI/4;
224 for(k=j+1; k<DIM; k++)
228 comm_vol += nw[i]*nw[j]*nw[k]*M_PI/6;
239 static bool inhomogeneous_z(const t_inputrec *ir)
241 return ((EEL_PME(ir->coulombtype) || ir->coulombtype==eelEWALD) &&
242 ir->ePBC==epbcXYZ && ir->ewald_geometry==eewg3DC);
245 static float comm_cost_est(gmx_domdec_t *dd,real limit,real cutoff,
246 matrix box,gmx_ddbox_t *ddbox,
247 int natoms,t_inputrec *ir,
249 int npme_tot,ivec nc)
252 int i,j,k,nk,overlap;
254 float comm_vol,comm_vol_xf,comm_pme,cost_pbcdx;
255 /* This is the cost of a pbc_dx call relative to the cost
256 * of communicating the coordinate and force of an atom.
257 * This will be machine dependent.
258 * These factors are for x86 with SMP or Infiniband.
260 float pbcdx_rect_fac = 0.1;
261 float pbcdx_tric_fac = 0.2;
263 /* Check the DD algorithm restrictions */
264 if ((ir->ePBC == epbcXY && ir->nwall < 2 && nc[ZZ] > 1) ||
265 (ir->ePBC == epbcSCREW && (nc[XX] == 1 || nc[YY] > 1 || nc[ZZ] > 1)))
270 if (inhomogeneous_z(ir) && nc[ZZ] > 1)
275 /* Check if the triclinic requirements are met */
278 for(j=i+1; j<ddbox->npbcdim; j++)
280 if (box[j][i] != 0 || ir->deform[j][i] != 0 ||
281 (ir->epc != epcNO && ir->compress[j][i] != 0))
283 if (nc[j] > 1 && nc[i] == 1)
293 bt[i] = ddbox->box_size[i]*ddbox->skew_fac[i];
295 /* Without PBC there are no cell size limits with 2 cells */
296 if (!(i >= ddbox->npbcdim && nc[i] <= 2) && bt[i] < nc[i]*limit)
304 /* The following choices should match those
305 * in init_domain_decomposition in domdec.c.
307 if (nc[XX] == 1 && nc[YY] > 1)
312 else if (nc[YY] == 1)
319 /* Will we use 1D or 2D PME decomposition? */
320 npme[XX] = (npme_tot % nc[XX] == 0) ? nc[XX] : npme_tot;
321 npme[YY] = npme_tot/npme[XX];
325 /* When two dimensions are (nearly) equal, use more cells
326 * for the smallest index, so the decomposition does not
327 * depend sensitively on the rounding of the box elements.
331 for(j=i+1; j<DIM; j++)
333 /* Check if the box size is nearly identical,
334 * in that case we prefer nx > ny and ny > nz.
336 if (fabs(bt[j] - bt[i]) < 0.01*bt[i] && nc[j] > nc[i])
338 /* The XX/YY check is a bit compact. If nc[YY]==npme[YY]
339 * this means the swapped nc has nc[XX]==npme[XX],
340 * and we can also swap X and Y for PME.
342 /* Check if dimension i and j are equivalent for PME.
343 * For x/y: if nc[YY]!=npme[YY], we can not swap x/y
344 * For y/z: we can not have PME decomposition in z
347 !((i == XX && j == YY && nc[YY] != npme[YY]) ||
348 (i == YY && j == ZZ && npme[YY] > 1)))
356 /* This function determines only half of the communication cost.
357 * All PP, PME and PP-PME communication is symmetric
358 * and the "back"-communication cost is identical to the forward cost.
361 comm_vol = comm_box_frac(nc,cutoff,ddbox);
366 /* Determine the largest volume for PME x/f redistribution */
367 if (nc[i] % npme[i] != 0)
371 comm_vol_xf = (npme[i]==2 ? 1.0/3.0 : 0.5);
375 comm_vol_xf = 1.0 - lcd(nc[i],npme[i])/(double)npme[i];
377 comm_pme += 3*natoms*comm_vol_xf;
380 /* Grid overlap communication */
383 nk = (i==0 ? ir->nkx : ir->nky);
384 overlap = (nk % npme[i] == 0 ? ir->pme_order-1 : ir->pme_order);
385 comm_pme += npme[i]*overlap*ir->nkx*ir->nky*ir->nkz/nk;
389 /* PME FFT communication volume.
390 * This only takes the communication into account and not imbalance
391 * in the calculation. But the imbalance in communication and calculation
392 * are similar and therefore these formulas also prefer load balance
393 * in the FFT and pme_solve calculation.
395 comm_pme += (npme[YY] - 1)*npme[YY]*div_up(ir->nky,npme[YY])*div_up(ir->nkz,npme[YY])*ir->nkx;
396 comm_pme += (npme[XX] - 1)*npme[XX]*div_up(ir->nkx,npme[XX])*div_up(ir->nky,npme[XX])*ir->nkz;
398 /* Add cost of pbc_dx for bondeds */
400 if ((nc[XX] == 1 || nc[YY] == 1) || (nc[ZZ] == 1 && ir->ePBC != epbcXY))
402 if ((ddbox->tric_dir[XX] && nc[XX] == 1) ||
403 (ddbox->tric_dir[YY] && nc[YY] == 1))
405 cost_pbcdx = pbcdxr*pbcdx_tric_fac;
409 cost_pbcdx = pbcdxr*pbcdx_rect_fac;
416 "nc %2d %2d %2d %2d %2d vol pp %6.4f pbcdx %6.4f pme %9.3e tot %9.3e\n",
417 nc[XX],nc[YY],nc[ZZ],npme[XX],npme[YY],
418 comm_vol,cost_pbcdx,comm_pme,
419 3*natoms*(comm_vol + cost_pbcdx) + comm_pme);
422 return 3*natoms*(comm_vol + cost_pbcdx) + comm_pme;
425 static void assign_factors(gmx_domdec_t *dd,
426 real limit,real cutoff,
427 matrix box,gmx_ddbox_t *ddbox,
428 int natoms,t_inputrec *ir,
429 float pbcdxr,int npme,
430 int ndiv,int *div,int *mdiv,ivec ir_try,ivec opt)
437 ce = comm_cost_est(dd,limit,cutoff,box,ddbox,
438 natoms,ir,pbcdxr,npme,ir_try);
439 if (ce >= 0 && (opt[XX] == 0 ||
440 ce < comm_cost_est(dd,limit,cutoff,box,ddbox,
444 copy_ivec(ir_try,opt);
450 for(x=mdiv[0]; x>=0; x--)
454 ir_try[XX] *= div[0];
456 for(y=mdiv[0]-x; y>=0; y--)
460 ir_try[YY] *= div[0];
462 for(i=0; i<mdiv[0]-x-y; i++)
464 ir_try[ZZ] *= div[0];
468 assign_factors(dd,limit,cutoff,box,ddbox,natoms,ir,pbcdxr,npme,
469 ndiv-1,div+1,mdiv+1,ir_try,opt);
471 for(i=0; i<mdiv[0]-x-y; i++)
473 ir_try[ZZ] /= div[0];
477 ir_try[YY] /= div[0];
482 ir_try[XX] /= div[0];
487 static real optimize_ncells(FILE *fplog,
488 int nnodes_tot,int npme_only,
489 bool bDynLoadBal,real dlb_scale,
490 gmx_mtop_t *mtop,matrix box,gmx_ddbox_t *ddbox,
493 real cellsize_limit,real cutoff,
494 bool bInterCGBondeds,bool bInterCGMultiBody,
497 int npp,npme,ndiv,*div,*mdiv,d,nmax;
503 limit = cellsize_limit;
509 npp = nnodes_tot - npme_only;
510 if (EEL_PME(ir->coulombtype))
512 npme = (npme_only > 0 ? npme_only : npp);
521 /* For Ewald exclusions pbc_dx is not called */
523 (IR_EXCL_FORCES(*ir) && !EEL_FULL(ir->coulombtype));
524 pbcdxr = (double)n_bonded_dx(mtop,bExcl_pbcdx)/(double)mtop->natoms;
528 /* Every molecule is a single charge group: no pbc required */
531 /* Add a margin for DLB and/or pressure scaling */
534 if (dlb_scale >= 1.0)
536 gmx_fatal(FARGS,"The value for option -dds should be smaller than 1");
540 fprintf(fplog,"Scaling the initial minimum size with 1/%g (option -dds) = %g\n",dlb_scale,1/dlb_scale);
544 else if (ir->epc != epcNO)
548 fprintf(fplog,"To account for pressure scaling, scaling the initial minimum size with %g\n",DD_GRID_MARGIN_PRES_SCALE);
549 limit *= DD_GRID_MARGIN_PRES_SCALE;
555 fprintf(fplog,"Optimizing the DD grid for %d cells with a minimum initial size of %.3f nm\n",npp,limit);
557 if (inhomogeneous_z(ir))
559 fprintf(fplog,"Ewald_geometry=%s: assuming inhomogeneous particle distribution in z, will not decompose in z.\n",eewg_names[ir->ewald_geometry]);
564 fprintf(fplog,"The maximum allowed number of cells is:");
567 nmax = (int)(ddbox->box_size[d]*ddbox->skew_fac[d]/limit);
568 if (d >= ddbox->npbcdim && nmax < 2)
572 if (d == ZZ && inhomogeneous_z(ir))
576 fprintf(fplog," %c %d",'X' + d,nmax);
584 fprintf(debug,"Average nr of pbc_dx calls per atom %.2f\n",pbcdxr);
587 /* Decompose npp in factors */
588 ndiv = factorize(npp,&div,&mdiv);
594 assign_factors(dd,limit,cutoff,box,ddbox,mtop->natoms,ir,pbcdxr,
595 npme,ndiv,div,mdiv,itry,nc);
603 real dd_choose_grid(FILE *fplog,
604 t_commrec *cr,gmx_domdec_t *dd,t_inputrec *ir,
605 gmx_mtop_t *mtop,matrix box,gmx_ddbox_t *ddbox,
606 bool bDynLoadBal,real dlb_scale,
607 real cellsize_limit,real cutoff_dd,
608 bool bInterCGBondeds,bool bInterCGMultiBody)
615 if (EEL_PME(ir->coulombtype))
617 if (cr->npmenodes >= 0)
619 if (cr->nnodes <= 2 && cr->npmenodes > 0)
622 "Can not have separate PME nodes with 2 or less nodes");
627 if (cr->nnodes <= 10)
633 cr->npmenodes = guess_npme(fplog,mtop,ir,box,cr->nnodes);
638 fprintf(fplog,"Using %d separate PME nodes\n",cr->npmenodes);
643 if (cr->npmenodes < 0)
649 limit = optimize_ncells(fplog,cr->nnodes,cr->npmenodes,
650 bDynLoadBal,dlb_scale,
651 mtop,box,ddbox,ir,dd,
652 cellsize_limit,cutoff_dd,
653 bInterCGBondeds,bInterCGMultiBody,
660 /* Communicate the information set by the master to all nodes */
661 gmx_bcast(sizeof(dd->nc),dd->nc,cr);
662 if (EEL_PME(ir->coulombtype))
664 gmx_bcast(sizeof(ir->nkx),&ir->nkx,cr);
665 gmx_bcast(sizeof(ir->nky),&ir->nky,cr);
666 gmx_bcast(sizeof(cr->npmenodes),&cr->npmenodes,cr);