ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/mpi_implementation/mpiSimulation_module.F90
Revision: 215
Committed: Thu Dec 19 21:59:51 2002 UTC (22 years, 4 months ago) by chuckv
File size: 15422 byte(s)
Log Message:
+ added lennard-jones force module and corresponding class.
+ created forceFactory directory.

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     !! @version $Id: mpiSimulation_module.F90,v 1.1 2002-12-19 21:59:51 chuckv Exp $, $Date: 2002-12-19 21:59:51 $, $Name: not supported by cvs2svn $, $Revision: 1.1 $
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 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     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     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     subroutine setupSimParallel(thisComponentPlan,status)
104     ! Passed Arguments
105     ! integer, intent(inout) :: nDim !! Number of dimensions
106     !! mpiComponentPlan struct from C
107     type (mpiComponentPlan), intent(inout) :: thisComponentPlan
108     integer, intent(out) :: status
109     integer, intnet(out) :: localStatus
110    
111     status = 0
112     if (componentPlanSet) then
113     return
114     endif
115    
116     componentPlanSet = .true.
117    
118     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    
130     !! initialize gather and scatter plans used in this simulation
131     call plan_gather_scatter(1,thisComponentPlan%nComponentsRow,&
132     thisComponentPlan,row_comm,plan_row)
133     call plan_gather_scatter(nDim,thisComponentPlan%nComponentsRow,&
134     thisComponentPlan,row_comm,plan_row3d)
135     call plan_gather_scatter(1,thisComponentPlan%nComponentsColumn,&
136     thisComponentPlan,col_comm,plan_col)
137     call plan_gather_scatter(nDim,thisComponentPlan%nComponentsColumn,&
138     thisComponentPlan,col_comm,plan_col3d)
139    
140    
141    
142     end subroutine setupSimParallel
143    
144     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     !! 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     subroutine plan_gather_scatter( nDim,nComponents,thisComponentPlan, &
303     thisComm, this_plan,status)
304     integer, intent(in) :: nDim !! Number of dimensions for gather scatter plan
305     integer, intent(in) :: nComponents
306     type (mpiComponentPlan), intent(in), target :: thisComponentPlan
307     type (gs_plan), intent(out) :: this_plan !! MPI Component Plan
308     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     integer :: ierror
313     integer :: i,junk
314    
315     if (present(status)) status = 0
316    
317    
318     !! 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    
325     ! Set this plan size for displs array.
326     this_plan%gsPlanSize = nDim * nComponents
327    
328     ! 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    
340     call mpi_comm_size(this_plan%myPlanComm,this_plan%planNprocs,mpi_err)
341    
342     if (mpi_err /= 0) then
343     if (present(status)) status = -1
344     return
345     end if
346    
347     !! 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     if (ierror /= 0) then
351     if (present(status)) status = -1
352     return
353     end if
354    
355     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     end if
360    
361     !! 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    
365     if (mpi_err /= 0) then
366     if (present(status)) status = -1
367     return
368     end if
369    
370    
371     !! figure out the total number of particles in this plan
372     this_plan%globalPlanSize = sum(this_plan%counts)
373    
374    
375     !! initialize plan displacements.
376     this_plan%displs(0) = 0
377     do i = 1, this_plan%planNprocs - 1,1
378     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    
387    
388     this_plan%gsComponentPlan => null()
389     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     noffset = this_plan%displs(this_plan%myPlanRank)
406    
407     call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_integer, &
408     rbuffer,this_plan%counts,this_plan%displs,mpi_integer, &
409     this_plan%myPlanComm, mpi_err)
410    
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     call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_double_precision, &
429     rbuffer,this_plan%counts,this_plan%displs,mpi_double_precision, &
430     this_plan%myPlanComm, mpi_err)
431    
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     call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_double_precision, &
453     rbuffer,this_plan%counts,this_plan%displs,mpi_double_precision, &
454     this_plan%myPlanComm, mpi_err)
455    
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     mpi_double_precision, MPI_SUM, this_plan%myPlanComm, mpi_err)
474    
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     mpi_double_precision, MPI_SUM, this_plan%myPlanComm, mpi_err)
492    
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