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

# User Rev Content
1 chuckv 215 !! 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 230 !! @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 chuckv 215
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 chuckv 230 public :: getNcol
27     public :: getNrow
28    
29    
30 chuckv 215 !! 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 chuckv 230
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 chuckv 215
517 chuckv 230
518    
519 chuckv 215 #endif
520     end module mpiSimulation
521