merge_index_arrays Subroutine

   SUBROUTINE merge_index_arrays(new_row_i, new_col_i, new_blk_p, new_size, &
                                 old_row_i, old_col_i, old_blk_p, old_size, &
                                 add_ip, add_size, new_blk_d, old_blk_d, &
                                 added_size_offset, added_sizes, added_size, added_nblks)
      !! Merges two indices
      !!
      !! Added sizes
      !! added_size_offset and added_sizes can be optionally
      !! specified. This is meant for cases where the added blocks may
      !! be duplicates of existing blocks. In this way it is possible
      !! to recalculate new block pointers to avoid wasted space.
      !! @note Used in local multiply
      !! Assumes they are both pre-sorted

      INTEGER, INTENT(IN)                                :: new_size
         !! size of merged index
      INTEGER, DIMENSION(new_size), INTENT(OUT)          :: new_blk_p, new_col_i, new_row_i
         !! merged result
         !! merged result
         !! merged result
      INTEGER, INTENT(IN)                                :: old_size
         !! size of current index
      INTEGER, DIMENSION(old_size), INTENT(IN)           :: old_blk_p, old_col_i, old_row_i
         !! current index
         !! current index
         !! current index
      INTEGER, INTENT(IN)                                :: add_size
         !! size of index to add into the current index
      INTEGER, DIMENSION(3, add_size), INTENT(IN)        :: add_ip
         !! index to add into the current index
      INTEGER, DIMENSION(new_size), INTENT(OUT), &
         OPTIONAL                                        :: new_blk_d
      INTEGER, DIMENSION(old_size), INTENT(IN), OPTIONAL :: old_blk_d
      INTEGER, INTENT(IN), OPTIONAL                      :: added_size_offset
         !! specify base of added sizes
      INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: added_sizes
         !! specify sizes of added blocks
      INTEGER, INTENT(OUT), OPTIONAL                     :: added_size, added_nblks
         !! counts number of sizes of added blocks
         !! actual number of new elements

      INTEGER                                            :: add_blk, bp, i, merge_from_whom, &
                                                            new_blk, old_blk
      LOGICAL                                            :: multidata

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

      bp = 0
      multidata = PRESENT(old_blk_d) .AND. PRESENT(new_blk_d)
      IF (old_size + add_size .NE. new_size) &
         DBCSR_WARN("Mismatch of new and old size")
      IF (PRESENT(added_size_offset) .NEQV. PRESENT(added_sizes)) &
         DBCSR_ABORT("Must specify a set of arguments")
      IF (PRESENT(added_sizes) .NEQV. PRESENT(added_size)) &
         DBCSR_ABORT("Must specify a set of arguments")
      IF (debug_mod) THEN
         WRITE (*, *) " Old array", old_size
         DO i = 1, old_size
            WRITE (*, '(I7,2X,I7,2X,I7)') old_row_i(i), old_col_i(i), old_blk_p(i)
         END DO
         WRITE (*, *) " Add array", add_size
         DO i = 1, add_size
            WRITE (*, '(I7,2X,I7,2X,I7)') add_ip(1:3, i)
         END DO
      END IF
      IF (PRESENT(added_nblks)) added_nblks = 0
      IF (PRESENT(added_size)) THEN
         added_size = 0
         bp = added_size_offset
      END IF
      IF (add_size .GT. 0) THEN
         old_blk = 1
         add_blk = 1
         new_blk = 1
         IF (old_size .EQ. 0) THEN
            new_row_i(1:add_size) = add_ip(1, 1:add_size)
            new_col_i(1:add_size) = add_ip(2, 1:add_size)
            new_blk_p(1:add_size) = add_ip(3, 1:add_size)
            !IF (multidata) new_blk_d(1:add_size) = add_ip(4, 1:add_size)
            IF (PRESENT(added_nblks)) added_nblks = add_size
            IF (PRESENT(added_size)) added_size = SUM(added_sizes)
         ELSE
            DO WHILE (new_blk .LE. new_size)
               merge_from_whom = 0
               IF (old_blk .LE. old_size .AND. add_blk .LE. add_size) THEN
                  IF (add_ip(1, add_blk) .EQ. old_row_i(old_blk) &
                      .AND. add_ip(2, add_blk) .EQ. old_col_i(old_blk)) THEN
                     IF (debug_mod) THEN
                        WRITE (*, *) "Duplicate block! addblk", &
                           add_blk, "oldblk", old_blk
                     END IF
                  END IF
                  ! Rows come first
                  IF (add_ip(1, add_blk) .LT. old_row_i(old_blk)) THEN
                     merge_from_whom = 2
                  ELSEIF (add_ip(1, add_blk) .GT. old_row_i(old_blk)) THEN
                     merge_from_whom = 1
                  ELSE ! Same rows, so now come the columns
                     IF (add_ip(2, add_blk) .LT. old_col_i(old_blk)) THEN
                        ! Merges from the add array
                        merge_from_whom = 2
                     ELSEIF (add_ip(2, add_blk) .GT. old_col_i(old_blk)) THEN
                        ! Merges from the old array
                        merge_from_whom = 1
                     ELSE
                        ! Merge from old array and skip one in the new array
                        IF (debug_mod) THEN
                           WRITE (*, *) "Duplicate, keeping old", &
                              add_ip(1, add_blk), add_ip(2, add_blk)
                        END IF
                        merge_from_whom = 1
                        add_blk = add_blk + 1
                     END IF
                  END IF
               ELSE
                  IF (add_blk .LE. add_size) THEN
                     ! Merges from the add array
                     merge_from_whom = 2
                  ELSEIF (old_blk .LE. old_size) THEN
                     ! Merges from the old array
                     merge_from_whom = 1
                  ELSE
                     ! Hmmm, nothing to merge...
                     merge_from_whom = 0
                     !WRITE(*,*)"Ran out of data to merge"
                  END IF
               END IF
               SELECT CASE (merge_from_whom)
               CASE (2)
                  ! Merges from the add array
                  new_row_i(new_blk) = add_ip(1, add_blk)
                  new_col_i(new_blk) = add_ip(2, add_blk)
                  new_blk_p(new_blk) = add_ip(3, add_blk)
                  !IF (multidata) new_blk_d(new_blk) = add_ip(4, add_blk)
                  IF (PRESENT(added_nblks)) added_nblks = added_nblks + 1
                  IF (PRESENT(added_sizes)) THEN
                     new_blk_p(new_blk) = bp
                     bp = bp + added_sizes(add_blk)
                     added_size = added_size + added_sizes(add_blk)
                  END IF
                  add_blk = add_blk + 1
               CASE (1)
                  ! Merges from the old array
                  new_row_i(new_blk) = old_row_i(old_blk)
                  new_col_i(new_blk) = old_col_i(old_blk)
                  new_blk_p(new_blk) = old_blk_p(old_blk)
                  IF (multidata) new_blk_p(new_blk) = old_blk_d(old_blk)
                  old_blk = old_blk + 1
               CASE DEFAULT
                  !WRITE(*,*)"Nothing to merge"
               END SELECT
               new_blk = new_blk + 1
            END DO
         END IF
      ELSE
         new_row_i(1:old_size) = old_row_i(1:old_size)
         new_col_i(1:old_size) = old_col_i(1:old_size)
         new_blk_p(1:old_size) = old_blk_p(1:old_size)
         IF (multidata) new_blk_d(1:old_size) = old_blk_d(1:old_size)
      END IF
      IF (debug_mod) THEN
         WRITE (*, *) " New array"
         DO i = 1, new_size
            WRITE (*, '(4(2X,I7))') new_row_i(i), new_col_i(i), new_blk_p(i)
         END DO
      END IF
   END SUBROUTINE merge_index_arrays