ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/mpi_implementation/mpiSimulation_module.F90
Revision: 259
Committed: Thu Jan 30 22:29:58 2003 UTC (22 years, 3 months ago) by chuckv
File size: 19459 byte(s)
Log Message:
MPI is Brooooooken...

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.6 2003-01-30 22:29:58 chuckv Exp $, $Date: 2003-01-30 22:29:58 $, $Name: not supported by cvs2svn $, $Revision: 1.6 $
9
10
11
12
13 module mpiSimulation
14 use definitions
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 public :: isMPISimSet
29
30
31 !! PUBLIC Subroutines contained in MPI module
32 public :: mpi_bcast
33 public :: mpi_allreduce
34 public :: mpi_reduce
35 public :: mpi_send
36 public :: mpi_recv
37 public :: mpi_get_processor_name
38 public :: mpi_finalize
39
40 !! PUBLIC mpi variables
41 public :: mpi_comm_world
42 public :: mpi_character
43 public :: mpi_integer
44 public :: mpi_double_precision
45 public :: mpi_sum
46 public :: mpi_max
47 public :: mpi_status_size
48 public :: mpi_any_source
49
50 !! Safety logical to prevent access to ComponetPlan until
51 !! set by C++.
52 logical :: ComponentPlanSet = .false.
53
54
55 !! generic mpi error declaration.
56 integer,public :: mpi_err
57
58
59
60 !! Include mpiComponentPlan type. mpiComponentPlan is a
61 !! dual header file for both c and fortran.
62 #define __FORTRAN90
63 #include "mpiComponentPlan.h"
64
65
66
67 !! Tags used during force loop for parallel simulation
68 integer, allocatable, dimension(:) :: tagLocal
69 integer, public, allocatable, dimension(:) :: tagRow
70 integer ,public, allocatable, dimension(:) :: tagColumn
71
72 !! Logical set true if mpiSimulation has been initialized
73 logical :: isSimSet = .false.
74
75
76 type (mpiComponentPlan) :: mpiSim
77
78 !! gs_plan contains plans for gather and scatter routines
79 type, public :: gs_plan
80 private
81 type (mpiComponentPlan), pointer :: gsComponentPlan => NULL()
82 integer :: gsPlanSize !! size of this plan (nDim*nComponents)
83 integer :: globalPlanSize !! size of all components in plan
84 integer, dimension(:), pointer :: displs !! Displacements array for mpi indexed from 0.
85 integer, dimension(:), pointer :: counts !! Counts array for mpi indexed from 0.
86 integer :: myPlanComm !! My communicator for this plan
87 integer :: myPlanRank !! My rank in this plan
88 integer :: planNprocs !! Number of processors in this plan
89 end type gs_plan
90
91 ! plans for different decompositions
92 type (gs_plan), public :: plan_row
93 type (gs_plan), public :: plan_row3d
94 type (gs_plan), public :: plan_col
95 type (gs_plan), public :: plan_col3d
96
97 type (mpiComponentPlan), pointer :: simComponentPlan
98
99 ! interface for gather scatter routines
100 !! Generic interface for gather.
101 !! Gathers an local array into row or column array
102 !! Interface provided for integer, double and double
103 !! rank 2 arrays.
104 interface gather
105 module procedure gather_integer
106 module procedure gather_double
107 module procedure gather_double_2d
108 end interface
109
110 !! Generic interface for scatter.
111 !! Scatters a row or column array, adding componets
112 !! and reducing them to a local nComponent array.
113 !! Interface provided for double and double rank=2 arrays.
114 interface scatter
115 module procedure scatter_double
116 module procedure scatter_double_2d
117 end interface
118
119
120
121 contains
122
123 !! Sets up mpiComponentPlan with structure passed from C++.
124 subroutine setupSimParallel(thisComponentPlan,ntags,tags,status)
125 ! Passed Arguments
126 ! integer, intent(inout) :: nDim !! Number of dimensions
127 !! mpiComponentPlan struct from C
128 type (mpiComponentPlan), intent(inout) :: thisComponentPlan
129 !! Number of tags passed, nlocal
130 integer, intent(in) :: ntags
131 !! Result status, 0 = normal, -1 = error
132 integer, intent(out) :: status
133 integer :: localStatus
134 !! Global reference tag for local particles
135 integer, dimension(ntags),intent(inout) :: tags
136
137
138 status = 0
139 if (componentPlanSet) then
140 return
141 endif
142 componentPlanSet = .true.
143
144 !! copy c component plan to fortran
145 mpiSim = thisComponentPlan
146 write(*,*) "Seting up simParallel"
147 if (mpiSim%myNode == 0) then
148 write(*,*) "SetupSimParallel: writing component plan"
149
150 write(*,*) "nMolGlobal: ", mpiSim%nMolGlobal
151 write(*,*) "nAtomsGlobal: ", mpiSim%nAtomsGlobal
152 write(*,*) "nBondGlobal: ", mpiSim%nBondsGlobal
153 write(*,*) "nTorsionsGlobal: ", mpiSim%nTorsionsGlobal
154 write(*,*) "nSRIGlobal: ", mpiSim%nSRIGlobal
155 write(*,*) "myMolStart: ", mpiSim%myMolStart
156 write(*,*) "myMolEnd: ", mpiSim%myMolEnd
157 write(*,*) "myMol: ", mpiSim%myMol
158 write(*,*) "myNlocal: ", mpiSim%myNlocal
159 write(*,*) "myNode: ", mpiSim%myNode
160 write(*,*) "numberProcessors: ", mpiSim%numberProcessors
161 write(*,*) "rowComm: ", mpiSim%rowComm
162 write(*,*) "columnComm: ", mpiSim%columnComm
163 write(*,*) "numberRows: ", mpiSim%numberRows
164 write(*,*) "numberColumns: ", mpiSim%numberColumns
165 write(*,*) "nComponentsRow: ", mpiSim%nComponentsRow
166 write(*,*) "nComponentsColumn: ", mpiSim%nComponentsColumn
167 write(*,*) "rowIndex: ", mpiSim%rowIndex
168 write(*,*) "columnIndex: ", mpiSim%columnIndex
169 endif
170
171
172 call make_Force_Grid(mpiSim,localStatus)
173 if (localStatus /= 0) then
174 write(default_error,*) "Error creating force grid"
175 status = -1
176 return
177 endif
178
179 call updateGridComponents(mpiSim,localStatus)
180 if (localStatus /= 0) then
181 write(default_error,*) "Error updating grid components"
182 status = -1
183 return
184 endif
185
186
187 !! initialize gather and scatter plans used in this simulation
188 call plan_gather_scatter(1,mpiSim%nComponentsRow,&
189 mpiSim,mpiSim%rowComm,plan_row)
190 call plan_gather_scatter(nDim,mpiSim%nComponentsRow,&
191 mpiSim,mpiSim%rowComm,plan_row3d)
192 call plan_gather_scatter(1,mpiSim%nComponentsColumn,&
193 mpiSim,mpiSim%columnComm,plan_col)
194 call plan_gather_scatter(nDim,mpiSim%nComponentsColumn,&
195 mpiSim,mpiSim%columnComm,plan_col3d)
196
197 ! Initialize tags
198 call setTags(tags,localStatus)
199 if (localStatus /= 0) then
200 status = -1
201 return
202 endif
203 isSimSet = .true.
204 end subroutine setupSimParallel
205
206 subroutine replanSimParallel(thisComponentPlan,status)
207 ! Passed Arguments
208 !! mpiComponentPlan struct from C
209 type (mpiComponentPlan), intent(inout) :: thisComponentPlan
210 integer, intent(out) :: status
211 integer :: localStatus
212 integer :: mpierror
213 status = 0
214
215 call updateGridComponents(thisComponentPlan,localStatus)
216 if (localStatus /= 0) then
217 status = -1
218 return
219 endif
220
221 !! Unplan Gather Scatter plans
222 call unplan_gather_scatter(plan_row)
223 call unplan_gather_scatter(plan_row3d)
224 call unplan_gather_scatter(plan_col)
225 call unplan_gather_scatter(plan_col3d)
226
227
228 !! initialize gather and scatter plans used in this simulation
229 call plan_gather_scatter(1,thisComponentPlan%nComponentsRow,&
230 thisComponentPlan,thisComponentPlan%rowComm,plan_row)
231 call plan_gather_scatter(nDim,thisComponentPlan%nComponentsRow,&
232 thisComponentPlan,thisComponentPlan%rowComm,plan_row3d)
233 call plan_gather_scatter(1,thisComponentPlan%nComponentsColumn,&
234 thisComponentPlan,thisComponentPlan%columnComm,plan_col)
235 call plan_gather_scatter(nDim,thisComponentPlan%nComponentsColumn,&
236 thisComponentPlan,thisComponentPlan%rowComm,plan_col3d)
237
238
239
240 end subroutine replanSimParallel
241
242 !! Updates number of row and column components for long range forces.
243 subroutine updateGridComponents(thisComponentPlan,status)
244 type (mpiComponentPlan) :: thisComponentPlan !! mpiComponentPlan
245
246 !! Status return
247 !! - 0 Success
248 !! - -1 Failure
249 integer, intent(out) :: status
250 integer :: nComponentsLocal
251 integer :: nComponentsRow = 0
252 integer :: nComponentsColumn = 0
253 integer :: mpiErrors
254
255 status = 0
256 if (.not. componentPlanSet) return
257
258 if (thisComponentPlan%myNlocal == 0 ) then
259 status = -1
260 return
261 endif
262
263 nComponentsLocal = thisComponentPlan%myNlocal
264
265 call mpi_allreduce(nComponentsLocal,nComponentsRow,1,mpi_integer,&
266 mpi_sum,thisComponentPlan%rowComm,mpiErrors)
267 if (mpiErrors /= 0) then
268 status = -1
269 return
270 endif
271
272 call mpi_allreduce(nComponentsLocal,nComponentsColumn,1,mpi_integer, &
273 mpi_sum,thisComponentPlan%columnComm,mpiErrors)
274 if (mpiErrors /= 0) then
275 status = -1
276 return
277 endif
278
279 thisComponentPlan%nComponentsRow = nComponentsRow
280 thisComponentPlan%nComponentsColumn = nComponentsColumn
281
282
283 end subroutine updateGridComponents
284
285
286 !! Creates a square force decomposition of processors into row and column
287 !! communicators.
288 subroutine make_Force_Grid(thisComponentPlan,status)
289 type (mpiComponentPlan) :: thisComponentPlan
290 integer, intent(out) :: status !! status returns -1 if error
291 integer :: nColumnsMax !! Maximum number of columns
292 integer :: nWorldProcessors !! Total number of processors in World comm.
293 integer :: rowIndex !! Row for this processor.
294 integer :: columnIndex !! Column for this processor.
295 integer :: nRows !! Total number of rows.
296 integer :: nColumns !! Total number of columns.
297 integer :: mpiErrors !! MPI errors.
298 integer :: rowCommunicator !! MPI row communicator.
299 integer :: columnCommunicator
300 integer :: myWorldRank
301 integer :: i
302
303
304 if (.not. ComponentPlanSet) return
305 status = 0
306
307 !! We make a dangerous assumption here that if numberProcessors is
308 !! zero, then we need to get the information from MPI.
309 if (thisComponentPlan%numberProcessors == 0 ) then
310 call mpi_comm_size( MPI_COMM_WORLD, nWorldProcessors,mpiErrors)
311 if ( mpiErrors /= 0 ) then
312 status = -1
313 return
314 endif
315 call mpi_comm_rank( MPI_COMM_WORLD,myWorldRank,mpiErrors)
316 if ( mpiErrors /= 0 ) then
317 status = -1
318 return
319 endif
320
321 else
322 nWorldProcessors = thisComponentPlan%numberProcessors
323 myWorldRank = thisComponentPlan%myNode
324 endif
325
326
327
328
329 nColumnsMax = nint(sqrt(real(nWorldProcessors,kind=dp)))
330
331 do i = 1, nColumnsMax
332 if (mod(nWorldProcessors,i) == 0) nColumns = i
333 end do
334
335 nRows = nWorldProcessors/nColumns
336
337 rowIndex = myWorldRank/nColumns
338
339
340
341 call mpi_comm_split(mpi_comm_world,rowIndex,0,rowCommunicator,mpiErrors)
342 if ( mpiErrors /= 0 ) then
343 write(default_error,*) "MPI comm split failed at row communicator"
344 status = -1
345 return
346 endif
347
348 columnIndex = mod(myWorldRank,nColumns)
349 call mpi_comm_split(mpi_comm_world,columnIndex,0,columnCommunicator,mpiErrors)
350 if ( mpiErrors /= 0 ) then
351 write(default_error,*) "MPI comm split faild at columnCommunicator"
352 status = -1
353 return
354 endif
355
356 ! Set appropriate components of thisComponentPlan
357 thisComponentPlan%rowComm = rowCommunicator
358 thisComponentPlan%columnComm = columnCommunicator
359 thisComponentPlan%rowIndex = rowIndex
360 thisComponentPlan%columnIndex = columnIndex
361 thisComponentPlan%numberRows = nRows
362 thisComponentPlan%numberColumns = nColumns
363
364
365 end subroutine make_Force_Grid
366
367
368 !! initalizes a gather scatter plan
369 subroutine plan_gather_scatter( nDim,nComponents,thisComponentPlan, &
370 thisComm, this_plan,status)
371 integer, intent(in) :: nDim !! Number of dimensions for gather scatter plan
372 integer, intent(in) :: nComponents
373 type (mpiComponentPlan), intent(in), target :: thisComponentPlan
374 type (gs_plan), intent(out) :: this_plan !! MPI Component Plan
375 integer, intent(in) :: thisComm !! MPI communicator for this plan
376
377 integer :: arraySize !! size to allocate plan for
378 integer, intent(out), optional :: status
379 integer :: ierror
380 integer :: i,junk
381
382 if (present(status)) status = 0
383
384
385 !! Set gsComponetPlan pointer
386 !! to the componet plan we want to use for this gather scatter plan.
387 !! WARNING this could be dangerous since thisComponentPlan was origionally
388 !! allocated in C++ and there is a significant difference between c and
389 !! f95 pointers....
390 this_plan%gsComponentPlan => thisComponentPlan
391
392 ! Set this plan size for displs array.
393 this_plan%gsPlanSize = nDim * nComponents
394
395 ! Duplicate communicator for this plan
396 call mpi_comm_dup(thisComm,this_plan%myPlanComm,mpi_err)
397 if (mpi_err /= 0) then
398 if (present(status)) status = -1
399 return
400 end if
401 call mpi_comm_rank(this_plan%myPlanComm,this_plan%myPlanRank,mpi_err)
402 if (mpi_err /= 0) then
403 if (present(status)) status = -1
404 return
405 end if
406
407 call mpi_comm_size(this_plan%myPlanComm,this_plan%planNprocs,mpi_err)
408
409 if (mpi_err /= 0) then
410 if (present(status)) status = -1
411 return
412 end if
413
414 !! counts and displacements arrays are indexed from 0 to be compatable
415 !! with MPI arrays.
416 allocate (this_plan%counts(0:this_plan%planNprocs-1),STAT=ierror)
417 if (ierror /= 0) then
418 if (present(status)) status = -1
419 return
420 end if
421
422 allocate (this_plan%displs(0:this_plan%planNprocs-1),STAT=ierror)
423 if (ierror /= 0) then
424 if (present(status)) status = -1
425 return
426 end if
427
428 !! gather all the local sizes into a size # processors array.
429 call mpi_allgather(this_plan%gsPlanSize,1,mpi_integer,this_plan%counts, &
430 1,mpi_integer,thisComm,mpi_err)
431
432 if (mpi_err /= 0) then
433 if (present(status)) status = -1
434 return
435 end if
436
437
438 !! figure out the total number of particles in this plan
439 this_plan%globalPlanSize = sum(this_plan%counts)
440
441
442 !! initialize plan displacements.
443 this_plan%displs(0) = 0
444 do i = 1, this_plan%planNprocs - 1,1
445 this_plan%displs(i) = this_plan%displs(i-1) + this_plan%counts(i-1)
446 end do
447
448 end subroutine plan_gather_scatter
449
450
451 subroutine unplan_gather_scatter(this_plan)
452 type (gs_plan), intent(inout) :: this_plan
453
454
455 this_plan%gsComponentPlan => null()
456 call mpi_comm_free(this_plan%myPlanComm,mpi_err)
457
458 deallocate(this_plan%counts)
459 deallocate(this_plan%displs)
460
461 end subroutine unplan_gather_scatter
462
463 subroutine gather_integer( sbuffer, rbuffer, this_plan, status)
464
465 type (gs_plan), intent(in) :: this_plan
466 integer, dimension(:), intent(in) :: sbuffer
467 integer, dimension(:), intent(in) :: rbuffer
468 integer :: noffset
469 integer, intent(out), optional :: status
470 integer :: i
471
472
473 if (present(status)) status = 0
474 noffset = this_plan%displs(this_plan%myPlanRank)
475
476 call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_integer, &
477 rbuffer,this_plan%counts,this_plan%displs,mpi_integer, &
478 this_plan%myPlanComm, mpi_err)
479
480 if (mpi_err /= 0) then
481 if (present(status)) status = -1
482 endif
483
484 end subroutine gather_integer
485
486 subroutine gather_double( sbuffer, rbuffer, this_plan, status)
487
488 type (gs_plan), intent(in) :: this_plan
489 real( kind = DP ), dimension(:), intent(in) :: sbuffer
490 real( kind = DP ), dimension(:), intent(in) :: rbuffer
491 integer :: noffset
492 integer, intent(out), optional :: status
493
494 if (present(status)) status = 0
495 noffset = this_plan%displs(this_plan%myPlanRank)
496
497 call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_double_precision, &
498 rbuffer,this_plan%counts,this_plan%displs,mpi_double_precision, &
499 this_plan%myPlanComm, mpi_err)
500
501 if (mpi_err /= 0) then
502 if (present(status)) status = -1
503 endif
504
505 end subroutine gather_double
506
507 subroutine gather_double_2d( sbuffer, rbuffer, this_plan, status)
508
509 type (gs_plan), intent(in) :: this_plan
510 real( kind = DP ), dimension(:,:), intent(in) :: sbuffer
511 real( kind = DP ), dimension(:,:), intent(in) :: rbuffer
512 integer :: noffset,i,ierror
513 integer, intent(out), optional :: status
514
515 external mpi_allgatherv
516
517 if (present(status)) status = 0
518
519 ! noffset = this_plan%displs(this_plan%me)
520
521 call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_double_precision, &
522 rbuffer,this_plan%counts,this_plan%displs,mpi_double_precision, &
523 this_plan%myPlanComm, mpi_err)
524
525 if (mpi_err /= 0) then
526 if (present(status)) status = -1
527 endif
528
529 end subroutine gather_double_2d
530
531 subroutine scatter_double( sbuffer, rbuffer, this_plan, status)
532
533 type (gs_plan), intent(in) :: this_plan
534 real( kind = DP ), dimension(:), intent(in) :: sbuffer
535 real( kind = DP ), dimension(:), intent(in) :: rbuffer
536 integer, intent(out), optional :: status
537 external mpi_reduce_scatter
538
539 if (present(status)) status = 0
540
541 call mpi_reduce_scatter(sbuffer,rbuffer, this_plan%counts, &
542 mpi_double_precision, MPI_SUM, this_plan%myPlanComm, mpi_err)
543
544 if (mpi_err /= 0) then
545 if (present(status)) status = -1
546 endif
547
548 end subroutine scatter_double
549
550 subroutine scatter_double_2d( sbuffer, rbuffer, this_plan, status)
551
552 type (gs_plan), intent(in) :: this_plan
553 real( kind = DP ), dimension(:,:), intent(in) :: sbuffer
554 real( kind = DP ), dimension(:,:), intent(in) :: rbuffer
555 integer, intent(out), optional :: status
556 external mpi_reduce_scatter
557
558 if (present(status)) status = 0
559 call mpi_reduce_scatter(sbuffer,rbuffer, this_plan%counts, &
560 mpi_double_precision, MPI_SUM, this_plan%myPlanComm, mpi_err)
561
562 if (mpi_err /= 0) then
563 if (present(status)) status = -1
564 endif
565
566 end subroutine scatter_double_2d
567
568
569 subroutine setTags(tags,status)
570 integer, dimension(:) :: tags
571 integer :: status
572
573 integer :: alloc_stat
574
575 integer :: ncol
576 integer :: nrow
577
578 status = 0
579 ! allocate row arrays
580 nrow = getNrow(plan_row)
581 ncol = getNcol(plan_col)
582
583 if (.not. allocated(tagRow)) then
584 allocate(tagRow(nrow),STAT=alloc_stat)
585 if (alloc_stat /= 0 ) then
586 status = -1
587 return
588 endif
589 else
590 deallocate(tagRow)
591 allocate(tagRow(nrow),STAT=alloc_stat)
592 if (alloc_stat /= 0 ) then
593 status = -1
594 return
595 endif
596
597 endif
598 ! allocate column arrays
599 if (.not. allocated(tagColumn)) then
600 allocate(tagColumn(ncol),STAT=alloc_stat)
601 if (alloc_stat /= 0 ) then
602 status = -1
603 return
604 endif
605 else
606 deallocate(tagColumn)
607 allocate(tagColumn(ncol),STAT=alloc_stat)
608 if (alloc_stat /= 0 ) then
609 status = -1
610 return
611 endif
612 endif
613
614 call gather(tags,tagRow,plan_row)
615 call gather(tags,tagColumn,plan_col)
616
617
618 end subroutine setTags
619
620 pure function getNcol(thisplan) result(ncol)
621 type (gs_plan), intent(in) :: thisplan
622 integer :: ncol
623 ncol = thisplan%gsComponentPlan%nComponentsColumn
624 end function getNcol
625
626 pure function getNrow(thisplan) result(ncol)
627 type (gs_plan), intent(in) :: thisplan
628 integer :: ncol
629 ncol = thisplan%gsComponentPlan%nComponentsrow
630 end function getNrow
631
632 function isMPISimSet() result(isthisSimSet)
633 logical :: isthisSimSet
634 if (isSimSet) then
635 isthisSimSet = .true.
636 else
637 isthisSimSet = .false.
638 endif
639 end function isMPISimSet
640
641
642
643
644 end module mpiSimulation
645