source: CONFIG_DEVT/IPSLCM6.5_work_ENSEMBLES/oasis3-mct/BLD/build/lib/mctdir/mct/m_GlobalSegMap.F90 @ 5725

Last change on this file since 5725 was 5725, checked in by aclsce, 3 years ago

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