dbcsr_reserve_blocks Subroutine

   SUBROUTINE dbcsr_reserve_blocks(matrix, rows, columns, blk_pointers)
      !! Inserts block reservations into a matrix, avoiding the work matrix.
      !!
      !! Data
      !! No data can be specified; instead, space is reserved and zeroed. To
      !! add data, call dbcsr_put_block afterwards.
      !!
      !! Reserving existing blocks
      !! Duplicates are not added, but allocations may be greater than
      !! the minimum necessary.
      !!
      !! blk_pointers
      !! When blk_pointers is passed, the newly added blocks use these pointers.
      !! No data is cleared in this case

      TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix
         !! Matrix into which blocks should be added.
      INTEGER, DIMENSION(:), INTENT(IN)                  :: rows, columns
         !! Rows of the blocks to add
         !! Columns of the blocks to add
      INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: blk_pointers
         !! block pointers to use for new blocks

      CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_reserve_blocks'

      INTEGER :: blk, blk_p, data_size_new, data_size_old, handle, nblkrows, nblks_actual_added, &
                 nblks_added, nblks_new, nblks_old, new_data_sizes, nze
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: add_blkp, add_cols, add_rows, &
                                                            added_sizes, new_blk_p, new_col_i, &
                                                            new_row_i, old_row_i
      INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: added_blk_info

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

      CALL timeset(routineN, handle)
      IF (SIZE(rows) /= SIZE(columns)) &
         DBCSR_ABORT("Size of rows and columns array must match.")
      IF (PRESENT(blk_pointers)) THEN
         IF (SIZE(rows) /= SIZE(blk_pointers)) &
            DBCSR_ABORT("Size of rows and block pointecs arrays must match.")
         data_size_old = 0
      ELSE
         ! Get current data size
         data_size_old = dbcsr_data_get_size_referenced(matrix%data_area)
      END IF
      ! Ensures that the rows and columns are sorted.
      nblks_added = SIZE(rows)
      ALLOCATE (add_rows(nblks_added))
      add_rows(:) = rows(:)
      ALLOCATE (add_cols(nblks_added))
      add_cols(:) = columns(:)
      IF (PRESENT(blk_pointers)) THEN
         ALLOCATE (add_blkp(nblks_added))
         add_blkp(:) = blk_pointers(:)
         CALL dbcsr_sort_indices(nblks_added, add_rows, add_cols, &
                                 blk_p=add_blkp)
      ELSE
         CALL dbcsr_sort_indices(nblks_added, add_rows, add_cols)
      END IF
      nblks_old = dbcsr_get_num_blocks(matrix)
      nblkrows = dbcsr_nblkrows_total(matrix)
      IF (SIZE(rows) .GT. 0 .AND. nblkrows .LE. 0) &
         DBCSR_ABORT("Can not add blocks to matrix with no rows.")
      ! Adjust the index.
      ! Get the old row indices
      ALLOCATE (old_row_i(nblks_old))
      CALL dbcsr_expand_row_index(matrix%row_p, old_row_i, &
                                  nblkrows, nblks_old)
      ! Calculate new block pointers. Possibly high estimates.
      new_data_sizes = 0
      blk_p = data_size_old + 1 ! New blocks start at the end of the old
      ALLOCATE (added_blk_info(3, nblks_added))
      ALLOCATE (added_sizes(nblks_added))
      DO blk = 1, nblks_added
         IF (PRESENT(blk_pointers)) THEN
            blk_p = add_blkp(blk)
         END IF
         added_blk_info(1:3, blk) = (/add_rows(blk), add_cols(blk), blk_p/)
         nze = dbcsr_blk_row_size(matrix, add_rows(blk)) &
               *dbcsr_blk_column_size(matrix, add_cols(blk))
         added_sizes(blk) = nze
         blk_p = blk_p + nze
      END DO
      DEALLOCATE (add_rows)
      DEALLOCATE (add_cols)
      IF (PRESENT(blk_pointers)) DEALLOCATE (add_blkp)
      !
      nblks_new = nblks_old + nblks_added ! Possibly high estimate
      ALLOCATE (new_row_i(nblks_new))
      ALLOCATE (new_col_i(nblks_new))
      ALLOCATE (new_blk_p(nblks_new))
      ! Merge the two indices
      IF (PRESENT(blk_pointers)) THEN
         CALL merge_index_arrays(new_row_i, new_col_i, new_blk_p, nblks_new, &
                                 old_row_i, matrix%col_i, matrix%blk_p, nblks_old, &
                                 added_blk_info, nblks_added, added_nblks=nblks_actual_added)
         data_size_new = 0
      ELSE
         CALL merge_index_arrays(new_row_i, new_col_i, new_blk_p, nblks_new, &
                                 old_row_i, matrix%col_i, matrix%blk_p, nblks_old, &
                                 added_blk_info, nblks_added, added_nblks=nblks_actual_added, &
                                 added_sizes=added_sizes, added_size_offset=data_size_old + 1, &
                                 added_size=data_size_new)
      END IF
      nblks_new = nblks_actual_added + nblks_old
      ! Free some memory
      DEALLOCATE (added_blk_info)
      DEALLOCATE (added_sizes)
      DEALLOCATE (old_row_i)
      ! We can skip this if no block was actually added.
      IF (nblks_actual_added .GT. 0) THEN
         ! Write the new index
         matrix%nblks = nblks_new
         matrix%nze = matrix%nze + data_size_new
         matrix%index(dbcsr_slot_nblks) = matrix%nblks
         matrix%index(dbcsr_slot_nze) = matrix%index(dbcsr_slot_nze)
         CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_col_i)
         CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_blk_p)
         CALL dbcsr_addto_index_array(matrix, dbcsr_slot_col_i, &
                                      new_col_i(1:nblks_new), &
                                      extra=nblks_new)
         CALL dbcsr_addto_index_array(matrix, dbcsr_slot_blk_p, &
                                      new_blk_p(1:nblks_new))
         CALL dbcsr_make_dbcsr_index(matrix%row_p, new_row_i(1:nblks_new), &
                                     nblkrows, nblks_new)
         IF (.NOT. PRESENT(blk_pointers)) THEN
            ! Resize data area to fit the new blocks.
            CALL dbcsr_data_ensure_size(matrix%data_area, &
                                        data_size=matrix%nze)
            ! Zero the new data blocks.
            CALL dbcsr_data_clear(matrix%data_area, &
                                  lb=data_size_old + 1, ub=matrix%nze)
         END IF
      END IF
      CALL timestop(handle)
   END SUBROUTINE dbcsr_reserve_blocks