ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/mpi_implementation/mpiSimulation_module.F90
Revision: 230
Committed: Thu Jan 9 19:40:38 2003 UTC (22 years, 3 months ago) by chuckv
File size: 15803 byte(s)
Log Message:
More changes to lj_FF.

File Contents

# Content
1 !! MPI support for long range forces using force decomposition
2 !! on a square grid of processors.
3 !! Corresponds to mpiSimunation.cpp for C++
4 !! mpi_module also contains a private interface for mpi f90 routines.
5 !!
6 !! @author Charles F. Vardeman II
7 !! @author Matthew Meineke
8 !! @version $Id: mpiSimulation_module.F90,v 1.2 2003-01-09 19:40:38 chuckv Exp $, $Date: 2003-01-09 19:40:38 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
9
10
11
12
13 module mpiSimulation
14 #ifdef MPI
15 use mpi
16 implicit none
17 PRIVATE
18
19
20 !! PUBLIC Subroutines contained in this module
21 !! gather and scatter are a generic interface
22 !! to gather and scatter routines
23 public :: gather, scatter
24 public :: setupSimParallel
25 public :: replanSimParallel
26 public :: getNcol
27 public :: getNrow
28
29
30 !! PUBLIC Subroutines contained in MPI module
31 public :: mpi_bcast
32 public :: mpi_allreduce
33 public :: mpi_reduce
34 public :: mpi_send
35 public :: mpi_recv
36 public :: mpi_get_processor_name
37 public :: mpi_finalize
38
39 !! PUBLIC mpi variables
40 public :: mpi_comm_world
41 public :: mpi_character
42 public :: mpi_integer
43 public :: mpi_double_precision
44 public :: mpi_sum
45 public :: mpi_max
46 public :: mpi_status_size
47 public :: mpi_any_source
48
49 !! Safety logical to prevent access to ComponetPlan until
50 !! set by C++.
51 logical :: ComponentPlanSet = .false.
52
53
54 !! generic mpi error declaration.
55 integer,public :: mpi_err
56
57 !! Include mpiComponentPlan type. mpiComponentPlan is a
58 !! dual header file for both c and fortran.
59 #define __FORTRAN90
60 #include "mpiComponentPlan.h"
61
62 !! gs_plan contains plans for gather and scatter routines
63 type, public :: gs_plan
64 private
65 type (mpiComponentPlan), pointer :: gsComponentPlan => NULL()
66 integer :: gsPlanSize !! size of this plan (nDim*nComponents)
67 integer :: globalPlanSize !! size of all components in plan
68 integer, dimension(:), pointer :: displs !! Displacements array for mpi indexed from 0.
69 integer, dimension(:), pointer :: counts !! Counts array for mpi indexed from 0.
70 integer :: myPlanComm !! My communicator for this plan
71 integer :: myPlanRank !! My rank in this plan
72 integer :: planNprocs !! Number of processors in this plan
73 end type gs_plan
74
75 ! plans for different decompositions
76 type (gs_plan), public :: plan_row
77 type (gs_plan), public :: plan_row3d
78 type (gs_plan), public :: plan_col
79 type (gs_plan), public :: plan_col3d
80
81 type (mpiComponentPlan), pointer :: simComponentPlan
82
83 ! interface for gather scatter routines
84 !! Generic interface for gather.
85 !! Gathers an local array into row or column array
86 !! Interface provided for integer, double and double
87 !! rank 2 arrays.
88 interface gather
89 module procedure gather_integer
90 module procedure gather_double
91 module procedure gather_double_2d
92 end interface
93
94 !! Generic interface for scatter.
95 !! Scatters a row or column array, adding componets
96 !! and reducing them to a local nComponent array.
97 !! Interface provided for double and double rank=2 arrays.
98 interface scatter
99 module procedure scatter_double
100 module procedure scatter_double_2d
101 end interface
102
103
104 contains
105
106 !! Sets up mpiComponentPlan with structure passed from C++.
107 subroutine setupSimParallel(thisComponentPlan,status)
108 ! Passed Arguments
109 ! integer, intent(inout) :: nDim !! Number of dimensions
110 !! mpiComponentPlan struct from C
111 type (mpiComponentPlan), intent(inout) :: thisComponentPlan
112 integer, intent(out) :: status
113 integer, intnet(out) :: localStatus
114
115 status = 0
116 if (componentPlanSet) then
117 return
118 endif
119
120 componentPlanSet = .true.
121
122 call make_Force_Grid(thisComponentPlan,localStatus)
123 if (localStatus /= 0) then
124 status = -1
125 return
126 endif
127
128 call updateGridComponents(thisComponentPlan,localStatus)
129 if (localStatus /= 0) then
130 status = -1
131 return
132 endif
133
134 !! initialize gather and scatter plans used in this simulation
135 call plan_gather_scatter(1,thisComponentPlan%nComponentsRow,&
136 thisComponentPlan,row_comm,plan_row)
137 call plan_gather_scatter(nDim,thisComponentPlan%nComponentsRow,&
138 thisComponentPlan,row_comm,plan_row3d)
139 call plan_gather_scatter(1,thisComponentPlan%nComponentsColumn,&
140 thisComponentPlan,col_comm,plan_col)
141 call plan_gather_scatter(nDim,thisComponentPlan%nComponentsColumn,&
142 thisComponentPlan,col_comm,plan_col3d)
143
144
145
146 end subroutine setupSimParallel
147
148 subroutine replanSimParallel(thisComponentPlan,status)
149 ! Passed Arguments
150 !! mpiComponentPlan struct from C
151 type (mpiComponentPlan), intent(inout) :: thisComponentPlan
152 integer, intent(out) :: status
153 integer, intnet(out) :: localStatus
154
155 status = 0
156
157 call updateGridComponents(thisComponentPlan,localStatus)
158 if (localStatus /= 0) then
159 status = -1
160 return
161 endif
162
163 !! Unplan Gather Scatter plans
164 call unplan_gather_scatter(plan_row)
165 call unplan_gather_scatter(plan_row3d)
166 call unplan_gather_scatter(plan_col)
167 call unplan_gather_scatter(plan_col3d)
168
169
170 !! initialize gather and scatter plans used in this simulation
171 call plan_gather_scatter(1,thisComponentPlan%nComponentsRow,&
172 thisComponentPlan,row_comm,plan_row)
173 call plan_gather_scatter(nDim,thisComponentPlan%nComponentsRow,&
174 thisComponentPlan,row_comm,plan_row3d)
175 call plan_gather_scatter(1,thisComponentPlan%nComponentsColumn,&
176 thisComponentPlan,col_comm,plan_col)
177 call plan_gather_scatter(nDim,thisComponentPlan%nComponentsColumn,&
178 thisComponentPlan,col_comm,plan_col3d)
179
180
181
182 end subroutine replanSimParallel
183
184 !! Updates number of row and column components for long range forces.
185 subroutine updateGridComponents(thisComponentPlan,status)
186 type (mpiComponentPlan) :: thisComponentPlan !! mpiComponentPlan
187
188 !! Status return
189 !! - 0 Success
190 !! - -1 Failure
191 integer, intent(out) :: status
192 integer :: nComponentsLocal
193 integer :: nComponentsRow = 0
194 integer :: nComponensColumn = 0
195 integer :: mpiErrors
196
197 status = 0
198 if (.not. componentPlanSet) return
199
200 if (thisComponentPlan%myNlocal == 0 ) then
201 status = -1
202 return
203 endif
204
205 nComponentsLocal = thisComponentPlan%myNlocal
206
207 call mpi_allreduce(nComponentsLocal,nComponentsRow,1,mpi_integer,&
208 mpi_sum,thisComponentPlan%rowComm,mpiError)
209 if (mpiErrors /= 0) then
210 status = -1
211 return
212 endif
213
214 call mpi_allreduce(nComponentsLocal,nComponentsColumn,1,mpi_integer, &
215 mpi_sum,thisComponentPlan%columnComm,mpiError)
216 if (mpiErrors /= 0) then
217 status = -1
218 return
219 endif
220
221 thisComponentPlan%nComponentsRow = nComponentsRow
222 thisComponentPlan%nComponentsColumn = nComponentsColumn
223
224
225 end subroutine updateGridComponents
226
227
228 !! Creates a square force decomposition of processors into row and column
229 !! communicators.
230 subroutine make_Force_Grid(thisComponentPlan,status)
231 type (mpiComponentPlan) :: thisComponentPlan
232 integer, intent(out) :: status !! status returns -1 if error
233 integer :: nColumnsMax !! Maximum number of columns
234 integer :: nWorldProcessors !! Total number of processors in World comm.
235 integer :: rowIndex !! Row for this processor.
236 integer :: columnIndex !! Column for this processor.
237 integer :: nRows !! Total number of rows.
238 integer :: nColumns !! Total number of columns.
239 integer :: mpiErrors !! MPI errors.
240 integer :: rowCommunicator !! MPI row communicator.
241 integer :: columnCommunicator
242 integer :: myWorldRank
243 integer :: i
244
245
246 if (.not. ComponentPlanSet) return
247 status = 0
248
249 !! We make a dangerous assumption here that if numberProcessors is
250 !! zero, then we need to get the information from MPI.
251 if (thisComponentPlan%numberProcessors == 0 ) then
252 call mpi_comm_size( MPI_COMM_WORLD, nWorldProcessors,mpiErrors)
253 if ( mpiErrors /= 0 ) then
254 status = -1
255 return
256 endif
257 call mpi_comm_rank( MPI_COMM_WORLD,myWorldRank,mpiErrors)
258 if ( mpiErrors /= 0 ) then
259 status = -1
260 return
261 endif
262
263 else
264 nWorldProcessors = thisComponentPlan%numberProcessors
265 myWorldRank = thisComponentPlan%myRank
266 endif
267
268
269
270
271 nColumnsMax = nint(sqrt(real(nWorldProcessors,kind=dp)))
272
273 do i = 1, nColumnsMax
274 if (mod(nWorldProcessors,i) == 0) nColumns = i
275 end do
276
277 nRows = nWorldProcessors/nColumns
278
279 rowIndex = myWorldRank/nColumns
280 call mpi_comm_split(mpi_comm_world,rowIndex,0,rowCommunicator,mpiError)
281 if ( mpiErrors /= 0 ) then
282 status = -1
283 return
284 endif
285
286 columnIndex = mod(myWorldRank,nColumns)
287 call mpi_comm_split(mpi_comm_world,columnIndex,0,columnCommunicator,mpiError)
288 if ( mpiErrors /= 0 ) then
289 status = -1
290 return
291 endif
292
293 ! Set appropriate components of thisComponentPlan
294 thisComponentPlan%rowComm = rowCommunicator
295 thisComponentPlan%columnComm = columnCommunicator
296 thisComponentPlan%rowIndex = rowIndex
297 thisComponentPlan%columnIndex = columnIndex
298 thisComponentPlan%numberRows = nRows
299 thisComponentPlan%numberColumns = nColumns
300
301
302 end subroutine make_Force_Grid
303
304
305 !! initalizes a gather scatter plan
306 subroutine plan_gather_scatter( nDim,nComponents,thisComponentPlan, &
307 thisComm, this_plan,status)
308 integer, intent(in) :: nDim !! Number of dimensions for gather scatter plan
309 integer, intent(in) :: nComponents
310 type (mpiComponentPlan), intent(in), target :: thisComponentPlan
311 type (gs_plan), intent(out) :: this_plan !! MPI Component Plan
312 integer, intent(in) :: thisComm !! MPI communicator for this plan
313
314 integer :: arraySize !! size to allocate plan for
315 integer, intent(out), optional :: status
316 integer :: ierror
317 integer :: i,junk
318
319 if (present(status)) status = 0
320
321
322 !! Set gsComponetPlan pointer
323 !! to the componet plan we want to use for this gather scatter plan.
324 !! WARNING this could be dangerous since thisComponentPlan was origionally
325 !! allocated in C++ and there is a significant difference between c and
326 !! f95 pointers....
327 gsComponentPlan => thisComponetPlan
328
329 ! Set this plan size for displs array.
330 this_plan%gsPlanSize = nDim * nComponents
331
332 ! Duplicate communicator for this plan
333 call mpi_comm_dup(thisComm,this_plan%myPlanComm,mpi_err)
334 if (mpi_err /= 0) then
335 if (present(status)) status = -1
336 return
337 end if
338 call mpi_comm_rank(this_plan%myPlanComm,this_plan%myPlanRank,mpi_err)
339 if (mpi_err /= 0) then
340 if (present(status)) status = -1
341 return
342 end if
343
344 call mpi_comm_size(this_plan%myPlanComm,this_plan%planNprocs,mpi_err)
345
346 if (mpi_err /= 0) then
347 if (present(status)) status = -1
348 return
349 end if
350
351 !! counts and displacements arrays are indexed from 0 to be compatable
352 !! with MPI arrays.
353 allocate (this_plan%counts(0:this_plan%planNprocs-1),STAT=ierror)
354 if (ierror /= 0) then
355 if (present(status)) status = -1
356 return
357 end if
358
359 allocate (this_plan%displs(0:this_plan%planNprocs-1),STAT=ierror)
360 if (ierror /= 0) then
361 if (present(status)) status = -1
362 return
363 end if
364
365 !! gather all the local sizes into a size # processors array.
366 call mpi_allgather(gs_plan%gsPlanSize,1,mpi_integer,this_plan%counts, &
367 1,mpi_integer,comm,mpi_err)
368
369 if (mpi_err /= 0) then
370 if (present(status)) status = -1
371 return
372 end if
373
374
375 !! figure out the total number of particles in this plan
376 this_plan%globalPlanSize = sum(this_plan%counts)
377
378
379 !! initialize plan displacements.
380 this_plan%displs(0) = 0
381 do i = 1, this_plan%planNprocs - 1,1
382 this_plan%displs(i) = this_plan%displs(i-1) + this_plan%counts(i-1)
383 end do
384
385 end subroutine plan_gather_scatter
386
387
388 subroutine unplan_gather_scatter(this_plan)
389 type (gs_plan), intent(inout) :: this_plan
390
391
392 this_plan%gsComponentPlan => null()
393 call mpi_comm_free(this_plan%comm,mpi_err)
394
395 deallocate(this_plan%counts)
396 deallocate(this_plan%displs)
397
398 end subroutine unplan_gather_scatter
399
400 subroutine gather_integer( sbuffer, rbuffer, this_plan, status)
401
402 type (gs_plan), intent(in) :: this_plan
403 integer, dimension(:), intent(in) :: sbuffer
404 integer, dimension(:), intent(in) :: rbuffer
405 integer :: noffset
406 integer, intent(out), optional :: status
407
408 if (present(status)) status = 0
409 noffset = this_plan%displs(this_plan%myPlanRank)
410
411 call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_integer, &
412 rbuffer,this_plan%counts,this_plan%displs,mpi_integer, &
413 this_plan%myPlanComm, mpi_err)
414
415 if (mpi_err /= 0) then
416 if (present(status)) status = -1
417 endif
418
419 end subroutine gather_integer
420
421 subroutine gather_double( sbuffer, rbuffer, this_plan, status)
422
423 type (gs_plan), intent(in) :: this_plan
424 real( kind = DP ), dimension(:), intent(in) :: sbuffer
425 real( kind = DP ), dimension(:), intent(in) :: rbuffer
426 integer :: noffset
427 integer, intent(out), optional :: status
428
429 if (present(status)) status = 0
430 noffset = this_plan%displs(this_plan%me)
431
432 call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_double_precision, &
433 rbuffer,this_plan%counts,this_plan%displs,mpi_double_precision, &
434 this_plan%myPlanComm, mpi_err)
435
436 if (mpi_err /= 0) then
437 if (present(status)) status = -1
438 endif
439
440 end subroutine gather_double
441
442 subroutine gather_double_2d( sbuffer, rbuffer, this_plan, status)
443
444 type (gs_plan), intent(in) :: this_plan
445 real( kind = DP ), dimension(:,:), intent(in) :: sbuffer
446 real( kind = DP ), dimension(:,:), intent(in) :: rbuffer
447 integer :: noffset,i,ierror
448 integer, intent(out), optional :: status
449
450 external mpi_allgatherv
451
452 if (present(status)) status = 0
453
454 ! noffset = this_plan%displs(this_plan%me)
455
456 call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_double_precision, &
457 rbuffer,this_plan%counts,this_plan%displs,mpi_double_precision, &
458 this_plan%myPlanComm, mpi_err)
459
460 if (mpi_err /= 0) then
461 if (present(status)) status = -1
462 endif
463
464 end subroutine gather_double_2d
465
466 subroutine scatter_double( sbuffer, rbuffer, this_plan, status)
467
468 type (gs_plan), intent(in) :: this_plan
469 real( kind = DP ), dimension(:), intent(in) :: sbuffer
470 real( kind = DP ), dimension(:), intent(in) :: rbuffer
471 integer, intent(out), optional :: status
472 external mpi_reduce_scatter
473
474 if (present(status)) status = 0
475
476 call mpi_reduce_scatter(sbuffer,rbuffer, this_plan%counts, &
477 mpi_double_precision, MPI_SUM, this_plan%myPlanComm, mpi_err)
478
479 if (mpi_err /= 0) then
480 if (present(status)) status = -1
481 endif
482
483 end subroutine scatter_double
484
485 subroutine scatter_double_2d( sbuffer, rbuffer, this_plan, status)
486
487 type (gs_plan), intent(in) :: this_plan
488 real( kind = DP ), dimension(:,:), intent(in) :: sbuffer
489 real( kind = DP ), dimension(:,:), intent(in) :: rbuffer
490 integer, intent(out), optional :: status
491 external mpi_reduce_scatter
492
493 if (present(status)) status = 0
494 call mpi_reduce_scatter(sbuffer,rbuffer, this_plan%counts, &
495 mpi_double_precision, MPI_SUM, this_plan%myPlanComm, mpi_err)
496
497 if (mpi_err /= 0) then
498 if (present(status)) status = -1
499 endif
500
501 end subroutine scatter_double_2d
502
503
504 function getNcol(thisplan) result(ncol)
505 type (gsPlan) :: thisplan
506 integer :: ncol
507 ncol = thisplan%gsComponentPlan%nComponentsCol
508 end function getNcol
509
510 function getNrow(thisplan) result(ncol)
511 type (gsPlan) :: thisplan
512 integer :: ncol
513 ncol = thisplan%gsComponentPlan%nComponentsrow
514 end function getNrow
515
516
517
518
519 #endif
520 end module mpiSimulation
521