source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/lib/mct/mct/m_GlobalSegMap.F90 @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 4 years ago
  • Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from ​https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818
File size: 74.3 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2!    Math and Computer Science Division, Argonne National Laboratory   !
3!-----------------------------------------------------------------------
4! CVS m_GlobalSegMap.F90,v 1.56 2009-03-17 16:51:49 jacob Exp
5! CVS MCT_2_8_0
6!BOP -------------------------------------------------------------------
7!
8! !MODULE: m_GlobalSegMap - a nontrivial 1-D decomposition of an array.
9!
10! !DESCRIPTION:
11! Consider the problem of the 1-dimensional decomposition of an array
12! across multiple processes.  If each process owns only one contiguous
13! segment, then the {\tt GlobalMap} (see {\tt m\_GlobalMap} or details)
14! is sufficient to describe the decomposition.  If, however, each 
15! process owns multiple, non-adjacent segments of the array, a more
16! sophisticated approach is needed.   The {\tt GlobalSegMap} data type
17! allows one to describe a one-dimensional decomposition of an array
18! with each process owning multiple, non-adjacent segments of the array.
19!
20! In the current implementation of the {\tt GlobalSegMap}, there is no
21! santity check to guarantee that
22!$${\tt GlobalSegMap\%gsize} = \sum_{{\tt i}=1}^{\tt ngseg}
23! {\tt GlobalSegMap\%length(i)} . $$
24! The reason we have not implemented such a check is to allow the user
25! to use the {\tt GlobalSegMap} type to support decompositions of both
26! {\em haloed} and {\em masked} data.
27!
28! !INTERFACE:
29
30 module m_GlobalSegMap
31
32      implicit none
33
34      private   ! except
35
36! !PUBLIC MEMBER FUNCTIONS:
37
38      public :: GlobalSegMap    ! The class data structure
39      public :: init            ! Create
40      public :: clean           ! Destroy
41      public :: comp_id         ! Return component ID number
42      public :: gsize           ! Return global vector size (excl. halos)
43      public :: GlobalStorage   ! Return total number of points in map,
44                                ! including halo points (if present).
45      public :: ProcessStorage  ! Return local storage on a given process.
46      public :: OrderedPoints   ! Return grid points of a given process in
47                                ! MCT-assumed order.
48      public :: lsize           ! Return local--that is, on-process--storage
49                                ! size (incl. halos)
50      public :: ngseg           ! Return global number of segments
51      public :: nlseg           ! Return local number of segments
52      public :: max_nlseg       ! Return max local number of segments
53      public :: active_pes      ! Return number of pes with at least 1
54                                ! datum, and if requested, a list of them.
55      public :: peLocs          ! Given an input list of point indices,
56                                ! return its (unique) process ID.
57      public :: haloed          ! Is the input GlobalSegMap haloed?
58      public :: rank            ! Rank which process owns a datum
59      public :: Sort            ! compute index permutation to re-order
60                                ! GlobalSegMap%start, GlobalSegMap%length,
61                                ! and GlobalSegMap%pe_loc
62      public :: Permute         ! apply index permutation to re-order
63                                ! GlobalSegMap%start, GlobalSegMap%length,
64                                ! and GlobalSegMap%pe_loc
65      public :: SortPermute     ! compute index permutation and apply it to
66                                ! re-order the GlobalSegMap components
67                                ! GlobalSegMap%start, GlobalSegMap%length,
68                                ! and GlobalSegMap%pe_loc
69      public :: increasing      ! Are the indices for each pe strictly
70                                ! increasing?
71      public :: copy            ! Copy the gsmap
72
73! !PUBLIC TYPES:
74
75    type GlobalSegMap
76#ifdef SEQUENCE
77      sequence
78#endif
79      integer :: comp_id                        ! Component ID number
80      integer :: ngseg                          ! No. of Global segments
81      integer :: gsize                          ! No. of Global elements
82      integer,dimension(:),pointer :: start     ! global seg. start index
83      integer,dimension(:),pointer :: length    ! segment lengths
84      integer,dimension(:),pointer :: pe_loc    ! PE locations
85    end type GlobalSegMap
86
87    interface init ; module procedure   &
88        initd_, &       ! initialize from all PEs
89        initr_, &       ! initialize from the root
90        initp_, &       ! initialize in parallel from replicated arrays
91        initp1_, &      ! initialize in parallel from 1 replicated array
92        initp0_, &      ! null constructor using replicated data
93        init_index_     ! initialize from local index arrays
94    end interface
95
96    interface clean   ; module procedure clean_   ; end interface
97    interface comp_id ; module procedure comp_id_ ; end interface
98    interface gsize   ; module procedure gsize_   ; end interface
99    interface GlobalStorage ; module procedure &
100       GlobalStorage_
101    end interface
102    interface ProcessStorage ; module procedure &
103       ProcessStorage_
104    end interface
105    interface OrderedPoints ; module procedure &
106       OrderedPoints_
107    end interface
108    interface lsize ; module procedure lsize_ ; end interface
109    interface ngseg ; module procedure ngseg_ ; end interface
110    interface nlseg ; module procedure nlseg_ ; end interface
111    interface max_nlseg ; module procedure max_nlseg_ ; end interface
112    interface active_pes ; module procedure active_pes_ ; end interface
113    interface peLocs ; module procedure peLocs_ ; end interface
114    interface haloed ; module procedure haloed_ ; end interface
115    interface rank  ; module procedure &
116        rank1_ , &      ! single rank case
117        rankm_          ! degenerate (multiple) ranks for halo case
118    end interface
119    interface Sort ; module procedure Sort_ ; end interface
120    interface Permute ; module procedure &
121        PermuteInPlace_ 
122    end interface
123    interface SortPermute ; module procedure &
124        SortPermuteInPlace_ 
125    end interface
126    interface increasing ; module procedure increasing_ ; end interface
127    interface copy ; module procedure copy_ ; end interface
128
129
130! !REVISION HISTORY:
131!       28Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
132!       26Jan01 - J.W. Larson <larson@mcs.anl.gov> - replaced the component
133!                 GlobalSegMap%comm with GlobalSegMap%comp_id.
134!       06Feb01 - J.W. Larson <larson@mcs.anl.gov> - removed the
135!                 GlobalSegMap%lsize component.  Also, added the
136!                 GlobalStorage query function.
137!       24Feb01 - J.W. Larson <larson@mcs.anl.gov> - Added the replicated
138!                 initialization routines initp_() and initp1().
139!       25Feb01 - J.W. Larson <larson@mcs.anl.gov> - Added the routine
140!                 ProcessStorage_().
141!       18Apr01 - J.W. Larson <larson@mcs.anl.gov> - Added the routine
142!                 peLocs().
143!       26Apr01 - R. Jacob <jacob@mcs.anl.gov> - Added the routine
144!                 OrderedPoints_().
145!       03Aug01 - E. Ong <eong@mcs.anl.gov> - In initd_, call initr_
146!                 with actual shaped arguments on non-root processes to satisfy
147!                 F90 standard. See comments in initd.         
148!       18Oct01 - J.W. Larson <larson@mcs.anl.gov> - Added the routine
149!                 bcast(), and also cleaned up prologues.
150!EOP ___________________________________________________________________
151
152  character(len=*),parameter :: myname='m_GlobalSegMap'
153
154 contains
155
156!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
157!    Math and Computer Science Division, Argonne National Laboratory   !
158!BOP -------------------------------------------------------------------
159!
160! !IROUTINE: initd_ - define the map from distributed data
161!
162! !DESCRIPTION:
163! This routine takes the {\em scattered} input {\tt INTEGER} arrays
164! {\tt start}, {\tt length}, and {\tt pe\_loc}, gathers these data to
165! the {\tt root} process, and from them creates a {\em global} set of
166! segment information for the output {\tt GlobalSegMap} argument
167! {\tt GSMap}.  The input {\tt INTEGER} arguments {\tt comp\_id},
168! {\tt gsize} provide the {\tt GlobalSegMap} component ID number and
169! global grid size, respectively.  The input argument {\tt my\_comm} is
170! the F90 {\tt INTEGER} handle for the MPI communicator.  If the input
171! arrays are overdimensioned, optional argument {\em numel} can be
172! used to specify how many elements should be used.
173!
174!
175! !INTERFACE:
176
177 subroutine initd_(GSMap, start, length, root, my_comm, &
178                   comp_id, pe_loc, gsize, numel)
179
180!
181! !USES:
182!
183      use m_mpif90
184      use m_die
185      use m_stdio
186      use m_FcComms, only : fc_gather_int, fc_gatherv_int
187
188      implicit none
189
190! !INPUT PARAMETERS:
191
192      integer,dimension(:),intent(in) :: start          ! segment local start
193                                                        ! indices
194      integer,dimension(:),intent(in) :: length         ! segment local lengths
195      integer,intent(in)              :: root           ! root on my_com
196      integer,intent(in)              :: my_comm        ! local communicatior
197      integer,intent(in)              :: comp_id        ! component model ID
198      integer,dimension(:), pointer, optional :: pe_loc ! process location
199      integer,intent(in), optional    :: gsize          ! global vector size
200                                                        ! (optional).  It can
201                                                        ! be computed by this
202                                                        ! routine if no haloing
203                                                        ! is assumed.
204      integer,intent(in), optional    :: numel          ! specify number of elements
205                                                        ! to use in start, length
206
207! !OUTPUT PARAMETERS:
208
209      type(GlobalSegMap),intent(out)  :: GSMap   ! Output GlobalSegMap
210
211! !REVISION HISTORY:
212!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
213!       14Nov00 - J.W. Larson <larson@mcs.anl.gov> - final working version
214!       09Jan01 - J.W. Larson <larson@mcs.anl.gov> - repaired:  a subtle
215!                 bug concerning the usage of the argument pe_loc (result
216!                 was the new pointer variable my_pe_loc); a mistake in
217!                 the tag arguments to MPI_IRECV; a bug in the declaration
218!                 of the array status used by MPI_WAITALL.
219!       26Jan01 - J.W. Larson <larson@mcs.anl.gov> - replaced optional
220!                 argument gsm_comm with required argument comp_id.
221!       23Sep02 - Add optional argument numel to allow start, length
222!                 arrays to be overdimensioned.
223!       31Jan09 - P.H. Worley <worleyph@ornl.gov> - replaced irecv/send/waitall
224!                 logic with calls to flow controlled gather routines
225!EOP ___________________________________________________________________
226
227  character(len=*),parameter :: myname_=myname//'::initd_'
228  integer :: nPEs, myID, ier, l, i
229  integer :: ngseg  ! number of global segments
230  integer :: nlseg  ! number of local segments
231  integer :: nlseg_tmp(1) ! workaround for explicit interface expecting an array
232
233        ! arrays allocated on the root to which data are gathered
234  integer, dimension(:), allocatable :: root_start, root_length, root_pe_loc
235        ! arrays allocated on the root to coordinate gathering of
236        ! data and non-blocking receives by the root
237  integer, dimension(:), allocatable :: counts, displs
238        ! data and non-blocking receives by the root
239  integer, dimension(:), pointer :: my_pe_loc
240
241        ! Determine local process ID:
242
243  call MP_COMM_RANK(my_comm, myID, ier)
244
245  if(ier /= 0) call MP_perr_die(myname_,'MP_comm_rank()',ier)
246
247
248        ! Check consistency of sizes of input arrays:
249
250  if(size(length) /= size(start)) then
251     ier = -1
252     call die(myname_,'length/start array size mismatch',ier)
253  endif
254
255  if(present(pe_loc)) then
256     if(size(pe_loc) /= size(start)) then
257        ier = -1
258        call die(myname_,'pe_loc/start array size mismatch',ier)
259     endif
260  endif
261
262        ! Store in the variable nlseg the local size
263        ! array start(:)
264
265  if(present(numel)) then
266    nlseg=numel
267  else
268    nlseg = size(start)
269  endif
270
271        ! If the argument pe_loc is not present, then we are
272        ! initializing the GlobalSegMap on the communicator
273        ! my_comm.  We will need pe_loc to be allocated and
274        ! with local size given by the input value of nlseg,
275        ! and then initialize it with the local process id myID.
276
277  if(present(pe_loc)) then
278     my_pe_loc => pe_loc
279  else
280     allocate(my_pe_loc(nlseg), stat=ier)
281     if(ier /= 0) call die(myname_,'allocate(my_pe_loc)',ier)
282     my_pe_loc = myID
283  endif
284
285  call MPI_COMM_SIZE(my_comm, npes, ier)
286  if(ier /= 0) call MP_perr_die(myname_,'MPI_COMM_SIZE()',ier)
287
288        ! Allocate an array of displacements (displs) and counts
289        ! to hold the local values of nlseg on the root
290
291  if(myID == root) then
292     allocate(counts(0:npes-1), displs(0:npes-1), stat=ier)
293     if (ier /= 0) then 
294        call die(myname_, 'allocate(counts,...',ier)
295     endif
296  else
297     allocate(counts(1), displs(1), stat=ier)
298     if (ier /= 0) then 
299        call die(myname_, 'allocate(counts,...',ier)
300     endif
301  endif
302
303        ! Send local number of segments to the root.
304
305  nlseg_tmp(1) = nlseg
306  call fc_gather_int(nlseg_tmp, 1, MP_INTEGER, counts, 1, MP_INTEGER, &
307                     root, my_comm)
308
309        ! On the root compute the value of ngseg, along with
310        ! the entries of counts and displs.
311
312  if(myID == root) then
313     ngseg = 0
314     do i=0,npes-1
315        ngseg = ngseg + counts(i)
316        if(i == 0) then
317           displs(i) = 0
318        else
319           displs(i) = displs(i-1) + counts(i-1)
320        endif
321     end do
322  endif
323
324        ! Now only the root has the correct value of ngseg.
325
326        ! On the root, allocate memory for the arrays root_start,
327        ! and root_length.  If the argument pe_loc is present,
328        ! allocate root_pe_loc, too.
329
330        ! Non-root processes call initr_ with root_start, root_length,
331        ! and root_pe_loc, although these arguments are not used in the
332        ! subroutine. Since these correspond to dummy shaped array arguments
333        ! in initr_, the Fortran 90 standard dictates that the actual
334        ! arguments must contain complete shape information. Therefore,
335        ! these array arguments must be allocated on all processes.
336 
337  if(myID == root) then
338
339     allocate(root_start(ngseg), root_length(ngseg), &
340              root_pe_loc(ngseg), stat=ier)
341     if (ier /= 0) then
342        call die(myname_, 'allocate(root_start...',ier)
343     endif
344
345  else
346
347     allocate(root_start(1), root_length(1), &
348              root_pe_loc(1), stat=ier)
349     if (ier /= 0) then
350        call die(myname_, 'allocate((non)root_start...',ier)
351     endif
352
353  endif
354
355        ! Now, each process sends its values of start(:) to fill in
356        ! the appropriate portion of root_start(:y) on the root.
357
358  call fc_gatherv_int(start, nlseg, MP_INTEGER, &
359                      root_start, counts, displs, MP_INTEGER, &
360                      root, my_comm)
361
362        ! Next, each process sends its values of length(:) to fill in
363        ! the appropriate portion of root_length(:) on the root.
364
365  call fc_gatherv_int(length, nlseg, MP_INTEGER, &
366                      root_length, counts, displs, MP_INTEGER, &
367                      root, my_comm)
368
369        ! Finally, if the argument pe_loc is present, each process sends
370        ! its values of pe_loc(:) to fill in the appropriate portion of
371        ! root_pe_loc(:) on the root. 
372   
373  call fc_gatherv_int(my_pe_loc, nlseg, MP_INTEGER, &
374                      root_pe_loc, counts, displs, MP_INTEGER, &
375                      root, my_comm)
376
377  call MPI_BARRIER(my_comm, ier)
378  if(ier /= 0) call MP_perr_die(myname_,'MPI_BARRIER my_pe_loc',ier)
379
380        ! Now, we have everything on the root needed to call initr_().
381
382  if(present(gsize)) then
383     call initr_(GSMap, ngseg, root_start, root_length, &
384                 root_pe_loc, root, my_comm, comp_id, gsize)
385  else
386     call initr_(GSMap, ngseg, root_start, root_length, &
387                 root_pe_loc, root, my_comm, comp_id)
388  endif
389
390
391        ! Clean up the array pe_loc(:) if it was allocated
392
393  if(present(pe_loc)) then
394     nullify(my_pe_loc) 
395  else
396     deallocate(my_pe_loc, stat=ier)
397     if(ier /= 0) call die(myname_, 'deallocate(my_pe_loc)', ier)
398  endif
399
400        ! Clean up the arrays root_start(:), et cetera...
401
402  deallocate(root_start, root_length, root_pe_loc, stat=ier)
403  if(ier /= 0) then
404     call die(myname_, 'deallocate(root_start,...)', ier)
405  endif
406
407        ! Clean up the arrays counts(:) and displs(:)
408
409   deallocate(counts, displs, stat=ier)
410   if(ier /= 0) then
411     call die(myname_, 'deallocate(counts,...)', ier)
412   endif
413
414 end subroutine initd_
415
416!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
417!    Math and Computer Science Division, Argonne National Laboratory   !
418!BOP -------------------------------------------------------------------
419!
420! !IROUTINE: initr_ initialize the map from the root
421!
422! !DESCRIPTION:
423! This routine takes the input {\tt INTEGER} arrays {\tt start},
424! {\tt length}, and {\tt pe\_loc} (all valid only on the {\tt root}
425! process), and from them creates a {\em global} set of segment
426! information for the output {\tt GlobalSegMap} argument
427! {\tt GSMap}.  The input {\tt INTEGER} arguments {\tt ngseg},
428! {\tt comp\_id}, {\tt gsize} (again, valid only on the {\tt root}
429! process) provide the {\tt GlobalSegMap} global segment count, component
430! ID number, and global grid size, respectively.  The input argument
431! {\tt my\_comm} is the F90 {\tt INTEGER} handle for the MPI communicator.
432!
433! !INTERFACE:
434
435 subroutine initr_(GSMap, ngseg, start, length, pe_loc, root,  &
436                   my_comm, comp_id, gsize)
437!
438! !USES:
439!
440      use m_mpif90
441      use m_die
442      use m_stdio
443 
444      implicit none
445
446! !INPUT PARAMETERS:
447
448      integer, intent(in)             :: ngseg   ! no. of global segments
449      integer,dimension(:),intent(in) :: start   ! segment local start index
450      integer,dimension(:),intent(in) :: length  ! the distributed sizes
451      integer,dimension(:),intent(in) :: pe_loc  ! process location
452      integer,intent(in)              :: root    ! root on my_com
453      integer,intent(in)              :: my_comm ! local communicatior
454      integer,intent(in)              :: comp_id ! component id number
455      integer,intent(in), optional    :: gsize   ! global vector size
456                                                 ! (optional).  It can
457                                                 ! be computed by this
458                                                 ! routine if no haloing
459                                                 ! is assumed.
460
461! !OUTPUT PARAMETERS:
462
463      type(GlobalSegMap),intent(out)  :: GSMap   ! Output GlobalSegMap
464
465! !REVISION HISTORY:
466!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
467!       09Nov00 - J.W. Larson <larson@mcs.anl.gov> - final working version
468!       10Jan01 - J.W. Larson <larson@mcs.anl.gov> - minor bug fix
469!       12Jan01 - J.W. Larson <larson@mcs.anl.gov> - minor bug fix regarding
470!                                                    disparities in ngseg on
471!                                                    the root and other
472!                                                    processes
473!       26Jan01 - J.W. Larson <larson@mcs.anl.gov> - replaced optional
474!                 argument gsm_comm with required argument comp_id.
475!EOP ___________________________________________________________________
476
477  character(len=*),parameter :: myname_=myname//'::initr_'
478  integer :: myID,ier,l,i
479
480        ! Determine the local process ID myID:
481
482  call MPI_COMM_RANK(my_comm, myID, ier)
483  if(ier/=0) call MP_perr_die(myname_,'MPI_COMM_RANK()',ier)
484
485        ! Argument checking:  check to make sure the arrays
486        ! start, length, and pe_loc each have ngseg elements.
487        ! If not, stop with an error.  This is done on the
488        ! root process since it owns the initialization data.
489
490  if(myID == root) then
491    if( size(start(:)) /= ngseg ) then
492      write(stderr,'(2a,2(a,i4))') myname_,     &
493        ': _root_ argument error',              &
494        ', size(start) =',size(start),  &
495        ', ngseg =',ngseg
496      call die(myname_)
497    endif
498    if( size(length(:)) /= ngseg ) then
499      write(stderr,'(2a,2(a,i4))') myname_,     &
500        ': _root_ argument error',              &
501        ', size(length) =',size(length),        &
502        ', ngseg =',ngseg
503      call die(myname_)
504    endif
505    if( size(pe_loc(:)) /= ngseg ) then
506      write(stderr,'(2a,2(a,i4))') myname_,     &
507        ': _root_ argument error',              &
508        ', size(pe_loc) =',size(pe_loc),        &
509        ', ngseg =',ngseg
510      call die(myname_)
511    endif
512  endif
513
514        ! Initialize GSMap%ngseg and GSMap%comp_id on the root:
515
516  if(myID == root) then
517     GSMap%ngseg = ngseg
518     GSMap%comp_id = comp_id
519  endif
520
521        ! Broadcast the value of GSMap%ngseg
522
523  call MPI_BCAST(GSMap%ngseg, 1, MP_INTEGER, root, my_comm, ier)
524  if(ier/=0) call MP_perr_die(myname_,'MPI_BCAST(GSmap%ngseg)',ier)
525
526        ! Broadcast the value of GSMap%comp_id
527
528  call MPI_BCAST(GSMap%comp_id, 1, MP_INTEGER, root, my_comm, ier)
529  if(ier/=0) call MP_perr_die(myname_,'MPI_BCAST(GSmap%comp_id)',ier)
530
531        ! Allocate the components GSMap%start(:), GSMap%length(:),
532        ! and GSMap%pe_loc(:)
533
534  allocate(GSMap%start(GSMap%ngseg), GSMap%length(GSMap%ngseg), &
535           GSMap%pe_loc(GSMap%ngseg), stat = ier)
536  if(ier/=0) call die(myname_,'allocate(GSmap%start(:),...',ier)
537
538#ifdef MALL_ON
539        call mall_ci(size(transfer(GSMap%start,(/1/))),myname_)
540        call mall_ci(size(transfer(GSMap%length,(/1/))),myname_)
541        call mall_ci(size(transfer(GSMap%pe_loc,(/1/))),myname_)
542#endif
543
544        ! On the root process, initialize GSMap%start(:), GSMap%length(:),
545        ! and GSMap%pe_loc(:) with the data contained in start(:),
546        ! length(:) and pe_loc(:), respectively
547
548  if(myID == root) then
549     GSMap%start(1:GSMap%ngseg) = start(1:GSMap%ngseg)
550     GSMap%length(1:GSMap%ngseg) = length(1:GSMap%ngseg)
551     GSMap%pe_loc(1:GSMap%ngseg) = pe_loc(1:GSMap%ngseg)
552  endif
553
554        ! Broadcast the root values of GSMap%start(:), GSMap%length(:),
555        ! and GSMap%pe_loc(:)
556
557  call MPI_BCAST(GSMap%start, GSMap%ngseg, MP_INTEGER, root, my_comm, ier)
558  if(ier/=0) call MP_perr_die(myname_,'MPI_BCAST(GSMap%start)',ier)
559
560  call MPI_BCAST(GSMap%length, GSMap%ngseg, MP_INTEGER, root, my_comm, ier)
561  if(ier/=0) call MP_perr_die(myname_,'MPI_BCAST(GSMap%length)',ier)
562
563  call MPI_BCAST(GSMap%pe_loc, GSMap%ngseg, MP_INTEGER, root, my_comm, ier)
564  if(ier/=0) call MP_perr_die(myname_,'MPI_BCAST(GSMap%pe_loc)',ier)
565
566        ! If the argument gsize is present, use the root value to
567        ! set GSMap%gsize and broadcast it.  If it is not present,
568        ! this will be computed by summing the entries of GSM%length(:).
569        ! Again, note that if one is storing halo points, the sum will
570        ! produce a result larger than the actual global vector.  If
571        ! halo points are to be used in the mapping we advise strongly
572        ! that the user specify the value gsize as an argument.
573
574  if(present(gsize)) then
575     if(myID == root) then
576        GSMap%gsize = gsize
577     endif
578     call MPI_BCAST(GSMap%gsize, 1, MP_INTEGER, root, my_comm, ier)
579     if(ier/=0) call MP_perr_die(myname_, 'MPI_BCAST(GSMap%gsize)', ier)
580  else
581     GSMap%gsize = 0
582     do i=1,GSMap%ngseg
583        GSMap%gsize = GSMap%gsize + GSMap%length(i)
584     end do
585  endif
586
587 end subroutine initr_
588
589!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
590!    Math and Computer Science Division, Argonne National Laboratory   !
591!BOP -------------------------------------------------------------------
592!
593! !IROUTINE: initp_ - define the map from replicated data.
594!
595! !DESCRIPTION:
596!
597! The routine {\tt initp\_()} takes the input {\em replicated} arguments
598! {\tt comp\_id}, {\tt ngseg}, {\tt gsize}, {\tt start(:)},
599! {\tt length(:)}, and {\tt pe\_loc(:)}, and uses them to initialize an
600! output {\tt GlobalSegMap} {\tt GSMap}.  This routine operates on the
601! assumption that these data are replicated across the communicator on
602! which the {\tt GlobalSegMap} is being created.
603!
604! !INTERFACE:
605
606 subroutine initp_(GSMap, comp_id, ngseg, gsize, start, length, pe_loc)
607
608!
609! !USES:
610!
611      use m_mpif90
612      use m_die, only : die
613      use m_stdio
614
615      implicit none
616
617! !INPUT PARAMETERS:
618
619      integer,intent(in)              :: comp_id ! component model ID
620      integer,intent(in)              :: ngseg   ! global number of segments
621      integer,intent(in)              :: gsize   ! global vector size
622      integer,dimension(:),intent(in) :: start   ! segment local start index
623      integer,dimension(:),intent(in) :: length  ! the distributed sizes
624      integer,dimension(:),intent(in) :: pe_loc  ! process location
625
626! !OUTPUT PARAMETERS:
627
628      type(GlobalSegMap),intent(out)  :: GSMap   ! Output GlobalSegMap
629
630! !REVISION HISTORY:
631!       24Feb01 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
632!EOP ___________________________________________________________________
633
634  character(len=*),parameter :: myname_=myname//'::initp_'
635  integer :: ierr, n
636
637       ! Argument Checks -- Is comp_id positive?
638
639  if(comp_id <= 0) then
640     call die(myname_,'non-positive value of comp_id',comp_id)
641  endif
642
643       ! Is gsize positive?
644
645  if(gsize <= 0) then
646     call die(myname_,'non-positive value of gsize',gsize)
647  endif
648
649
650       ! Is ngseg positive?
651
652  if(ngseg <= 0) then
653     call die(myname_,'non-positive value of ngseg',ngseg)
654  endif
655
656       ! Are the arrays start(:), length(:), and pe_loc(:) the
657       !correct size?
658
659  if(size(start) /= ngseg) then
660     call die(myname_,'start(:)/ngseg size mismatch',ngseg)
661  endif
662  if (size(length) /= ngseg) then
663     call die(myname_,'length(:)/ngseg size mismatch',ngseg)
664  endif
665  if (size(pe_loc) /= ngseg) then
666     call die(myname_,'pe_loc(:)/ngseg size mismatch',ngseg)
667  endif
668
669       ! Allocate index and location arrays for GSMap:
670
671  allocate(GSMap%start(ngseg), GSMap%length(ngseg), GSMap%pe_loc(ngseg), &
672           stat = ierr)
673  if (ierr /= 0) then
674     call die(myname_,'allocate(GSMap%start...',ngseg)
675  endif
676       
677       ! Assign the components of GSMap:
678
679  GSMap%comp_id = comp_id
680  GSMap%ngseg = ngseg
681  GSMap%gsize = gsize
682
683  do n=1,ngseg
684     GSMap%start(n) = start(n)
685     GSMap%length(n) = length(n)
686     GSMap%pe_loc(n) = pe_loc(n)
687  end do
688
689  end subroutine initp_     
690
691!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
692!    Math and Computer Science Division, Argonne National Laboratory   !
693!BOP -------------------------------------------------------------------
694!
695! !IROUTINE: initp1_ - define the map from replicated data using 1 array.
696!
697! !DESCRIPTION:
698!
699! The routine {\tt initp1\_()} takes the input {\em replicated} arguments
700! {\tt comp\_id}, {\tt ngseg}, {\tt gsize}, and {\tt all\_arrays(:)},
701! and uses them to initialize an output {\tt GlobalSegMap} {\tt GSMap}. 
702! This routine operates on the assumption that these data are replicated
703! across the communicator on which the {\tt GlobalSegMap} is being created.
704! The input array {\tt all\_arrays(:)} should be of length {\tt 2 * ngseg},
705! and is packed so that
706! $$ {\tt all\_arrays(1:ngseg)} = {\tt GSMap\%start(1:ngseg)} $$
707! $$ {\tt all\_arrays(ngseg+1:2*ngseg)} = {\tt GSMap\%length(1:ngseg)} $$
708! $$ {\tt all\_arrays(2*ngseg+1:3*ngseg)} = {\tt GSMap\%pe\_loc(1:ngseg)} .$$
709!
710! !INTERFACE:
711
712 subroutine initp1_(GSMap, comp_id, ngseg, gsize, all_arrays)
713
714!
715! !USES:
716!
717      use m_mpif90
718      use m_die, only : die
719      use m_stdio
720
721      implicit none
722
723! !INPUT PARAMETERS:
724
725      integer,intent(in)              :: comp_id    ! component model ID
726      integer,intent(in)              :: ngseg      ! global no. of segments
727      integer,intent(in)              :: gsize      ! global vector size
728      integer,dimension(:),intent(in) :: all_arrays ! packed array of length
729                                                    ! 3*ngseg containing (in
730                                                    ! this order):  start(:),
731                                                    ! length(:), and pe_loc(:)
732
733! !OUTPUT PARAMETERS:
734
735      type(GlobalSegMap),intent(out)  :: GSMap      ! Output GlobalSegMap
736
737! !REVISION HISTORY:
738!       24Feb01 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
739!EOP ___________________________________________________________________
740
741  character(len=*),parameter :: myname_=myname//'::initp1_'
742  integer :: ierr, n
743
744       ! Argument Checks -- Is comp_id positive?
745
746  if(comp_id <= 0) then
747     call die(myname_,'non-positive value of comp_id',comp_id)
748  endif
749
750       ! Is gsize positive?
751
752  if(gsize <= 0) then
753     call die(myname_,'non-positive value of gsize',gsize)
754  endif
755
756
757       ! Is ngseg positive?
758
759  if(ngseg <= 0) then
760     call die(myname_,'non-positive value of ngseg',ngseg)
761  endif
762
763       ! Is the array all_arrays(:) the right length?
764
765  if(size(all_arrays) /= 3*ngseg) then
766     call die(myname_,'all_arrays(:)/3*ngseg size mismatch',ngseg)
767  endif
768
769       ! Allocate index and location arrays for GSMap:
770
771  allocate(GSMap%start(ngseg), GSMap%length(ngseg), GSMap%pe_loc(ngseg), &
772           stat = ierr)
773  if (ierr /= 0) then
774     call die(myname_,'allocate(GSMap%start...',ngseg)
775  endif
776       
777       ! Assign the components of GSMap:
778
779  GSMap%comp_id = comp_id
780  GSMap%ngseg = ngseg
781  GSMap%gsize = gsize
782
783  do n=1,ngseg
784     GSMap%start(n) = all_arrays(n)
785     GSMap%length(n) = all_arrays(ngseg + n)
786     GSMap%pe_loc(n) = all_arrays(2*ngseg + n)
787  end do
788
789  end subroutine initp1_     
790
791!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
792!    Math and Computer Science Division, Argonne National Laboratory   !
793!BOP -------------------------------------------------------------------
794!
795! !IROUTINE: initp0_ - Null Constructor Using Replicated Data
796!
797! !DESCRIPTION:
798!
799! The routine {\tt initp0\_()} takes the input {\em replicated} arguments
800! {\tt comp\_id}, {\tt ngseg}, {\tt gsize}, and uses them perform null
801! construction of the output {\tt GlobalSegMap} {\tt GSMap}.  This is a
802! null constructor in the sense that we are not filling in the segment
803! information arrays.  This routine operates on the assumption that these
804! data are replicated across the communicator on which the
805! {\tt GlobalSegMap} is being created.
806!
807! !INTERFACE:
808
809 subroutine initp0_(GSMap, comp_id, ngseg, gsize)
810
811!
812! !USES:
813!
814      use m_die, only : die
815      use m_stdio
816
817      implicit none
818
819! !INPUT PARAMETERS:
820
821      integer,intent(in)              :: comp_id ! component model ID
822      integer,intent(in)              :: ngseg   ! global number of segments
823      integer,intent(in)              :: gsize   ! global vector size
824
825! !OUTPUT PARAMETERS:
826
827      type(GlobalSegMap),intent(out)  :: GSMap   ! Output GlobalSegMap
828
829! !REVISION HISTORY:
830! 13Aug03 - J.W. Larson <larson@mcs.anl.gov> - Initial version.
831!EOP ___________________________________________________________________
832
833  character(len=*),parameter :: myname_=myname//'::initp0_'
834
835  integer :: ierr
836
837  nullify(GSMap%start)
838  nullify(GSMap%length)
839  nullify(GSMap%pe_loc)
840
841  GSMap%comp_id = comp_id
842  GSMap%ngseg = ngseg
843  GSMap%gsize = gsize
844
845  allocate(GSMap%start(ngseg), GSMap%length(ngseg), GSMap%pe_loc(ngseg), &
846           stat=ierr)
847  if(ierr /= 0) then
848     write(stderr,'(3a,i8)') myname_, &
849          ':: FATAL--allocate of segment information storage space failed.', &
850          '  ierr = ',ierr
851     call die(myname_)
852  endif
853
854 end subroutine initp0_
855
856
857
858!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
859!    Math and Computer Science Division, Argonne National Laboratory   !
860!BOP -------------------------------------------------------------------
861!
862! !IROUTINE: init_index_ - initialize GSM from local index arrays
863!
864! !DESCRIPTION:
865!
866! The routine {\tt init\_index\_()} takes a local array of indices
867! {\tt lindx} and uses them to create a {\tt GlobalSegMap}. 
868! {\tt lindx} is parsed to determine the lengths of the runs, and
869! then a call is made to {\tt initd\_}.  The optional argument
870! {\tt lsize} can be used if only the first {\tt lsize} number
871! of elements of {\tt lindx} are valid.  The optional argument
872! {\tt gsize} is used to specify the global number of unique points
873! if this can not be determined from the collective {\tt lindx}.
874!
875!
876! !INTERFACE:
877
878    subroutine init_index_(GSMap, lindx, my_comm, comp_id, lsize, gsize)
879
880!
881! !USES:
882!
883
884!  use m_GlobalSegMap,only: GlobalSegMap
885!  use m_GlobalSegMap,only: MCT_GSMap_init => init
886
887!  use shr_sys_mod
888
889  use m_die
890  implicit none
891
892! !INPUT PARAMETERS:
893
894     integer , dimension(:),intent(in) :: lindx   ! index buffer
895     integer , intent(in) :: my_comm         ! mpi communicator group (mine)
896     integer , intent(in) :: comp_id         ! component id (mine)
897
898     integer , intent(in),optional :: lsize  ! size of index buffer
899     integer , intent(in),optional :: gsize  ! global vector size
900
901! !OUTPUT PARAMETERS:
902
903     type(GlobalSegMap),intent(out) :: GSMap ! Output GlobalSegMap
904     
905
906! !REVISION HISTORY:
907!       30Jul02 - T. Craig - initial version in cpl6.
908!       17Nov05 - R. Loy <rloy@mcs.anl.gov> - install into MCT
909!       18Nov05 - R. Loy <rloy@mcs.anl.gov> - make lsize optional
910!       25Jul06 - R. Loy <rloy@mcs.anl.gov> - error check on lindex/alloc/dealloc
911!EOP ___________________________________________________________________
912
913
914     !--- local ---
915
916     character(len=*),parameter :: myname_=myname//'::init_index_'
917
918     integer             :: i,j,k,n      ! generic indicies
919     integer             :: nseg         ! counts number of segments for GSMap
920     integer,allocatable :: start(:)     ! used to init GSMap
921     integer,allocatable :: count(:)     ! used to init GSMap
922     integer,parameter   :: pid0=0       ! mpi process id for root pe
923     integer,parameter   :: debug=0      !
924
925     integer rank,ierr
926     integer mysize
927
928
929     if (present(lsize)) then
930       mysize=lsize
931     else
932       mysize=size(lindx)
933     endif
934
935     if (mysize<0) call die(myname_, &
936        'lindx size is negative (you may have run out of points)')
937
938!!
939!! Special case if this processor doesn't have any data indices
940!!
941   if (mysize==0) then
942     allocate(start(0),count(0),stat=ierr)
943     if(ierr/=0) call die(myname_,'allocate(start,count)',ierr)
944   
945     nseg=0
946   else
947
948     call MPI_COMM_RANK(my_comm,rank, ierr)
949
950     ! compute segment's start indicies and length counts
951
952     ! first pass - count how many runs of consecutive numbers
953
954     nseg=1
955     do n = 2,mysize
956        i = lindx(n-1)
957        j = lindx(n)
958        if ( j-i /= 1) nseg=nseg+1
959     end do
960
961     allocate(start(nseg),count(nseg),stat=ierr)
962     if(ierr/=0) call die(myname_,'allocate(start,count)',ierr)
963
964     ! second pass - determine how long each run is
965
966     nseg = 1
967     start(nseg) = lindx(1)
968     count(nseg) = 1
969     do n = 2,mysize
970        i = lindx(n-1)
971        j = lindx(n)
972        if ( j-i /= 1) then
973            nseg = nseg+1
974            start(nseg) = lindx(n)
975            count(nseg) = 1
976         else
977            count(nseg) = count(nseg)+1
978         end if
979     end do
980
981   endif ! if mysize==0
982
983
984     if (debug.ne.0) then
985        write(6,*) rank,'init_index: SIZE ',nseg
986
987        do n=1,nseg
988          write(6,*) rank,'init_index: START,COUNT ',start(n),count(n)
989        end do
990     endif
991
992
993      if (present(gsize)) then
994        call initd_( GSMap, start, count, pid0, my_comm,  &
995                     comp_id, gsize=gsize)
996      else
997        call initd_( GSMap, start, count, pid0, my_comm,  &
998                     comp_id)
999      endif
1000
1001
1002      deallocate(start, count, stat=ierr)
1003      if(ierr/=0) call warn(myname_,'deallocate(start,count)',ierr)
1004     
1005
1006   end subroutine init_index_
1007
1008
1009
1010!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1011!    Math and Computer Science Division, Argonne National Laboratory   !
1012!BOP -------------------------------------------------------------------
1013!
1014! !IROUTINE: clean_ - clean the map
1015!
1016! !DESCRIPTION:
1017! This routine deallocates the array components of the {\tt GlobalSegMap}
1018! argument {\tt GSMap}: {\tt GSMap\%start}, {\tt GSMap\%length}, and
1019! {\tt GSMap\%pe\_loc}.  It also zeroes out the values of the integer
1020! components {\tt GSMap\%ngseg}, {\tt GSMap\%comp\_id}, and
1021! {\tt GSMap\%gsize}.
1022!
1023! !INTERFACE:
1024
1025    subroutine clean_(GSMap,stat)
1026!
1027! !USES:
1028!
1029      use m_die
1030
1031      implicit none
1032 
1033! !INPUT/OUTPUT PARAMETERS:
1034
1035      type(GlobalSegMap), intent(inout) :: GSMap
1036      integer, optional,  intent(out)   :: stat
1037
1038! !REVISION HISTORY:
1039!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
1040!       01Mar02 - E.T. Ong <eong@mcs.anl.gov> - added stat argument.
1041!                 Removed dies to prevent crashing.
1042!EOP ___________________________________________________________________
1043
1044  character(len=*),parameter :: myname_=myname//'::clean_'
1045  integer :: ier
1046
1047#ifdef MALL_ON
1048 
1049  if( (associated(GSMap%start) .and. associated(GSMap%length)) &
1050       .and. associated(GSMap%pe_loc) )
1051     call mall_co(size(transfer(GSMap%start,(/1/))),myname_)
1052     call mall_co(size(transfer(GSMap%length,(/1/))),myname_)
1053     call mall_co(size(transfer(GSMap%pe_loc,(/1/))),myname_)
1054  endif
1055
1056#endif
1057
1058  deallocate(GSMap%start, GSMap%length, GSMap%pe_loc, stat=ier)
1059
1060  if(present(stat)) then
1061     stat=ier
1062  else
1063     if(ier /= 0) call warn(myname_,'deallocate(GSMap%start,...)',ier)
1064  endif
1065
1066  GSMap%ngseg = 0
1067  GSMap%comp_id  = 0
1068  GSMap%gsize = 0
1069
1070 end subroutine clean_
1071
1072!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1073!    Math and Computer Science Division, Argonne National Laboratory   !
1074!BOP -------------------------------------------------------------------
1075!
1076! !IROUTINE: ngseg_ - Return the global number of segments from the map
1077!
1078! !DESCRIPTION:
1079! The function {\tt ngseg\_()} returns the global number of vector
1080! segments in the {\tt GlobalSegMap} argument {\tt GSMap}.  This is
1081! merely the value of {\tt GSMap\%ngseg}.
1082!
1083! !INTERFACE:
1084
1085 integer function ngseg_(GSMap)
1086
1087      implicit none
1088
1089! !INPUT PARAMETERS:
1090
1091      type(GlobalSegMap),intent(in) :: GSMap
1092
1093! !REVISION HISTORY:
1094!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
1095!EOP ___________________________________________________________________
1096
1097  character(len=*),parameter :: myname_=myname//'::ngseg_'
1098
1099  ngseg_=GSMap%ngseg
1100
1101 end function ngseg_
1102
1103!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1104!    Math and Computer Science Division, Argonne National Laboratory   !
1105!BOP -------------------------------------------------------------------
1106!
1107! !IROUTINE: nlseg_ - Return the local number of segments from the map
1108!
1109! !DESCRIPTION:
1110! The function {\tt nlseg\_()} returns the number of vector segments
1111! in the {\tt GlobalSegMap} argument {\tt GSMap} that reside on the
1112! process specified by the input argument {\tt pID}.  This is the
1113! number of entries {\tt GSMap\%pe\_loc} whose value equals {\tt pID}.
1114!
1115! !INTERFACE:
1116
1117 integer function nlseg_(GSMap, pID)
1118
1119      implicit none
1120
1121! !INPUT PARAMETERS:
1122
1123      type(GlobalSegMap),intent(in) :: GSMap
1124      integer,           intent(in) :: pID
1125
1126! !REVISION HISTORY:
1127!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
1128!       14Jun01 - J.W. Larson <larson@mcs.anl.gov> - Bug fix in lower
1129!                 limit of loop over elements of GSMap%pe_loc(:).  The
1130!                 original code had this lower limit set to 0, which
1131!                 was out-of-bounds (but uncaught).  The correct lower
1132!                 index is 1.  This bug was discovered by Everest Ong.
1133!EOP ___________________________________________________________________
1134
1135  character(len=*),parameter :: myname_=myname//'::nlseg_'
1136  integer :: i, nlocseg
1137
1138        ! Initialize the number of segments residing on pID, nlocseg
1139
1140  nlocseg = 0
1141
1142        ! Compute the number of segments residing on pID, nlocseg
1143
1144  do i=1,GSMap%ngseg
1145     if(GSMap%pe_loc(i) == pID) then
1146        nlocseg = nlocseg + 1
1147     endif
1148  end do
1149
1150        ! Return the total
1151
1152  nlseg_ = nlocseg
1153
1154 end function nlseg_
1155
1156
1157
1158!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1159!    Math and Computer Science Division, Argonne National Laboratory   !
1160!BOP -------------------------------------------------------------------
1161!
1162! !IROUTINE: max_nlseg_ - Return the max number of segments over all procs
1163!
1164! !DESCRIPTION:
1165! The function {\tt max\_nlseg\_()} returns the maximum number
1166! over all processors of the vector
1167! segments in the {\tt GlobalSegMap} argument {\tt gsap}
1168! E.g. max\_p(nlseg(gsmap,p)) but computed more efficiently
1169!
1170! !INTERFACE:
1171
1172        integer function max_nlseg_(gsmap)
1173
1174! !USES:
1175
1176        use m_MCTWorld,      only :ThisMCTWorld
1177        use m_mpif90
1178        use m_die
1179
1180        use m_stdio    ! rml
1181
1182        implicit none
1183
1184! !INPUT PARAMETERS:
1185
1186        type(GlobalSegMap), intent(in) :: gsmap
1187
1188
1189! !REVISION HISTORY:
1190!       17Jan07 - R. Loy <rloy@mcs.anl.gov> - initial prototype
1191!EOP ___________________________________________________________________
1192
1193
1194
1195! Local variables
1196
1197        character(len=*),parameter :: myname_=myname//'::max_local_segs'
1198
1199        integer i
1200        integer this_comp_id
1201        integer nprocs
1202
1203        integer, allocatable::  segcount(:)  ! segments on proc i
1204        integer ier
1205
1206        integer this_ngseg
1207        integer segment_pe
1208        integer max_segcount
1209
1210
1211! Start of routine
1212
1213        this_comp_id = comp_id(gsmap)
1214        nprocs=ThisMCTWorld%nprocspid(this_comp_id)
1215
1216        allocate( segcount(nprocs), stat=ier )
1217        if (ier/=0) call die(myname_,'allocate segcount')
1218
1219        segcount=0
1220
1221        this_ngseg=ngseg(gsmap)
1222
1223        do i=1,this_ngseg
1224
1225          segment_pe = gsmap%pe_loc(i) + 1     ! want value 1..nprocs
1226
1227          if (segment_pe < 1 .OR. segment_pe > nprocs) then
1228            call die(myname_,'bad segment location',segment_pe)
1229          endif
1230
1231          segcount(segment_pe) = segcount(segment_pe) + 1
1232        enddo
1233
1234        max_segcount=0
1235        do i=1,nprocs
1236          max_segcount= max( max_segcount, segcount(i) )
1237        enddo
1238
1239        deallocate(segcount, stat=ier)
1240        if (ier/=0) call die(myname_,'deallocate segcount')
1241
1242
1243        max_nlseg_=max_segcount
1244
1245        end function max_nlseg_
1246
1247
1248!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1249!    Math and Computer Science Division, Argonne National Laboratory   !
1250!BOP -------------------------------------------------------------------
1251!
1252! !IROUTINE: comp_id_ - Return the commponent ID from the GlobalSegMap.
1253!
1254! !DESCRIPTION:
1255! The function {\tt comp\_id\_()} returns component ID number stored in
1256! {\tt GSMap\%comp\_id}.
1257!
1258! !INTERFACE:
1259
1260 integer function comp_id_(GSMap)
1261
1262! !USES:
1263
1264      use m_die,only: die
1265      use m_stdio, only :stderr
1266
1267      implicit none
1268
1269! !INPUT PARAMETERS:
1270
1271      type(GlobalSegMap),intent(in) :: GSMap
1272
1273! !REVISION HISTORY:
1274!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
1275!       26Jan01 - J.W. Larson <larson@mcs.anl.gov> - renamed comp_id_
1276!                 to fit within MCT_World component ID context.
1277!       01May01 - R.L. Jacob  <jacob@mcs.anl.gov> - make sure GSMap
1278!                 is defined.
1279!EOP ___________________________________________________________________
1280
1281  character(len=*),parameter :: myname_=myname//'::comp_id_'
1282
1283  if(.not.associated(GSMap%start) ) then
1284    write(stderr,'(2a)') myname_, &
1285    ' MCTERROR:  GSMap argument not initialized...exiting'
1286    call die(myname_)
1287  endif
1288
1289  comp_id_ = GSMap%comp_id
1290
1291 end function comp_id_
1292
1293!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1294!    Math and Computer Science Division, Argonne National Laboratory   !
1295!BOP -------------------------------------------------------------------
1296!
1297! !IROUTINE: gsize_ - Return the global vector size from the GlobalSegMap.
1298!
1299! !DESCRIPTION:
1300! The function {\tt gsize\_()} takes the input {\tt GlobalSegMap}
1301! arguement {\tt GSMap} and returns the global vector length stored
1302! in {\tt GlobalSegMap\%gsize}.
1303!
1304! !INTERFACE:
1305
1306 integer function gsize_(GSMap)
1307
1308      implicit none
1309
1310! !INPUT PARAMETERS:
1311
1312      type(GlobalSegMap),intent(in) :: GSMap
1313
1314! !REVISION HISTORY:
1315!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
1316!EOP ___________________________________________________________________
1317
1318  character(len=*),parameter :: myname_=myname//'::gsize_'
1319
1320  gsize_=GSMap%gsize
1321
1322 end function gsize_
1323
1324!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1325!    Math and Computer Science Division, Argonne National Laboratory   !
1326!BOP -------------------------------------------------------------------
1327!
1328! !IROUTINE: GlobalStorage_ - Return global storage space required.
1329!
1330! !DESCRIPTION:
1331! The function {\tt GlobalStorage\_()} takes the input {\tt GlobalSegMap}
1332! arguement {\tt GSMap} and returns the global storage space required
1333! ({\em i.e.}, the vector length) to hold all the data specified by
1334! {\tt GSMap}.
1335!
1336! {\bf N.B.:  } If {\tt GSMap} contains halo or masked points, the value
1337! by {\tt GlobalStorage\_()} may differ from {\tt GSMap\%gsize}.
1338!
1339! !INTERFACE:
1340
1341 integer function GlobalStorage_(GSMap)
1342
1343      implicit none
1344
1345! !INPUT PARAMETERS:
1346
1347      type(GlobalSegMap),intent(in) :: GSMap
1348
1349! !REVISION HISTORY:
1350!       06Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version
1351!EOP ___________________________________________________________________
1352
1353  character(len=*),parameter :: myname_=myname//'::GlobalStorage_'
1354
1355  integer :: global_storage, ngseg, n
1356
1357      ! Return global number of segments:
1358
1359  ngseg = ngseg_(GSMap)
1360
1361      ! Initialize global_storage (the total number of points in the
1362      ! GlobalSegMap:
1363
1364  global_storage = 0
1365
1366      ! Add up the number of points present in the GlobalSegMap:
1367
1368  do n=1,ngseg
1369     global_storage = global_storage + GSMap%length(n)
1370  end do
1371
1372  GlobalStorage_ = global_storage
1373
1374 end function GlobalStorage_
1375
1376!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1377!    Math and Computer Science Division, Argonne National Laboratory   !
1378!BOP -------------------------------------------------------------------
1379!
1380! !IROUTINE: ProcessStorage_ - Number of points on a given process.
1381!
1382! !DESCRIPTION:
1383! The function {\tt ProcessStorage\_()} takes the input {\tt GlobalSegMap}
1384! arguement {\tt GSMap} and returns the storage space required by process
1385! {\tt PEno} ({\em i.e.}, the vector length) to hold all the data specified
1386! by {\tt GSMap}.
1387!
1388! !INTERFACE:
1389
1390 integer function ProcessStorage_(GSMap, PEno)
1391
1392      implicit none
1393
1394! !INPUT PARAMETERS:
1395
1396      type(GlobalSegMap),intent(in) :: GSMap
1397      integer,           intent(in) :: PEno
1398
1399! !REVISION HISTORY:
1400!       06Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version
1401!EOP ___________________________________________________________________
1402
1403  character(len=*),parameter :: myname_=myname//'::ProcessStorage_'
1404
1405  integer :: pe_storage, ngseg, n
1406
1407      ! Return global number of segments:
1408
1409  ngseg = ngseg_(GSMap)
1410
1411      ! Initialize pe_storage (the total number of points on process
1412      ! PEno in the GlobalSegMap):
1413
1414  pe_storage = 0
1415
1416      ! Add up the number of points on process PEno in the GlobalSegMap:
1417
1418  do n=1,ngseg
1419     if(GSMap%pe_loc(n) == PEno) then
1420        pe_storage = pe_storage + GSMap%length(n)
1421     endif
1422  end do
1423
1424  ProcessStorage_ = pe_storage
1425
1426 end function ProcessStorage_
1427
1428!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1429!    Math and Computer Science Division, Argonne National Laboratory   !
1430!BOP -------------------------------------------------------------------
1431!
1432! !IROUTINE: OrderedPoints_ - The grid points on a given process
1433!                               returned in the assumed MCT order.
1434!
1435! !DESCRIPTION:
1436! The function {\tt OrderedPoints\_()} takes the input {\tt GlobalSegMap}
1437! arguement {\tt GSMap} and returns a vector of the points owned by
1438! {\tt PEno}.  {\tt Points} is allocated here.  The calling process
1439! is responsible for deallocating the space.
1440!
1441! !INTERFACE:
1442
1443    subroutine OrderedPoints_(GSMap, PEno, Points)
1444
1445!
1446! !USES:
1447!
1448      use m_die,only: die
1449
1450      implicit none
1451
1452 ! !INPUT PARAMETERS:
1453
1454      type(GlobalSegMap), intent(in) :: GSMap   ! input GlobalSegMap
1455      integer,            intent(in) :: PEno    ! input process number
1456      integer,dimension(:),pointer   :: Points  ! the vector of points
1457
1458! !REVISION HISTORY:
1459!       25Apr01 - R. Jacob <jacob@mcs.anl.gov> - initial prototype
1460!EOP ___________________________________________________________________
1461
1462  character(len=*),parameter :: myname_=myname//'::OrderedPoints_'
1463  integer :: nlsegs,mysize,ier,i,j,k
1464  integer,dimension(:),allocatable :: mystarts,mylengths
1465
1466  nlsegs = nlseg(GSMap,PEno)
1467  mysize=ProcessStorage(GSMap,PEno)
1468
1469  allocate(mystarts(nlsegs),mylengths(nlsegs), &
1470      Points(mysize),stat=ier)
1471  if(ier/=0) call die(myname_,'allocate(mystarts,..)',ier)
1472
1473! pull out the starts and lengths that PEno owns in the order
1474! they appear in the GSMap.
1475  j=1
1476  do i=1,GSMap%ngseg
1477    if(GSMap%pe_loc(i)==PEno) then
1478      mystarts(j)=GSMap%start(i)
1479      mylengths(j)=GSMap%length(i)
1480      j=j+1
1481    endif
1482  enddo
1483
1484! now recalculate the values of the grid point numbers
1485! based on the starts and lengths
1486! form one long vector which is all local GSMap points
1487  i=1
1488  do j=1,nlsegs
1489    do k=1,mylengths(j)
1490     Points(i)=mystarts(j)+k-1
1491     i=i+1
1492    enddo
1493  enddo
1494
1495  deallocate(mystarts,mylengths, stat=ier)
1496  if(ier/=0) call die(myname_,'deallocate(mystarts,..)',ier)
1497
1498 end subroutine OrderedPoints_
1499
1500!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1501!    Math and Computer Science Division, Argonne National Laboratory   !
1502!BOP -------------------------------------------------------------------
1503!
1504! !IROUTINE: lsize_ - find the local storage size from the map
1505!
1506! !DESCRIPTION:
1507! This function returns the number of points owned by the local process,
1508! as defined by the input {\tt GlobalSegMap} argument {\tt GSMap}.  The
1509! local process ID is determined through use of the input {\tt INTEGER}
1510! argument {\tt comm}, which is the Fortran handle for the MPI
1511! communicator.
1512!
1513! !INTERFACE:
1514
1515 integer function lsize_(GSMap, comm)
1516!
1517! !USES:
1518!
1519      use m_mpif90
1520      use m_die ,          only : MP_perr_die
1521
1522      implicit none
1523
1524! !INPUT PARAMETERS:
1525
1526      type(GlobalSegMap), intent(in) :: GSMap
1527      integer,            intent(in) :: comm
1528
1529
1530! !REVISION HISTORY:
1531!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
1532!       06Feb01 - J.W. Larson <larson@mcs.anl.gov> - Computed directly
1533!                 from the GlobalSegMap, rather than returning a hard-
1534!                 wired local attribute. This required the addition of
1535!                 the communicator argument.
1536!EOP ___________________________________________________________________
1537
1538  character(len=*),parameter :: myname_=myname//'::lsize_'
1539  integer :: ierr, local_size, myID, n, ngseg
1540
1541        ! Determine local rank myID:
1542
1543  call MP_COMM_RANK(comm, myID, ierr)
1544  if(ierr /= 0) call MP_perr_die(myname_,'MP_COMM_RANK',ierr)
1545
1546        ! Determine global number of segments:
1547
1548  ngseg = ngseg_(GSMap)
1549
1550        ! Compute the local size of the distributed vector by summing
1551        ! the entries of GSMap%length(:) whose corresponding values in
1552        ! GSMap%pe_loc(:) equal the local process ID.  This automatically
1553        ! takes into account haloing (if present).
1554
1555  local_size = 0
1556
1557  do n=1,ngseg
1558     if(GSMap%pe_loc(n) == myID) then
1559        local_size = local_size + GSMap%length(n)
1560     endif
1561  end do
1562
1563  lsize_ = local_size
1564
1565 end function lsize_
1566
1567!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1568!    Math and Computer Science Division, Argonne National Laboratory   !
1569!BOP -------------------------------------------------------------------
1570!
1571! !IROUTINE: rank1_ - rank which process owns a datum with given global
1572! index.
1573!
1574! !DESCRIPTION:
1575! This routine assumes that there is one process that owns the datum with
1576! a given global index.  It should not be used when the input
1577! {\tt GlobalSegMap} argument {\tt GSMap} has been built to incorporate
1578! halo points.
1579!
1580! !INTERFACE:
1581
1582    subroutine rank1_(GSMap, i_g, rank)
1583
1584      implicit none
1585
1586! !INPUT PARAMETERS:
1587
1588      type(GlobalSegMap), intent(in)  :: GSMap   ! input GlobalSegMap
1589      integer,            intent(in)  :: i_g     ! a global index
1590
1591! !OUTPUT PARAMETERS:
1592
1593      integer,            intent(out) :: rank    ! the pe on which this
1594                                                 ! element resides
1595! !REVISION HISTORY:
1596!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
1597!EOP ___________________________________________________________________
1598
1599  character(len=*),parameter :: myname_=myname//'::rank1_'
1600  integer :: i,ilc,ile
1601
1602        ! Initially, set the rank to -1 (invalid).
1603  rank=-1       
1604
1605  do i=1,size(GSMap%start)
1606    ilc = GSMap%start(i)
1607    ile = ilc + GSMap%length(i) - 1
1608
1609                ! If i_g in [ilc,ile].  Note that i_g := [1:..]
1610
1611    if(ilc <= i_g .and. i_g <= ile) then
1612      rank = GSMap%pe_loc(i)
1613      return
1614    endif
1615  end do
1616
1617 end subroutine rank1_
1618
1619!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1620!    Math and Computer Science Division, Argonne National Laboratory   !
1621!BOP -------------------------------------------------------------------
1622!
1623! !IROUTINE: rankm_ - rank which processes own a datum with given global
1624! index.
1625!
1626! !DESCRIPTION:
1627! This routine assumes that there may be more than one process that owns
1628! the datum with a given global index.  This routine should  be used when
1629! the input {\tt GlobalSegMap} argument {\tt GSMap} has been built to
1630! incorporate ! halo points.  {\em Nota Bene}:  The output array {\tt rank}
1631! is allocated in this routine and must be deallocated by the routine calling
1632! {\tt rankm\_()}.  Failure to do so could result in a memory leak.
1633!
1634! !INTERFACE:
1635
1636    subroutine rankm_(GSMap, i_g, num_loc, rank)
1637
1638      implicit none
1639
1640! !INPUT PARAMETERS:
1641
1642      type(GlobalSegMap), intent(in)  :: GSMap   ! input GlobalSegMap
1643      integer,            intent(in)  :: i_g     ! a global index
1644
1645! !OUTPUT PARAMETERS:
1646
1647      integer,            intent(out) :: num_loc ! the number of processes
1648                                                 ! which own element i_g
1649      integer, dimension(:), pointer  :: rank    ! the process(es) on which
1650                                                 ! element i_g resides
1651! !REVISION HISTORY:
1652!       29Sep00 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
1653!EOP ___________________________________________________________________
1654
1655  character(len=*),parameter :: myname_=myname//'::rankm_'
1656  integer :: i, ilc, ile, ier, n
1657
1658        ! First sweep:  determine the number of processes num_loc
1659        ! that own the given datum:
1660
1661  num_loc = 0
1662
1663  do i=1,size(GSMap%start)
1664
1665    ilc = GSMap%start(i)
1666    ile = ilc + GSMap%length(i) - 1
1667
1668                ! If i_g in [ilc,ile].  Note that i_g := [1:..]
1669
1670    if(ilc <= i_g .and. i_g <= ile) then
1671      num_loc = num_loc + 1
1672    endif
1673
1674  end do
1675
1676  if(num_loc == 0) then
1677
1678        ! If i_g is nowhere to be found in GSMap, set num_loc to
1679        ! unity and return a null value for rank
1680
1681     num_loc = 1
1682     allocate(rank(num_loc), stat=ier)
1683     rank = -1  ! null value
1684     return
1685
1686  else
1687        ! Allocate output array rank(1:num_loc)
1688
1689     allocate(rank(num_loc), stat=ier)
1690
1691        ! Second sweep:  fill in the entries to rank(:)
1692
1693     n = 0 ! counter
1694
1695     do i=1,size(GSMap%start)
1696
1697        ilc = GSMap%start(i)
1698        ile = ilc + GSMap%length(i) - 1
1699
1700        ! If i_g in [ilc,ile].  Note that i_g := [1:..]
1701
1702        if(ilc <= i_g .and. i_g <= ile) then
1703           n = n + 1
1704           rank(n) = GSMap%pe_loc(i)
1705        endif
1706
1707     end do
1708
1709  endif
1710
1711 end subroutine rankm_
1712
1713!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1714!    Math and Computer Science Division, Argonne National Laboratory   !
1715!BOP -------------------------------------------------------------------
1716!
1717! !IROUTINE: active_pes_ - number of processes that own data.
1718! index.
1719!
1720! !DESCRIPTION:
1721! This routine scans the pe location list of the input {\tt GlobalSegMap}
1722! {\tt GSMap\%pe\_loc(:)}, and counts the number of pe locations that
1723! own at least one datum.  This value is returned in the {\tt INTEGER}
1724! argument {\tt n\_active}.  If the optional {\tt INTEGER} array argument
1725! {\tt list} is included in the call, a sorted list (in ascending order) of
1726! the active processes will be returned.
1727!
1728! {\bf N.B.:} If {\tt active\_pes\_()} is invoked with the optional argument
1729! {\tt pe\_list} included, this routine will allocate and return this array.
1730! The user must deallocate this array once it is no longer needed.  Failure
1731! to do so will result in a memory leak.
1732!
1733! !INTERFACE:
1734
1735    subroutine active_pes_(GSMap, n_active, pe_list)
1736!
1737! !USES:
1738!
1739      use m_die ,          only : die
1740      use m_SortingTools , only : IndexSet
1741      use m_SortingTools , only : IndexSort
1742      use m_SortingTools , only : Permute
1743
1744      implicit none
1745
1746! !INPUT PARAMETERS:
1747
1748      type(GlobalSegMap),    intent(in)        :: GSMap 
1749
1750! !OUTPUT PARAMETERS:
1751
1752      integer,               intent(out)       :: n_active
1753      integer, dimension(:), pointer, optional :: pe_list
1754
1755! !REVISION HISTORY:
1756!       03Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version.
1757!EOP ___________________________________________________________________
1758
1759  character(len=*),parameter :: myname_=myname//'::active_pes_'
1760
1761  integer :: count, i, n, ngseg, ierr
1762  logical :: new
1763  integer, dimension(:), allocatable :: temp_list
1764  integer, dimension(:), allocatable :: perm
1765
1766        ! retrieve total number of segments in the map:
1767
1768  ngseg = ngseg_(GSMap)
1769
1770        ! allocate workspace to tally process id list:
1771
1772  allocate(temp_list(ngseg), stat=ierr)
1773  if(ierr /= 0) call die(myname_,'allocate(temp_list...',ierr)
1774 
1775        ! initialize temp_list to -1 (which can never be a process id)
1776
1777  temp_list = -1
1778
1779        ! initialize the distinct active process count:
1780
1781  count = 0
1782
1783        ! scan entries of GSMap%pe_loc to count active processes:
1784
1785  do n=1,ngseg
1786     if(GSMap%pe_loc(n) >= 0) then ! a legitimate pe_location
1787
1788        ! assume initially that GSMap%pe_loc(n) is a process id previously
1789        ! not encountered
1790
1791        new = .true.
1792
1793        ! test this proposition against the growing list of distinct
1794        ! process ids stored in temp_list(:)
1795
1796        do i=1, count
1797           if(GSMap%pe_loc(n) == temp_list(i)) new = .false.
1798        end do
1799
1800        ! If GSMap%pe_loc(n) represents a previously unencountered
1801        ! process id, increment the count, and add this id to the list
1802
1803        if(new) then
1804           count = count + 1
1805           temp_list(count) = GSMap%pe_loc(n)
1806        endif
1807
1808     else  ! a negative entry in GSMap%pe_loc(n)
1809        ierr = 2
1810        call die(myname_,'negative value of GSMap%pe_loc',ierr)
1811     endif
1812  end do
1813
1814        ! If the argument pe_list is present, we must allocate this
1815        ! array, fill it, and sort it
1816
1817  if(present(pe_list)) then
1818
1819        ! allocate pe_list and permutation array perm
1820
1821     allocate(pe_list(count), perm(count), stat=ierr)
1822     if (ierr /= 0) then
1823        call die(myname_,'allocate(pe_list...',ierr)
1824     endif
1825
1826     do n=1,count
1827        pe_list(n) = temp_list(n)
1828     end do
1829
1830        ! sorting and permutation...
1831
1832     call IndexSet(perm)
1833     call IndexSort(count, perm, pe_list, descend=.false.)
1834     call Permute(pe_list, perm, count)
1835
1836        ! deallocate permutation array...
1837
1838     deallocate(perm, stat=ierr)
1839     if (ierr /= 0) then
1840        call die(myname_,'deallocate(perm)',ierr)
1841     endif
1842
1843  endif ! if(present(pe_list))...
1844
1845        ! deallocate work array temp_list...
1846
1847  deallocate(temp_list, stat=ierr)
1848  if (ierr /= 0) then
1849     call die(myname_,'deallocate(temp_list)',ierr)
1850  endif
1851
1852        ! finally, store the active process count in output variable
1853        ! n_active:
1854
1855  n_active = count
1856
1857 end subroutine active_pes_
1858
1859!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1860!    Math and Computer Science Division, Argonne National Laboratory   !
1861!BOP -------------------------------------------------------------------
1862!
1863! !IROUTINE: peLocs_ - process ID locations for distributed points.
1864! index.
1865!
1866! !DESCRIPTION:
1867! This routine takes an input {\tt INTEGER} array of point indices
1868! {\tt points(:)}, compares them with an input {\tt GlobalSegMap}
1869! {\tt pointGSMap}, and returns the {\em unique} process ID location
1870! for each point.  Note the emphasize on unique.  The assumption here
1871! (which is tested) is that {\tt pointGSMap} is not haloed.  The process
1872! ID locations for the points is returned in the array {\tt pe\_locs(:)}.
1873!
1874! {\bf N.B.:} The test of {\tt pointGSMap} for halo points, and the
1875! subsequent search for the process ID for each point is very slow.  This
1876! first version of the routine is serial.  A parallel version of this
1877! routine will need to be developed.
1878!
1879! !INTERFACE:
1880
1881 subroutine peLocs_(pointGSMap, npoints, points, pe_locs)
1882!
1883! !USES:
1884!
1885      use m_die ,          only : die
1886
1887      implicit none
1888
1889! !INPUT PARAMETERS:
1890
1891      type(GlobalSegMap),    intent(in)   :: pointGSMap 
1892      integer,               intent(in)   :: npoints
1893      integer, dimension(:), intent(in)   :: points
1894
1895! !OUTPUT PARAMETERS:
1896
1897      integer, dimension(:), intent(out)  :: pe_locs
1898
1899! !REVISION HISTORY:
1900!       18Apr01 - J.W. Larson <larson@mcs.anl.gov> - initial version.
1901!EOP ___________________________________________________________________
1902
1903  character(len=*),parameter :: myname_=myname//'::peLocs_'
1904  integer :: ierr
1905  integer :: iseg, ngseg, ipoint
1906  integer :: lower_index, upper_index
1907
1908! Input argument checks:
1909
1910  if(size(points) < npoints) then
1911     ierr = size(points)
1912     call die(myname_,'input points list array too small',ierr)
1913  endif
1914
1915  if(size(pe_locs) < npoints) then
1916     ierr = size(pe_locs)
1917     call die(myname_,'output pe_locs array too small',ierr)
1918  endif
1919
1920  if(haloed_(pointGSMap)) then
1921     ierr = 1
1922     call die(myname_,'input pointGSMap haloed--not valid',ierr)
1923  endif
1924
1925! Brute-force indexing...no assumptions regarding sorting of points(:)
1926! or pointGSMap%start(:)
1927
1928! Number of segments in pointGSMap:
1929
1930  ngseg = ngseg_(pointGSMap)
1931
1932  do ipoint=1,npoints ! loop over points
1933
1934     do iseg=1,ngseg  ! loop over segments
1935
1936        lower_index = pointGSMap%start(iseg)
1937        upper_index = lower_index + pointGSMap%length(iseg) - 1
1938
1939        if((points(ipoint) >= lower_index) .and. &
1940             (points(ipoint) <= upper_index)) then
1941           pe_locs(ipoint) = pointGSMap%pe_loc(iseg)
1942        endif
1943     
1944     end do ! do iseg=1, ngseg
1945  end do ! do ipoint=1,npoints
1946
1947 end subroutine peLocs_
1948
1949!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1950!    Math and Computer Science Division, Argonne National Laboratory   !
1951!BOP -------------------------------------------------------------------
1952!
1953! !IROUTINE: haloed_ - test GlobalSegMap for presence of halo points.
1954! index.
1955!
1956! !DESCRIPTION:
1957! This {\tt LOGICAL} function tests the input {\tt GlobalSegMap}
1958! {\tt GSMap} for the presence of halo points.  Halo points are points
1959! that appear in more than one segment of a {\tt GlobalSegMap}.  If
1960! {\em any} halo point is found, the function {\tt haloed\_()} returns
1961! immediately with value {\tt .TRUE.}  If, after an exhaustive search
1962! of the map has been completed, no halo points are found, the function
1963! {\tt haloed\_()} returns with value {\tt .FALSE.}
1964!
1965! The search algorithm is:
1966!
1967! \begin{enumerate}
1968! \item Extract the segment start and length information from
1969! {\tt GSMap\%start} and {\tt GSMap\%length} into the temporary
1970! arrays {\tt start(:)} and {\tt length(:)}.
1971! \item Sort these arrays in {\em ascending order} keyed by {\tt start}.
1972! \item Scan the arrays {\tt start} and{\tt length}.  A halo point is
1973! present if for at least one value of the index
1974! $1 \leq {\tt n} \leq {\tt GSMap\%ngseg}$
1975! $${\tt start(n)} + {\tt length(n)} - 1 \geq {\tt start(n+1)}$$.
1976! \end{enumerate}
1977!
1978! {\bf N.B.:} Beware that the search for halo points is potentially
1979! expensive. 
1980!
1981! !INTERFACE:
1982
1983    logical function haloed_(GSMap)
1984!
1985! !USES:
1986!
1987      use m_die ,          only : die
1988      use m_SortingTools , only : IndexSet
1989      use m_SortingTools , only : IndexSort
1990      use m_SortingTools , only : Permute
1991
1992      implicit none
1993
1994 ! !INPUT PARAMETERS:
1995
1996     type(GlobalSegMap), intent(in)           :: GSMap 
1997
1998! !REVISION HISTORY:
1999!       08Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version.
2000!       26Apr01 - J.W. Larson <larson@mcs.anl.gov> - Bug fix.
2001!EOP ___________________________________________________________________
2002
2003  character(len=*),parameter :: myname_=myname//'::haloed_'
2004
2005! Error Flag
2006
2007  integer :: ierr
2008
2009! Loop index and storage for number of segments in GSMap
2010
2011  integer :: n, ngseg
2012
2013! Temporary storage for GSMap%start, GSMap%length, and index
2014! permutation array:
2015
2016  integer, dimension(:), allocatable :: start, length, perm
2017
2018! Logical flag indicating segment overlap
2019
2020  logical :: overlap
2021
2022       ! How many segments in GSMap?
2023
2024  ngseg = ngseg_(GSMap)
2025
2026       ! allocate temporary arrays:
2027
2028  allocate(start(ngseg), length(ngseg), perm(ngseg), stat=ierr)
2029  if (ierr /= 0) then
2030     call die(myname_,'allocate(start...',ierr)
2031  endif
2032
2033       ! Fill the temporary arrays start(:) and length(:)
2034
2035  do n=1,ngseg
2036     start(n) = GSMap%start(n)
2037     length(n) = GSMap%length(n)
2038  end do
2039
2040       ! Initialize the index permutation array:
2041
2042  call IndexSet(perm)
2043
2044       ! Create the index permutation that will order the data so the
2045       ! entries of start(:) appear in ascending order:
2046
2047  call IndexSort(ngseg, perm, start, descend=.false.)
2048
2049       ! Permute the data so the entries of start(:) are now in
2050       ! ascending order:
2051
2052  call Permute(start,perm,ngseg)
2053
2054       ! Apply this same permutation to length(:)
2055
2056  call Permute(length,perm,ngseg)
2057
2058       ! Set LOGICAL flag indicating segment overlap to .FALSE.
2059
2060  overlap = .FALSE.
2061
2062       ! Now, scan the segments, looking for overlapping segments.  Upon
2063       ! discovery of the first overlapping pair of segments, set the
2064       ! flag overlap to .TRUE. and exit.
2065
2066  n = 0
2067
2068  SCAN_LOOP: do
2069     n = n + 1
2070     if(n == ngseg) EXIT ! we are finished, and there were no halo pts.
2071     if((start(n) + length(n) - 1) >= start(n+1)) then ! found overlap
2072        overlap = .TRUE.
2073        EXIT
2074     endif
2075  end do SCAN_LOOP
2076
2077       ! Clean up allocated memory:
2078
2079  deallocate(start, length, perm, stat=ierr)
2080  if (ierr /= 0) then
2081     call die(myname_,'deallocate(start...',ierr)
2082  endif
2083
2084       ! Assign function return value:
2085
2086 haloed_ = overlap
2087
2088 end function haloed_
2089
2090!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2091!    Math and Computer Science Division, Argonne National Laboratory   !
2092!BOP -------------------------------------------------------------------
2093!
2094! !IROUTINE: Sort_ - generate index permutation for GlobalSegMap.
2095!
2096! !DESCRIPTION:
2097! {\tt Sort\_()} uses the supplied keys {\tt key1} and {\tt key2} to
2098! generate a permutation {\tt perm} that will put the entries of the
2099! components {\tt GlobalSegMap\%start}, {\tt GlobalSegMap\%length} and
2100! {\tt GlobalSegMap\%pe\_loc} in {\em ascending} lexicographic order.
2101!
2102! {\bf N.B.:} {\tt Sort\_()} returns an allocated array {\tt perm(:)}.  It
2103! the user must deallocate this array once it is no longer needed.  Failure
2104! to do so could create a memory leak.
2105!
2106! !INTERFACE:
2107
2108    subroutine Sort_(GSMap, key1, key2, perm)
2109!
2110! !USES:
2111!
2112      use m_die ,          only : die
2113      use m_SortingTools , only : IndexSet
2114      use m_SortingTools , only : IndexSort
2115
2116      implicit none
2117
2118! !INPUT PARAMETERS:
2119
2120      type(GlobalSegMap),    intent(in)           :: GSMap ! input GlobalSegMap
2121      integer, dimension(:), intent(in)           :: key1  ! first sort key
2122      integer, dimension(:), intent(in), optional :: key2  ! second sort key
2123
2124! !OUTPUT PARAMETERS:
2125
2126      integer, dimension(:), pointer :: perm    ! output index permutation
2127
2128! !REVISION HISTORY:
2129!       02Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version
2130!EOP ___________________________________________________________________
2131
2132  character(len=*),parameter :: myname_=myname//'::Sort_'
2133
2134  integer :: ierr, length
2135
2136  length = ngseg_(GSMap)
2137
2138        ! Argument checking.  are key1 and key2 (if supplied) the
2139        ! same length as the components of GSMap?  If not, stop with
2140        ! an error.
2141
2142  ierr = 0
2143
2144  if(size(key1) /= length) then
2145     ierr = 1
2146     call die(myname_,'key1 GSMap size mismatch',ierr)
2147  endif
2148
2149  if(present(key2)) then
2150     if(size(key2) /= length) then
2151        ierr = 2
2152        call die(myname_,'key2 GSMap size mismatch',ierr)
2153     endif
2154     if(size(key1) /= size(key2)) then
2155        ierr = 3
2156        call die(myname_,'key1 key2 size mismatch',ierr)
2157     endif
2158  endif
2159
2160        ! allocate space for permutation array perm(:)
2161
2162  allocate(perm(length), stat=ierr)
2163  if(ierr /= 0) call die(myname_,'allocate(perm)',ierr)
2164
2165        ! Initialize perm(i)=i, for i=1,length
2166
2167  call IndexSet(perm)
2168 
2169        ! Index permutation is achieved by successive calls to IndexSort(),
2170        ! with the keys supplied one at a time in the order reversed from
2171        ! the desired sort order.
2172
2173  if(present(key2)) then
2174     call IndexSort(length, perm, key2, descend=.false.)
2175  endif
2176
2177  call IndexSort(length, perm, key1, descend=.false.)
2178
2179        ! Yes, it is that simple.  The desired index permutation is now
2180        ! stored in perm(:)
2181
2182 end subroutine Sort_
2183
2184!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2185!    Math and Computer Science Division, Argonne National Laboratory   !
2186!BOP -------------------------------------------------------------------
2187!
2188! !IROUTINE: PermuteInPlace_ - apply index permutation to GlobalSegMap.
2189!
2190! !DESCRIPTION:
2191! {\tt PermuteInPlace\_()} uses a supplied index permutation {\tt perm}
2192! to re-order {\tt GlobalSegMap\%start}, {\tt GlobalSegMap\%length} and
2193! {\tt GlobalSegMap\%pe\_loc}.
2194!
2195! !INTERFACE:
2196
2197    subroutine PermuteInPlace_(GSMap, perm)
2198!
2199! !USES:
2200!
2201      use m_die ,          only : die
2202      use m_SortingTools , only : Permute
2203
2204      implicit none
2205
2206! !INPUT PARAMETERS:
2207
2208      integer, dimension(:), intent(in) :: perm
2209
2210! !INPUT/OUTPUT PARAMETERS:
2211
2212      type(GlobalSegMap), intent(inout) :: GSMap
2213
2214! !REVISION HISTORY:
2215!       02Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version.
2216!EOP ___________________________________________________________________
2217
2218  character(len=*),parameter :: myname_=myname//'::PermuteInPlace_'
2219
2220  integer :: length, ierr
2221
2222  length = ngseg_(GSMap)
2223
2224        ! Argument checking.  Do the components of GSMap
2225        ! (e.g. GSMap%start) have the same length as the
2226        ! permutation array perm?  If not, stop with an error.
2227
2228  ierr = 0
2229
2230  if(size(perm) /= length) then
2231     ierr = 1
2232     call die(myname_,'perm GSMap size mismatch',ierr)
2233  endif
2234
2235        ! In-place index permutation using perm(:) :
2236
2237  call Permute(GSMap%start,perm,length)
2238  call Permute(GSMap%length,perm,length)
2239  call Permute(GSMap%pe_loc,perm,length)
2240
2241        ! Now, the components of GSMap are ordered according to
2242        ! perm(:).
2243
2244 end subroutine PermuteInPlace_
2245
2246!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2247!    Math and Computer Science Division, Argonne National Laboratory   !
2248!BOP -------------------------------------------------------------------
2249!
2250! !IROUTINE: SortPermuteInPlace_ - Sort in-place GlobalSegMap components.
2251!
2252! !DESCRIPTION:
2253! {\tt SortPermuteInPlace\_()} uses a the supplied key(s) to generate
2254! and apply an index permutation that will place the {\tt GlobalSegMap}
2255! components {\tt GlobalSegMap\%start}, {\tt GlobalSegMap\%length} and
2256! {\tt GlobalSegMap\%pe\_loc} in lexicographic order.
2257!
2258! !INTERFACE:
2259
2260    subroutine SortPermuteInPlace_(GSMap, key1, key2)
2261!
2262! !USES:
2263!
2264      use m_die ,          only : die
2265
2266      implicit none
2267
2268! !INPUT PARAMETERS:
2269
2270      integer, dimension(:), intent(in)           :: key1
2271      integer, dimension(:), intent(in), optional :: key2
2272
2273! !INPUT/OUTPUT PARAMETERS:
2274
2275      type(GlobalSegMap),    intent(inout)        :: GSMap
2276
2277! !REVISION HISTORY:
2278!       02Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version.
2279!EOP ___________________________________________________________________
2280
2281  character(len=*),parameter :: myname_=myname//'::SortPermuteInPlace_'
2282
2283  integer :: length, ierr
2284  integer, dimension(:), pointer :: perm
2285
2286   length = ngseg_(GSMap)
2287
2288        ! Argument checking.  are key1 and key2 (if supplied) the
2289        ! same length as the components of GSMap?  If not, stop with
2290        ! an error.
2291  ierr = 0
2292  if(size(key1) /= length) then
2293     ierr = 1
2294     call die(myname_,'key1 GSMap size mismatch',ierr)
2295  endif
2296
2297  if(present(key2)) then
2298     if(size(key2) /= length) then
2299        ierr = 2
2300        call die(myname_,'key2 GSMap size mismatch',ierr)
2301     endif
2302     if(size(key1) /= size(key2)) then
2303        ierr = 3
2304        call die(myname_,'key1 key2 size mismatch',ierr)
2305     endif
2306  endif
2307
2308        ! Generate desired index permutation:     
2309
2310  if(present(key2)) then
2311     call Sort_(GSMap, key1, key2, perm)
2312  else
2313     call Sort_(GSMap, key1=key1, perm=perm)
2314  endif
2315
2316        ! Apply index permutation:     
2317
2318  call PermuteInPlace_(GSMap, perm)
2319
2320        ! Now the components of GSMap have been re-ordered.
2321        ! Deallocate the index permutation array perm(:)
2322
2323  deallocate(perm, stat=ierr)
2324  if(ierr /= 0) call die(myname_,'deallocate(perm...)',ierr)
2325
2326 end subroutine SortPermuteInPlace_
2327
2328
2329!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2330!    Math and Computer Science Division, Argonne National Laboratory   !
2331!BOP -------------------------------------------------------------------
2332!
2333! !IROUTINE: increasing_ - Return .TRUE. if GSMap has increasing indices
2334!
2335! !DESCRIPTION:
2336! The function {\tt increasing\_()} returns .TRUE. if each proc's
2337! indices in the {\tt GlobalSegMap} argument {\tt GSMap} have
2338! strictly increasing indices.  I.e. the proc's segments have indices
2339! in ascending order and are non-overlapping.
2340!
2341! !INTERFACE:
2342
2343 logical function increasing_(gsmap)
2344
2345! !USES:
2346      use m_MCTWorld, only: ThisMCTWorld
2347      use m_die
2348
2349      implicit none
2350
2351! !INPUT PARAMETERS:
2352
2353      type(GlobalSegMap),intent(in) :: gsmap
2354
2355! !REVISION HISTORY:
2356!       06Jun07 - R. Loy <rloy@mcs.anl.gov> - initial version
2357!EOP ___________________________________________________________________
2358
2359  character(len=*),parameter :: myname_=myname//'::increasing_'
2360
2361  integer comp_id
2362  integer nprocs
2363  integer i
2364  integer this_ngseg
2365  integer ier
2366  integer, allocatable:: last_index(:)
2367  integer pe_loc
2368
2369  comp_id = gsmap%comp_id
2370  nprocs=ThisMCTWorld%nprocspid(comp_id)
2371
2372  allocate( last_index(nprocs), stat=ier )
2373  if (ier/=0) call die(myname_,'allocate last_index')
2374
2375  last_index= -1
2376  increasing_ = .TRUE.
2377  this_ngseg=ngseg(gsmap)
2378
2379  iloop: do i=1,this_ngseg
2380    pe_loc=gsmap%pe_loc(i)+1  ! want value 1..nprocs
2381    if (gsmap%start(i) <= last_index(pe_loc)) then
2382      increasing_ = .FALSE.
2383      exit iloop
2384    endif
2385    last_index(pe_loc)=gsmap%start(i)+gsmap%length(i)-1
2386  enddo iloop
2387
2388  deallocate( last_index, stat=ier )
2389  if (ier/=0) call die(myname_,'deallocate last_index')
2390
2391 end function increasing_
2392
2393
2394!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2395!    Math and Computer Science Division, Argonne National Laboratory   !
2396!BOP -------------------------------------------------------------------
2397!
2398! !IROUTINE: copy_ - Copy the gsmap to a new gsmap
2399!
2400! !DESCRIPTION:
2401! Make a copy of a gsmap.
2402! Note this is a deep copy of all arrays.
2403!
2404! !INTERFACE:
2405
2406  subroutine copy_(src,dest)
2407
2408! !USES:
2409      use m_MCTWorld, only: ThisMCTWorld
2410      use m_die
2411
2412      implicit none
2413
2414! !INPUT PARAMETERS:
2415
2416      type(GlobalSegMap),intent(in) :: src
2417
2418! !OUTPUT PARAMETERS:
2419
2420      type(GlobalSegMap),intent(out) :: dest
2421
2422
2423! !REVISION HISTORY:
2424!       27Jul07 - R. Loy <rloy@mcs.anl.gov> - initial version
2425!EOP ___________________________________________________________________
2426
2427
2428    call initp_( dest, src%comp_id, src%ngseg, src%gsize, &
2429                 src%start, src%length, src%pe_loc )
2430
2431  end subroutine copy_
2432
2433 end module m_GlobalSegMap
2434
Note: See TracBrowser for help on using the repository browser.