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

!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!    Math and Computer Science Division, Argonne National Laboratory   !
!-----------------------------------------------------------------------
! CVS m_SparseMatrixDecomp.F90,v 1.15 2008-05-12 01:59:32 jacob Exp
! CVS MCT_2_8_0 
!BOP -------------------------------------------------------------------
!
! !MODULE: m_SparseMatrixDecomp -- Parallel sparse matrix decomposition.
!
! !DESCRIPTION:
! The {\tt SparseMatrix} datatype provides sparse matrix storage for 
! the parallel matrix-vector multiplication ${\bf y} = {\bf M} {\bf x}$.
! This module provides services to create decompositions for the 
! {\tt SparseMatrix}.  The matrix decompositions available are row 
! and column decompositions.  They are generated by invoking the 
! appropriate routine in this module, and passing the corresponding 
! {\em vector} decomposition.  For a row (column) decomposition, one 
! invokes the routine {\tt ByRow()} ({\tt ByColumn()}), passing the 
! domain decomposition for the vector {\bf y} ({\bf x}).
!
! !INTERFACE:

 module m_SparseMatrixDecomp

      private   ! except

! !PUBLIC MEMBER FUNCTIONS:
!
      public :: ByColumn
      public :: ByRow


    interface ByColumn ; module procedure &
         ByColumnGSMap_
    end interface

    interface ByRow ; module procedure &
         ByRowGSMap_
    end interface

! !REVISION HISTORY:
! 13Apr01 - J.W. Larson <larson@mcs.anl.gov> - initial prototype
!           and API specifications.
! 03Aug01 - E. Ong <eong@mcs.anl.gov> - in ByRowGSMap and ByColumnGSMap,
!           call GlobalSegMap_init on non-root processes with actual 
!           shaped arguments to satisfy Fortran 90 standard. See
!           comments in ByRowGSMap/ByColumnGSMap.
!EOP ___________________________________________________________________

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

 contains

!-------------------------------------------------------------------------
!     Math + Computer Science Division / Argonne National Laboratory     !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE:  ByColumnGSMap_ - Generate Row-based GlobalSegMap for SparseMatrix
! 
! !INTERFACE:

 subroutine ByColumnGSMap_(xGSMap, sMat, sMGSMap, root, comm)
!
! !USES:
!
   use m_die,  only: MP_perr_die,die

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

   use m_AttrVect, only: AttrVect
   use m_AttrVect, only: AttrVect_init => init
   use m_AttrVect, only: AttrVect_zero => zero
   use m_AttrVect, only: AttrVect_lsize => lsize
   use m_AttrVect, only: AttrVect_indexIA => indexIA
   use m_AttrVect, only: AttrVect_copy => copy
   use m_AttrVect, only: AttrVect_clean => clean
   
   use m_AttrVectComms, only: AttrVect_scatter => scatter
   use m_AttrVectComms, only: AttrVect_gather => gather

   use m_GlobalMap, only : GlobalMap
   use m_GlobalMap, only : GlobalMap_init => init
   use m_GlobalMap, only : GlobalMap_clean => clean

   use m_GlobalSegMap, only: GlobalSegMap
   use m_GlobalSegMap, only: GlobalSegMap_init => init
   use m_GlobalSegMap, only: GlobalSegMap_peLocs => peLocs
   use m_GlobalSegMap, only: GlobalSegMap_comp_id => comp_id

   use m_SparseMatrix, only: SparseMatrix
   use m_SparseMatrix, only: SparseMatrix_lsize => lsize
   use m_SparseMatrix, only: SparseMatrix_SortPermute => SortPermute

   implicit none

! !INPUT PARAMETERS: 
!
   type(GlobalSegMap), intent(in)    :: xGSMap
   integer,            intent(in)    :: root
   integer,            intent(in)    :: comm

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

! !OUTPUT PARAMETERS:
!
   type(GlobalSegMap), intent(out) :: sMGSMap

! !DESCRIPTION: This routine is invoked from all processes on the 
! communicator {\tt comm} to create from an input {\tt SparseMatrix}
! {\tt sMat} (valid only on the {\tt root} process) and an input 
! {\bf x}-vector decomposition described by the {\tt GlobalSegMap} 
! argument {\tt xGSMap} (valid at least on the {\tt root}) to create 
! an output {\tt GlobalSegMap} decomposition of the matrix elements 
! {\tt sMGSMap}, which is valid on all processes on the communicator.
! This matrix {\tt GlobalSegMap} describes the corresponding column 
! decomposition of {\tt sMat}.
!
! {\bf N.B.}:  The argument {\tt sMat} is returned sorted in lexicographic 
! order by column and row.
!
! !REVISION HISTORY: 
!
! 13Apr01 - J.W. Larson <larson@mcs.anl.gov> - initial API spec.
! 26Apr01 - R.L. Jacob <jacob@mcs.anl.gov> - add use statements for
!           GlobalSegMap_init and GSMap_peLocs.
!           Add gsize argument required to GSMap_peLocs.
!           Add underscore to ComputeSegments call so it matches
!           the subroutine decleration.
!           change attribute on starts,lengths, and pe_locs to
!           pointer to match GSMap_init.
!           add use m_die statement
! 26Apr01 - J.W. Larson <larson@mcs.anl.gov> - fixed major logic bug
!           that had all processes executing some operations that 
!           should only occur on the root.
! 09Jul03 - E.T. Ong <eong@mcs.anl.gov> - call pe_locs in parallel. 
!           reduce the serial sort from gcol:grow to just gcol.
!EOP
!-------------------------------------------------------------------------

  character(len=*),parameter :: myname_=myname//'ByColumnGSMap_'
! Process ID number
  integer :: myID, mySIZE
! Attributes for the output GlobalSegMap
  integer :: gsize, comp_id, ngseg
! Temporary array for identifying each matrix element column and 
! process ID destination
  type(AttrVect) :: gcol
  type(AttrVect) :: dist_gcol
  type(AttrVect) :: element_pe_locs 
  type(AttrVect) :: dist_element_pe_locs
! Index variables for the AttrVects
  integer :: dist_gsize
  integer :: gcol_index
  integer :: element_pe_locs_index
! Temporary array for initializing GlobalMap Decomposition
  integer,dimension(:), allocatable :: counts
! GlobalMap for setting up decomposition to call pe_locs
  type(GlobalMap) :: dist_GMap
! Temporary arrays for matrix GlobalSegMap attributes
  integer, dimension(:), pointer :: starts, lengths, pe_locs
! List storage for sorting keys
  type(List) :: sort_keys
! Error flag
  integer :: ierr
! Loop index
  integer :: i

       ! Determine process id number myID

  call MPI_COMM_RANK(comm, myID, ierr)
  if(ierr /= 0) then
     call MP_perr_die(myname_,'call MPI_COMM_RANK(...',ierr)
  endif

      ! Determine the number of processors in communicator

  call MPI_COMM_SIZE(comm, mySIZE, ierr)
  if(ierr /= 0) then
     call MP_perr_die(myname_,'call MPI_COMM_SIZE(...',ierr)
  endif

       ! Allocate space for GlobalMap length information

  allocate(counts(0:mySIZE-1),stat=ierr)
  if(ierr/=0) call die(myname_,"allocate(counts)",ierr)

       ! First step:  a lot of prep work on the root only:

  if(myID == root) then

       ! Sort the matrix entries in sMat by column.  
       ! First, create the key list...

     call List_init(sort_keys,'gcol')

       ! Now perform the sort/permute...

     call SparseMatrix_SortPermute(sMat, sort_keys)

     call List_clean(sort_keys)

       ! The global size of matrix GlobalSegMap is the number nonzero
       ! elements in sMat.  

     gsize = SparseMatrix_lsize(sMat)

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

     call AttrVect_init(aV=gcol, iList="gcol", lsize=gsize)

       ! Extract global column information and place in array gCol

     call AttrVect_copy(aVin=sMat%data, aVout=gcol, iList="gcol")

       ! Setup GlobalMap decomposition lengths:

     do i=0,mySIZE-1
        counts(i) = gsize/mySIZE
     enddo
     counts(mySIZE-1) = counts(mySIZE-1) + mod(gsize,mySIZE)

  endif

       ! Initialize GlobalMap so that we can scatter the global row
       ! information. The GlobalMap will inherit the component ID
       ! from xGSMap

  comp_id = GlobalSegMap_comp_id(xGSMap)

  call GlobalMap_init(GMap=dist_GMap, comp_id=comp_id, lns=counts, &
                      root=root, comm=comm)

  call AttrVect_scatter(iV=gcol, oV=dist_gcol, GMap=dist_GMap, &
                        root=root, comm=comm)

       ! Similarly, we want to scatter the element_pe_locs using the 
       ! same decomposition
  
  dist_gsize = AttrVect_lsize(dist_gcol)

  call AttrVect_init(aV=dist_element_pe_locs, iList="element_pe_locs", &
                     lsize=dist_gsize)
  call AttrVect_zero(dist_element_pe_locs)
     
       ! Compute process ID destination for each matrix element,
       ! and store in the AttrVect element_pe_locs

  gcol_index = AttrVect_indexIA(dist_gcol,"gcol", dieWith=myname_)
  element_pe_locs_index = AttrVect_indexIA(dist_element_pe_locs, &
                            "element_pe_locs", dieWith=myname_)

  call GlobalSegMap_peLocs(xGSMap, dist_gsize,      &
          dist_gcol%iAttr(gcol_index,1:dist_gsize), &
          dist_element_pe_locs%iAttr(element_pe_locs_index,1:dist_gsize))

  call AttrVect_gather(iV=dist_element_pe_locs, oV=element_pe_locs, &
                       GMap=dist_GMap, root=root, comm=comm)

       ! Back to the root operations

  if(myID == root) then

       ! Sanity check: Is the globalsize of sMat the same as the 
       ! gathered size of element_pe_locs?
                  
     if(gsize /= AttrVect_lsize(element_pe_locs)) then
        call die(myname_,"gsize /= AttrVect_lsize(element_pe_locs) &
                 & on root process")
     endif
       
       ! Using the entries of gCol and element_pe_locs, build the
       ! output GlobalSegMap attribute arrays starts(:), lengths(:),
       ! and pe_locs(:)

     gcol_index = AttrVect_indexIA(gcol,"gcol", dieWith=myname_)
     element_pe_locs_index = AttrVect_indexIA(element_pe_locs, &
                                "element_pe_locs", dieWith=myname_)

     call ComputeSegments_(element_pe_locs%iAttr(element_pe_locs_index, &
                                                              1:gsize), &
                           gcol%iAttr(gcol_index,1:gsize), &
                           gsize, ngseg, starts, lengths, pe_locs)
       ! Clean up on the root

     call AttrVect_clean(gcol)
     call AttrVect_clean(element_pe_locs)

  endif ! if(myID == root)

       ! Non-root processes call GlobalSegMap_init with root_start, 
       ! root_length, and root_pe_loc, although these arguments are 
       ! not used in the subroutine. Since these correspond to dummy 
       ! shaped array arguments in initr_, the Fortran 90 standard 
       ! dictates that the actual arguments must contain complete shape 
       ! information. Therefore, these array arguments must be 
       ! allocated on all processes.

  if(myID /= root) then
     allocate(starts(0),lengths(0),pe_locs(0),stat=ierr)
     if(ierr /= 0) then
        call die(myname_,'non-root allocate(starts...',ierr)
     endif
  endif

       ! Using this local data on the root, create the SparseMatrix 
       ! GlobalSegMap sMGSMap (which will be valid on all processes
       ! on the communicator:

  call GlobalSegMap_init(sMGSMap, ngseg, starts, lengths, pe_locs, &
                         root, comm, comp_id, gsize)

       ! Clean up 

  call GlobalMap_clean(dist_GMap)
  call AttrVect_clean(dist_gcol)
  call AttrVect_clean(dist_element_pe_locs)

  deallocate(starts, lengths, pe_locs, counts, stat=ierr)
  if(ierr /= 0) then
     call die(myname_,'deallocate(starts...',ierr)
  endif


 end subroutine ByColumnGSMap_

!-------------------------------------------------------------------------
!     Math + Computer Science Division / Argonne National Laboratory     !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE:  ByRowGSMap_ - Generate Row-based GlobalSegMap for SparseMatrix
! 
! !INTERFACE:

 subroutine ByRowGSMap_(yGSMap, sMat, sMGSMap, root, comm)
!
! !USES:
!

   use m_die,  only: MP_perr_die,die

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

   use m_AttrVect, only: AttrVect
   use m_AttrVect, only: AttrVect_init => init
   use m_AttrVect, only: AttrVect_lsize => lsize
   use m_AttrVect, only: AttrVect_indexIA => indexIA
   use m_AttrVect, only: AttrVect_copy => copy
   use m_AttrVect, only: AttrVect_clean => clean
   use m_AttrVect, only: AttrVect_zero => zero
   
   use m_AttrVectComms, only: AttrVect_scatter => scatter
   use m_AttrVectComms, only: AttrVect_gather => gather

   use m_GlobalMap, only : GlobalMap
   use m_GlobalMap, only : GlobalMap_init => init
   use m_GlobalMap, only : GlobalMap_clean => clean

   use m_GlobalSegMap, only: GlobalSegMap
   use m_GlobalSegMap, only: GlobalSegMap_init => init
   use m_GlobalSegMap, only: GlobalSegMap_peLocs => peLocs
   use m_GlobalSegMap, only: GlobalSegMap_comp_id => comp_id

   use m_SparseMatrix, only: SparseMatrix
   use m_SparseMatrix, only: SparseMatrix_lsize => lsize
   use m_SparseMatrix, only: SparseMatrix_SortPermute => SortPermute

   implicit none

! !INPUT PARAMETERS: 
!
   type(GlobalSegMap), intent(in)    :: yGSMap
   integer,            intent(in)    :: root
   integer,            intent(in)    :: comm

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

! !OUTPUT PARAMETERS:
!
   type(GlobalSegMap), intent(out) :: sMGSMap

! !DESCRIPTION: This routine is invoked from all processes on the 
! communicator {\tt comm} to create from an input {\tt SparseMatrix}
! {\tt sMat} (valid only on the {\tt root} process) and an input 
! {\bf y}-vector decomposition described by the {\tt GlobalSegMap} 
! argument {\tt yGSMap} (valid at least on the {\tt root}) to create 
! an output {\tt GlobalSegMap} decomposition of the matrix elements 
! {\tt sMGSMap}, which is valid on all processes on the communicator.
! This matrix {\tt GlobalSegMap} describes the corresponding row 
! decomposition of {\tt sMat}.
!
! {\bf N.B.}:  The argument {\tt sMat} is returned sorted in lexicographic 
! order by row and column.
!
! !REVISION HISTORY: 
!
! 13Apr01 - J.W. Larson <larson@mcs.anl.gov> - initial API spec.
! 26Apr01 - R.L. Jacob <jacob@mcs.anl.gov> - add use statements for
!           GlobalSegMap_init and GSMap_peLocs.
!           Add gsize argument required to GSMap_peLocs.
!           Add underscore to ComputeSegments call so it matches
!           the subroutine decleration.
!           change attribute on starts,lengths, and pe_locs to
!            pointer to match GSMap_init.
! 26Apr01 - J.W. Larson <larson@mcs.anl.gov> - fixed major logic bug
!           that had all processes executing some operations that 
!           should only occur on the root.
! 09Jun03 - E.T. Ong <eong@mcs.anl.gov> - call peLocs in parallel.
!           reduce the serial sort from grow:gcol to just grow.
!EOP
!-------------------------------------------------------------------------

  character(len=*),parameter :: myname_=myname//'ByRowGSMap_'
! Process ID number and communicator size
  integer :: myID, mySIZE
! Attributes for the output GlobalSegMap
  integer :: gsize, comp_id, ngseg
! Temporary array for identifying each matrix element row and 
! process ID destination
  type(AttrVect) :: grow
  type(AttrVect) :: dist_grow
  type(AttrVect) :: element_pe_locs 
  type(AttrVect) :: dist_element_pe_locs
! Index variables for AttrVects
  integer :: dist_gsize
  integer :: grow_index
  integer :: element_pe_locs_index
! Temporary array for initializing GlobalMap Decomposition
  integer,dimension(:), allocatable :: counts
! GlobalMap for setting up decomposition to call pe_locs
  type(GlobalMap) :: dist_GMap
! Temporary arrays for matrix GlobalSegMap attributes
  integer, dimension(:), pointer :: starts, lengths, pe_locs
! List storage for sorting keys
  type(List) :: sort_keys
! Error flag
  integer :: ierr
! Loop index
  integer :: i

       ! Determine process id number myID

  call MPI_COMM_RANK(comm, myID, ierr)
  if(ierr /= 0) then
     call MP_perr_die(myname_,'call MPI_COMM_RANK(...',ierr)
  endif

       ! Determine the number of processors in communicator

  call MPI_COMM_SIZE(comm, mySIZE, ierr)
  if(ierr /= 0) then
     call MP_perr_die(myname_,'call MPI_COMM_SIZE(...',ierr)
  endif

       ! Allocate space for GlobalMap length information

  allocate(counts(0:mySIZE-1),stat=ierr)
  if(ierr/=0) call die(myname_,"allocate(counts)",ierr)

       ! First step:  a lot of prep work on the root only:

  if(myID == root) then

       ! Sort the matrix entries in sMat by row. 
       ! First, create the key list...

     call List_init(sort_keys,'grow')

       ! Now perform the sort/permute...

     call SparseMatrix_SortPermute(sMat, sort_keys)

     call List_clean(sort_keys)

       ! The global size of matrix GlobalSegMap is the number of rows.  

     gsize = SparseMatrix_lsize(sMat)

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

     call AttrVect_init(aV=grow, iList="grow", lsize=gsize)

       ! Extract global row information and place in AttrVect grow

     call AttrVect_copy(aVin=sMat%data, aVout=grow, iList="grow")

       ! Setup GlobalMap decomposition lengths: 
       ! Give any extra points to the last process

     do i=0,mySIZE-1
        counts(i) = gsize/mySIZE
     enddo
     counts(mySIZE-1) = counts(mySIZE-1) + mod(gsize,mySIZE)

  endif
  
       ! Initialize GlobalMap and scatter the global row information.
       ! The GlobalMap will inherit the component ID from yGSMap

  comp_id = GlobalSegMap_comp_id(yGSMap)
     
  call GlobalMap_init(GMap=dist_GMap, comp_id=comp_id, lns=counts, &
                      root=root, comm=comm)

  call AttrVect_scatter(iV=grow, oV=dist_grow, GMap=dist_GMap, &
                        root=root, comm=comm)

       ! Similarly, we want to scatter the element_pe_locs using the 
       ! same decomposition
  
  dist_gsize = AttrVect_lsize(dist_grow)

  call AttrVect_init(aV=dist_element_pe_locs, iList="element_pe_locs", &
                     lsize=dist_gsize)
  call AttrVect_zero(dist_element_pe_locs)
     
       ! Compute process ID destination for each matrix element,
       ! and store in the AttrVect element_pe_locs

  grow_index = AttrVect_indexIA(dist_grow,"grow", dieWith=myname_)
  element_pe_locs_index = AttrVect_indexIA(dist_element_pe_locs, &
                            "element_pe_locs", dieWith=myname_)

  call GlobalSegMap_peLocs(yGSMap, dist_gsize,      &
          dist_grow%iAttr(grow_index,1:dist_gsize), &
          dist_element_pe_locs%iAttr(element_pe_locs_index,1:dist_gsize))

       ! Gather element_pe_locs on root so that we can call compute_segments

  call AttrVect_gather(iV=dist_element_pe_locs, oV=element_pe_locs, &
                       GMap=dist_GMap, root=root, comm=comm)

       ! Back to the root operations

  if(myID == root) then

       ! Sanity check: Is the globalsize of sMat the same as the 
       ! gathered size of element_pe_locs?
                  
     if(gsize /= AttrVect_lsize(element_pe_locs)) then
        call die(myname_,"gsize /= AttrVect_lsize(element_pe_locs) &
                 & on root process")
     endif

       ! Using the entries of grow and element_pe_locs, build the
       ! output GlobalSegMap attribute arrays starts(:), lengths(:),
       ! and pe_locs(:)

     grow_index = AttrVect_indexIA(grow,"grow", dieWith=myname_)
     element_pe_locs_index = AttrVect_indexIA(element_pe_locs, &
                                "element_pe_locs", dieWith=myname_)

     call ComputeSegments_(element_pe_locs%iAttr(element_pe_locs_index, &
                                                              1:gsize), &
                           grow%iAttr(grow_index,1:gsize), &
                           gsize, ngseg, starts, lengths, pe_locs)

       ! Clean up on the root

     call AttrVect_clean(grow)
     call AttrVect_clean(element_pe_locs)

  endif ! if(myID == root)

       ! Non-root processes call GlobalSegMap_init with root_start, 
       ! root_length, and root_pe_loc, although these arguments are 
       ! not used in the subroutine. Since these correspond to dummy 
       ! shaped array arguments in initr_, the Fortran 90 standard 
       ! dictates that the actual arguments must contain complete shape 
       ! information. Therefore, these array arguments must be 
       ! allocated on all processes.

  if(myID /= root) then
     allocate(starts(0),lengths(0),pe_locs(0),stat=ierr)
     if(ierr /= 0) then
	call die(myname_,'non-root allocate(starts...',ierr)
     endif
  endif

       ! Using this local data on the root, create the SparseMatrix 
       ! GlobalSegMap sMGSMap (which will be valid on all processes
       ! on the communicator. The GlobalSegMap will inherit the 
       ! component ID from yGSMap

  call GlobalSegMap_init(sMGSMap, ngseg, starts, lengths, pe_locs, &
                         root, comm, comp_id, gsize)

       ! Clean up:

  call GlobalMap_clean(dist_GMap)
  call AttrVect_clean(dist_grow)
  call AttrVect_clean(dist_element_pe_locs)

  deallocate(starts, lengths, pe_locs, counts, stat=ierr)
  if(ierr /= 0) then
     call die(myname_,'deallocate(starts...',ierr)
  endif


 end subroutine ByRowGSMap_

!-------------------------------------------------------------------------
!     Math + Computer Science Division / Argonne National Laboratory     !
!-------------------------------------------------------------------------
!BOP
!
! !IROUTINE:  ComputeSegments_ - Create segments from list data.
! 
! !INTERFACE:

 subroutine ComputeSegments_(element_pe_locs, elements, num_elements, &
                             nsegs, seg_starts, seg_lengths, seg_pe_locs)
!
! !USES:
!

   use m_die,  only: die

   implicit none

! !INPUT PARAMETERS: 
!
   integer, dimension(:), intent(in)  :: element_pe_locs
   integer, dimension(:), intent(in)  :: elements
   integer,               intent(in)  :: num_elements

! !OUTPUT PARAMETERS:
!
   integer,               intent(out) :: nsegs
   integer, dimension(:), pointer     :: seg_starts
   integer, dimension(:), pointer     :: seg_lengths
   integer, dimension(:), pointer     :: seg_pe_locs

! !DESCRIPTION: This routine examins an input list of {\tt num\_elements} 
! process ID locations stored in the array {\tt element\_pe\_locs}, counts
! the number of contiguous segments {\tt nsegs}, and returns the segment 
! start index, length, and process ID location in the arrays {\tt seg\_starts(:)},
! {\tt seg\_lengths(:)}, and {\tt seg\_pe\_locs(:)}, respectively.
!
! {\bf N.B.}:  The argument {\tt sMat} is returned sorted in lexicographic 
! order by row and column.
!
! !REVISION HISTORY: 
!
! 18Apr01 - J.W. Larson <larson@mcs.anl.gov> - initial version.
! 28Aug01 - M.J. Zavislak <zavislak@mcs.anl.gov>
!               Changed first sanity check to get size(element_pe_locs)
!               instead of size(elements)
!EOP
!-------------------------------------------------------------------------
  character(len=*),parameter :: myname_=myname//'ComputeSegments_'

  integer :: i, ierr, iseg

       ! Input argument sanity checks:

  if(size(element_pe_locs) < num_elements) then
     call die(myname_,'input argument array element_pe_locs too small', &
	  num_elements-size(element_pe_locs))
  endif

  if(size(elements) < num_elements) then
     call die(myname_,'input argument array elements too small', &
	  num_elements-size(elements))
  endif

       ! First pass:  how many segments?

  do i=1,num_elements

     if(i == 1) then ! bootstrap segment count

	nsegs = 1

     else ! usual point/segment processing

       ! New segment?  If so, increment nsegs.

	if((elements(i) > elements(i-1) + 1) .or. &
	     (element_pe_locs(i) /= element_pe_locs(i-1))) then ! new segment
	   nsegs = nsegs + 1
	endif

     endif ! if(i == 1) block

  end do  ! do i=1,num_elements

  allocate(seg_starts(nsegs), seg_lengths(nsegs), seg_pe_locs(nsegs), &
           stat=ierr)

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

       ! Second pass:  fill in segment data.

       ! NOTE: Structure of this loop was changed from a for loop
       ! to avoid a faulty vectorization on the SUPER-UX compiler

  i=1
  ASSIGN_LOOP: do

     if(i == 1) then  ! bootstrap first segment info.

        iseg = 1
        seg_starts(iseg) = 1
        seg_lengths(iseg) = 1
        seg_pe_locs(iseg) = element_pe_locs(iseg)

     else ! do usual point/segment processing

	! New segment?  This happens if 1) elements(i) > elements(i-1) + 1, or
	! 2) element_pe_locs(i) /= element_pe_locs(i-1).
     
	if((elements(i) > elements(i-1) + 1) .or. &
	     (element_pe_locs(i) /= element_pe_locs(i-1))) then ! new segment

	   ! Initialize new segment
	   iseg = iseg + 1
	   seg_starts(iseg) = i
	   seg_lengths(iseg) = 1
	   seg_pe_locs(iseg) = element_pe_locs(i)

	else

	   ! Increment current segment length
	   seg_lengths(iseg) = seg_lengths(iseg) + 1

	endif ! If new segment block

     endif ! if(i == 1) block

     ! Prepare index i for the next loop around; 
     if(i>=num_elements) EXIT
     i = i + 1

  end do ASSIGN_LOOP 

  if(iseg /= nsegs) then
     call die(myname_,'segment number difference',iseg-nsegs)
  endif

 end subroutine ComputeSegments_

 end module m_SparseMatrixDecomp