Update copyright statements and change license to LGPL
[alexxy/gromacs.git] / src / gmxlib / network.c
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5  * Copyright (c) 2001-2004, The GROMACS development team,
6  * check out http://www.gromacs.org for more information.
7  * Copyright (c) 2012, by the GROMACS development team, led by
8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9  * others, as listed in the AUTHORS file in the top-level source
10  * directory and at http://www.gromacs.org.
11  *
12  * GROMACS is free software; you can redistribute it and/or
13  * modify it under the terms of the GNU Lesser General Public License
14  * as published by the Free Software Foundation; either version 2.1
15  * of the License, or (at your option) any later version.
16  *
17  * GROMACS is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
20  * Lesser General Public License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public
23  * License along with GROMACS; if not, see
24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
26  *
27  * If you want to redistribute modifications to GROMACS, please
28  * consider that scientific software is very special. Version
29  * control is crucial - bugs must be traceable. We will be happy to
30  * consider code for inclusion in the official distribution, but
31  * derived work must not be called official GROMACS. Details are found
32  * in the README & COPYING files - if they are missing, get the
33  * official version at http://www.gromacs.org.
34  *
35  * To help us fund GROMACS development, we humbly ask that you cite
36  * the research papers on the package. Check out http://www.gromacs.org.
37  */
38 #ifdef HAVE_CONFIG_H
39 #include <config.h>
40 #endif
41
42 #include <string.h>
43 #include "gmx_fatal.h"
44 #include "main.h"
45 #include "smalloc.h"
46 #include "network.h"
47 #include "copyrite.h"
48 #include "statutil.h"
49 #include <ctype.h>
50 #include "macros.h"
51
52 #ifdef GMX_LIB_MPI
53 #include <mpi.h>
54 #endif
55
56 #ifdef GMX_THREAD_MPI
57 #include "tmpi.h"
58 #endif
59
60 #include "mpelogging.h"
61
62 /* The source code in this file should be thread-safe. 
63       Please keep it that way. */
64
65 gmx_bool gmx_mpi_initialized(void)
66 {
67   int n;
68 #ifndef GMX_MPI
69   return 0;
70 #else
71   MPI_Initialized(&n);
72   
73   return n;
74 #endif
75 }
76
77 int gmx_setup(int *argc,char **argv,int *nnodes)
78 {
79 #ifndef GMX_MPI
80   gmx_call("gmx_setup");
81   return 0;
82 #else
83   char   buf[256];
84   int    resultlen;               /* actual length of node name      */
85   int    i,flag;
86   int  mpi_num_nodes;
87   int  mpi_my_rank;
88   char mpi_hostname[MPI_MAX_PROCESSOR_NAME];
89
90   /* Call the MPI routines */
91 #ifdef GMX_LIB_MPI
92 #ifdef GMX_FAHCORE
93   (void) fah_MPI_Init(argc,&argv);
94 #else
95   (void) MPI_Init(argc,&argv);
96 #endif
97 #endif
98   (void) MPI_Comm_size( MPI_COMM_WORLD, &mpi_num_nodes );
99   (void) MPI_Comm_rank( MPI_COMM_WORLD, &mpi_my_rank );
100   (void) MPI_Get_processor_name( mpi_hostname, &resultlen );
101
102
103 #ifdef USE_MPE
104   /* MPE logging routines. Get event IDs from MPE: */
105   /* General events */
106   ev_timestep1               = MPE_Log_get_event_number( );
107   ev_timestep2               = MPE_Log_get_event_number( );
108   ev_force_start             = MPE_Log_get_event_number( );
109   ev_force_finish            = MPE_Log_get_event_number( );
110   ev_do_fnbf_start           = MPE_Log_get_event_number( );
111   ev_do_fnbf_finish          = MPE_Log_get_event_number( );
112   ev_ns_start                = MPE_Log_get_event_number( );
113   ev_ns_finish               = MPE_Log_get_event_number( );
114   ev_calc_bonds_start        = MPE_Log_get_event_number( );
115   ev_calc_bonds_finish       = MPE_Log_get_event_number( );
116   ev_global_stat_start       = MPE_Log_get_event_number( );
117   ev_global_stat_finish      = MPE_Log_get_event_number( );
118   ev_virial_start            = MPE_Log_get_event_number( );
119   ev_virial_finish           = MPE_Log_get_event_number( );
120   
121   /* Shift related events */
122   ev_shift_start             = MPE_Log_get_event_number( );
123   ev_shift_finish            = MPE_Log_get_event_number( );
124   ev_unshift_start           = MPE_Log_get_event_number( );
125   ev_unshift_finish          = MPE_Log_get_event_number( );
126   ev_mk_mshift_start         = MPE_Log_get_event_number( );
127   ev_mk_mshift_finish        = MPE_Log_get_event_number( );
128   
129   /* PME related events */
130   ev_pme_start               = MPE_Log_get_event_number( );
131   ev_pme_finish              = MPE_Log_get_event_number( );
132   ev_spread_on_grid_start    = MPE_Log_get_event_number( );
133   ev_spread_on_grid_finish   = MPE_Log_get_event_number( );
134   ev_sum_qgrid_start         = MPE_Log_get_event_number( );
135   ev_sum_qgrid_finish        = MPE_Log_get_event_number( );
136   ev_gmxfft3d_start          = MPE_Log_get_event_number( );
137   ev_gmxfft3d_finish         = MPE_Log_get_event_number( );
138   ev_solve_pme_start         = MPE_Log_get_event_number( );
139   ev_solve_pme_finish        = MPE_Log_get_event_number( );
140   ev_gather_f_bsplines_start = MPE_Log_get_event_number( );
141   ev_gather_f_bsplines_finish= MPE_Log_get_event_number( );
142   ev_reduce_start            = MPE_Log_get_event_number( );
143   ev_reduce_finish           = MPE_Log_get_event_number( );
144   ev_rscatter_start          = MPE_Log_get_event_number( );
145   ev_rscatter_finish         = MPE_Log_get_event_number( );
146   ev_alltoall_start          = MPE_Log_get_event_number( );
147   ev_alltoall_finish         = MPE_Log_get_event_number( );
148   ev_pmeredist_start         = MPE_Log_get_event_number( );
149   ev_pmeredist_finish        = MPE_Log_get_event_number( );
150   ev_init_pme_start          = MPE_Log_get_event_number( );      
151   ev_init_pme_finish         = MPE_Log_get_event_number( );
152   ev_send_coordinates_start  = MPE_Log_get_event_number( );
153   ev_send_coordinates_finish = MPE_Log_get_event_number( );
154   ev_update_fr_start         = MPE_Log_get_event_number( );
155   ev_update_fr_finish        = MPE_Log_get_event_number( );
156   ev_clear_rvecs_start       = MPE_Log_get_event_number( );
157   ev_clear_rvecs_finish      = MPE_Log_get_event_number( ); 
158   ev_update_start            = MPE_Log_get_event_number( ); 
159   ev_update_finish           = MPE_Log_get_event_number( ); 
160   ev_output_start            = MPE_Log_get_event_number( ); 
161   ev_output_finish           = MPE_Log_get_event_number( ); 
162   ev_sum_lrforces_start      = MPE_Log_get_event_number( ); 
163   ev_sum_lrforces_finish     = MPE_Log_get_event_number( ); 
164   ev_sort_start              = MPE_Log_get_event_number( );
165   ev_sort_finish             = MPE_Log_get_event_number( );
166   ev_sum_qgrid_start         = MPE_Log_get_event_number( );
167   ev_sum_qgrid_finish        = MPE_Log_get_event_number( );
168   
169   /* Essential dynamics related events */
170   ev_edsam_start             = MPE_Log_get_event_number( );
171   ev_edsam_finish            = MPE_Log_get_event_number( );
172   ev_get_coords_start        = MPE_Log_get_event_number( );
173   ev_get_coords_finish       = MPE_Log_get_event_number( );
174   ev_ed_apply_cons_start     = MPE_Log_get_event_number( );
175   ev_ed_apply_cons_finish    = MPE_Log_get_event_number( );
176   ev_fit_to_reference_start  = MPE_Log_get_event_number( );
177   ev_fit_to_reference_finish = MPE_Log_get_event_number( );
178   
179   /* describe events: */
180   if ( mpi_my_rank == 0 ) 
181   {
182     /* General events */
183     MPE_Describe_state(ev_timestep1,               ev_timestep2,                "timestep START",  "magenta" );
184     MPE_Describe_state(ev_force_start,             ev_force_finish,             "force",           "cornflower blue" );
185     MPE_Describe_state(ev_do_fnbf_start,           ev_do_fnbf_finish,           "do_fnbf",         "navy" );
186     MPE_Describe_state(ev_ns_start,                ev_ns_finish,                "neighbor search", "tomato" );
187     MPE_Describe_state(ev_calc_bonds_start,        ev_calc_bonds_finish,        "bonded forces",   "slate blue" );
188     MPE_Describe_state(ev_global_stat_start,       ev_global_stat_finish,       "global stat",     "firebrick3");
189     MPE_Describe_state(ev_update_fr_start,         ev_update_fr_finish,         "update forcerec", "goldenrod");
190     MPE_Describe_state(ev_clear_rvecs_start,       ev_clear_rvecs_finish,       "clear rvecs",     "bisque");
191     MPE_Describe_state(ev_update_start,            ev_update_finish,            "update",          "cornsilk");
192     MPE_Describe_state(ev_output_start,            ev_output_finish,            "output",          "black");
193     MPE_Describe_state(ev_virial_start,            ev_virial_finish,            "calc_virial",     "thistle4");
194     
195     /* PME related events */
196     MPE_Describe_state(ev_pme_start,               ev_pme_finish,               "doing PME",       "grey" );
197     MPE_Describe_state(ev_spread_on_grid_start,    ev_spread_on_grid_finish,    "spread",          "dark orange" );   
198     MPE_Describe_state(ev_sum_qgrid_start,         ev_sum_qgrid_finish,         "sum qgrid",       "slate blue");
199     MPE_Describe_state(ev_gmxfft3d_start,          ev_gmxfft3d_finish,          "fft3d",           "snow2" );   
200     MPE_Describe_state(ev_solve_pme_start,         ev_solve_pme_finish,         "solve PME",       "indian red" );   
201     MPE_Describe_state(ev_gather_f_bsplines_start, ev_gather_f_bsplines_finish, "bsplines",        "light sea green" );   
202     MPE_Describe_state(ev_reduce_start,            ev_reduce_finish,            "reduce",          "cyan1" );
203     MPE_Describe_state(ev_rscatter_start,          ev_rscatter_finish,          "rscatter",        "cyan3" );
204     MPE_Describe_state(ev_alltoall_start,          ev_alltoall_finish,          "alltoall",        "LightCyan4" );
205     MPE_Describe_state(ev_pmeredist_start,         ev_pmeredist_finish,         "pmeredist",       "thistle" );
206     MPE_Describe_state(ev_init_pme_start,          ev_init_pme_finish,          "init PME",        "snow4");
207     MPE_Describe_state(ev_send_coordinates_start,  ev_send_coordinates_finish,  "send_coordinates","blue");
208     MPE_Describe_state(ev_sum_lrforces_start,      ev_sum_lrforces_finish,      "sum_LRforces",    "lime green");
209     MPE_Describe_state(ev_sort_start,              ev_sort_finish,              "sort pme atoms",  "brown");
210     MPE_Describe_state(ev_sum_qgrid_start,         ev_sum_qgrid_finish,         "sum charge grid", "medium orchid");
211     
212     /* Shift related events */
213     MPE_Describe_state(ev_shift_start,             ev_shift_finish,             "shift",           "orange");
214     MPE_Describe_state(ev_unshift_start,           ev_unshift_finish,           "unshift",         "dark orange");    
215     MPE_Describe_state(ev_mk_mshift_start,         ev_mk_mshift_finish,         "mk_mshift",       "maroon");
216         
217     /* Essential dynamics related events */
218     MPE_Describe_state(ev_edsam_start,             ev_edsam_finish,             "EDSAM",           "deep sky blue");
219     MPE_Describe_state(ev_get_coords_start,        ev_get_coords_finish,        "ED get coords",   "steel blue");
220     MPE_Describe_state(ev_ed_apply_cons_start,     ev_ed_apply_cons_finish,     "ED apply constr", "forest green");
221     MPE_Describe_state(ev_fit_to_reference_start,  ev_fit_to_reference_finish,  "ED fit to ref",   "lavender");
222        
223   }
224   MPE_Init_log();
225 #endif
226  
227 #ifdef GMX_LIB_MPI 
228   fprintf(stderr,"NNODES=%d, MYRANK=%d, HOSTNAME=%s\n",
229           mpi_num_nodes,mpi_my_rank,mpi_hostname);
230 #endif
231   
232   *nnodes=mpi_num_nodes;
233   
234   return mpi_my_rank;
235 #endif
236 }
237
238 int  gmx_node_num(void)
239 {
240 #ifndef GMX_MPI
241   return 1;
242 #else
243   int i;
244   (void) MPI_Comm_size(MPI_COMM_WORLD, &i);
245   return i;
246 #endif
247 }
248
249 int gmx_node_rank(void)
250 {
251 #ifndef GMX_MPI
252   return 0;
253 #else
254   int i;
255   (void) MPI_Comm_rank(MPI_COMM_WORLD, &i);
256   return i;
257 #endif
258 }
259
260
261 int gmx_hostname_num()
262 {
263 #ifndef GMX_MPI
264   return 0;
265 #else
266 #ifdef GMX_THREAD_MPI
267   /* thread-MPI currently puts the thread number in the process name,
268    * we might want to change this, as this is inconsistent with what
269    * most MPI implementations would do when running on a single node.
270    */
271   return 0;
272 #else
273   int  resultlen,hostnum,i,j;
274   char mpi_hostname[MPI_MAX_PROCESSOR_NAME],hostnum_str[MPI_MAX_PROCESSOR_NAME];
275
276   MPI_Get_processor_name(mpi_hostname,&resultlen);
277   /* This procedure can only differentiate nodes with host names
278    * that end on unique numbers.
279    */
280   i = 0;
281   j = 0;
282   /* Only parse the host name up to the first dot */
283   while(i < resultlen && mpi_hostname[i] != '.') {
284     if (isdigit(mpi_hostname[i])) {
285       hostnum_str[j++] = mpi_hostname[i];
286     }
287     i++;
288   }
289   hostnum_str[j] = '\0';
290   if (j == 0) {
291     hostnum = 0;
292   } else {
293     /* Use only the last 9 decimals, so we don't overflow an int */
294     hostnum = strtol(hostnum_str + max(0,j-9), NULL, 10);
295   }
296
297   if (debug) {
298     fprintf(debug,"In gmx_setup_nodecomm: hostname '%s', hostnum %d\n",
299         mpi_hostname,hostnum);
300   }
301   return hostnum;
302 #endif
303 #endif
304 }
305
306 void gmx_setup_nodecomm(FILE *fplog,t_commrec *cr)
307 {
308     gmx_nodecomm_t *nc;
309     int  n,rank,hostnum,ng,ni;
310
311     /* Many MPI implementations do not optimize MPI_Allreduce
312      * (and probably also other global communication calls)
313      * for multi-core nodes connected by a network.
314      * We can optimize such communication by using one MPI call
315      * within each node and one between the nodes.
316      * For MVAPICH2 and Intel MPI this reduces the time for
317      * the global_stat communication by 25%
318      * for 2x2-core 3 GHz Woodcrest connected by mixed DDR/SDR Infiniband.
319      * B. Hess, November 2007
320      */
321
322     nc = &cr->nc;
323
324     nc->bUse = FALSE;
325 #ifndef GMX_THREAD_MPI
326 #ifdef GMX_MPI
327     MPI_Comm_size(cr->mpi_comm_mygroup,&n);
328     MPI_Comm_rank(cr->mpi_comm_mygroup,&rank);
329
330     hostnum = gmx_hostname_num();
331
332     if (debug)
333     {
334         fprintf(debug,"In gmx_setup_nodecomm: splitting communicator of size %d\n",n);
335     }
336
337
338     /* The intra-node communicator, split on node number */
339     MPI_Comm_split(cr->mpi_comm_mygroup,hostnum,rank,&nc->comm_intra);
340     MPI_Comm_rank(nc->comm_intra,&nc->rank_intra);
341     if (debug)
342     {
343         fprintf(debug,"In gmx_setup_nodecomm: node rank %d rank_intra %d\n",
344                 rank,nc->rank_intra);
345     }
346     /* The inter-node communicator, split on rank_intra.
347      * We actually only need the one for rank=0,
348      * but it is easier to create them all.
349      */
350     MPI_Comm_split(cr->mpi_comm_mygroup,nc->rank_intra,rank,&nc->comm_inter);
351     /* Check if this really created two step communication */
352     MPI_Comm_size(nc->comm_inter,&ng);
353     MPI_Comm_size(nc->comm_intra,&ni);
354     if (debug)
355     {
356         fprintf(debug,"In gmx_setup_nodecomm: groups %d, my group size %d\n",
357                 ng,ni);
358     }
359
360     if (getenv("GMX_NO_NODECOMM") == NULL &&
361         ((ng > 1 && ng < n) || (ni > 1 && ni < n)))
362     {
363         nc->bUse = TRUE;
364         if (fplog)
365         {
366             fprintf(fplog,"Using two step summing over %d groups of on average %.1f processes\n\n",
367                     ng,(real)n/(real)ng);
368         }
369         if (nc->rank_intra > 0)
370         {
371             MPI_Comm_free(&nc->comm_inter);
372         }
373     }
374     else
375     {
376         /* One group or all processes in a separate group, use normal summing */
377         MPI_Comm_free(&nc->comm_inter);
378         MPI_Comm_free(&nc->comm_intra);
379         if (debug)
380         {
381             fprintf(debug,"In gmx_setup_nodecomm: not unsing separate inter- and intra-node communicators.\n");
382         }
383     }
384 #endif
385 #else
386     /* tMPI runs only on a single node so just use the nodeid */
387     nc->rank_intra = cr->nodeid;
388 #endif
389 }
390
391 void gmx_init_intra_counters(t_commrec *cr)
392 {
393     /* counters for PP+PME and PP-only processes on my node */
394     int nnodes, nnodes_pp, id_mynode=-1, id_mynode_group=-1, nproc_mynode, nproc_mynode_pp;
395 #if defined GMX_MPI && !defined GMX_THREAD_MPI
396     int i, mynum, *num, *num_s, *num_pp, *num_pp_s;
397 #endif
398
399     nnodes    = cr->nnodes;
400     nnodes_pp = nnodes - cr->npmenodes;
401
402 #if defined GMX_MPI && !defined GMX_THREAD_MPI
403     /* We have MPI and can expect to have different compute nodes */
404     mynum = gmx_hostname_num();
405
406     /* We can't rely on MPI_IN_PLACE, so we need send and receive buffers */
407     snew(num,   nnodes);
408     snew(num_s, nnodes);
409     snew(num_pp,   nnodes_pp);
410     snew(num_pp_s, nnodes_pp);
411
412     num_s[cr->sim_nodeid] = mynum;
413     if (cr->duty & DUTY_PP)
414     {
415         num_pp_s[cr->nodeid] = mynum;
416     }
417
418     MPI_Allreduce(num_s, num, nnodes, MPI_INT, MPI_SUM, cr->mpi_comm_mysim);
419     MPI_Allreduce(num_pp_s, num_pp, nnodes_pp, MPI_INT, MPI_SUM, cr->mpi_comm_mygroup);
420
421     id_mynode       = 0;
422     id_mynode_group = 0;
423     nproc_mynode    = 0;
424     nproc_mynode_pp = 0;
425     for(i=0; i<nnodes; i++)
426     {
427         if (num[i] == mynum)
428         {
429             nproc_mynode++;
430             if (i < cr->sim_nodeid)
431             {
432                 id_mynode++;
433             }
434             if (i < cr->nodeid)
435             {
436                 id_mynode_group++;
437             }
438         }
439     }
440     for(i=0; i<nnodes_pp; i++)
441     {
442         if (num_pp[i] == mynum)
443         {
444             nproc_mynode_pp++;
445         }
446     }
447     sfree(num);
448     sfree(num_s);
449     sfree(num_pp);
450     sfree(num_pp_s);
451 #else
452     /* Serial or thread-MPI code, we are running within a node */
453     id_mynode       = cr->sim_nodeid;
454     id_mynode_group = cr->nodeid;
455     nproc_mynode    = cr->nnodes;
456     nproc_mynode_pp = cr->nnodes - cr->npmenodes;
457 #endif
458
459     if (debug)
460     {
461         char sbuf[STRLEN];
462         if (cr->duty & DUTY_PP && cr->duty & DUTY_PME)
463         {
464             sprintf(sbuf, "PP+PME");
465         }
466         else
467         {
468             sprintf(sbuf, "%s", cr->duty & DUTY_PP ? "PP" : "PME");
469         }
470         fprintf(debug, "On %3s node %d: nodeid_intra=%d, nodeid_group_intra=%d, "
471                 "nnodes_intra=%d, nnodes_pp_intra=%d\n", sbuf, cr->sim_nodeid,
472                 id_mynode, id_mynode_group, nproc_mynode, nproc_mynode_pp);
473     }
474
475     cr->nodeid_intra        = id_mynode;
476     cr->nodeid_group_intra  = id_mynode_group;
477     cr->nnodes_intra        = nproc_mynode;
478     cr->nnodes_pp_intra     = nproc_mynode_pp;
479 }
480
481
482 void gmx_barrier(const t_commrec *cr)
483 {
484 #ifndef GMX_MPI
485   gmx_call("gmx_barrier");
486 #else
487   MPI_Barrier(cr->mpi_comm_mygroup);
488 #endif
489 }
490
491 void gmx_abort(int noderank,int nnodes,int errorno)
492 {
493 #ifndef GMX_MPI
494   gmx_call("gmx_abort");
495 #else
496 #ifdef GMX_THREAD_MPI
497   fprintf(stderr,"Halting program %s\n",ShortProgram());
498   thanx(stderr);
499   exit(1);
500 #else
501   if (nnodes > 1)
502   {
503       fprintf(stderr,"Halting parallel program %s on CPU %d out of %d\n",
504               ShortProgram(),noderank,nnodes);
505   }
506   else
507   {
508       fprintf(stderr,"Halting program %s\n",ShortProgram());
509   }
510
511   thanx(stderr);
512   MPI_Abort(MPI_COMM_WORLD,errorno);
513   exit(1);
514 #endif
515 #endif
516 }
517
518 void gmx_bcast(int nbytes,void *b,const t_commrec *cr)
519 {
520 #ifndef GMX_MPI
521   gmx_call("gmx_bast");
522 #else
523   MPI_Bcast(b,nbytes,MPI_BYTE,MASTERRANK(cr),cr->mpi_comm_mygroup);
524 #endif
525 }
526
527 void gmx_bcast_sim(int nbytes,void *b,const t_commrec *cr)
528 {
529 #ifndef GMX_MPI
530   gmx_call("gmx_bast");
531 #else
532   MPI_Bcast(b,nbytes,MPI_BYTE,MASTERRANK(cr),cr->mpi_comm_mysim);
533 #endif
534 }
535
536 void gmx_sumd(int nr,double r[],const t_commrec *cr)
537 {
538 #ifndef GMX_MPI
539     gmx_call("gmx_sumd");
540 #else
541 #if defined(MPI_IN_PLACE_EXISTS) || defined(GMX_THREAD_MPI)
542     if (cr->nc.bUse) {
543         if (cr->nc.rank_intra == 0)
544         {
545             /* Use two step summing. */
546             MPI_Reduce(MPI_IN_PLACE,r,nr,MPI_DOUBLE,MPI_SUM,0,
547                        cr->nc.comm_intra);
548             /* Sum the roots of the internal (intra) buffers. */
549             MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_DOUBLE,MPI_SUM,
550                           cr->nc.comm_inter);
551         }
552         else
553         {
554             /* This is here because of the silly MPI specification
555                 that MPI_IN_PLACE should be put in sendbuf instead of recvbuf */
556             MPI_Reduce(r,NULL,nr,MPI_DOUBLE,MPI_SUM,0,cr->nc.comm_intra);
557         }
558         MPI_Bcast(r,nr,MPI_DOUBLE,0,cr->nc.comm_intra);
559     } 
560     else 
561     {
562         MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_DOUBLE,MPI_SUM, 
563                       cr->mpi_comm_mygroup);
564     }
565 #else
566     int i;
567
568     if (nr > cr->mpb->dbuf_alloc) {
569         cr->mpb->dbuf_alloc = nr;
570         srenew(cr->mpb->dbuf,cr->mpb->dbuf_alloc);
571     }
572     if (cr->nc.bUse) {
573         /* Use two step summing */
574         MPI_Allreduce(r,cr->mpb->dbuf,nr,MPI_DOUBLE,MPI_SUM,cr->nc.comm_intra);
575         if (cr->nc.rank_intra == 0) {
576             /* Sum with the buffers reversed */
577             MPI_Allreduce(cr->mpb->dbuf,r,nr,MPI_DOUBLE,MPI_SUM, 
578                           cr->nc.comm_inter);
579         }
580         MPI_Bcast(r,nr,MPI_DOUBLE,0,cr->nc.comm_intra);
581     } else {
582         MPI_Allreduce(r,cr->mpb->dbuf,nr,MPI_DOUBLE,MPI_SUM,
583                       cr->mpi_comm_mygroup);
584         for(i=0; i<nr; i++)
585             r[i] = cr->mpb->dbuf[i];
586     }
587 #endif
588 #endif
589 }
590
591 void gmx_sumf(int nr,float r[],const t_commrec *cr)
592 {
593 #ifndef GMX_MPI
594     gmx_call("gmx_sumf");
595 #else
596 #if defined(MPI_IN_PLACE_EXISTS) || defined(GMX_THREAD_MPI)
597     if (cr->nc.bUse) {
598         /* Use two step summing.  */
599         if (cr->nc.rank_intra == 0)
600         {
601             MPI_Reduce(MPI_IN_PLACE,r,nr,MPI_FLOAT,MPI_SUM,0,
602                        cr->nc.comm_intra);
603             /* Sum the roots of the internal (intra) buffers */
604             MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_FLOAT,MPI_SUM,
605                           cr->nc.comm_inter);
606         }
607         else
608         {
609             /* This is here because of the silly MPI specification
610                 that MPI_IN_PLACE should be put in sendbuf instead of recvbuf */
611             MPI_Reduce(r,NULL,nr,MPI_FLOAT,MPI_SUM,0,cr->nc.comm_intra);
612         }
613         MPI_Bcast(r,nr,MPI_FLOAT,0,cr->nc.comm_intra);
614     } 
615     else 
616     {
617         MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_FLOAT,MPI_SUM,cr->mpi_comm_mygroup);
618     }
619 #else
620     int i;
621
622     if (nr > cr->mpb->fbuf_alloc) {
623         cr->mpb->fbuf_alloc = nr;
624         srenew(cr->mpb->fbuf,cr->mpb->fbuf_alloc);
625     }
626     if (cr->nc.bUse) {
627         /* Use two step summing */
628         MPI_Allreduce(r,cr->mpb->fbuf,nr,MPI_FLOAT,MPI_SUM,cr->nc.comm_intra);
629         if (cr->nc.rank_intra == 0) {
630             /* Sum with the buffers reversed */
631             MPI_Allreduce(cr->mpb->fbuf,r,nr,MPI_FLOAT,MPI_SUM, 
632                           cr->nc.comm_inter);
633         }
634         MPI_Bcast(r,nr,MPI_FLOAT,0,cr->nc.comm_intra);
635     } else {
636         MPI_Allreduce(r,cr->mpb->fbuf,nr,MPI_FLOAT,MPI_SUM,
637                       cr->mpi_comm_mygroup);
638         for(i=0; i<nr; i++)
639             r[i] = cr->mpb->fbuf[i];
640     }
641 #endif
642 #endif
643 }
644
645 void gmx_sumi(int nr,int r[],const t_commrec *cr)
646 {
647 #ifndef GMX_MPI
648     gmx_call("gmx_sumi");
649 #else
650 #if defined(MPI_IN_PLACE_EXISTS) || defined(GMX_THREAD_MPI)
651     if (cr->nc.bUse) {
652         /* Use two step summing */
653         if (cr->nc.rank_intra == 0) 
654         {
655             MPI_Reduce(MPI_IN_PLACE,r,nr,MPI_INT,MPI_SUM,0,cr->nc.comm_intra);
656             /* Sum with the buffers reversed */
657             MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_INT,MPI_SUM,cr->nc.comm_inter);
658         }
659         else
660         {
661             /* This is here because of the silly MPI specification
662                 that MPI_IN_PLACE should be put in sendbuf instead of recvbuf */
663             MPI_Reduce(r,NULL,nr,MPI_INT,MPI_SUM,0,cr->nc.comm_intra);
664         }
665         MPI_Bcast(r,nr,MPI_INT,0,cr->nc.comm_intra);
666     } 
667     else 
668     {
669         MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_INT,MPI_SUM,cr->mpi_comm_mygroup);
670     }
671 #else
672     int i;
673
674     if (nr > cr->mpb->ibuf_alloc) {
675         cr->mpb->ibuf_alloc = nr;
676         srenew(cr->mpb->ibuf,cr->mpb->ibuf_alloc);
677     }
678     if (cr->nc.bUse) {
679         /* Use two step summing */
680         MPI_Allreduce(r,cr->mpb->ibuf,nr,MPI_INT,MPI_SUM,cr->nc.comm_intra);
681         if (cr->nc.rank_intra == 0) {
682             /* Sum with the buffers reversed */
683             MPI_Allreduce(cr->mpb->ibuf,r,nr,MPI_INT,MPI_SUM,cr->nc.comm_inter);
684         }
685         MPI_Bcast(r,nr,MPI_INT,0,cr->nc.comm_intra);
686     } else {
687         MPI_Allreduce(r,cr->mpb->ibuf,nr,MPI_INT,MPI_SUM,cr->mpi_comm_mygroup);
688         for(i=0; i<nr; i++)
689             r[i] = cr->mpb->ibuf[i];
690     }
691 #endif
692 #endif
693 }
694
695 void gmx_sumli(int nr,gmx_large_int_t r[],const t_commrec *cr)
696 {
697 #ifndef GMX_MPI
698     gmx_call("gmx_sumli");
699 #else
700 #if defined(MPI_IN_PLACE_EXISTS) || defined(GMX_THREAD_MPI)
701     if (cr->nc.bUse) {
702         /* Use two step summing */
703         if (cr->nc.rank_intra == 0) 
704         {
705             MPI_Reduce(MPI_IN_PLACE,r,nr,GMX_MPI_LARGE_INT,MPI_SUM,0,
706                        cr->nc.comm_intra);
707             /* Sum with the buffers reversed */
708             MPI_Allreduce(MPI_IN_PLACE,r,nr,GMX_MPI_LARGE_INT,MPI_SUM,
709                           cr->nc.comm_inter);
710         }
711         else
712         {
713             /* This is here because of the silly MPI specification
714                 that MPI_IN_PLACE should be put in sendbuf instead of recvbuf */
715             MPI_Reduce(r,NULL,nr,GMX_MPI_LARGE_INT,MPI_SUM,0,cr->nc.comm_intra);
716         }
717         MPI_Bcast(r,nr,GMX_MPI_LARGE_INT,0,cr->nc.comm_intra);
718     } 
719     else 
720     {
721         MPI_Allreduce(MPI_IN_PLACE,r,nr,GMX_MPI_LARGE_INT,MPI_SUM,cr->mpi_comm_mygroup);
722     }
723 #else
724     int i;
725
726     if (nr > cr->mpb->libuf_alloc) {
727         cr->mpb->libuf_alloc = nr;
728         srenew(cr->mpb->libuf,cr->mpb->libuf_alloc);
729     }
730     if (cr->nc.bUse) {
731         /* Use two step summing */
732         MPI_Allreduce(r,cr->mpb->libuf,nr,GMX_MPI_LARGE_INT,MPI_SUM,
733                       cr->nc.comm_intra);
734         if (cr->nc.rank_intra == 0) {
735             /* Sum with the buffers reversed */
736             MPI_Allreduce(cr->mpb->libuf,r,nr,GMX_MPI_LARGE_INT,MPI_SUM,
737                           cr->nc.comm_inter);
738         }
739         MPI_Bcast(r,nr,GMX_MPI_LARGE_INT,0,cr->nc.comm_intra);
740     } else {
741         MPI_Allreduce(r,cr->mpb->libuf,nr,GMX_MPI_LARGE_INT,MPI_SUM,
742                       cr->mpi_comm_mygroup);
743         for(i=0; i<nr; i++)
744             r[i] = cr->mpb->libuf[i];
745     }
746 #endif
747 #endif
748 }
749
750
751
752 #ifdef GMX_MPI
753 void gmx_sumd_comm(int nr,double r[],MPI_Comm mpi_comm)
754 {
755 #if defined(MPI_IN_PLACE_EXISTS) || defined(GMX_THREAD_MPI)
756     MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_DOUBLE,MPI_SUM,mpi_comm);
757 #else
758     /* this function is only used in code that is not performance critical,
759        (during setup, when comm_rec is not the appropriate communication  
760        structure), so this isn't as bad as it looks. */
761     double *buf;
762     int i;
763
764     snew(buf, nr);
765     MPI_Allreduce(r,buf,nr,MPI_DOUBLE,MPI_SUM,mpi_comm);
766     for(i=0; i<nr; i++)
767         r[i] = buf[i];
768     sfree(buf);
769 #endif
770 }
771 #endif
772
773 #ifdef GMX_MPI
774 void gmx_sumf_comm(int nr,float r[],MPI_Comm mpi_comm)
775 {
776 #if defined(MPI_IN_PLACE_EXISTS) || defined(GMX_THREAD_MPI)
777     MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_FLOAT,MPI_SUM,mpi_comm);
778 #else
779     /* this function is only used in code that is not performance critical,
780        (during setup, when comm_rec is not the appropriate communication  
781        structure), so this isn't as bad as it looks. */
782     float *buf;
783     int i;
784
785     snew(buf, nr);
786     MPI_Allreduce(r,buf,nr,MPI_FLOAT,MPI_SUM,mpi_comm);
787     for(i=0; i<nr; i++)
788         r[i] = buf[i];
789     sfree(buf);
790 #endif
791 }
792 #endif
793
794 void gmx_sumd_sim(int nr,double r[],const gmx_multisim_t *ms)
795 {
796 #ifndef GMX_MPI
797   gmx_call("gmx_sumd_sim");
798 #else
799   gmx_sumd_comm(nr,r,ms->mpi_comm_masters);
800 #endif
801 }
802
803 void gmx_sumf_sim(int nr,float r[],const gmx_multisim_t *ms)
804 {
805 #ifndef GMX_MPI
806   gmx_call("gmx_sumf_sim");
807 #else
808   gmx_sumf_comm(nr,r,ms->mpi_comm_masters);
809 #endif
810 }
811
812 void gmx_sumi_sim(int nr,int r[], const gmx_multisim_t *ms)
813 {
814 #ifndef GMX_MPI
815     gmx_call("gmx_sumi_sim");
816 #else
817 #if defined(MPI_IN_PLACE_EXISTS) || defined(GMX_THREAD_MPI)
818     MPI_Allreduce(MPI_IN_PLACE,r,nr,MPI_INT,MPI_SUM,ms->mpi_comm_masters);
819 #else
820     /* this is thread-unsafe, but it will do for now: */
821     int i;
822
823     if (nr > ms->mpb->ibuf_alloc) {
824         ms->mpb->ibuf_alloc = nr;
825         srenew(ms->mpb->ibuf,ms->mpb->ibuf_alloc);
826     }
827     MPI_Allreduce(r,ms->mpb->ibuf,nr,MPI_INT,MPI_SUM,ms->mpi_comm_masters);
828     for(i=0; i<nr; i++)
829         r[i] = ms->mpb->ibuf[i];
830 #endif
831 #endif
832 }
833
834 void gmx_sumli_sim(int nr,gmx_large_int_t r[], const gmx_multisim_t *ms)
835 {
836 #ifndef GMX_MPI
837     gmx_call("gmx_sumli_sim");
838 #else
839 #if defined(MPI_IN_PLACE_EXISTS) || defined(GMX_THREAD_MPI)
840     MPI_Allreduce(MPI_IN_PLACE,r,nr,GMX_MPI_LARGE_INT,MPI_SUM,
841                   ms->mpi_comm_masters);
842 #else
843     /* this is thread-unsafe, but it will do for now: */
844     int i;
845
846     if (nr > ms->mpb->libuf_alloc) {
847         ms->mpb->libuf_alloc = nr;
848         srenew(ms->mpb->libuf,ms->mpb->libuf_alloc);
849     }
850     MPI_Allreduce(r,ms->mpb->libuf,nr,GMX_MPI_LARGE_INT,MPI_SUM,
851                   ms->mpi_comm_masters);
852     for(i=0; i<nr; i++)
853         r[i] = ms->mpb->libuf[i];
854 #endif
855 #endif
856 }
857
858
859 void gmx_finalize_par(void)
860 {
861 #ifndef GMX_MPI
862     /* Compiled without MPI, no MPI finalizing needed */
863     return;
864 #else
865     int initialized,finalized;
866     int ret;
867
868     MPI_Initialized(&initialized);
869     if (!initialized)
870     {
871         return;
872     }
873     /* just as a check; we don't want to finalize twice */
874     MPI_Finalized(&finalized);
875     if (finalized)
876     {
877         return;
878     }
879
880   /* We sync the processes here to try to avoid problems
881    * with buggy MPI implementations that could cause
882    * unfinished processes to terminate.
883    */
884   MPI_Barrier(MPI_COMM_WORLD);
885
886   /*
887   if (DOMAINDECOMP(cr)) {
888     if (cr->npmenodes > 0 || cr->dd->bCartesian) 
889       MPI_Comm_free(&cr->mpi_comm_mygroup);
890     if (cr->dd->bCartesian)
891       MPI_Comm_free(&cr->mpi_comm_mysim);
892   }
893   */
894
895   /* Apparently certain mpich implementations cause problems
896    * with MPI_Finalize. In that case comment out MPI_Finalize.
897    */
898   if (debug)
899     fprintf(debug,"Will call MPI_Finalize now\n");
900
901   ret = MPI_Finalize();
902   if (debug)
903     fprintf(debug,"Return code from MPI_Finalize = %d\n",ret);
904 #endif
905 }
906