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