src/gromacs/gmxana/gmx_tune_pme.c

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