ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/mpiSimulation.F90
Revision: 211
Committed: Fri Dec 13 19:15:57 2002 UTC (21 years, 8 months ago) by chuckv
File size: 15419 byte(s)
Log Message:
Added replan of gather-scatter routines for future load balancing.

File Contents

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