src/tools/gmx_tune_pme.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2008, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  *
  31  * And Hey:
  32  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  33  */
  34 #ifdef HAVE_CONFIG_H
  35 #include <config.h>
  36 #endif
  37
  38
  39 #include <time.h>
  40 #ifdef HAVE_SYS_TIME_H
  41 #include <sys/time.h>
  42 #endif
  43
  44
  45
  46 #include "statutil.h"
  47 #include "typedefs.h"
  48 #include "smalloc.h"
  49 #include "vec.h"
  50 #include "copyrite.h"
  51 #include "statutil.h"
  52 #include "tpxio.h"
  53 #include "string2.h"
  54 #include "readinp.h"
  55 #include "calcgrid.h"
  56 #include "checkpoint.h"
  57 #include "gmx_ana.h"
  58 #include "names.h"
  59 #include "perf_est.h"
  60
  61
  62
  63 /* Enum for situations that can occur during log file parsing, the
  64  * corresponding string entries can be found in do_the_tests() in
  65  * const char* ParseLog[] */
  66 enum {
  67     eParselogOK,
  68     eParselogNotFound,
  69     eParselogNoPerfData,
  70     eParselogTerm,
  71     eParselogResetProblem,
  72     eParselogNoDDGrid,
  73     eParselogTPXVersion,
  74     eParselogNotParallel,
  75     eParselogFatal,
  76     eParselogNr
  77 };
  78
  79
  80 typedef struct
  81 {
  82     int  nPMEnodes;       /* number of PME-only nodes used in this test */
  83     int  nx, ny, nz;      /* DD grid */
  84     int  guessPME;        /* if nPMEnodes == -1, this is the guessed number of PME nodes */
  85     double *Gcycles;      /* This can contain more than one value if doing multiple tests */
  86     double Gcycles_Av;
  87     float *ns_per_day;
  88     float ns_per_day_Av;
  89     float *PME_f_load;    /* PME mesh/force load average*/
  90     float PME_f_load_Av;  /* Average average ;) ... */
  91     char *mdrun_cmd_line; /* Mdrun command line used for this test */
  92 } t_perf;
  93
  94
  95 typedef struct
  96 {
  97     int  nr_inputfiles;         /* The number of tpr and mdp input files */
  98     gmx_large_int_t orig_sim_steps;  /* Number of steps to be done in the real simulation */
  99     real *rcoulomb;             /* The coulomb radii [0...nr_inputfiles] */
 100     real *rvdw;                 /* The vdW radii */
 101     real *rlist;                /* Neighbourlist cutoff radius */
 102     real *rlistlong;
 103     int  *nkx, *nky, *nkz;
 104     real *fsx, *fsy, *fsz;      /* Fourierspacing in x,y,z dimension */
 105 } t_inputinfo;
 106
 107
 108 static void sep_line(FILE *fp)
 109 {
 110     fprintf(fp, "\n------------------------------------------------------------\n");
 111 }
 112
 113
 114 /* Wrapper for system calls */
 115 static int gmx_system_call(char *command)
 116 {
 117 #ifdef GMX_NO_SYSTEM
 118     gmx_fatal(FARGS,"No calls to system(3) supported on this platform. Attempted to call:\n'%s'\n",command);
 119 #else
 120     return ( system(command) );
 121 #endif
 122 }
 123
 124
 125 /* Check if string starts with substring */
 126 static gmx_bool str_starts(const char *string, const char *substring)
 127 {
 128     return ( strncmp(string, substring, strlen(substring)) == 0);
 129 }
 130
 131
 132 static void cleandata(t_perf *perfdata, int test_nr)
 133 {
 134     perfdata->Gcycles[test_nr]    = 0.0;
 135     perfdata->ns_per_day[test_nr] = 0.0;
 136     perfdata->PME_f_load[test_nr] = 0.0;
 137
 138     return;
 139 }
 140
 141
 142 static gmx_bool is_equal(real a, real b)
 143 {
 144     real diff, eps=1.0e-7;
 145
 146
 147     diff = a - b;
 148
 149     if (diff < 0.0) diff = -diff;
 150
 151     if (diff < eps)
 152         return TRUE;
 153     else
 154         return FALSE;
 155 }
 156
 157
 158 static void finalize(const char *fn_out)
 159 {
 160     char buf[STRLEN];
 161     FILE *fp;
 162
 163
 164     fp = fopen(fn_out,"r");
 165     fprintf(stdout,"\n\n");
 166
 167     while( fgets(buf,STRLEN-1,fp) != NULL )
 168     {
 169         fprintf(stdout,"%s",buf);
 170     }
 171     fclose(fp);
 172     fprintf(stdout,"\n\n");
 173 }
 174
 175
 176 enum {eFoundNothing, eFoundDDStr, eFoundAccountingStr, eFoundCycleStr};
 177
 178 static int parse_logfile(const char *logfile, const char *errfile,
 179         t_perf *perfdata, int test_nr, int presteps, gmx_large_int_t cpt_steps,
 180         int nnodes)
 181 {
 182     FILE  *fp;
 183     char  line[STRLEN], dumstring[STRLEN], dumstring2[STRLEN];
 184     const char matchstrdd[]="Domain decomposition grid";
 185     const char matchstrcr[]="resetting all time and cycle counters";
 186     const char matchstrbal[]="Average PME mesh/force load:";
 187     const char matchstring[]="R E A L   C Y C L E   A N D   T I M E   A C C O U N T I N G";
 188     const char errSIG[]="signal, stopping at the next";
 189     int   iFound;
 190     int   procs;
 191     float  dum1,dum2,dum3;
 192     int   npme;
 193     gmx_large_int_t resetsteps=-1;
 194     gmx_bool  bFoundResetStr = FALSE;
 195     gmx_bool  bResetChecked  = FALSE;
 196
 197
 198     if (!gmx_fexist(logfile))
 199     {
 200         fprintf(stderr, "WARNING: Could not find logfile %s.\n", logfile);
 201         cleandata(perfdata, test_nr);
 202         return eParselogNotFound;
 203     }
 204
 205     fp = fopen(logfile, "r");
 206     perfdata->PME_f_load[test_nr] = -1.0;
 207     perfdata->guessPME            = -1;
 208
 209     iFound = eFoundNothing;
 210     if (1 == nnodes)
 211         iFound = eFoundDDStr; /* Skip some case statements */
 212
 213     while (fgets(line, STRLEN, fp) != NULL)
 214     {
 215         /* Remove leading spaces */
 216         ltrim(line);
 217
 218         /* Check for TERM and INT signals from user: */
 219         if ( strstr(line, errSIG) != NULL )
 220         {
 221             fclose(fp);
 222             cleandata(perfdata, test_nr);
 223             return eParselogTerm;
 224         }
 225
 226         /* Check whether cycle resetting  worked */
 227         if (presteps > 0 && !bFoundResetStr)
 228         {
 229             if (strstr(line, matchstrcr) != NULL)
 230             {
 231                 sprintf(dumstring, "Step %s", gmx_large_int_pfmt);
 232                 sscanf(line, dumstring, &resetsteps);
 233                 bFoundResetStr = TRUE;
 234                 if (resetsteps == presteps+cpt_steps)
 235                 {
 236                     bResetChecked = TRUE;
 237                 }
 238                 else
 239                 {
 240                     sprintf(dumstring , gmx_large_int_pfmt, resetsteps);
 241                     sprintf(dumstring2, gmx_large_int_pfmt, presteps+cpt_steps);
 242                     fprintf(stderr, "WARNING: Time step counters were reset at step %s,\n"
 243                                     "         though they were supposed to be reset at step %s!\n",
 244                             dumstring, dumstring2);
 245                 }
 246             }
 247         }
 248
 249         /* Look for strings that appear in a certain order in the log file: */
 250         switch(iFound)
 251         {
 252             case eFoundNothing:
 253                 /* Look for domain decomp grid and separate PME nodes: */
 254                 if (str_starts(line, matchstrdd))
 255                 {
 256                     sscanf(line, "Domain decomposition grid %d x %d x %d, separate PME nodes %d",
 257                             &(perfdata->nx), &(perfdata->ny), &(perfdata->nz), &npme);
 258                     if (perfdata->nPMEnodes == -1)
 259                         perfdata->guessPME = npme;
 260                     else if (perfdata->nPMEnodes != npme)
 261                         gmx_fatal(FARGS, "PME nodes from command line and output file are not identical");
 262                     iFound = eFoundDDStr;
 263                 }
 264                 /* Catch a few errors that might have occured: */
 265                 else if (str_starts(line, "There is no domain decomposition for"))
 266                 {
 267                     fclose(fp);
 268                     return eParselogNoDDGrid;
 269                 }
 270                 else if (str_starts(line, "reading tpx file"))
 271                 {
 272                     fclose(fp);
 273                     return eParselogTPXVersion;
 274                 }
 275                 else if (str_starts(line, "The -dd or -npme option request a parallel simulation"))
 276                 {
 277                     fclose(fp);
 278                     return eParselogNotParallel;
 279                 }
 280                 break;
 281             case eFoundDDStr:
 282                 /* Look for PME mesh/force balance (not necessarily present, though) */
 283                 if (str_starts(line, matchstrbal))
 284                     sscanf(&line[strlen(matchstrbal)], "%f", &(perfdata->PME_f_load[test_nr]));
 285                 /* Look for matchstring */
 286                 if (str_starts(line, matchstring))
 287                     iFound = eFoundAccountingStr;
 288                 break;
 289             case eFoundAccountingStr:
 290                 /* Already found matchstring - look for cycle data */
 291                 if (str_starts(line, "Total  "))
 292                 {
 293                     sscanf(line,"Total %d %lf",&procs,&(perfdata->Gcycles[test_nr]));
 294                     iFound = eFoundCycleStr;
 295                 }
 296                 break;
 297             case eFoundCycleStr:
 298                 /* Already found cycle data - look for remaining performance info and return */
 299                 if (str_starts(line, "Performance:"))
 300                 {
 301                     sscanf(line,"%s %f %f %f %f", dumstring, &dum1, &dum2, &(perfdata->ns_per_day[test_nr]), &dum3);
 302                     fclose(fp);
 303                     if (bResetChecked || presteps == 0)
 304                         return eParselogOK;
 305                     else
 306                         return eParselogResetProblem;
 307                 }
 308                 break;
 309         }
 310     } /* while */
 311
 312     /* Close the log file */
 313     fclose(fp);
 314
 315     /* Check why there is no performance data in the log file.
 316      * Did a fatal errors occur? */
 317     if (gmx_fexist(errfile))
 318     {
 319         fp = fopen(errfile, "r");
 320         while (fgets(line, STRLEN, fp) != NULL)
 321         {
 322             if ( str_starts(line, "Fatal error:") )
 323             {
 324                 if (fgets(line, STRLEN, fp) != NULL)
 325                     fprintf(stderr, "\nWARNING: An error occured during this benchmark:\n"
 326                                     "%s\n", line);
 327                 fclose(fp);
 328                 cleandata(perfdata, test_nr);
 329                 return eParselogFatal;
 330             }
 331         }
 332         fclose(fp);
 333     }
 334     else
 335     {
 336         fprintf(stderr, "WARNING: Could not find stderr file %s.\n", errfile);
 337     }
 338
 339     /* Giving up ... we could not find out why there is no performance data in
 340      * the log file. */
 341     fprintf(stdout, "No performance data in log file.\n");
 342     cleandata(perfdata, test_nr);
 343
 344     return eParselogNoPerfData;
 345 }
 346
 347
 348 static gmx_bool analyze_data(
 349         FILE        *fp,
 350         const char  *fn,
 351         t_perf      **perfdata,
 352         int         nnodes,
 353         int         ntprs,
 354         int         ntests,
 355         int         nrepeats,
 356         t_inputinfo *info,
 357         int         *index_tpr,    /* OUT: Nr of mdp file with best settings */
 358         int         *npme_optimal) /* OUT: Optimal number of PME nodes */
 359 {
 360     int  i,j,k;
 361     int line=0, line_win=-1;
 362     int  k_win=-1, i_win=-1, winPME;
 363     double s=0.0;  /* standard deviation */
 364     t_perf *pd;
 365     char strbuf[STRLEN];
 366     char str_PME_f_load[13];
 367     gmx_bool bCanUseOrigTPR;
 368     gmx_bool bRefinedCoul, bRefinedVdW, bRefinedGrid;
 369
 370
 371     if (nrepeats > 1)
 372     {
 373         sep_line(fp);
 374         fprintf(fp, "Summary of successful runs:\n");
 375         fprintf(fp, "Line tpr PME nodes  Gcycles Av.     Std.dev.       ns/day        PME/f");
 376         if (nnodes > 1)
 377             fprintf(fp, "    DD grid");
 378         fprintf(fp, "\n");
 379     }
 380
 381
 382     for (k=0; k<ntprs; k++)
 383     {
 384         for (i=0; i<ntests; i++)
 385         {
 386             /* Select the right dataset: */
 387             pd = &(perfdata[k][i]);
 388
 389             pd->Gcycles_Av    = 0.0;
 390             pd->PME_f_load_Av = 0.0;
 391             pd->ns_per_day_Av = 0.0;
 392
 393             if (pd->nPMEnodes == -1)
 394                 sprintf(strbuf, "(%3d)", pd->guessPME);
 395             else
 396                 sprintf(strbuf, "     ");
 397
 398             /* Get the average run time of a setting */
 399             for (j=0; j<nrepeats; j++)
 400             {
 401                 pd->Gcycles_Av    += pd->Gcycles[j];
 402                 pd->PME_f_load_Av += pd->PME_f_load[j];
 403             }
 404             pd->Gcycles_Av    /= nrepeats;
 405             pd->PME_f_load_Av /= nrepeats;
 406
 407             for (j=0; j<nrepeats; j++)
 408             {
 409                 if (pd->ns_per_day[j] > 0.0)
 410                     pd->ns_per_day_Av += pd->ns_per_day[j];
 411                 else
 412                 {
 413                     /* Somehow the performance number was not aquired for this run,
 414                      * therefor set the average to some negative value: */
 415                     pd->ns_per_day_Av = -1.0f*nrepeats;
 416                     break;
 417                 }
 418             }
 419             pd->ns_per_day_Av /= nrepeats;
 420
 421             /* Nicer output: */
 422             if (pd->PME_f_load_Av > 0.0)
 423                 sprintf(str_PME_f_load, "%12.3f", pd->PME_f_load_Av);
 424             else
 425                 sprintf(str_PME_f_load, "%s", "         -  ");
 426
 427
 428             /* We assume we had a successful run if both averages are positive */
 429             if (pd->Gcycles_Av > 0.0 && pd->ns_per_day_Av > 0.0)
 430             {
 431                 /* Output statistics if repeats were done */
 432                 if (nrepeats > 1)
 433                 {
 434                     /* Calculate the standard deviation */
 435                     s = 0.0;
 436                     for (j=0; j<nrepeats; j++)
 437                         s += pow( pd->Gcycles[j] - pd->Gcycles_Av, 2 );
 438                     s /= (nrepeats - 1);
 439                     s = sqrt(s);
 440
 441                     fprintf(fp, "%4d %3d %4d%s %12.3f %12.3f %12.3f %s",
 442                             line, k, pd->nPMEnodes, strbuf, pd->Gcycles_Av, s,
 443                             pd->ns_per_day_Av, str_PME_f_load);
 444                     if (nnodes > 1)
 445                         fprintf(fp, "  %3d %3d %3d", pd->nx, pd->ny, pd->nz);
 446                     fprintf(fp, "\n");
 447                 }
 448                 /* Store the index of the best run found so far in 'winner': */
 449                 if ( (k_win == -1) || (pd->Gcycles_Av < perfdata[k_win][i_win].Gcycles_Av) )
 450                 {
 451                     k_win = k;
 452                     i_win = i;
 453                     line_win = line;
 454                 }
 455                 line++;
 456             }
 457         }
 458     }
 459
 460     if (k_win == -1)
 461         gmx_fatal(FARGS, "None of the runs was successful! Check %s for problems.", fn);
 462
 463     sep_line(fp);
 464
 465     winPME = perfdata[k_win][i_win].nPMEnodes;
 466
 467     if (1 == ntests)
 468     {
 469         /* We stuck to a fixed number of PME-only nodes */
 470         sprintf(strbuf, "settings No. %d", k_win);
 471     }
 472     else
 473     {
 474         /* We have optimized the number of PME-only nodes */
 475         if (winPME == -1)
 476             sprintf(strbuf, "%s", "the automatic number of PME nodes");
 477         else
 478             sprintf(strbuf, "%d PME nodes", winPME);
 479     }
 480     fprintf(fp, "Best performance was achieved with %s", strbuf);
 481     if ((nrepeats > 1) && (ntests > 1))
 482         fprintf(fp, " (see line %d)", line_win);
 483     fprintf(fp, "\n");
 484
 485     /* Only mention settings if they were modified: */
 486     bRefinedCoul = !is_equal(info->rcoulomb[k_win], info->rcoulomb[0]);
 487     bRefinedVdW  = !is_equal(info->rvdw[k_win]    , info->rvdw[0]    );
 488     bRefinedGrid = !(info->nkx[k_win] == info->nkx[0] &&
 489                      info->nky[k_win] == info->nky[0] &&
 490                      info->nkz[k_win] == info->nkz[0]);
 491
 492     if (bRefinedCoul || bRefinedVdW || bRefinedGrid)
 493     {
 494         fprintf(fp, "Optimized PME settings:\n");
 495         bCanUseOrigTPR = FALSE;
 496     }
 497     else
 498     {
 499         bCanUseOrigTPR = TRUE;
 500     }
 501
 502     if (bRefinedCoul)
 503         fprintf(fp, "   New Coulomb radius: %f nm (was %f nm)\n", info->rcoulomb[k_win], info->rcoulomb[0]);
 504
 505     if (bRefinedVdW)
 506         fprintf(fp, "   New Van der Waals radius: %f nm (was %f nm)\n", info->rvdw[k_win], info->rvdw[0]);
 507
 508     if (bRefinedGrid)
 509         fprintf(fp, "   New Fourier grid xyz: %d %d %d (was %d %d %d)\n", info->nkx[k_win], info->nky[k_win], info->nkz[k_win],
 510                                                                           info->nkx[0], info->nky[0], info->nkz[0]);
 511
 512     if (bCanUseOrigTPR && ntprs > 1)
 513         fprintf(fp, "and original PME settings.\n");
 514
 515     fflush(fp);
 516
 517     /* Return the index of the mdp file that showed the highest performance
 518      * and the optimal number of PME nodes */
 519     *index_tpr    = k_win;
 520     *npme_optimal = winPME;
 521
 522     return bCanUseOrigTPR;
 523 }
 524
 525
 526 /* Get the commands we need to set up the runs from environment variables */
 527 static void get_program_paths(gmx_bool bThreads, char *cmd_mpirun[], char cmd_np[],
 528                               char *cmd_mdrun[], int repeats)
 529 {
 530     char *command=NULL;
 531     char *cp;
 532     char *cp2;
 533     char line[STRLEN];
 534     FILE *fp;
 535     const char def_mpirun[] = "mpirun";
 536     const char def_mdrun[]  = "mdrun";
 537     const char filename[]   = "benchtest.log";
 538     const char match_mpi[]  = "NNODES=";
 539     const char match_mdrun[]= "Program: ";
 540     const char empty_mpirun[] = "";
 541     gmx_bool  bMdrun = FALSE;
 542     gmx_bool  bMPI   = FALSE;
 543
 544
 545     /* Get the commands we need to set up the runs from environment variables */
 546     if (!bThreads)
 547     {
 548         if ( (cp = getenv("MPIRUN")) != NULL)
 549             *cmd_mpirun = strdup(cp);
 550         else
 551             *cmd_mpirun = strdup(def_mpirun);
 552     }
 553     else
 554     {
 555         *cmd_mpirun = strdup(empty_mpirun);
 556     }
 557
 558     if ( (cp = getenv("MDRUN" )) != NULL )
 559         *cmd_mdrun  = strdup(cp);
 560     else
 561         *cmd_mdrun  = strdup(def_mdrun);
 562
 563
 564     /* If no simulations have to be performed, we are done here */
 565     if (repeats <= 0)
 566         return;
 567
 568     /* Run a small test to see whether mpirun + mdrun work  */
 569     fprintf(stdout, "Making sure that mdrun can be executed. ");
 570     if (bThreads)
 571     {
 572         snew(command, strlen(*cmd_mdrun) + strlen(cmd_np) + strlen(filename) + 50);
 573         sprintf(command, "%s%s-version -maxh 0.001 1> %s 2>&1", *cmd_mdrun, cmd_np, filename);
 574     }
 575     else
 576     {
 577         snew(command, strlen(*cmd_mpirun) + strlen(cmd_np) + strlen(*cmd_mdrun) + strlen(filename) + 50);
 578         sprintf(command, "%s%s%s -version -maxh 0.001 1> %s 2>&1", *cmd_mpirun, cmd_np, *cmd_mdrun, filename);
 579     }
 580     fprintf(stdout, "Trying '%s' ... ", command);
 581     make_backup(filename);
 582     gmx_system_call(command);
 583
 584     /* Check if we find the characteristic string in the output: */
 585     if (!gmx_fexist(filename))
 586         gmx_fatal(FARGS, "Output from test run could not be found.");
 587
 588     fp = fopen(filename, "r");
 589     /* We need to scan the whole output file, since sometimes the queuing system
 590      * also writes stuff to stdout/err */
 591     while ( !feof(fp) )
 592     {
 593         cp2=fgets(line, STRLEN, fp);
 594         if (cp2!=NULL)
 595         {
 596             if ( str_starts(line, match_mdrun) )
 597                 bMdrun = TRUE;
 598             if ( str_starts(line, match_mpi) )
 599                 bMPI = TRUE;
 600         }
 601     }
 602     fclose(fp);
 603
 604     if (bThreads)
 605     {
 606         if (bMPI)
 607         {
 608             gmx_fatal(FARGS, "Need a threaded version of mdrun. This one\n"
 609                     "(%s)\n"
 610                     "seems to have been compiled with MPI instead.",
 611                     *cmd_mdrun);
 612         }
 613     }
 614     else
 615     {
 616         if (bMdrun && !bMPI)
 617         {
 618             gmx_fatal(FARGS, "Need an MPI-enabled version of mdrun. This one\n"
 619                     "(%s)\n"
 620                     "seems to have been compiled without MPI support.",
 621                     *cmd_mdrun);
 622         }
 623     }
 624
 625     if (!bMdrun)
 626     {
 627         gmx_fatal(FARGS, "Cannot execute mdrun. Please check %s for problems!",
 628                 filename);
 629     }
 630
 631     fprintf(stdout, "passed.\n");
 632
 633     /* Clean up ... */
 634     remove(filename);
 635     sfree(command);
 636 }
 637
 638
 639 static void launch_simulation(
 640         gmx_bool bLaunch,       /* Should the simulation be launched? */
 641         FILE *fp,               /* General log file */
 642         gmx_bool bThreads,      /* whether to use threads */
 643         char *cmd_mpirun,       /* Command for mpirun */
 644         char *cmd_np,           /* Switch for -np or -nt or empty */
 645         char *cmd_mdrun,        /* Command for mdrun */
 646         char *args_for_mdrun,   /* Arguments for mdrun */
 647         const char *simulation_tpr,   /* This tpr will be simulated */
 648         int  nnodes,            /* Number of nodes to run on */
 649         int  nPMEnodes)         /* Number of PME nodes to use */
 650 {
 651     char  *command;
 652
 653
 654     /* Make enough space for the system call command,
 655      * (100 extra chars for -npme ... etc. options should suffice): */
 656     snew(command, strlen(cmd_mpirun)+strlen(cmd_mdrun)+strlen(cmd_np)+strlen(args_for_mdrun)+strlen(simulation_tpr)+100);
 657
 658     /* Note that the -passall options requires args_for_mdrun to be at the end
 659      * of the command line string */
 660     if (bThreads)
 661     {
 662         sprintf(command, "%s%s-npme %d -s %s %s",
 663                 cmd_mdrun, cmd_np, nPMEnodes, simulation_tpr, args_for_mdrun);
 664     }
 665     else
 666     {
 667         sprintf(command, "%s%s%s -npme %d -s %s %s",
 668                 cmd_mpirun, cmd_np, cmd_mdrun, nPMEnodes, simulation_tpr, args_for_mdrun);
 669     }
 670
 671     fprintf(fp, "%s this command line to launch the simulation:\n\n%s", bLaunch? "Using":"Please use", command);
 672     sep_line(fp);
 673     fflush(fp);
 674
 675     /* Now the real thing! */
 676     if (bLaunch)
 677     {
 678         fprintf(stdout, "\nLaunching simulation with best parameters now.\nExecuting '%s'", command);
 679         sep_line(stdout);
 680         fflush(stdout);
 681         gmx_system_call(command);
 682     }
 683 }
 684
 685
 686 static void modify_PMEsettings(
 687         gmx_large_int_t simsteps,  /* Set this value as number of time steps */
 688         const char *fn_best_tpr,   /* tpr file with the best performance */
 689         const char *fn_sim_tpr)    /* name of tpr file to be launched */
 690 {
 691     t_inputrec   *ir;
 692     t_state      state;
 693     gmx_mtop_t   mtop;
 694     char         buf[200];
 695
 696     snew(ir,1);
 697     read_tpx_state(fn_best_tpr,ir,&state,NULL,&mtop);
 698
 699     /* Set nsteps to the right value */
 700     ir->nsteps = simsteps;
 701
 702     /* Write the tpr file which will be launched */
 703     sprintf(buf, "Writing optimized simulation file %s with nsteps=%s.\n", fn_sim_tpr, gmx_large_int_pfmt);
 704     fprintf(stdout,buf,ir->nsteps);
 705     fflush(stdout);
 706     write_tpx_state(fn_sim_tpr,ir,&state,&mtop);
 707
 708     sfree(ir);
 709 }
 710
 711
 712 #define EPME_SWITCHED(e) ((e) == eelPMESWITCH || (e) == eelPMEUSERSWITCH)
 713
 714 /* Make additional TPR files with more computational load for the
 715  * direct space processors: */
 716 static void make_benchmark_tprs(
 717         const char *fn_sim_tpr,     /* READ : User-provided tpr file                 */
 718         char *fn_bench_tprs[],      /* WRITE: Names of benchmark tpr files           */
 719         gmx_large_int_t benchsteps, /* Number of time steps for benchmark runs       */
 720         gmx_large_int_t statesteps, /* Step counter in checkpoint file               */
 721         real rmin,                  /* Minimal Coulomb radius                        */
 722         real rmax,                  /* Maximal Coulomb radius                        */
 723         int *ntprs,                 /* No. of TPRs to write, each with a different
 724                                        rcoulomb and fourierspacing                   */
 725         t_inputinfo *info,          /* Contains information about mdp file options   */
 726         FILE *fp)                   /* Write the output here                         */
 727 {
 728     int          i,j,d;
 729     t_inputrec   *ir;
 730     t_state      state;
 731     gmx_mtop_t   mtop;
 732     real         nlist_buffer;      /* Thickness of the buffer regions for PME-switch potentials */
 733     char         buf[200];
 734     rvec         box_size;
 735     gmx_bool     bNote = FALSE;
 736     real         add;               /* Add this to rcoul for the next test    */
 737     real         fac = 1.0;         /* Scaling factor for Coulomb radius      */
 738     real         fourierspacing;    /* Basic fourierspacing from tpr          */
 739
 740
 741     sprintf(buf, "Making benchmark tpr file%s with %s time step%s",
 742             *ntprs > 1? "s":"", gmx_large_int_pfmt, benchsteps>1?"s":"");
 743     fprintf(stdout, buf, benchsteps);
 744     if (statesteps > 0)
 745     {
 746         sprintf(buf, " (adding %s steps from checkpoint file)", gmx_large_int_pfmt);
 747         fprintf(stdout, buf, statesteps);
 748         benchsteps += statesteps;
 749     }
 750     fprintf(stdout, ".\n");
 751
 752
 753     snew(ir,1);
 754     read_tpx_state(fn_sim_tpr,ir,&state,NULL,&mtop);
 755
 756     /* Check if some kind of PME was chosen */
 757     if (EEL_PME(ir->coulombtype) == FALSE)
 758         gmx_fatal(FARGS, "Can only do optimizations for simulations with %s electrostatics.",
 759                 EELTYPE(eelPME));
 760
 761     /* Check if rcoulomb == rlist, which is necessary for plain PME. */
 762     if (  (eelPME == ir->coulombtype) && !(ir->rcoulomb == ir->rlist) )
 763     {
 764         gmx_fatal(FARGS, "%s requires rcoulomb (%f) to be equal to rlist (%f).",
 765                 EELTYPE(eelPME), ir->rcoulomb, ir->rlist);
 766     }
 767     /* For other PME types, rcoulomb is allowed to be smaller than rlist */
 768     else if (ir->rcoulomb > ir->rlist)
 769     {
 770         gmx_fatal(FARGS, "%s requires rcoulomb (%f) to be equal to or smaller than rlist (%f)",
 771                 EELTYPE(ir->coulombtype), ir->rcoulomb, ir->rlist);
 772     }
 773
 774     /* Reduce the number of steps for the benchmarks */
 775     info->orig_sim_steps = ir->nsteps;
 776     ir->nsteps           = benchsteps;
 777
 778     /* For PME-switch potentials, keep the radial distance of the buffer region */
 779     nlist_buffer   = ir->rlist - ir->rcoulomb;
 780
 781     /* Determine length of triclinic box vectors */
 782     for(d=0; d<DIM; d++)
 783     {
 784         box_size[d] = 0;
 785         for(i=0;i<DIM;i++)
 786             box_size[d] += state.box[d][i]*state.box[d][i];
 787         box_size[d] = sqrt(box_size[d]);
 788     }
 789
 790     /* Reconstruct fourierspacing per dimension from the number of grid points and box size */
 791     info->fsx[0] = box_size[XX]/ir->nkx;
 792     info->fsy[0] = box_size[YY]/ir->nky;
 793     info->fsz[0] = box_size[ZZ]/ir->nkz;
 794
 795     /* Reconstruct the fourierspacing from the number of PME grid points found in the tpr */
 796     fourierspacing = max(box_size[ZZ]/ir->nkz, max(box_size[XX]/ir->nkx, box_size[YY]/ir->nky));
 797
 798     fprintf(stdout, "Calculating PME grid points on the basis of a fourierspacing of %f nm\n", fourierspacing);
 799
 800     /* For performance comparisons the number of particles is useful to have */
 801     fprintf(fp, "   Number of particles  : %d\n", mtop.natoms);
 802
 803     /* Print information about settings of which some are potentially modified: */
 804     fprintf(fp, "   Coulomb type         : %s\n", EELTYPE(ir->coulombtype));
 805     fprintf(fp, "   Grid spacing x y z   : %f %f %f\n",
 806             box_size[XX]/ir->nkx, box_size[YY]/ir->nky, box_size[ZZ]/ir->nkz);
 807     fprintf(fp, "   Van der Waals type   : %s\n", EVDWTYPE(ir->vdwtype));
 808     if (EVDW_SWITCHED(ir->vdwtype))
 809         fprintf(fp, "   rvdw_switch          : %f nm\n", ir->rvdw_switch);
 810     if (EPME_SWITCHED(ir->coulombtype))
 811         fprintf(fp, "   rlist                : %f nm\n", ir->rlist);
 812     if (ir->rlistlong != max_cutoff(ir->rvdw,ir->rcoulomb))
 813         fprintf(fp, "   rlistlong            : %f nm\n", ir->rlistlong);
 814
 815     /* Print a descriptive line about the tpr settings tested */
 816     fprintf(fp, "\nWill try these real/reciprocal workload settings:\n");
 817     fprintf(fp, " No.   scaling  rcoulomb");
 818     fprintf(fp, "  nkx  nky  nkz");
 819     fprintf(fp, "   spacing");
 820     if (evdwCUT == ir->vdwtype)
 821         fprintf(fp, "      rvdw");
 822     if (EPME_SWITCHED(ir->coulombtype))
 823         fprintf(fp, "     rlist");
 824     if ( ir->rlistlong != max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb)) )
 825         fprintf(fp, " rlistlong");
 826     fprintf(fp, "  tpr file\n");
 827
 828     /* Loop to create the requested number of tpr input files */
 829     for (j = 0; j < *ntprs; j++)
 830     {
 831         /* The first .tpr is the provided one, just need to modify nsteps,
 832          * so skip the following block */
 833         if (j != 0)
 834         {
 835             /* Determine which Coulomb radii rc to use in the benchmarks */
 836             add = (rmax-rmin)/(*ntprs-1);
 837             if (is_equal(rmin,info->rcoulomb[0]))
 838             {
 839                 ir->rcoulomb = rmin + j*add;
 840             }
 841             else if (is_equal(rmax,info->rcoulomb[0]))
 842             {
 843                 ir->rcoulomb = rmin + (j-1)*add;
 844             }
 845             else
 846             {
 847                 /* rmin != rcoul != rmax, ergo test between rmin and rmax */
 848                 add = (rmax-rmin)/(*ntprs-2);
 849                 ir->rcoulomb = rmin + (j-1)*add;
 850             }
 851
 852             /* Determine the scaling factor fac */
 853             fac = ir->rcoulomb/info->rcoulomb[0];
 854
 855             /* Scale the Fourier grid spacing */
 856             ir->nkx = ir->nky = ir->nkz = 0;
 857             calc_grid(NULL,state.box,fourierspacing*fac,&ir->nkx,&ir->nky,&ir->nkz);
 858
 859             /* Adjust other radii since various conditions neet to be fulfilled */
 860             if (eelPME == ir->coulombtype)
 861             {
 862                 /* plain PME, rcoulomb must be equal to rlist */
 863                 ir->rlist = ir->rcoulomb;
 864             }
 865             else
 866             {
 867                 /* rlist must be >= rcoulomb, we keep the size of the buffer region */
 868                 ir->rlist = ir->rcoulomb + nlist_buffer;
 869             }
 870
 871             if (evdwCUT == ir->vdwtype)
 872             {
 873                 /* For vdw cutoff, rvdw >= rlist */
 874                 ir->rvdw = max(info->rvdw[0], ir->rlist);
 875             }
 876
 877             ir->rlistlong = max_cutoff(ir->rlist,max_cutoff(ir->rvdw,ir->rcoulomb));
 878
 879         } /* end of "if (j != 0)" */
 880
 881         /* for j==0: Save the original settings
 882          * for j >0: Save modified radii and Fourier grids */
 883         info->rcoulomb[j]  = ir->rcoulomb;
 884         info->rvdw[j]      = ir->rvdw;
 885         info->nkx[j]       = ir->nkx;
 886         info->nky[j]       = ir->nky;
 887         info->nkz[j]       = ir->nkz;
 888         info->rlist[j]     = ir->rlist;
 889         info->rlistlong[j] = ir->rlistlong;
 890         info->fsx[j]       = fac*fourierspacing;
 891         info->fsy[j]       = fac*fourierspacing;
 892         info->fsz[j]       = fac*fourierspacing;
 893
 894         /* Write the benchmark tpr file */
 895         strncpy(fn_bench_tprs[j],fn_sim_tpr,strlen(fn_sim_tpr)-strlen(".tpr"));
 896         sprintf(buf, "_bench%.2d.tpr", j);
 897         strcat(fn_bench_tprs[j], buf);
 898         fprintf(stdout,"Writing benchmark tpr %s with nsteps=", fn_bench_tprs[j]);
 899         fprintf(stdout, gmx_large_int_pfmt, ir->nsteps);
 900         if (j > 0)
 901             fprintf(stdout,", scaling factor %f\n", fac);
 902         else
 903             fprintf(stdout,", unmodified settings\n");
 904
 905         write_tpx_state(fn_bench_tprs[j],ir,&state,&mtop);
 906
 907         /* Write information about modified tpr settings to log file */
 908         fprintf(fp, "%4d%10f%10f", j, fac, ir->rcoulomb);
 909         fprintf(fp, "%5d%5d%5d", ir->nkx, ir->nky, ir->nkz);
 910         fprintf(fp, " %9f ", info->fsx[j]);
 911         if (evdwCUT == ir->vdwtype)
 912             fprintf(fp, "%10f", ir->rvdw);
 913         if (EPME_SWITCHED(ir->coulombtype))
 914             fprintf(fp, "%10f", ir->rlist);
 915         if ( info->rlistlong[0] != max_cutoff(info->rlist[0],max_cutoff(info->rvdw[0],info->rcoulomb[0])) )
 916             fprintf(fp, "%10f", ir->rlistlong);
 917         fprintf(fp, "  %-14s\n",fn_bench_tprs[j]);
 918
 919         /* Make it clear to the user that some additional settings were modified */
 920         if (   !is_equal(ir->rvdw     , info->rvdw[0])
 921             || !is_equal(ir->rlistlong, info->rlistlong[0]) )
 922         {
 923             bNote = TRUE;
 924         }
 925     }
 926     if (bNote)
 927         fprintf(fp, "\nNote that in addition to the Coulomb radius and the Fourier grid\n"
 928                     "other input settings were also changed (see table above).\n"
 929                     "Please check if the modified settings are appropriate.\n");
 930     fflush(stdout);
 931     fflush(fp);
 932     sfree(ir);
 933 }
 934
 935
 936 /* Rename the files we want to keep to some meaningful filename and
 937  * delete the rest */
 938 static void cleanup(const t_filenm *fnm, int nfile, int k, int nnodes,
 939                     int nPMEnodes, int nr, gmx_bool bKeepStderr)
 940 {
 941     char numstring[STRLEN];
 942     char newfilename[STRLEN];
 943     const char *fn=NULL;
 944     int i;
 945     const char *opt;
 946
 947
 948     fprintf(stdout, "Cleaning up, deleting benchmark temp files ...\n");
 949
 950     for (i=0; i<nfile; i++)
 951     {
 952         opt = (char *)fnm[i].opt;
 953         if ( strcmp(opt, "-p") == 0 )
 954         {
 955             /* do nothing; keep this file */
 956             ;
 957         }
 958         else if (strcmp(opt, "-bg") == 0)
 959         {
 960             /* Give the log file a nice name so one can later see which parameters were used */
 961             numstring[0] = '\0';
 962             if (nr > 0)
 963                 sprintf(numstring, "_%d", nr);
 964             sprintf(newfilename, "%s_no%d_np%d_npme%d%s", opt2fn("-bg",nfile,fnm), k, nnodes, nPMEnodes, numstring);
 965             if (gmx_fexist(opt2fn("-bg",nfile,fnm)))
 966             {
 967                 fprintf(stdout, "renaming log file to %s\n", newfilename);
 968                 make_backup(newfilename);
 969                 rename(opt2fn("-bg",nfile,fnm), newfilename);
 970             }
 971         }
 972         else if (strcmp(opt, "-err") == 0)
 973         {
 974             /* This file contains the output of stderr. We want to keep it in
 975              * cases where there have been problems. */
 976             fn = opt2fn(opt, nfile, fnm);
 977             numstring[0] = '\0';
 978             if (nr > 0)
 979                 sprintf(numstring, "_%d", nr);
 980             sprintf(newfilename, "%s_no%d_np%d_npme%d%s", fn, k, nnodes, nPMEnodes, numstring);
 981             if (gmx_fexist(fn))
 982             {
 983                 if (bKeepStderr)
 984                 {
 985                     fprintf(stdout, "Saving stderr output in %s\n", newfilename);
 986                     make_backup(newfilename);
 987                     rename(fn, newfilename);
 988                 }
 989                 else
 990                 {
 991                     fprintf(stdout, "Deleting %s\n", fn);
 992                     remove(fn);
 993                 }
 994             }
 995         }
 996         /* Delete the files which are created for each benchmark run: (options -b*) */
 997         else if ( (0 == strncmp(opt, "-b", 2)) && (opt2bSet(opt,nfile,fnm) || !is_optional(&fnm[i])) )
 998         {
 999             fn = opt2fn(opt, nfile, fnm);
1000             if (gmx_fexist(fn))
1001             {
1002                 fprintf(stdout, "Deleting %s\n", fn);
1003                 remove(fn);
1004             }
1005         }
1006     }
1007 }
1008
1009
1010 /* Returns the largest common factor of n1 and n2 */
1011 static int largest_common_factor(int n1, int n2)
1012 {
1013     int factor, nmax;
1014
1015     nmax = min(n1, n2);
1016     for (factor=nmax; factor > 0; factor--)
1017     {
1018         if ( 0==(n1 % factor) && 0==(n2 % factor) )
1019         {
1020             return(factor);
1021         }
1022     }
1023     return 0; /* one for the compiler */
1024 }
1025
1026 enum {eNpmeAuto, eNpmeAll, eNpmeReduced, eNpmeSubset, eNpmeNr};
1027
1028 /* Create a list of numbers of PME nodes to test */
1029 static void make_npme_list(
1030         const char *npmevalues_opt,  /* Make a complete list with all
1031                            * possibilities or a short list that keeps only
1032                            * reasonable numbers of PME nodes                  */
1033         int *nentries,    /* Number of entries we put in the nPMEnodes list   */
1034         int *nPMEnodes[], /* Each entry contains the value for -npme          */
1035         int nnodes,       /* Total number of nodes to do the tests on         */
1036         int minPMEnodes,  /* Minimum number of PME nodes                      */
1037         int maxPMEnodes)  /* Maximum number of PME nodes                      */
1038 {
1039     int i,npme,npp;
1040     int min_factor=1;     /* We request that npp and npme have this minimal
1041                            * largest common factor (depends on npp)           */
1042     int nlistmax;         /* Max. list size                                   */
1043     int nlist;            /* Actual number of entries in list                 */
1044     int eNPME=0;
1045
1046
1047     /* Do we need to check all possible values for -npme or is a reduced list enough? */
1048     if ( 0 == strcmp(npmevalues_opt, "all") )
1049     {
1050         eNPME = eNpmeAll;
1051     }
1052     else if ( 0 == strcmp(npmevalues_opt, "subset") )
1053     {
1054         eNPME = eNpmeSubset;
1055     }
1056     else /* "auto" or "range" */
1057     {
1058         if (nnodes <= 64)
1059             eNPME = eNpmeAll;
1060         else if (nnodes < 128)
1061             eNPME = eNpmeReduced;
1062         else
1063             eNPME = eNpmeSubset;
1064     }
1065
1066     /* Calculate how many entries we could possibly have (in case of -npme all) */
1067     if (nnodes > 2)
1068     {
1069         nlistmax = maxPMEnodes - minPMEnodes + 3;
1070         if (0 == minPMEnodes)
1071             nlistmax--;
1072     }
1073     else
1074         nlistmax = 1;
1075
1076     /* Now make the actual list which is at most of size nlist */
1077     snew(*nPMEnodes, nlistmax);
1078     nlist = 0; /* start counting again, now the real entries in the list */
1079     for (i = 0; i < nlistmax - 2; i++)
1080     {
1081         npme = maxPMEnodes - i;
1082         npp  = nnodes-npme;
1083         switch (eNPME)
1084         {
1085             case eNpmeAll:
1086                 min_factor = 1;
1087                 break;
1088             case eNpmeReduced:
1089                 min_factor = 2;
1090                 break;
1091             case eNpmeSubset:
1092                 /* For 2d PME we want a common largest factor of at least the cube
1093                  * root of the number of PP nodes */
1094                 min_factor = (int) pow(npp, 1.0/3.0);
1095                 break;
1096             default:
1097                 gmx_fatal(FARGS, "Unknown option for eNPME in make_npme_list");
1098                 break;
1099         }
1100         if (largest_common_factor(npp, npme) >= min_factor)
1101         {
1102             (*nPMEnodes)[nlist] = npme;
1103             nlist++;
1104         }
1105     }
1106     /* We always test 0 PME nodes and the automatic number */
1107     *nentries = nlist + 2;
1108     (*nPMEnodes)[nlist  ] =  0;
1109     (*nPMEnodes)[nlist+1] = -1;
1110
1111     fprintf(stderr, "Will try the following %d different values for -npme:\n", *nentries);
1112     for (i=0; i<*nentries-1; i++)
1113         fprintf(stderr, "%d, ", (*nPMEnodes)[i]);
1114     fprintf(stderr, "and %d (auto).\n", (*nPMEnodes)[*nentries-1]);
1115 }
1116
1117
1118 /* Allocate memory to store the performance data */
1119 static void init_perfdata(t_perf *perfdata[], int ntprs, int datasets, int repeats)
1120 {
1121     int i, j, k;
1122
1123
1124     for (k=0; k<ntprs; k++)
1125     {
1126         snew(perfdata[k], datasets);
1127         for (i=0; i<datasets; i++)
1128         {
1129             for (j=0; j<repeats; j++)
1130             {
1131                 snew(perfdata[k][i].Gcycles   , repeats);
1132                 snew(perfdata[k][i].ns_per_day, repeats);
1133                 snew(perfdata[k][i].PME_f_load, repeats);
1134             }
1135         }
1136     }
1137 }
1138
1139
1140 /* Check for errors on mdrun -h */
1141 static void make_sure_it_runs(char *mdrun_cmd_line, int length, FILE *fp)
1142 {
1143     char *command, *msg;
1144     int  ret;
1145
1146
1147     snew(command, length +  15);
1148     snew(msg    , length + 500);
1149
1150     fprintf(stdout, "Making shure the benchmarks can be executed ...\n");
1151     sprintf(command, "%s-h -quiet", mdrun_cmd_line);
1152     ret = gmx_system_call(command);
1153
1154     if (0 != ret)
1155     {
1156         /* To prevent confusion, do not again issue a gmx_fatal here since we already
1157          * get the error message from mdrun itself */
1158         sprintf(msg,    "Cannot run the benchmark simulations! Please check the error message of\n"
1159                         "mdrun for the source of the problem. Did you provide a command line\n"
1160                         "argument that neither g_tune_pme nor mdrun understands? Offending command:\n"
1161                         "\n%s\n\n", command);
1162
1163         fprintf(stderr, "%s", msg);
1164         sep_line(fp);
1165         fprintf(fp    , "%s", msg);
1166
1167         exit(ret);
1168     }
1169
1170     sfree(command);
1171     sfree(msg    );
1172 }
1173
1174
1175 static void do_the_tests(
1176         FILE *fp,                   /* General g_tune_pme output file         */
1177         char **tpr_names,           /* Filenames of the input files to test   */
1178         int maxPMEnodes,            /* Max fraction of nodes to use for PME   */
1179         int minPMEnodes,            /* Min fraction of nodes to use for PME   */
1180         int npme_fixed,             /* If >= -1, test fixed number of PME
1181                                      * nodes only                             */
1182         const char *npmevalues_opt, /* Which -npme values should be tested    */
1183         t_perf **perfdata,          /* Here the performace data is stored     */
1184         int *pmeentries,            /* Entries in the nPMEnodes list          */
1185         int repeats,                /* Repeat each test this often            */
1186         int nnodes,                 /* Total number of nodes = nPP + nPME     */
1187         int nr_tprs,                /* Total number of tpr files to test      */
1188         gmx_bool bThreads,          /* Threads or MPI?                        */
1189         char *cmd_mpirun,           /* mpirun command string                  */
1190         char *cmd_np,               /* "-np", "-n", whatever mpirun needs     */
1191         char *cmd_mdrun,            /* mdrun command string                   */
1192         char *cmd_args_bench,       /* arguments for mdrun in a string        */
1193         const t_filenm *fnm,        /* List of filenames from command line    */
1194         int nfile,                  /* Number of files specified on the cmdl. */
1195         int sim_part,               /* For checkpointing                      */
1196         int presteps,               /* DLB equilibration steps, is checked    */
1197         gmx_large_int_t cpt_steps)  /* Time step counter in the checkpoint    */
1198 {
1199     int     i,nr,k,ret,count=0,totaltests;
1200     int     *nPMEnodes=NULL;
1201     t_perf  *pd=NULL;
1202     int     cmdline_length;
1203     char    *command, *cmd_stub;
1204     char    buf[STRLEN];
1205     gmx_bool bResetProblem=FALSE;
1206     gmx_bool bFirst=TRUE;
1207
1208
1209     /* This string array corresponds to the eParselog enum type at the start
1210      * of this file */
1211     const char* ParseLog[] = {"OK.",
1212                               "Logfile not found!",
1213                               "No timings, logfile truncated?",
1214                               "Run was terminated.",
1215                               "Counters were not reset properly.",
1216                               "No DD grid found for these settings.",
1217                               "TPX version conflict!",
1218                               "mdrun was not started in parallel!",
1219                               "An error occured." };
1220     char    str_PME_f_load[13];
1221
1222
1223     /* Allocate space for the mdrun command line. 100 extra characters should
1224        be more than enough for the -npme etcetera arguments */
1225     cmdline_length =  strlen(cmd_mpirun)
1226                     + strlen(cmd_np)
1227                     + strlen(cmd_mdrun)
1228                     + strlen(cmd_args_bench)
1229                     + strlen(tpr_names[0]) + 100;
1230     snew(command , cmdline_length);
1231     snew(cmd_stub, cmdline_length);
1232
1233     /* Construct the part of the command line that stays the same for all tests: */
1234     if (bThreads)
1235     {
1236         sprintf(cmd_stub, "%s%s", cmd_mdrun, cmd_np);
1237     }
1238     else
1239     {
1240         sprintf(cmd_stub, "%s%s%s ", cmd_mpirun, cmd_np, cmd_mdrun);
1241     }
1242
1243     /* Create a list of numbers of PME nodes to test */
1244     if (npme_fixed < -1)
1245     {
1246         make_npme_list(npmevalues_opt, pmeentries, &nPMEnodes,
1247                 nnodes, minPMEnodes, maxPMEnodes);
1248     }
1249     else
1250     {
1251         *pmeentries  = 1;
1252         snew(nPMEnodes, 1);
1253         nPMEnodes[0] = npme_fixed;
1254         fprintf(stderr, "Will use a fixed number of %d PME-only nodes.\n", nPMEnodes[0]);
1255     }
1256
1257     if (0 == repeats)
1258     {
1259         fprintf(fp, "\nNo benchmarks done since number of repeats (-r) is 0.\n");
1260         ffclose(fp);
1261         finalize(opt2fn("-p", nfile, fnm));
1262         exit(0);
1263     }
1264
1265     /* Allocate one dataset for each tpr input file: */
1266     init_perfdata(perfdata, nr_tprs, *pmeentries, repeats);
1267
1268     /*****************************************/
1269     /* Main loop over all tpr files to test: */
1270     /*****************************************/
1271     totaltests = nr_tprs*(*pmeentries)*repeats;
1272     for (k=0; k<nr_tprs;k++)
1273     {
1274         fprintf(fp, "\nIndividual timings for input file %d (%s):\n", k, tpr_names[k]);
1275         fprintf(fp, "PME nodes      Gcycles       ns/day        PME/f    Remark\n");
1276         /* Loop over various numbers of PME nodes: */
1277         for (i = 0; i < *pmeentries; i++)
1278         {
1279             pd = &perfdata[k][i];
1280
1281             /* Loop over the repeats for each scenario: */
1282             for (nr = 0; nr < repeats; nr++)
1283             {
1284                 pd->nPMEnodes = nPMEnodes[i];
1285
1286                 /* Add -npme and -s to the command line and save it. Note that
1287                  * the -passall (if set) options requires cmd_args_bench to be
1288                  * at the end of the command line string */
1289                 snew(pd->mdrun_cmd_line, cmdline_length);
1290                 sprintf(pd->mdrun_cmd_line, "%s-npme %d -s %s %s",
1291                         cmd_stub, pd->nPMEnodes, tpr_names[k], cmd_args_bench);
1292
1293                 /* To prevent that all benchmarks fail due to a show-stopper argument
1294                  * on the mdrun command line, we make a quick check with mdrun -h first */
1295                 if (bFirst)
1296                     make_sure_it_runs(pd->mdrun_cmd_line, cmdline_length, fp);
1297                 bFirst = FALSE;
1298
1299                 /* Do a benchmark simulation: */
1300                 if (repeats > 1)
1301                     sprintf(buf, ", pass %d/%d", nr+1, repeats);
1302                 else
1303                     buf[0]='\0';
1304                 fprintf(stdout, "\n=== Progress %2.0f%%, tpr %d/%d, run %d/%d%s:\n",
1305                         (100.0*count)/totaltests,
1306                         k+1, nr_tprs, i+1, *pmeentries, buf);
1307                 make_backup(opt2fn("-err",nfile,fnm));
1308                 sprintf(command, "%s 1> /dev/null 2>%s", pd->mdrun_cmd_line, opt2fn("-err",nfile,fnm));
1309                 fprintf(stdout, "%s\n", pd->mdrun_cmd_line);
1310                 gmx_system_call(command);
1311
1312                 /* Collect the performance data from the log file; also check stderr
1313                  * for fatal errors */
1314                 ret = parse_logfile(opt2fn("-bg",nfile,fnm), opt2fn("-err",nfile,fnm),
1315                         pd, nr, presteps, cpt_steps, nnodes);
1316                 if ((presteps > 0) && (ret == eParselogResetProblem))
1317                     bResetProblem = TRUE;
1318
1319                 if (-1 == pd->nPMEnodes)
1320                     sprintf(buf, "(%3d)", pd->guessPME);
1321                 else
1322                     sprintf(buf, "     ");
1323
1324                 /* Nicer output */
1325                 if (pd->PME_f_load[nr] > 0.0)
1326                     sprintf(str_PME_f_load, "%12.3f", pd->PME_f_load[nr]);
1327                 else
1328                     sprintf(str_PME_f_load, "%s", "         -  ");
1329
1330                 /* Write the data we got to disk */
1331                 fprintf(fp, "%4d%s %12.3f %12.3f %s    %s", pd->nPMEnodes,
1332                         buf, pd->Gcycles[nr], pd->ns_per_day[nr], str_PME_f_load, ParseLog[ret]);
1333                 if (! (ret==eParselogOK || ret==eParselogNoDDGrid || ret==eParselogNotFound) )
1334                     fprintf(fp, " Check %s file for problems.", ret==eParselogFatal? "err":"log");
1335                 fprintf(fp, "\n");
1336                 fflush(fp);
1337                 count++;
1338
1339                 /* Do some cleaning up and delete the files we do not need any more */
1340                 cleanup(fnm, nfile, k, nnodes, pd->nPMEnodes, nr, ret==eParselogFatal);
1341
1342                 /* If the first run with this number of processors already failed, do not try again: */
1343                 if (pd->Gcycles[0] <= 0.0 && repeats > 1)
1344                 {
1345                     fprintf(stdout, "Skipping remaining passes of unsuccessful setting, see log file for details.\n");
1346                     count += repeats-(nr+1);
1347                     break;
1348                 }
1349             } /* end of repeats loop */
1350         } /* end of -npme loop */
1351     } /* end of tpr file loop */
1352
1353     if (bResetProblem)
1354     {
1355         sep_line(fp);
1356         fprintf(fp, "WARNING: The cycle and time step counters could not be reset\n"
1357                     "properly. The reason could be that mpirun did not manage to\n"
1358                     "export the environment variable GMX_RESET_COUNTER. You might\n"
1359                     "have to give a special switch to mpirun for that.\n"
1360                     "Alternatively, you can manually set GMX_RESET_COUNTER to the\n"
1361                     "value normally provided by -presteps.");
1362         sep_line(fp);
1363     }
1364     sfree(command);
1365     sfree(cmd_stub);
1366 }
1367
1368
1369 static void check_input(
1370         int nnodes,
1371         int repeats,
1372         int *ntprs,
1373         real *rmin,
1374         real rcoulomb,
1375         real *rmax,
1376         real maxPMEfraction,
1377         real minPMEfraction,
1378         int  npme_fixed,
1379         gmx_large_int_t bench_nsteps,
1380         const t_filenm *fnm,
1381         int nfile,
1382         int sim_part,
1383         int presteps,
1384         int npargs,
1385         t_pargs *pa)
1386 {
1387     int old;
1388
1389
1390     /* Make sure the input file exists */
1391     if (!gmx_fexist(opt2fn("-s",nfile,fnm)))
1392         gmx_fatal(FARGS, "File %s not found.", opt2fn("-s",nfile,fnm));
1393
1394     /* Make sure that the checkpoint file is not overwritten during benchmarking */
1395     if ( (0 == strcmp(opt2fn("-cpi",nfile,fnm), opt2fn("-bcpo",nfile,fnm)) ) && (sim_part > 1) )
1396         gmx_fatal(FARGS, "Checkpoint input (-cpi) and benchmark checkpoint output (-bcpo) files must not be identical.\n"
1397                          "The checkpoint input file must not be overwritten during the benchmarks.\n");
1398
1399     /* Make sure that repeats is >= 0 (if == 0, only write tpr files) */
1400     if (repeats < 0)
1401         gmx_fatal(FARGS, "Number of repeats < 0!");
1402
1403     /* Check number of nodes */
1404     if (nnodes < 1)
1405         gmx_fatal(FARGS, "Number of nodes/threads must be a positive integer.");
1406
1407     /* Automatically choose -ntpr if not set */
1408     if (*ntprs < 1)
1409     {
1410         if (nnodes < 16)
1411             *ntprs = 1;
1412         else
1413         {
1414             *ntprs = 3;
1415             /* Set a reasonable scaling factor for rcoulomb */
1416             if (*rmax <= 0)
1417                 *rmax = rcoulomb * 1.2;
1418         }
1419         fprintf(stderr, "Will test %d tpr file%s.\n", *ntprs, *ntprs==1?"":"s");
1420     }
1421     else
1422     {
1423         if (1 == *ntprs)
1424             fprintf(stderr, "Note: Choose ntpr>1 to shift PME load between real and reciprocal space.\n");
1425     }
1426
1427     /* Make shure that rmin <= rcoulomb <= rmax */
1428     if (*rmin <= 0) *rmin = rcoulomb;
1429     if (*rmax <= 0) *rmax = rcoulomb;
1430     if ( !(*rmin <= *rmax) )
1431     {
1432         gmx_fatal(FARGS, "Please choose the Coulomb radii such that rmin <= rmax.\n"
1433                          "rmin = %g, rmax = %g, actual rcoul from .tpr file = %g\n", *rmin, *rmax, rcoulomb);
1434     }
1435     /* Add test scenarios if rmin or rmax were set */
1436     if (*ntprs <= 2)
1437     {
1438         if ( !is_equal(*rmin, rcoulomb) && (*ntprs == 1) )
1439         {
1440             (*ntprs)++;
1441             fprintf(stderr, "NOTE: Setting -rmin to %g changed -ntpr to %d\n",
1442                     *rmin, *ntprs);
1443         }
1444         if (!is_equal(*rmax, rcoulomb) && (*ntprs == 1) )
1445         {
1446             (*ntprs)++;
1447             fprintf(stderr, "NOTE: Setting -rmax to %g changed -ntpr to %d\n",
1448                     *rmax, *ntprs);
1449         }
1450     }
1451     old = *ntprs;
1452     /* If one of rmin, rmax is set, we need 2 tpr files at minimum */
1453     if ( !is_equal(*rmax, rcoulomb) || !is_equal(*rmin, rcoulomb) )
1454         *ntprs = max(*ntprs, 2);
1455
1456     /* If both rmin, rmax are set, we need 3 tpr files at minimum */
1457     if ( !is_equal(*rmax, rcoulomb) && !is_equal(*rmin, rcoulomb) )
1458         *ntprs = max(*ntprs, 3);
1459
1460     if (old != *ntprs)
1461     {
1462         fprintf(stderr, "NOTE: Your rmin, rmax setting changed -ntpr to %d\n", *ntprs);
1463     }
1464
1465     if (*ntprs > 1)
1466     {
1467         if (is_equal(*rmin,rcoulomb) && is_equal(rcoulomb,*rmax)) /* We have just a single rc */
1468         {
1469             fprintf(stderr, "WARNING: Resetting -ntpr to 1 since no Coulomb radius scaling is requested.\n"
1470                             "Please set rmin < rmax to test Coulomb radii in the [rmin, rmax] interval\n"
1471                             "with correspondingly adjusted PME grid settings\n");
1472             *ntprs = 1;
1473         }
1474     }
1475
1476     /* Check whether max and min fraction are within required values */
1477     if (maxPMEfraction > 0.5 || maxPMEfraction < 0)
1478         gmx_fatal(FARGS, "-max must be between 0 and 0.5");
1479     if (minPMEfraction > 0.5 || minPMEfraction < 0)
1480         gmx_fatal(FARGS, "-min must be between 0 and 0.5");
1481     if (maxPMEfraction < minPMEfraction)
1482         gmx_fatal(FARGS, "-max must be larger or equal to -min");
1483
1484     /* Check whether the number of steps - if it was set - has a reasonable value */
1485     if (bench_nsteps < 0)
1486         gmx_fatal(FARGS, "Number of steps must be positive.");
1487
1488     if (bench_nsteps > 10000 || bench_nsteps < 100)
1489     {
1490         fprintf(stderr, "WARNING: steps=");
1491         fprintf(stderr, gmx_large_int_pfmt, bench_nsteps);
1492         fprintf(stderr, ". Are you sure you want to perform so %s steps for each benchmark?\n", (bench_nsteps < 100)? "few" : "many");
1493     }
1494
1495     if (presteps < 0)
1496     {
1497         gmx_fatal(FARGS, "Cannot have a negative number of presteps.\n");
1498     }
1499
1500     /* Check for rcoulomb scaling if more than one .tpr file is tested */
1501     if (*ntprs > 1)
1502     {
1503         if (*rmin/rcoulomb < 0.75 || *rmax/rcoulomb > 1.25)
1504             fprintf(stderr, "WARNING: Applying extreme scaling factor. I hope you know what you are doing.\n");
1505     }
1506
1507     /* If a fixed number of PME nodes is set we do rcoulomb and PME gird tuning
1508      * only. We need to check whether the requested number of PME-only nodes
1509      * makes sense. */
1510     if (npme_fixed > -1)
1511     {
1512         /* No more than 50% of all nodes can be assigned as PME-only nodes. */
1513         if (2*npme_fixed > nnodes)
1514         {
1515             gmx_fatal(FARGS, "Cannot have more than %d PME-only nodes for a total of %d nodes (you chose %d).\n",
1516                              nnodes/2, nnodes, npme_fixed);
1517         }
1518         if ((npme_fixed > 0) && (5*npme_fixed < nnodes))
1519         {
1520             fprintf(stderr, "WARNING: Only %g percent of the nodes are assigned as PME-only nodes.\n",
1521                              100.0*((real)npme_fixed / (real)nnodes));
1522         }
1523         if (opt2parg_bSet("-min",npargs,pa) || opt2parg_bSet("-max",npargs,pa))
1524         {
1525             fprintf(stderr, "NOTE: The -min, -max, and -npme options have no effect when a\n"
1526                             "      fixed number of PME-only nodes is requested with -fix.\n");
1527         }
1528     }
1529 }
1530
1531
1532 /* Returns TRUE when "opt" is needed at launch time */
1533 static gmx_bool is_launch_file(char *opt, gmx_bool bSet)
1534 {
1535     /* Apart from the input .tpr we need all options that were set
1536      * on the command line and that do not start with -b */
1537     if (0 == strncmp(opt,"-b", 2) || 0 == strncmp(opt,"-s", 2))
1538         return FALSE;
1539
1540     if (bSet)
1541         return TRUE;
1542     else
1543         return FALSE;
1544 }
1545
1546
1547 /* Returns TRUE when "opt" defines a file which is needed for the benchmarks runs */
1548 static gmx_bool is_bench_file(char *opt, gmx_bool bSet, gmx_bool bOptional, gmx_bool bIsOutput)
1549 {
1550     /* Apart from the input .tpr, all files starting with "-b" are for
1551      * _b_enchmark files exclusively */
1552     if (0 == strncmp(opt,"-s", 2)) return FALSE;
1553     if (0 == strncmp(opt,"-b", 2) || 0 == strncmp(opt,"-s", 2))
1554     {
1555         if (!bOptional || bSet)
1556             return TRUE;
1557         else
1558             return FALSE;
1559     }
1560     else
1561     {
1562         if (bIsOutput)
1563             return FALSE;
1564         else
1565             if (bSet) /* These are additional input files like -cpi -ei */
1566                 return TRUE;
1567             else
1568                 return FALSE;
1569     }
1570 }
1571
1572
1573 /* Adds 'buf' to 'str' */
1574 static void add_to_string(char **str, char *buf)
1575 {
1576     int len;
1577
1578
1579     len = strlen(*str) + strlen(buf) + 1;
1580     srenew(*str, len);
1581     strcat(*str, buf);
1582 }
1583
1584
1585 /* Create the command line for the benchmark as well as for the real run */
1586 static void create_command_line_snippets(
1587         gmx_bool bThreads,
1588         gmx_bool bAppendFiles,
1589         gmx_bool bKeepAndNumCPT,
1590         gmx_bool bResetHWay,
1591         int      presteps,
1592         int      nfile,
1593         t_filenm fnm[],
1594         int      npargs,
1595         t_pargs  *pa,
1596         const char *procstring,      /* How to pass the number of processors to $MPIRUN */
1597         char     *cmd_np[],          /* Actual command line snippet, e.g. '-np <N>' */
1598         char     *cmd_args_bench[],  /* command line arguments for benchmark runs */
1599         char     *cmd_args_launch[], /* command line arguments for simulation run */
1600         char     extra_args[])       /* Add this to the end of the command line */
1601 {
1602     int        i;
1603     char       *opt;
1604     const char *name;
1605     char       strbuf[STRLEN];
1606
1607
1608     /* strlen needs at least '\0' as a string: */
1609     snew(*cmd_args_bench ,1);
1610     snew(*cmd_args_launch,1);
1611     *cmd_args_launch[0]='\0';
1612     *cmd_args_bench[0] ='\0';
1613
1614
1615     /*******************************************/
1616     /* 1. Process other command line arguments */
1617     /*******************************************/
1618     if (presteps > 0)
1619     {
1620         /* Add equilibration steps to benchmark options */
1621         sprintf(strbuf, "-resetstep %d ", presteps);
1622         add_to_string(cmd_args_bench, strbuf);
1623     }
1624     /* These switches take effect only at launch time */
1625     if (FALSE == bAppendFiles)
1626     {
1627         add_to_string(cmd_args_launch, "-noappend ");
1628     }
1629     if (bKeepAndNumCPT)
1630     {
1631         add_to_string(cmd_args_launch, "-cpnum ");
1632     }
1633     if (bResetHWay)
1634     {
1635         add_to_string(cmd_args_launch, "-resethway ");
1636     }
1637
1638     /********************/
1639     /* 2. Process files */
1640     /********************/
1641     for (i=0; i<nfile; i++)
1642     {
1643         opt  = (char *)fnm[i].opt;
1644         name = opt2fn(opt,nfile,fnm);
1645
1646         /* Strbuf contains the options, now let's sort out where we need that */
1647         sprintf(strbuf, "%s %s ", opt, name);
1648
1649         if ( is_bench_file(opt, opt2bSet(opt,nfile,fnm), is_optional(&fnm[i]), is_output(&fnm[i])) )
1650         {
1651             /* All options starting with -b* need the 'b' removed,
1652              * therefore overwrite strbuf */
1653             if (0 == strncmp(opt, "-b", 2))
1654                 sprintf(strbuf, "-%s %s ", &opt[2], name);
1655
1656             add_to_string(cmd_args_bench, strbuf);
1657         }
1658
1659         if ( is_launch_file(opt,opt2bSet(opt,nfile,fnm)) )
1660             add_to_string(cmd_args_launch, strbuf);
1661     }
1662
1663     add_to_string(cmd_args_bench , extra_args);
1664     add_to_string(cmd_args_launch, extra_args);
1665 }
1666
1667
1668 /* Set option opt */
1669 static void setopt(const char *opt,int nfile,t_filenm fnm[])
1670 {
1671   int i;
1672
1673   for(i=0; (i<nfile); i++)
1674     if (strcmp(opt,fnm[i].opt)==0)
1675       fnm[i].flag |= ffSET;
1676 }
1677
1678
1679 /* This routine inspects the tpr file and ...
1680  * 1. checks for output files that get triggered by a tpr option. These output
1681  *    files are marked as 'set' to allow for a proper cleanup after each
1682  *    tuning run.
1683  * 2. returns the PME:PP load ratio
1684  * 3. returns rcoulomb from the tpr */
1685 static float inspect_tpr(int nfile, t_filenm fnm[], real *rcoulomb)
1686 {
1687     gmx_bool     bPull;     /* Is pulling requested in .tpr file?             */
1688     gmx_bool     bTpi;      /* Is test particle insertion requested?          */
1689     gmx_bool     bFree;     /* Is a free energy simulation requested?         */
1690     gmx_bool     bNM;       /* Is a normal mode analysis requested?           */
1691     t_inputrec   ir;
1692     t_state      state;
1693     gmx_mtop_t   mtop;
1694
1695
1696     /* Check tpr file for options that trigger extra output files */
1697     read_tpx_state(opt2fn("-s",nfile,fnm),&ir,&state,NULL,&mtop);
1698     bPull = (epullNO != ir.ePull);
1699     bFree = (efepNO  != ir.efep );
1700     bNM   = (eiNM    == ir.eI   );
1701     bTpi  = EI_TPI(ir.eI);
1702
1703     /* Set these output files on the tuning command-line */
1704     if (bPull)
1705     {
1706         setopt("-pf"  , nfile, fnm);
1707         setopt("-px"  , nfile, fnm);
1708     }
1709     if (bFree)
1710     {
1711         setopt("-dhdl", nfile, fnm);
1712     }
1713     if (bTpi)
1714     {
1715         setopt("-tpi" , nfile, fnm);
1716         setopt("-tpid", nfile, fnm);
1717     }
1718     if (bNM)
1719     {
1720         setopt("-mtx" , nfile, fnm);
1721     }
1722
1723     *rcoulomb = ir.rcoulomb;
1724
1725     /* Return the estimate for the number of PME nodes */
1726     return pme_load_estimate(&mtop,&ir,state.box);
1727 }
1728
1729
1730 static void couple_files_options(int nfile, t_filenm fnm[])
1731 {
1732     int i;
1733     gmx_bool bSet,bBench;
1734     char *opt;
1735     char buf[20];
1736
1737
1738     for (i=0; i<nfile; i++)
1739     {
1740         opt  = (char *)fnm[i].opt;
1741         bSet = ((fnm[i].flag & ffSET) != 0);
1742         bBench = (0 == strncmp(opt,"-b", 2));
1743
1744         /* Check optional files */
1745         /* If e.g. -eo is set, then -beo also needs to be set */
1746         if (is_optional(&fnm[i]) && bSet && !bBench)
1747         {
1748             sprintf(buf, "-b%s", &opt[1]);
1749             setopt(buf,nfile,fnm);
1750         }
1751         /* If -beo is set, then -eo also needs to be! */
1752         if (is_optional(&fnm[i]) && bSet && bBench)
1753         {
1754             sprintf(buf, "-%s", &opt[2]);
1755             setopt(buf,nfile,fnm);
1756         }
1757     }
1758 }
1759
1760
1761 static double gettime()
1762 {
1763 #ifdef HAVE_GETTIMEOFDAY
1764     struct timeval t;
1765     double seconds;
1766
1767     gettimeofday(&t,NULL);
1768
1769     seconds = (double) t.tv_sec + 1e-6*(double)t.tv_usec;
1770
1771     return seconds;
1772 #else
1773     double  seconds;
1774
1775     seconds = time(NULL);
1776
1777     return seconds;
1778 #endif
1779 }
1780
1781
1782 #define BENCHSTEPS (1000)
1783
1784 int gmx_tune_pme(int argc,char *argv[])
1785 {
1786     const char *desc[] = {
1787             "For a given number [TT]-np[tt] or [TT]-nt[tt] of processors/threads, this program systematically",
1788             "times [TT]mdrun[tt] with various numbers of PME-only nodes and determines",
1789             "which setting is fastest. It will also test whether performance can",
1790             "be enhanced by shifting load from the reciprocal to the real space",
1791             "part of the Ewald sum. ",
1792             "Simply pass your [TT].tpr[tt] file to [TT]g_tune_pme[tt] together with other options",
1793             "for [TT]mdrun[tt] as needed.[PAR]",
1794             "Which executables are used can be set in the environment variables",
1795             "MPIRUN and MDRUN. If these are not present, 'mpirun' and 'mdrun'",
1796             "will be used as defaults. Note that for certain MPI frameworks you",
1797             "need to provide a machine- or hostfile. This can also be passed",
1798             "via the MPIRUN variable, e.g.[PAR]",
1799             "[TT]export MPIRUN=\"/usr/local/mpirun -machinefile hosts\"[tt][PAR]",
1800             "Please call [TT]g_tune_pme[tt] with the normal options you would pass to",
1801             "[TT]mdrun[tt] and add [TT]-np[tt] for the number of processors to perform the",
1802             "tests on, or [TT]-nt[tt] for the number of threads. You can also add [TT]-r[tt]",
1803             "to repeat each test several times to get better statistics. [PAR]",
1804             "[TT]g_tune_pme[tt] can test various real space / reciprocal space workloads",
1805             "for you. With [TT]-ntpr[tt] you control how many extra [TT].tpr[tt] files will be",
1806             "written with enlarged cutoffs and smaller Fourier grids respectively.",
1807             "Typically, the first test (number 0) will be with the settings from the input",
1808             "[TT].tpr[tt] file; the last test (number [TT]ntpr[tt]) will have the Coulomb cutoff",
1809             "specified by [TT]-rmax[tt] with a somwhat smaller PME grid at the same time. ",
1810             "In this last test, the Fourier spacing is multiplied with [TT]rmax[tt]/rcoulomb. ",
1811             "The remaining [TT].tpr[tt] files will have equally-spaced Coulomb radii (and Fourier "
1812             "spacings) between these extremes. [BB]Note[bb] that you can set [TT]-ntpr[tt] to 1",
1813             "if you just seek the optimal number of PME-only nodes; in that case",
1814             "your input [TT].tpr[tt] file will remain unchanged.[PAR]",
1815             "For the benchmark runs, the default of 1000 time steps should suffice for most",
1816             "MD systems. The dynamic load balancing needs about 100 time steps",
1817             "to adapt to local load imbalances, therefore the time step counters",
1818             "are by default reset after 100 steps. For large systems (>1M atoms), as well as ",
1819             "for a higher accuarcy of the measurements, you should set [TT]-resetstep[tt] to a higher value.",
1820             "From the 'DD' load imbalance entries in the md.log output file you",
1821             "can tell after how many steps the load is sufficiently balanced. Example call:[PAR]"
1822             "[TT]g_tune_pme -np 64 -s protein.tpr -launch[tt][PAR]",
1823             "After calling [TT]mdrun[tt] several times, detailed performance information",
1824             "is available in the output file [TT]perf.out.[tt] ",
1825             "[BB]Note[bb] that during the benchmarks, a couple of temporary files are written",
1826             "(options [TT]-b[tt]*), these will be automatically deleted after each test.[PAR]",
1827             "If you want the simulation to be started automatically with the",
1828             "optimized parameters, use the command line option [TT]-launch[tt].[PAR]",
1829     };
1830
1831     int        nnodes =1;
1832     int        repeats=2;
1833     int        pmeentries=0; /* How many values for -npme do we actually test for each tpr file */
1834     real       maxPMEfraction=0.50;
1835     real       minPMEfraction=0.25;
1836     int        maxPMEnodes, minPMEnodes;
1837     float      guessPMEratio;   /* guessed PME:PP ratio based on the tpr file */
1838     float      guessPMEnodes;
1839     int        npme_fixed=-2;             /* If >= -1, use only this number
1840                                            * of PME-only nodes                */
1841     int        ntprs=0;
1842     real       rmin=0.0,rmax=0.0;  /* min and max value for rcoulomb if scaling is requested */
1843     real       rcoulomb=-1.0;             /* Coulomb radius as set in .tpr file */
1844     gmx_large_int_t bench_nsteps=BENCHSTEPS;
1845     gmx_large_int_t new_sim_nsteps=-1;   /* -1 indicates: not set by the user */
1846     gmx_large_int_t cpt_steps=0;         /* Step counter in .cpt input file   */
1847     int        presteps=100;    /* Do a full cycle reset after presteps steps */
1848     gmx_bool   bOverwrite=FALSE, bKeepTPR;
1849     gmx_bool   bLaunch=FALSE;
1850     char       *ExtraArgs=NULL;
1851     char       **tpr_names=NULL;
1852     const char *simulation_tpr=NULL;
1853     int        best_npme, best_tpr;
1854     int        sim_part = 1;     /* For benchmarks with checkpoint files */
1855     char       bbuf[STRLEN];
1856
1857     /* Default program names if nothing else is found */
1858     char        *cmd_mpirun=NULL, *cmd_mdrun=NULL;
1859     char        *cmd_args_bench, *cmd_args_launch;
1860     char        *cmd_np=NULL;
1861
1862     t_perf      **perfdata=NULL;
1863     t_inputinfo *info;
1864     int         i;
1865     FILE        *fp;
1866     t_commrec   *cr;
1867
1868     /* Print out how long the tuning took */
1869     double      seconds;
1870
1871     static t_filenm fnm[] = {
1872       /* g_tune_pme */
1873       { efOUT, "-p",      "perf",     ffWRITE },
1874       { efLOG, "-err",    "bencherr", ffWRITE },
1875       { efTPX, "-so",     "tuned",    ffWRITE },
1876       /* mdrun: */
1877       { efTPX, NULL,      NULL,       ffREAD },
1878       { efTRN, "-o",      NULL,       ffWRITE },
1879       { efXTC, "-x",      NULL,       ffOPTWR },
1880       { efCPT, "-cpi",    NULL,       ffOPTRD },
1881       { efCPT, "-cpo",    NULL,       ffOPTWR },
1882       { efSTO, "-c",      "confout",  ffWRITE },
1883       { efEDR, "-e",      "ener",     ffWRITE },
1884       { efLOG, "-g",      "md",       ffWRITE },
1885       { efXVG, "-dhdl",   "dhdl",     ffOPTWR },
1886       { efXVG, "-field",  "field",    ffOPTWR },
1887       { efXVG, "-table",  "table",    ffOPTRD },
1888       { efXVG, "-tabletf", "tabletf",   ffOPTRD },
1889       { efXVG, "-tablep", "tablep",   ffOPTRD },
1890       { efXVG, "-tableb", "table",    ffOPTRD },
1891       { efTRX, "-rerun",  "rerun",    ffOPTRD },
1892       { efXVG, "-tpi",    "tpi",      ffOPTWR },
1893       { efXVG, "-tpid",   "tpidist",  ffOPTWR },
1894       { efEDI, "-ei",     "sam",      ffOPTRD },
1895       { efEDO, "-eo",     "sam",      ffOPTWR },
1896       { efGCT, "-j",      "wham",     ffOPTRD },
1897       { efGCT, "-jo",     "bam",      ffOPTWR },
1898       { efXVG, "-ffout",  "gct",      ffOPTWR },
1899       { efXVG, "-devout", "deviatie", ffOPTWR },
1900       { efXVG, "-runav",  "runaver",  ffOPTWR },
1901       { efXVG, "-px",     "pullx",    ffOPTWR },
1902       { efXVG, "-pf",     "pullf",    ffOPTWR },
1903       { efXVG, "-ro",     "rotation", ffOPTWR },
1904       { efLOG, "-ra",     "rotangles",ffOPTWR },
1905       { efLOG, "-rs",     "rotslabs", ffOPTWR },
1906       { efLOG, "-rt",     "rottorque",ffOPTWR },
1907       { efMTX, "-mtx",    "nm",       ffOPTWR },
1908       { efNDX, "-dn",     "dipole",   ffOPTWR },
1909       /* Output files that are deleted after each benchmark run */
1910       { efTRN, "-bo",     "bench",    ffWRITE },
1911       { efXTC, "-bx",     "bench",    ffWRITE },
1912       { efCPT, "-bcpo",   "bench",    ffWRITE },
1913       { efSTO, "-bc",     "bench",    ffWRITE },
1914       { efEDR, "-be",     "bench",    ffWRITE },
1915       { efLOG, "-bg",     "bench",    ffWRITE },
1916       { efEDO, "-beo",    "bench",    ffOPTWR },
1917       { efXVG, "-bdhdl",  "benchdhdl",ffOPTWR },
1918       { efXVG, "-bfield", "benchfld" ,ffOPTWR },
1919       { efXVG, "-btpi",   "benchtpi", ffOPTWR },
1920       { efXVG, "-btpid",  "benchtpid",ffOPTWR },
1921       { efGCT, "-bjo",    "bench",    ffOPTWR },
1922       { efXVG, "-bffout", "benchgct", ffOPTWR },
1923       { efXVG, "-bdevout","benchdev", ffOPTWR },
1924       { efXVG, "-brunav", "benchrnav",ffOPTWR },
1925       { efXVG, "-bpx",    "benchpx",  ffOPTWR },
1926       { efXVG, "-bpf",    "benchpf",  ffOPTWR },
1927       { efXVG, "-bro",    "benchrot", ffOPTWR },
1928       { efLOG, "-bra",    "benchrota",ffOPTWR },
1929       { efLOG, "-brs",    "benchrots",ffOPTWR },
1930       { efLOG, "-brt",    "benchrott",ffOPTWR },
1931       { efMTX, "-bmtx",   "benchn",   ffOPTWR },
1932       { efNDX, "-bdn",    "bench",    ffOPTWR }
1933     };
1934
1935     gmx_bool bThreads     = FALSE;
1936
1937     int  nthreads=1;
1938
1939     const char *procstring[] =
1940       { NULL, "-np", "-n", "none", NULL };
1941     const char *npmevalues_opt[] =
1942       { NULL, "auto", "all", "subset", NULL };
1943
1944     gmx_bool bAppendFiles=TRUE;
1945     gmx_bool bKeepAndNumCPT=FALSE;
1946     gmx_bool bResetCountersHalfWay=FALSE;
1947     gmx_bool bBenchmark=TRUE;
1948
1949     output_env_t oenv=NULL;
1950
1951     t_pargs pa[] = {
1952       /***********************/
1953       /* g_tune_pme options: */
1954       /***********************/
1955       { "-np",       FALSE, etINT,  {&nnodes},
1956         "Number of nodes to run the tests on (must be > 2 for separate PME nodes)" },
1957       { "-npstring", FALSE, etENUM, {procstring},
1958         "Specify the number of processors to [TT]$MPIRUN[tt] using this string"},
1959       { "-nt",       FALSE, etINT,  {&nthreads},
1960         "Number of threads to run the tests on (turns MPI & mpirun off)"},
1961       { "-r",        FALSE, etINT,  {&repeats},
1962         "Repeat each test this often" },
1963       { "-max",      FALSE, etREAL, {&maxPMEfraction},
1964         "Max fraction of PME nodes to test with" },
1965       { "-min",      FALSE, etREAL, {&minPMEfraction},
1966         "Min fraction of PME nodes to test with" },
1967       { "-npme",     FALSE, etENUM, {npmevalues_opt},
1968         "Within -min and -max, benchmark all possible values for [TT]-npme[tt], or just a reasonable subset. "
1969         "Auto neglects -min and -max and chooses reasonable values around a guess for npme derived from the .tpr"},
1970       { "-fix",      FALSE, etINT,  {&npme_fixed},
1971         "If >= -1, do not vary the number of PME-only nodes, instead use this fixed value and only vary rcoulomb and the PME grid spacing."},
1972       { "-rmax",     FALSE, etREAL, {&rmax},
1973         "If >0, maximal rcoulomb for -ntpr>1 (rcoulomb upscaling results in fourier grid downscaling)" },
1974       { "-rmin",     FALSE, etREAL, {&rmin},
1975         "If >0, minimal rcoulomb for -ntpr>1" },
1976       { "-ntpr",     FALSE, etINT,  {&ntprs},
1977         "Number of [TT].tpr[tt] files to benchmark. Create this many files with different rcoulomb scaling factors depending on -rmin and -rmax. "
1978         "If < 1, automatically choose the number of [TT].tpr[tt] files to test" },
1979       { "-steps",    FALSE, etGMX_LARGE_INT, {&bench_nsteps},
1980         "Take timings for this many steps in the benchmark runs" },
1981       { "-resetstep",FALSE, etINT,  {&presteps},
1982         "Let dlb equilibrate this many steps before timings are taken (reset cycle counters after this many steps)" },
1983       { "-simsteps", FALSE, etGMX_LARGE_INT, {&new_sim_nsteps},
1984         "If non-negative, perform this many steps in the real run (overwrites nsteps from [TT].tpr[tt], add [TT].cpt[tt] steps)" },
1985       { "-launch",   FALSE, etBOOL, {&bLaunch},
1986         "Launch the real simulation after optimization" },
1987       { "-bench",    FALSE, etBOOL, {&bBenchmark},
1988         "Run the benchmarks or just create the input [TT].tpr[tt] files?" },
1989       /******************/
1990       /* mdrun options: */
1991       /******************/
1992       /* We let g_tune_pme parse and understand these options, because we need to
1993        * prevent that they appear on the mdrun command line for the benchmarks */
1994       { "-append",   FALSE, etBOOL, {&bAppendFiles},
1995         "Append to previous output files when continuing from checkpoint instead of adding the simulation part number to all file names (for launch only)" },
1996       { "-cpnum",    FALSE, etBOOL, {&bKeepAndNumCPT},
1997         "Keep and number checkpoint files (launch only)" },
1998       { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
1999         "HIDDENReset the cycle counters after half the number of steps or halfway [TT]-maxh[tt] (launch only)" }
2000     };
2001
2002
2003 #define NFILE asize(fnm)
2004
2005     CopyRight(stderr,argv[0]);
2006
2007     seconds = gettime();
2008
2009     parse_common_args(&argc,argv,PCA_NOEXIT_ON_ARGS,
2010                       NFILE,fnm,asize(pa),pa,asize(desc),desc,
2011                       0,NULL,&oenv);
2012
2013     /* Store the remaining unparsed command line entries in a string which
2014      * is then attached to the mdrun command line */
2015     snew(ExtraArgs, 1);
2016     ExtraArgs[0] = '\0';
2017     for (i=1; i<argc; i++) /* argc will now be 1 if everything was understood */
2018     {
2019         add_to_string(&ExtraArgs, argv[i]);
2020         add_to_string(&ExtraArgs, " ");
2021     }
2022
2023     if (opt2parg_bSet("-nt",asize(pa),pa))
2024     {
2025         bThreads=TRUE;
2026         if (opt2parg_bSet("-npstring",asize(pa),pa))
2027             fprintf(stderr, "WARNING: -npstring has no effect when using threads.\n");
2028
2029         if (nnodes > 1)
2030             gmx_fatal(FARGS, "Can't run multi-threaded MPI simulation yet!");
2031         /* and now we just set this; a bit of an ugly hack*/
2032         nnodes=nthreads;
2033     }
2034     /* Check for PME:PP ratio and whether tpr triggers additional output files */
2035     guessPMEratio = inspect_tpr(NFILE,fnm,&rcoulomb);
2036
2037     /* Automatically set -beo options if -eo is set etc. */
2038     couple_files_options(NFILE,fnm);
2039
2040     /* Construct the command line arguments for benchmark runs
2041      * as well as for the simulation run */
2042     if (bThreads)
2043         sprintf(bbuf," -nt %d ", nthreads);
2044     else
2045         sprintf(bbuf," -np %d ", nnodes);
2046
2047     cmd_np = bbuf;
2048
2049     create_command_line_snippets(bThreads,bAppendFiles,bKeepAndNumCPT,bResetCountersHalfWay,presteps,
2050                                  NFILE,fnm,asize(pa),pa,procstring[0],
2051                                  &cmd_np, &cmd_args_bench, &cmd_args_launch,
2052                                  ExtraArgs);
2053
2054     /* Read in checkpoint file if requested */
2055     sim_part = 1;
2056     if(opt2bSet("-cpi",NFILE,fnm))
2057     {
2058         snew(cr,1);
2059         cr->duty=DUTY_PP; /* makes the following routine happy */
2060         read_checkpoint_simulation_part(opt2fn("-cpi",NFILE,fnm),
2061                     &sim_part,&cpt_steps,cr,
2062                     FALSE,NFILE,fnm,NULL,NULL);
2063         sfree(cr);
2064         sim_part++;
2065         /* sim_part will now be 1 if no checkpoint file was found */
2066         if (sim_part<=1)
2067             gmx_fatal(FARGS, "Checkpoint file %s not found!", opt2fn("-cpi",NFILE,fnm));
2068     }
2069
2070     /* Open performance output file and write header info */
2071     fp = ffopen(opt2fn("-p",NFILE,fnm),"w");
2072
2073     /* Make a quick consistency check of command line parameters */
2074     check_input(nnodes, repeats, &ntprs, &rmin, rcoulomb, &rmax,
2075                 maxPMEfraction, minPMEfraction, npme_fixed,
2076                 bench_nsteps, fnm, NFILE, sim_part, presteps,
2077                 asize(pa),pa);
2078
2079     /* Determine the maximum and minimum number of PME nodes to test,
2080      * the actual list of settings is build in do_the_tests(). */
2081     if ((nnodes > 2) && (npme_fixed < -1))
2082     {
2083         if (0 == strcmp(npmevalues_opt[0], "auto"))
2084         {
2085             /* Determine the npme range automatically based on the PME:PP load guess */
2086             if (guessPMEratio > 1.0)
2087             {
2088                 /* More PME than PP work, probably we do not need separate PME nodes at all! */
2089                 maxPMEnodes=nnodes/2;
2090                 minPMEnodes=nnodes/2;
2091             }
2092             else
2093             {
2094                 /* PME : PP load is in the range 0..1, let's test around the guess */
2095                 guessPMEnodes = nnodes/(1.0 + 1.0/guessPMEratio);
2096                 minPMEnodes = floor(0.7*guessPMEnodes);
2097                 maxPMEnodes =  ceil(1.6*guessPMEnodes);
2098                 maxPMEnodes = min(maxPMEnodes, nnodes/2);
2099             }
2100         }
2101         else
2102         {
2103             /* Determine the npme range based on user input */
2104             maxPMEnodes = floor(maxPMEfraction*nnodes);
2105             minPMEnodes = max(floor(minPMEfraction*nnodes), 0);
2106             fprintf(stdout, "Will try runs with %d ", minPMEnodes);
2107             if (maxPMEnodes != minPMEnodes)
2108                 fprintf(stdout, "- %d ", maxPMEnodes);
2109             fprintf(stdout, "PME-only nodes.\n  Note that the automatic number of PME-only nodes and no separate PME nodes are always tested.\n");
2110         }
2111     }
2112     else
2113     {
2114         maxPMEnodes = 0;
2115         minPMEnodes = 0;
2116     }
2117
2118     /* Get the commands we need to set up the runs from environment variables */
2119     get_program_paths(bThreads, &cmd_mpirun, cmd_np, &cmd_mdrun, repeats);
2120
2121     /* Print some header info to file */
2122     sep_line(fp);
2123     fprintf(fp, "\n      P E R F O R M A N C E   R E S U L T S\n");
2124     sep_line(fp);
2125     fprintf(fp, "%s for Gromacs %s\n", ShortProgram(),GromacsVersion());
2126     if (!bThreads)
2127     {
2128         fprintf(fp, "Number of nodes         : %d\n", nnodes);
2129         fprintf(fp, "The mpirun command is   : %s\n", cmd_mpirun);
2130         if ( strcmp(procstring[0], "none") != 0)
2131             fprintf(fp, "Passing # of nodes via  : %s\n", procstring[0]);
2132         else
2133             fprintf(fp, "Not setting number of nodes in system call\n");
2134     }
2135     else
2136         fprintf(fp, "Number of threads       : %d\n", nnodes);
2137
2138     fprintf(fp, "The mdrun  command is   : %s\n", cmd_mdrun);
2139     fprintf(fp, "mdrun args benchmarks   : %s\n", cmd_args_bench);
2140     fprintf(fp, "Benchmark steps         : ");
2141     fprintf(fp, gmx_large_int_pfmt, bench_nsteps);
2142     fprintf(fp, "\n");
2143     fprintf(fp, "dlb equilibration steps : %d\n", presteps);
2144     if (sim_part > 1)
2145     {
2146         fprintf(fp, "Checkpoint time step    : ");
2147         fprintf(fp, gmx_large_int_pfmt, cpt_steps);
2148         fprintf(fp, "\n");
2149     }
2150     fprintf(fp, "mdrun args at launchtime: %s\n", cmd_args_launch);
2151
2152     if (new_sim_nsteps >= 0)
2153     {
2154         bOverwrite = TRUE;
2155         fprintf(stderr, "Note: Simulation input file %s will have ", opt2fn("-so",NFILE,fnm));
2156         fprintf(stderr, gmx_large_int_pfmt, new_sim_nsteps+cpt_steps);
2157         fprintf(stderr, " steps.\n");
2158         fprintf(fp, "Simulation steps        : ");
2159         fprintf(fp, gmx_large_int_pfmt, new_sim_nsteps);
2160         fprintf(fp, "\n");
2161     }
2162     if (repeats > 1)
2163         fprintf(fp, "Repeats for each test   : %d\n", repeats);
2164
2165     if (npme_fixed >= -1)
2166     {
2167         fprintf(fp, "Fixing -npme at         : %d\n", npme_fixed);
2168     }
2169
2170     fprintf(fp, "Input file              : %s\n", opt2fn("-s",NFILE,fnm));
2171     fprintf(fp, "   PME/PP load estimate : %g\n", guessPMEratio);
2172
2173     /* Allocate memory for the inputinfo struct: */
2174     snew(info, 1);
2175     info->nr_inputfiles = ntprs;
2176     for (i=0; i<ntprs; i++)
2177     {
2178         snew(info->rcoulomb , ntprs);
2179         snew(info->rvdw     , ntprs);
2180         snew(info->rlist    , ntprs);
2181         snew(info->rlistlong, ntprs);
2182         snew(info->nkx      , ntprs);
2183         snew(info->nky      , ntprs);
2184         snew(info->nkz      , ntprs);
2185         snew(info->fsx      , ntprs);
2186         snew(info->fsy      , ntprs);
2187         snew(info->fsz      , ntprs);
2188     }
2189     /* Make alternative tpr files to test: */
2190     snew(tpr_names, ntprs);
2191     for (i=0; i<ntprs; i++)
2192         snew(tpr_names[i], STRLEN);
2193
2194     /* It can be that ntprs is reduced by make_benchmark_tprs if not enough
2195      * different grids could be found. */
2196     make_benchmark_tprs(opt2fn("-s",NFILE,fnm), tpr_names, bench_nsteps+presteps,
2197             cpt_steps, rmin, rmax, &ntprs, info, fp);
2198
2199     /********************************************************************************/
2200     /* Main loop over all scenarios we need to test: tpr files, PME nodes, repeats  */
2201     /********************************************************************************/
2202     snew(perfdata, ntprs);
2203     if (bBenchmark)
2204     {
2205         do_the_tests(fp, tpr_names, maxPMEnodes, minPMEnodes, npme_fixed, npmevalues_opt[0], perfdata, &pmeentries,
2206                 repeats, nnodes, ntprs, bThreads, cmd_mpirun, cmd_np, cmd_mdrun,
2207                 cmd_args_bench, fnm, NFILE, sim_part, presteps, cpt_steps);
2208
2209         fprintf(fp, "\nTuning took%8.1f minutes.\n", (gettime()-seconds)/60.0);
2210
2211         /* Analyse the results and give a suggestion for optimal settings: */
2212         bKeepTPR = analyze_data(fp, opt2fn("-p", NFILE, fnm), perfdata, nnodes, ntprs, pmeentries,
2213                 repeats, info, &best_tpr, &best_npme);
2214
2215         /* Take the best-performing tpr file and enlarge nsteps to original value */
2216         if ( bKeepTPR && !bOverwrite )
2217         {
2218             simulation_tpr = opt2fn("-s",NFILE,fnm);
2219         }
2220         else
2221         {
2222             simulation_tpr = opt2fn("-so",NFILE,fnm);
2223             modify_PMEsettings(bOverwrite? (new_sim_nsteps+cpt_steps) :
2224                     info->orig_sim_steps, tpr_names[best_tpr],
2225                     simulation_tpr);
2226         }
2227
2228         /* Let's get rid of the temporary benchmark input files */
2229         for (i=0; i<ntprs; i++)
2230         {
2231             fprintf(stdout, "Deleting temporary benchmark input file %s\n", tpr_names[i]);
2232             remove(tpr_names[i]);
2233         }
2234
2235         /* Now start the real simulation if the user requested it ... */
2236         launch_simulation(bLaunch, fp, bThreads, cmd_mpirun, cmd_np, cmd_mdrun,
2237                 cmd_args_launch, simulation_tpr, nnodes, best_npme);
2238     }
2239     ffclose(fp);
2240
2241     /* ... or simply print the performance results to screen: */
2242     if (!bLaunch)
2243         finalize(opt2fn("-p", NFILE, fnm));
2244
2245     return 0;
2246 }