ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/mpiSimulation.F90
Revision: 210
Committed: Fri Dec 13 18:06:58 2002 UTC (21 years, 8 months ago) by chuckv
File size: 14206 byte(s)
Log Message:
Finished updating mpiSimulation to work with c++.

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