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

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!-----------------------------------------------------------------------
! CVS m_GlobalToLocal.F90,v 1.16 2006-04-20 23:54:48 rloy Exp
! CVS MCT_2_8_0 
!BOP -------------------------------------------------------------------
!
! !MODULE: m_GlobalToLocal - Global to Local Index Translation
!
! !DESCRIPTION:
! This module contains routines for translating global array indices 
! into their local counterparts (that is, the indices into the local 
! data structure holding a given process' chunk of a distributed array).
! The MCT domain decomposition descriptors {\tt GlobalMap} and 
! {\tt GlobalSegMap} are both supported.   Indices can be translated 
! one-at-a-time using the {\tt GlobalToLocalIndex} routine or many 
! at once using the {\tt GlobalToLocalIndices} routine.
!
! This module also provides facilities for setting the local row and 
! column indices for a {\tt SparseMatrix} through the 
! {\tt GlobalToLocalMatrix} routines.
!
! !INTERFACE:

 module m_GlobalToLocal

! !USES:
! No external modules are used in the declaration section of this module.

      implicit none

      private   ! except

! !PUBLIC MEMBER FUNCTIONS:

      public :: GlobalToLocalIndex   ! Translate Global to Local index
                                     ! (i.e. recover local index for a
                                     ! point from its global index). 

      public :: GlobalToLocalIndices ! Translate Global to Local indices
                                     ! (i.e. recover local starts/lengths 
                                     ! of distributed data segments).
                                     
      public :: GlobalToLocalMatrix  ! Re-indexing of row or column
                                     ! indices for a SparseMatrix

    interface GlobalToLocalIndices ; module procedure   &
       GlobalSegMapToIndices_,  &   ! local arrays of starts/lengths
       GlobalSegMapToNavigator_, &  ! return local indices as Navigator
       GlobalSegMapToIndexArr_
    end interface

    interface GlobalToLocalIndex ; module procedure &
       GlobalSegMapToIndex_, &
       GlobalMapToIndex_
    end interface

    interface GlobalToLocalMatrix ; module procedure &
       GlobalSegMapToLocalMatrix_
    end interface


! !SEE ALSO:
!
! The MCT modules {\tt m\_GlobalMap} and {m\_GlobalSegMap} for more 
! information regarding MCT's domain decomposition descriptors.
!
! The MCT module {\tt m\_SparseMatrix} for more information regarding 
! the {\tt SparseMatrix} datatype.
!
! !REVISION HISTORY:
!  2Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname='MCT::m_GlobalToLocal'

 contains

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: GlobalSegMapToIndices_ - Return _local_ indices in arrays.
!
! !DESCRIPTION:  {\tt GlobalSegMapToIndices\_()} takes a user-supplied
! {\tt GlobalSegMap} data type {\tt GSMap}, which desribes a decomposition 
! on the input MPI communicator corresponding to the Fortran {\tt INTEGER} 
! handle {\tt comm} to translate the global directory of segment locations 
! into local indices for referencing the on-pe storage of the mapped 
! distributed data.
!
! {\bf N.B.:}  This routine returns two allocated arrays---{\tt start(:)} 
! and {\tt length(:)}---which must be deallocated once the user no longer
! needs them.  Failure to do this will create a memory leak.
!
! !INTERFACE:

 subroutine GlobalSegMapToIndices_(GSMap, comm, start, length)

!
! !USES:
!
      use m_mpif90
      use m_die,          only : MP_perr_die, die, warn
      use m_GlobalSegMap, only : GlobalSegMap
      use m_GlobalSegMap, only : GlobalSegMap_ngseg => ngseg
      use m_GlobalSegMap, only : GlobalSegMap_nlseg => nlseg

      implicit none

! !INPUT PARAMETERS:

      type(GlobalSegMap),   intent(in) :: GSMap ! Output GlobalSegMap
      integer,              intent(in) :: comm  ! communicator handle

! !OUTPUT PARAMETERS:

      integer,dimension(:), pointer :: start  ! local segment start indices
      integer,dimension(:), pointer :: length ! local segment sizes

! !REVISION HISTORY:
!  2Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::GlobalSegMapToIndices_'

  integer :: myID, ierr, ngseg, nlseg, n, count
 
          ! determine local process id myID

  call MP_COMM_RANK(comm, myID, ierr)
  if(ierr /= 0) call MP_perr_die(myname_,'MP_COMM_RANK',ierr)

          ! determine number of global segments ngseg:

  ngseg = GlobalSegMap_ngseg(GSMap)

          ! determine number of local segments on process myID nlseg:

  nlseg = GlobalSegMap_nlseg(GSMap, myID)

          ! allocate arrays start(:) and length(:) to store local
          ! segment information.

  allocate(start(nlseg), length(nlseg), stat=ierr)
  if(ierr /= 0) call die(myname_,'allocate(start...',ierr)

          ! Loop over GlobalSegMap%pe_loc(:) values to isolate
          ! global index values of local data.  Record number of
          ! matches in the INTEGER count.

  count = 0
  do n=1, ngseg
     if(GSMap%pe_loc(n) == myID) then
        count = count + 1
        if(count > nlseg) then
           ierr = 2
           call die(myname_,'too many pe matches',ierr)
        endif
        start(count) = GSMap%start(n)
        length(count) = GSMap%length(n)
     endif
  end do

  if(count < nlseg) then
     ierr = 3
     call die(myname_,'too few pe matches',ierr)
  endif

          ! translate global start indices to their local 
          ! values, based on their storage order and number
          ! of elements in each segment

  do n=1, count
     if(n == 1) then
        start(n) = 1
     else
        start(n) = start(n-1) + length(n-1)
     endif
  end do

 end subroutine GlobalSegMapToIndices_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: GlobalSegMapToIndex_ - Global to Local Index Translation
!
! !DESCRIPTION:  This {\tt INTEGER} query function takes a user-supplied
! {\tt GlobalSegMap} data type {\tt GSMap}, which desribes a decomposition 
! on the input MPI communicator corresponding to the Fortran {\tt INTEGER} 
! handle {\tt comm}, and the input global index value {\tt i\_g}, and 
! returns a positive local index value if the datum {\tt i\_g}.   If 
! the datum {\tt i\_g} is not stored on the local process ID, a value 
! of {\tt -1} is returned.
!
! !INTERFACE:


 integer function GlobalSegMapToIndex_(GSMap, i_g, comm)

!
! !USES:
!
      use m_mpif90
      use m_die,          only : MP_perr_die, die, warn
      use m_GlobalSegMap, only : GlobalSegMap
      use m_GlobalSegMap, only : GlobalSegMap_ngseg => ngseg
      use m_GlobalSegMap, only : GlobalSegMap_nlseg => nlseg

      implicit none

! !INPUT PARAMETERS:

      type(GlobalSegMap), intent(in)  :: GSMap ! Output GlobalSegMap
      integer,            intent(in)  :: i_g   ! global index
      integer,            intent(in)  :: comm  ! communicator handle

! !REVISION HISTORY:
!  2Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::GlobalSegMapToIndex_'

  integer :: myID
  integer :: count, ierr, ngseg, nlseg, n
  integer :: lower_bound, upper_bound
  integer :: local_start, local_index
  logical :: found

  ! Determine local process id myID:

  call MP_COMM_RANK(comm, myID, ierr)
  if(ierr /= 0) call MP_perr_die(myname_,'MP_COMM_RANK()',ierr)

  ! Extract the global number of segments in GSMap

  ngseg = GlobalSegMap_ngseg(GSMap)

  ! Extract the global number of segments in GSMap for myID

  nlseg = GlobalSegMap_nlseg(GSMap, myID)

  ! set the counter count, which records the number of times myID
  ! matches entries in GSMap%pe_loc(:)

  count = 0

  ! set local_start, which is the current local storage segment
  ! starting position

  local_start = 1

  ! set logical flag found to signify we havent found i_g:

  found = .false.

  n = 0

  SEARCH_LOOP: do 
     
     n = n+1
     if (n > ngseg) EXIT

     if(GSMap%pe_loc(n) == myID) then

  ! increment / check the pe_loc match counter

        count = count + 1
        if(count > nlseg) then
           ierr = 2
           call die(myname_,'too many pe matches',ierr)
        endif

  ! is i_g in this segment?

        lower_bound = GSMap%start(n)
        upper_bound = GSMap%start(n) + GSMap%length(n) - 1

        if((lower_bound <= i_g) .and. (i_g <= upper_bound)) then
           local_index = local_start + (i_g - GSMap%start(n))
           found = .true.
           EXIT
        else
           local_start = local_start + GSMap%length(n)
        endif

     endif
  end do SEARCH_LOOP

  ! We either found the local index, or have exhausted our options.

  if(found) then
     GlobalSegMapToIndex_ = local_index
  else
     GlobalSegMapToIndex_ = -1
  endif

 end function GlobalSegMapToIndex_



!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: GlobalSegMapToIndexArr_ - Global to Local Index Array Translation
!
! !DESCRIPTION:  Given a {\tt GlobalSegMap} data type {\tt GSMap}
! and MPI communicator corresponding to the Fortran {\tt INTEGER} 
! handle {\tt comm}, convert an array of global index values
! {\tt i\_global()} to an array of local index values {\tt i\_local()}.  If 
! the datum {\tt i\_global(j)} is not stored on the local process ID,
! then {\tt i\_local(j)} will be set to {\tt -1}/
!
! !INTERFACE:


subroutine GlobalSegMapToIndexArr_(GSMap, i_global, i_local, nindex, comm)

!
! !USES:
!
      use m_stdio
      use m_mpif90
      use m_die,          only : MP_perr_die, die, warn
      use m_GlobalSegMap, only : GlobalSegMap
      use m_GlobalSegMap, only : GlobalSegMap_ngseg => ngseg
      use m_GlobalSegMap, only : GlobalSegMap_nlseg => nlseg

      implicit none

! !INPUT PARAMETERS:

      type(GlobalSegMap), intent(in)  :: GSMap ! Output GlobalSegMap
      integer,            intent(in)  :: i_global(:)   ! global index
      integer,            intent(out) :: i_local(:)    ! local index
      integer,            intent(in)  :: nindex          ! size of i_global()
      integer,            intent(in)  :: comm  ! communicator handle

! !REVISION HISTORY:
!  12-apr-2006   R. Loy <rloy@mcs.anl.gov> - initial version
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::GlobalSegMapToIndexArr_'

  integer :: myID
  integer :: count, ierr, ngseg, nlseg
  integer,allocatable  :: mygs_lb(:),mygs_ub(:),mygs_len(:),mygs_lstart(:)

  integer :: i,j,n,startj

  ! Determine local process id myID:

  call MP_COMM_RANK(comm, myID, ierr)
  if(ierr /= 0) call MP_perr_die(myname_,'MP_COMM_RANK()',ierr)


  ngseg = GlobalSegMap_ngseg(GSMap)
  nlseg = GlobalSegMap_nlseg(GSMap, myID)

  if (nlseg <= 0) return;

  allocate( mygs_lb(nlseg), mygs_ub(nlseg), mygs_len(nlseg) )
  allocate( mygs_lstart(nlseg) )


!! 
!! determine the global segments on this processor 
!! just once, so the info be used repeatedly below
!!

  n = 0
  do i=1,ngseg
    if (GSMap%pe_loc(i) == myID ) then
      n=n+1
      mygs_lb(n)=GSMap%start(i)
      mygs_ub(n)=GSMap%start(i) + GSMap%length(i) -1
      mygs_len(n)=GSMap%length(i)
    endif
  enddo

  if (n .ne. nlseg) then
    write(stderr,*) myname_,"mismatch nlseg",n,nlseg
    call die(myname)
  endif

  mygs_lstart(1)=1
  do j=2,nlseg
    mygs_lstart(j)=mygs_lstart(j-1)+mygs_len(j-1)
  enddo


!! 
!! this loop is optimized for the case that the indices in iglobal()
!! are in the same order that they appear in the global segments,
!! which seems usually (always?) to be the case. 
!!
!! note that the j loop exit condition is only executed when the index
!! is not found in the current segment, which saves a factor of 2
!! since many consecutive indices are in the same segment.
!! 


  j=1
  do i=1,nindex

    i_local(i)= -1

    startj=j
    SEARCH_LOOP: do

      if ( (mygs_lb(j) <= i_global(i)) .and. &
           (i_global(i) <= mygs_ub(j))) then
        i_local(i) = mygs_lstart(j) + (i_global(i) - mygs_lb(j))
        EXIT SEARCH_LOOP
      else
        j=j+1
        if (j > nlseg) j=1                      ! wrap around
        if (j == startj) EXIT SEARCH_LOOP
      endif

    end do SEARCH_LOOP

  end do

!!!! this version vectorizes (outer loop)
!!!! performance for in-order input is slightly slower than the above
!!!! but performance on out-of-order input is probably much better
!!!! at the moment we are going on the assumption that caller is
!!!! likely providing in-order, so we won't use this version.
!!
!!  do i=1,nindex
!!
!!    i_local(i)= -1
!!
!!    SEARCH_LOOP: do j=1,nlseg
!!
!!      if ( (mygs_lb(j) <= i_global(i)) .and. &
!!           (i_global(i) <= mygs_ub(j))) then
!!        i_local(i) = mygs_lstart(j) + (i_global(i) - mygs_lb(j))
!!      endif
!!
!!    end do SEARCH_LOOP
!!
!!  end do


  deallocate( mygs_lb, mygs_ub, mygs_len, mygs_lstart )

 end subroutine GlobalSegMapToIndexArr_





!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: GlobalMapToIndex_ - Global to Local Index Translation
!
! !DESCRIPTION:  
! This {\tt INTEGER} query function takes as its input a user-supplied
! {\tt GlobalMap} data type {\tt GMap}, which desribes a decomposition 
! on the input MPI communicator corresponding to the Fortran {\tt INTEGER} 
! handle {\tt comm}, and the input global index value {\tt i\_g}, and 
! returns a positive local index value if the datum {\tt i\_g}.   If 
! the datum {\tt i\_g} is not stored on the local process ID, a value 
! of {\tt -1} is returned.
!
! !INTERFACE:


 integer function GlobalMapToIndex_(GMap, i_g, comm)

!
! !USES:
!
      use m_mpif90
      use m_die,          only : MP_perr_die, die, warn
      use m_GlobalMap,    only : GlobalMap

      implicit none

! !INPUT PARAMETERS:

      type(GlobalMap), intent(in)  :: GMap  ! Input GlobalMap
      integer,         intent(in)  :: i_g   ! global index
      integer,         intent(in)  :: comm  ! communicator handle

! !REVISION HISTORY:
!  2Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::GlobalMapToIndex_'

  integer :: myID
  integer :: count, ierr, ngseg, nlseg, n
  integer :: lower_bound, upper_bound
  integer :: local_start, local_index
  logical :: found

  ! Determine local process id myID:

  call MP_COMM_RANK(comm, myID, ierr)
  if(ierr /= 0) call MP_perr_die(myname_,'MP_COMM_RANK()',ierr)

  ! Initialize logical "point located" flag found as false

  found = .false.

  lower_bound = GMap%displs(myID) + 1
  upper_bound = GMap%displs(myID) + GMap%counts(myID) 

  if((lower_bound <= i_g) .and. (i_g <= upper_bound)) then
     found = .true.
     local_index = i_g - lower_bound + 1
  endif

  if(found) then
     GlobalMapToIndex_ = local_index
  else
     GlobalMapToIndex_ = -1
  endif

 end function GlobalMapToIndex_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: GlobalSegMapToNavigator_ - Return Navigator to Local Segments
!
! !DESCRIPTION:  
! This routine takes as its input takes a user-supplied
! {\tt GlobalSegMap} data type {\tt GSMap}, which desribes a decomposition 
! on the input MPI communicator corresponding to the Fortran {\tt INTEGER} 
! handle {\tt comm}, and returns the local segment start index and length 
! information for referencing the on-pe storage of the mapped distributed
! data.  These data are returned in the form of the output {\tt Navigator} 
! argument {Nav}.
!
! {\bf N.B.:}  This routine returns a {\tt Navigator} variable {\tt Nav},
! which must be deallocated once the user no longer needs it.  Failure to 
! do this will create a memory leak.
!
! !INTERFACE:

 subroutine GlobalSegMapToNavigator_(GSMap, comm, oNav)

!
! !USES:
!
      use m_mpif90
      use m_die,          only : MP_perr_die, die, warn
      use m_GlobalSegMap, only : GlobalSegMap
      use m_GlobalSegMap, only : GlobalSegMap_ngseg => ngseg
      use m_GlobalSegMap, only : GlobalSegMap_nlseg => nlseg
      use m_Navigator, only    : Navigator
      use m_Navigator, only    : Navigator_init => init

      implicit none

! !INPUT PARAMETERS:

      type(GlobalSegMap), intent(in)  :: GSMap ! Input GlobalSegMap
      integer,            intent(in)  :: comm  ! communicator handle

! !OUTPUT PARAMETERS:

      type(Navigator),    intent(out) :: oNav   ! Output Navigator

! !REVISION HISTORY:
!  2Feb01 - J.W. Larson <larson@mcs.anl.gov> - initial version
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::GlobalSegMapToNavigator_'

  integer :: myID, ierr, ngseg, nlseg, n, count
 
          ! determine local process id myID

  call MP_COMM_RANK(comm, myID, ierr)
  if(ierr /= 0) call MP_perr_die(myname_,'MP_COMM_RANK',ierr)

          ! determine number of global segments ngseg:

  ngseg = GlobalSegMap_ngseg(GSMap)

          ! determine number of local segments on process myID nlseg:

  nlseg = GlobalSegMap_nlseg(GSMap, myID)

          ! Allocate space for the Navigator oNav:

  call Navigator_init(oNav, nlseg, ierr)
  if(ierr /= 0) call die(myname_,'Navigator_init',ierr)

  call GlobalSegMapToIndices_(GSMap, comm, oNav%displs, oNav%counts)

 end subroutine GlobalSegMapToNavigator_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: GlobalSegMapToLocalMatrix_ - Set Local SparseMatrix Indices
!
! !DESCRIPTION:  
! This routine takes as its input a user-supplied {\tt GlobalSegMap} 
! domain decomposition {\tt GSMap}, which describes the decomposition of 
! either the rows or columns of the input/output {\tt SparseMatrix} 
! argument {\tt sMat} on the communicator associated with the {\tt INTEGER}
! handle {\tt comm}, and to translate the global row or column indices 
! of {\tt sMat} into their local counterparts.  The choice of either row
! or column is governed by the value of the input {\tt CHARACTER} 
! argument {\tt RCFlag}.  One sets this variable to either {\tt 'ROW'} or 
! {\tt 'row'} to specify row re-indexing (which are stored in
! {\tt sMat} and retrieved by indexing the attribute {\tt lrow}), and 
! {\tt 'COLUMN'} or {\tt 'column'} to specify column re-indexing (which 
! are stored in {\tt sMat} and retrieved by indexing the {\tt SparseMatrix} 
! attribute {\tt lcol}).
!
! !INTERFACE:

 subroutine GlobalSegMapToLocalMatrix_(sMat, GSMap, RCFlag, comm)

!
! !USES:
!
      use m_stdio
      use m_die,          only : die

      use m_SparseMatrix, only : SparseMatrix
      use m_SparseMatrix, only : SparseMatrix_indexIA => indexIA
      use m_SparseMatrix, only : SparseMatrix_lsize => lsize

      use m_GlobalSegMap, only : GlobalSegMap


      implicit none

! !INPUT PARAMETERS:

      type(GlobalSegMap), intent(in)    :: GSMap  ! Input GlobalSegMap
      character(len=*),   intent(in)    :: RCFlag ! 'row' or 'column'
      integer,            intent(in)    :: comm   ! communicator handle

! !INPUT/OUTPUT PARAMETERS:

      type(SparseMatrix), intent(inout) :: sMat

! !SEE ALSO:
! The MCT module m_SparseMatrix for more information about the 
! SparseMatrix type and its storage of global and local row-and
! column indices.
!
! !REVISION HISTORY:
!  3May01 - J.W. Larson <larson@mcs.anl.gov> - initial version, which
!           is _extremely_ slow, but safe.  This must be re-examined
!           later.
!EOP ___________________________________________________________________

  character(len=*),parameter :: myname_=myname//'::GlobalSegMapToLocalMatrix_'


  integer :: i, GlobalIndex, gindex, lindex, lsize

   integer, allocatable :: temp_gindex(:)  !! rml
   integer, allocatable :: temp_lindex(:)  !! rml


       ! What are we re-indexing, rows or columns?

  select case(RCFlag)
  case('ROW','row')
     gindex = SparseMatrix_indexIA(sMat, 'grow', dieWith=myname_)
     lindex = SparseMatrix_indexIA(sMat,'lrow', dieWith=myname_)
  case('COLUMN','column')
     gindex = SparseMatrix_indexIA(sMat,'gcol', dieWith=myname_)
     lindex = SparseMatrix_indexIA(sMat,'lcol', dieWith=myname_)
  case default
     write(stderr,'(3a)') myname_,":: unrecognized value of RCFLag ",RCFlag
     call die(myname)
  end select


       ! How many matrix elements are there?

  lsize = SparseMatrix_lsize(sMat)


  !! rml new code from here down - do the mapping all in one
  !! function call which has been tuned for speed

  allocate( temp_gindex(lsize) )
  allocate( temp_lindex(lsize) )


  do i=1,lsize
     temp_gindex(i) = sMat%data%iAttr(gindex,i)
  end do

  call GlobalSegMapToIndexArr_(GSMap, temp_gindex, temp_lindex, lsize, comm)  

  do i=1,lsize
    sMat%data%iAttr(lindex,i) = temp_lindex(i)
  end do


  deallocate(temp_gindex)  ! rml
  deallocate(temp_lindex)  ! rml


 end subroutine GlobalSegMapToLocalMatrix_

 end module m_GlobalToLocal