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