dbcsr_sort_indices Subroutine

   SUBROUTINE dbcsr_sort_indices(n, row_i, col_i, blk_p, blk_d)
      !! Sorts the rows & columns of a work matrix
      !!
      !! Description
      !! Sorts the row and column indices so that the rows monotonically
      !! increase and the columns monotonically increase within each row.
      !! Passing the blk_p array rearranges the block pointers accordingly.
      !! This must be done if they are pointing to valid data, otherwise
      !! they become invalid.

      INTEGER, INTENT(IN)                                :: n
         !! number of blocks (elements) to sort
      INTEGER, DIMENSION(1:), INTENT(INOUT)              :: row_i, col_i
         !! row indices
         !! column indices
      INTEGER, DIMENSION(1:), INTENT(INOUT), OPTIONAL    :: blk_p, blk_d
         !! block pointers
         !! data storage

      CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_sort_indices', &
                                     routineP = moduleN//':'//routineN
      INTEGER(KIND=int_8), PARAMETER                     :: lmask8 = 4294967295_int_8

      INTEGER                                            :: error_handle, i
      INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:)     :: sort_keys
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: buf, buf_d

!   ---------------------------------------------------------------------------

      IF (n .LE. 0) RETURN
      IF (SIZE(row_i) .EQ. 0) RETURN

      CALL timeset(routineN, error_handle)

      IF (SIZE(row_i) < n) DBCSR_ABORT('row_i too small')
      IF (SIZE(col_i) < n) DBCSR_ABORT('col_i too small')
      IF (PRESENT(blk_p)) THEN
         IF (SIZE(blk_p) < n) DBCSR_ABORT('blk_p too small')
         ALLOCATE (buf(n))
         buf(1:n) = blk_p(1:n)
      END IF
      IF (PRESENT(blk_d)) THEN
         ALLOCATE (buf_d(n))
         buf_d(1:n) = blk_d(1:n)
      END IF
      ! Create an ordering for both rows and columns. If the blk_p must
      ! be rearranged, then the col_i array will be used as a
      ! permutation vector.
      ALLOCATE (sort_keys(n))
      sort_keys(:) = IOR(ISHFT(INT(row_i(1:n), int_8), 32), INT(col_i(1:n), int_8))
      IF (PRESENT(blk_p)) col_i(1:n) = (/(i, i=1, n)/)
      ! Now do a nice quicksort.
      CALL sort(sort_keys, n, col_i)
      ! Since blk_d is usually not present we can have two loops that
      ! are essentially the same.
      IF (PRESENT(blk_p)) THEN
         DO i = 1, n
            blk_p(i) = buf(col_i(i))
         END DO
         DEALLOCATE (buf)
      END IF
      IF (PRESENT(blk_d)) THEN
         DO i = 1, n
            blk_d(i) = buf_d(col_i(i))
         END DO
         DEALLOCATE (buf_d)
      END IF
      DO i = 1, n
         col_i(i) = INT(IAND(sort_keys(i), lmask8), int_4)
         row_i(i) = INT(ISHFT(sort_keys(i), -32), int_4)
      END DO
      DEALLOCATE (sort_keys)
      IF (debug_mod .AND. PRESENT(blk_p)) &
         WRITE (*, *) routineP//' sort, blk_p =', blk_p
      CALL timestop(error_handle)

   END SUBROUTINE dbcsr_sort_indices