ecearth3/sources/oasis3-mct/lib/mct/mct/m_SparseMatrixToMaps.F90

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!-----------------------------------------------------------------------
! CVS m_SparseMatrixToMaps.F90,v 1.13 2007-05-11 02:09:48 rloy Exp
! CVS MCT_2_8_0 
!BOP -------------------------------------------------------------------
!
! !MODULE: m_SparseMatrixToMaps -- Maps from the Sparse Matrix 
!
! !DESCRIPTION:
! The {\tt SparseMatrix} provides consolidated (on one process) or 
! distributed sparse matrix storage for the operation
! ${\bf y} = {\bf M} {\bf x}$, where {\bf x} and {\bf y} are vectors, 
! and {\bf M} is a matrix.  In performing parallel matrix-vector 
! multiplication, one has numerous options regarding the decomposition 
! of the matrix {\bf M}, and the vectors {\bf y} and {\bf x}.
! This module provides services to generate mct mapping components---the
! {\tt GlobalMap} and {\tt GlobalSegMap} for the vectors {\bf y} and/or 
! {\bf x} based on the decomposition of the sparse matrix {\bf M}.
!
! !INTERFACE:

 module m_SparseMatrixToMaps
!
! !USES:
!
      use m_SparseMatrix, only : SparseMatrix

      implicit none

      private   ! except

      public :: SparseMatrixToXGlobalSegMap
      public :: SparseMatrixToYGlobalSegMap

    interface SparseMatrixToXGlobalSegMap ; module procedure &
	 SparseMatrixToXGlobalSegMap_
    end interface

    interface SparseMatrixToYGlobalSegMap ; module procedure &
	 SparseMatrixToYGlobalSegMap_
    end interface

! !REVISION HISTORY:
! 13Apr01 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
!           and API specifications.
!EOP ___________________________________________________________________

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

 contains

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: SparseMatrixToXGlobalSegMap_ - Generate X GlobalSegmap.
!
! !DESCRIPTION:  Given an input {\tt SparseMatrix} argument {\tt sMat},
! this routine generates an output {\tt GlobalSegMap} variable 
! {\tt xGSMap}, which describes the domain decomposition of the vector 
! {\bf x} in the distributed matrix-vector multiplication 
! $${\bf y} = {\bf M} {\bf x}.$$
!
! !INTERFACE:

 subroutine SparseMatrixToXGlobalSegMap_(sMat, xGSMap, root, comm, comp_id)
!
! !USES:
!
      use m_stdio, only : stderr
      use m_die,   only : die
      use m_mpif90

      use m_List, only : List
      use m_List, only : List_init => init
      use m_List, only : List_clean => clean

      use m_SparseMatrix, only : SparseMatrix
      use m_SparseMatrix, only : SparseMatrix_nCols => nCols
      use m_SparseMatrix, only : SparseMatrix_lsize => lsize
      use m_SparseMatrix, only : SparseMatrix_indexIA => indexIA
      use m_SparseMatrix, only : SparseMatrix_SortPermute => SortPermute

      use m_GlobalSegMap, only : GlobalSegMap
      use m_GlobalSegMap, only : GlobalSegMap_init => init

      implicit none

! !INPUT PARAMETERS: 
!
      integer,            intent(in)    :: root    ! communicator root
      integer,            intent(in)    :: comm    ! communicator handle
      integer,            intent(in)    :: comp_id ! component id

! !INPUT/OUTPUT PARAMETERS: 
!
      type(SparseMatrix), intent(inout) :: sMat    ! input SparseMatrix

! !OUTPUT PARAMETERS: 
!
      type(GlobalSegMap), intent(out)   :: xGSMap  ! segmented decomposition
                                                   ! for x
! !REVISION HISTORY:
! 13Apr01 - J.W. Larson <larson@mcs.anl.gov> - API specification.
! 25Apr01 - J.W. Larson <larson@mcs.anl.gov> - First version.
! 27Apr01 - J.W. Larson <larson@mcs.anl.gov> - Bug fix--intent of
!           argument sMat changed from (IN) to (INOUT)
! 27Apr01 - R.L. Jacob <jacob@mcs.anl.gov> - bug fix-- add use 
!           statement for SortPermute
! 01May01 - R.L. Jacob <jacob@mcs.anl.gov> - make comp_id an
!           input argument
!EOP ___________________________________________________________________

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

! SparseMatrix attributes:
  integer :: lsize
! GlobalSegMap input attributes:
  integer :: gsize, ngseg
  integer, dimension(:), pointer :: starts, lengths
! Temporary array for identifying each matrix element column and 
! process ID destination
  integer, dimension(:), allocatable :: gCol, element_pe_locs 
! Index to identify the gcol attribute in sMat:
  integer :: igCol
! Matrix element sorting keys list:
  type(List) :: sort_keys
! Loop index and error flag:
  integer :: i, ierr

       ! Determine he local number of matrix elements lsize

  lsize = SparseMatrix_lsize(sMat)

       ! The value of gsize is taken from the number of columns in sMat:

  gsize = SparseMatrix_nCols(sMat)

       ! Sort SparseMatrix entries by global column index gcol, then
       ! global row index.

       ! Create Sort keys list

  call List_init(sort_keys,'gcol:grow')

       ! Sort and permute the entries of sMat into lexicographic order
       ! by global column, then global row.

  call SparseMatrix_SortPermute(sMat, sort_keys)

       ! Clean up sort keys list

  call List_clean(sort_keys)

       ! Allocate storage space for matrix element column indices and
       ! process ID destinations

  allocate(gCol(lsize), stat=ierr)

  if(ierr /= 0) then
     call die(myname_,'allocate(gCol...',ierr)
  endif

       ! Extract global column information and place in array gCol

  igCol = SparseMatrix_indexIA(sMat, 'gcol', dieWith=myname_)

  do i=1, lsize
     gCol(i) = sMat%data%iAttr(igCol,i)
  end do

       ! Scan sorted entries of gCol to count segments (ngseg), and
       ! their starting indices and lengths (returned in the arrays
       ! starts(:) and lengths(:), respectively)

  call ComputeSegments_(gCol, lsize, ngseg, starts, lengths)

       ! Now we have sufficient data to call the GlobalSegMap
       ! initialization using distributed data:

  call GlobalSegMap_init(xGSMap, starts, lengths, root, comm, &
                         comp_id, gsize=gsize)
  
       ! clean up temporary arrays gCol(:), starts(:) and lengths(:),
       ! (the latter two were allocated in the call to the routine
       ! ComputeSegments_())

  deallocate(gCol, starts, lengths, stat=ierr)

  if(ierr /= 0) then
     call die(myname_,'deallocate(gCol...',ierr)
  endif

 end subroutine SparseMatrixToXGlobalSegMap_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: SparseMatrixToYGlobalSegMap_ - Generate Y GlobalSegmap.
!
! !DESCRIPTION:  Given an input {\tt SparseMatrix} argument {\tt sMat},
! this routine generates an output {\tt GlobalSegMap} variable 
! {\tt yGSMap}, which describes the domain decomposition of the vector 
! {\bf y} in the distributed matrix-vector multiplication 
! ${\bf y} = {\bf M} {\bf x}$.
!
! !INTERFACE:

 subroutine SparseMatrixToYGlobalSegMap_(sMat, yGSMap, root, comm, comp_id)
!
! !USES:
!
      use m_stdio, only : stderr
      use m_die,   only : die

      use m_List, only : List
      use m_List, only : List_init => init
      use m_List, only : List_clean => clean

      use m_SparseMatrix, only : SparseMatrix
      use m_SparseMatrix, only : SparseMatrix_nRows => nRows
      use m_SparseMatrix, only : SparseMatrix_lsize => lsize
      use m_SparseMatrix, only : SparseMatrix_indexIA => indexIA
      use m_SparseMatrix, only : SparseMatrix_SortPermute => SortPermute

      use m_GlobalSegMap, only : GlobalSegMap
      use m_GlobalSegMap, only : GlobalSegMap_init => init

      implicit none

! !INPUT PARAMETERS: 
!
      integer,            intent(in)    :: root    ! communicator root
      integer,            intent(in)    :: comm    ! communicator handle
      integer,            intent(in)    :: comp_id ! component id

! !INPUT/OUTPUT PARAMETERS: 
!
      type(SparseMatrix), intent(inout) :: sMat    ! input SparseMatrix

! !OUTPUT PARAMETERS: 
!
      type(GlobalSegMap), intent(out)   :: yGSMap  ! segmented decomposition
                                                   ! for y
! !REVISION HISTORY:
! 13Apr01 - J.W. Larson <larson@mcs.anl.gov> - API specification.
! 25Apr01 - J.W. Larson <larson@mcs.anl.gov> - initial code.
! 27Apr01 - J.W. Larson <larson@mcs.anl.gov> - Bug fix--intent of
!           argument sMat changed from (IN) to (INOUT)
! 27Apr01 - R.L. Jacob <jacob@mcs.anl.gov> - bug fix-- add use 
!           statement for SortPermute
! 01May01 - R.L. Jacob <jacob@mcs.anl.gov> - make comp_id an
!           input argument
! 07May02 - J.W. Larson <larson@mcs.anl.gov> - Changed interface to
!           make it consistent with SparseMatrixToXGlobalSegMap_().
!EOP ___________________________________________________________________

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

! SparseMatrix attributes:
  integer :: lsize
! GlobalSegMap input attributes:
  integer :: gsize, ngseg
  integer, dimension(:), pointer :: starts, lengths
! Temporary array for identifying each matrix element column and 
! process ID destination
  integer, dimension(:), allocatable :: gRow, element_pe_locs 
! Index to identify the gRow attribute in sMat:
  integer :: igRow
! Matrix element sorting keys list:
  type(List) :: sort_keys
! Loop index and error flag:
  integer :: i, ierr

       ! Determine he local number of matrix elements lsize

  lsize = SparseMatrix_lsize(sMat)

       ! The value of gsize is taken from the number of columns in sMat:

  gsize = SparseMatrix_nRows(sMat)

       ! Sort SparseMatrix entries by global column index grow, then
       ! global row index.

       ! Create Sort keys list

  call List_init(sort_keys,'grow:gcol')

       ! Sort and permute the entries of sMat into lexicographic order
       ! by global column, then global row.

  call SparseMatrix_SortPermute(sMat, sort_keys)

       ! Clean up sort keys list

  call List_clean(sort_keys)

       ! Allocate storage space for matrix element column indices and
       ! process ID destinations

  allocate(gRow(lsize), stat=ierr)

  if(ierr /= 0) then
     call die(myname_,'allocate(gRow...',ierr)
  endif

       ! Extract global column information and place in array gRow

  igRow = SparseMatrix_indexIA(sMat,'grow', dieWith=myname_)

  do i=1, lsize
     gRow(i) = sMat%data%iAttr(igRow,i)
  end do

       ! Scan sorted entries of gRow to count segments (ngseg), and
       ! their starting indices and lengths (returned in the arrays
       ! starts(:) and lengths(:), respectively)

  call ComputeSegments_(gRow, lsize, ngseg, starts, lengths)

       ! Now we have sufficient data to call the GlobalSegMap
       ! initialization using distributed data:

  call GlobalSegMap_init(yGSMap, starts, lengths, root, comm, &
                         comp_id, gsize=gsize)
  
       ! clean up temporary arrays gRow(:), starts(:) and lengths(:),
       ! (the latter two were allocated in the call to the routine
       ! ComputeSegments_())

  deallocate(gRow, starts, lengths, stat=ierr)

  if(ierr /= 0) then
     call die(myname_,'deallocate(gRow...',ierr)
  endif

 end subroutine SparseMatrixToYGlobalSegMap_

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!BOP -------------------------------------------------------------------
!
! !IROUTINE: CreateSegments_ - Generate segment information.
!
! !DESCRIPTION:  This routine examines an input {\tt INTEGER} list of
! numbers {\tt indices} (of length {\tt num\_indices}), determines the
! number of segments of consecutive numbers (or runs) {\tt nsegs}.  The
! starting indices for each run, and their lengths are returned in the 
! {\tt INTEGER} arrays {\tt starts(:)} and {\tt lengths(:)}, respectively.
!
! !INTERFACE:

 subroutine ComputeSegments_(indices, num_indices, nsegs, starts, lengths)

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

      implicit none
!
! !INPUT PARAMETERS: 
!

      integer, dimension(:), intent(in)  :: indices
      integer,               intent(in)  :: num_indices
!
! !OUTPUT PARAMETERS:
!
      integer,               intent(out) :: nsegs
      integer, dimension(:), pointer     :: starts
      integer, dimension(:), pointer     :: lengths


! !REVISION HISTORY:
! 19Apr01 - J.W. Larson <larson@mcs.anl.gov> - API specification.
! 25Apr01 - J.W. Larson <larson@mcs.anl.gov> - Initial code.
! 27Apr01 - J.W. Larson <larson@mcs.anl.gov> - Bug fix--error in
!           computation of segment starts/lengths.
! 27Nov01 - E.T. Ong <eong@mcs.anl.gov> - Bug fix--initialize
!           nsegs=0 in case num_indices=0.
!EOP ___________________________________________________________________

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

  integer :: i, ierr

       ! First pass:  count the segments

  nsegs = 0

  do i=1,num_indices

     if(i == 1) then ! bootstrap segment counting process

	nsegs = 1

     else

	if(indices(i) > indices(i-1) + 1) then ! new segment
	   nsegs = nsegs + 1
	endif

     endif ! if(i==1)

  end do ! do i=1, num_indices

       ! Allocate storage space for starts(:) and lengths(:)

  allocate(starts(nsegs), lengths(nsegs), stat=ierr)

  if(ierr /= 0) then
     call die(myname_,'allocate(starts...',ierr)
  endif

       ! Second pass:  compute segment start/length info

  do i=1,num_indices

     select case(i)
     case(1)  ! bootstrap segment counting process
	nsegs = 1
	starts(nsegs) = indices(i)
! rml patch
        lengths(nsegs) = 1
     case default

	if(i == num_indices) then ! last point
	   if(indices(i) > indices(i-1) + 1) then ! new segment with 1 pt.
       ! first, close the books on the penultimate segment:
	      lengths(nsegs) = indices(i-1) - starts(nsegs) + 1 
	      nsegs = nsegs + 1
	      starts(nsegs) = indices(i)
	      lengths(nsegs) = 1  ! (just one point)
	   else
	      lengths(nsegs) = indices(i) - starts(nsegs) + 1
	   endif
	else
	   if(indices(i) > indices(i-1) + 1) then ! new segment
	      lengths(nsegs) = indices(i-1) - starts(nsegs) + 1
	      nsegs = nsegs + 1
	      starts(nsegs) = indices(i)
	   endif
	endif

     end select ! select case(i)

  end do ! do i=1, num_indices

 end subroutine ComputeSegments_

 end  module m_SparseMatrixToMaps