src/gromacs/gmxana/gmx_tune_pme.c

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