dbcsr_crop_matrix Subroutine

public subroutine dbcsr_crop_matrix(matrix_b, matrix_a, full_row_bounds, full_column_bounds, shallow_data)

Crop and copies a matrix.

Arguments

Type IntentOptional Attributes Name
type(dbcsr_type), intent(inout) :: matrix_b

target DBCSR matrix

type(dbcsr_type), intent(in) :: matrix_a

source DBCSR matrix

integer, intent(in), optional, DIMENSION(2) :: full_row_bounds

rows to extract (array of size 2 holding the lower and upper inclusive bounds) columns to extract (array of size 2 holding the lower and upper inclusive bounds)

integer, intent(in), optional, DIMENSION(2) :: full_column_bounds

rows to extract (array of size 2 holding the lower and upper inclusive bounds) columns to extract (array of size 2 holding the lower and upper inclusive bounds)

logical, intent(in), optional :: shallow_data

Source Code

   SUBROUTINE dbcsr_crop_matrix(matrix_b, matrix_a, &
                                full_row_bounds, full_column_bounds, &
                                shallow_data)
      !! Crop and copies a matrix.

      TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_b
         !! target DBCSR matrix
      TYPE(dbcsr_type), INTENT(IN)                       :: matrix_a
         !! source DBCSR matrix
      INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL        :: full_row_bounds, full_column_bounds
         !! rows to extract (array of size 2 holding the lower and upper inclusive bounds)
         !! columns to extract (array of size 2 holding the lower and upper inclusive bounds)
      LOGICAL, INTENT(IN), OPTIONAL                      :: shallow_data

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

      INTEGER                                            :: col, f_col_f, f_row_f, handle, l_col_l, &
                                                            l_row_l, row
      INTEGER, DIMENSION(2)                              :: block_col_bounds, block_row_bounds
      LOGICAL                                            :: part_col, part_f_col, part_f_row, &
                                                            part_l_col, part_l_row, part_row, &
                                                            shallow, tr
      TYPE(dbcsr_data_obj)                               :: data_block
      TYPE(dbcsr_iterator)                               :: iter

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

      CALL timeset(routineN, handle)
      part_l_col = .FALSE.
      part_f_col = .FALSE.
      part_l_row = .FALSE.
      part_f_row = .FALSE.
      IF (PRESENT(shallow_data)) THEN
         shallow = shallow_data
      ELSE
         shallow = .FALSE.
      END IF
      block_row_bounds = 0
      block_col_bounds = 0
      part_col = .FALSE.
      part_row = .FALSE.
      !
      ! If row bounds are present, they must be converted to block
      ! addressing.
      IF (PRESENT(full_row_bounds)) THEN
         IF (SIZE(full_row_bounds) /= 2) &
            DBCSR_ABORT("Size of bounds specifier must be 2")
         IF (full_row_bounds(1) < 0) &
            DBCSR_ABORT("Invalid first row bound.")
         IF (full_row_bounds(2) > dbcsr_nfullrows_total(matrix_a)) &
            DBCSR_ABORT("Invalid last row bound.")
         IF (full_row_bounds(1) .EQ. 0) THEN
            block_row_bounds(1) = 1
         ELSE
            CALL find_block_of_element(full_row_bounds(1), block_row_bounds(1), &
                                       dbcsr_nblkrows_total(matrix_a), &
                                       dbcsr_row_block_offsets(matrix_a), &
                                       hint=0)
            part_f_row = array_get(dbcsr_row_block_offsets(matrix_a), block_row_bounds(1)) &
                         .NE. full_row_bounds(1)
         END IF
         f_row_f = -7
         IF (part_f_row) THEN
            ! Block offset of last cleared row
            f_row_f = full_row_bounds(1) - &
                      array_get(dbcsr_row_block_offsets(matrix_a), block_row_bounds(1))
         END IF
         IF (full_row_bounds(2) .EQ. 0) THEN
            block_row_bounds(2) = dbcsr_nblkrows_total(matrix_a)
         ELSE
            CALL find_block_of_element(full_row_bounds(2), block_row_bounds(2), &
                                       dbcsr_nblkrows_total(matrix_a), &
                                       dbcsr_row_block_offsets(matrix_a), &
                                       hint=0)
            part_l_row = array_get(dbcsr_row_block_offsets(matrix_a), block_row_bounds(2) + 1) - 1 &
                         .NE. full_row_bounds(2)
         END IF
         ! Block offset of first cleared row
         l_row_l = -7
         IF (part_l_row) THEN
            l_row_l = 2 + full_row_bounds(2) - &
                      array_get(dbcsr_row_block_offsets(matrix_a), block_row_bounds(2))
         END IF
         part_row = part_f_row .OR. part_l_row
      END IF
      !
      ! If column bounds are present, they must be converted to block
      ! addressing.
      IF (PRESENT(full_column_bounds)) THEN
         IF (SIZE(full_column_bounds) /= 2) &
            DBCSR_ABORT("Size of bounds specifier must be 2")
         IF (full_column_bounds(1) < 0) &
            DBCSR_ABORT("Invalid first column bound.")
         IF (full_column_bounds(2) > dbcsr_nfullcols_total(matrix_a)) &
            DBCSR_ABORT("Invalid last column bound.")
         IF (full_column_bounds(1) .EQ. 0) THEN
            block_col_bounds(1) = 1
         ELSE
            CALL find_block_of_element(full_column_bounds(1), block_col_bounds(1), &
                                       dbcsr_nblkcols_total(matrix_a), &
                                       dbcsr_col_block_offsets(matrix_a), &
                                       hint=0)
            part_f_col = array_get(dbcsr_col_block_offsets(matrix_a), block_col_bounds(1)) &
                         .NE. full_column_bounds(1)
         END IF
         f_col_f = -7
         IF (part_f_col) THEN
            ! Block offset of last cleared column
            f_col_f = full_column_bounds(1) - &
                      array_get(dbcsr_col_block_offsets(matrix_a), block_col_bounds(1))
         END IF
         IF (full_column_bounds(2) .EQ. 0) THEN
            block_col_bounds(2) = dbcsr_nblkcols_total(matrix_a)
         ELSE
            CALL find_block_of_element(full_column_bounds(2), block_col_bounds(2), &
                                       dbcsr_nblkcols_total(matrix_a), &
                                       dbcsr_col_block_offsets(matrix_a), &
                                       hint=0)
            part_l_col = array_get(dbcsr_col_block_offsets(matrix_a), block_col_bounds(2) + 1) - 1 &
                         .NE. full_column_bounds(2)
         END IF
         l_col_l = -7
         IF (part_l_col) THEN
            ! Block offset of first cleared column
            l_col_l = 2 + full_column_bounds(2) - &
                      array_get(dbcsr_col_block_offsets(matrix_a), block_col_bounds(2))
         END IF
         part_col = part_f_col .OR. part_l_col
      END IF
      !
      ! First copy the blocks then perform the intra-block zeroing.
      CALL dbcsr_copy_submatrix(matrix_b, matrix_a, &
                                block_row_bounds=block_row_bounds, &
                                block_column_bounds=block_col_bounds, &
                                shallow_data=shallow)
      IF (part_row .OR. part_col) THEN
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP          PRIVATE (data_block, iter, row, col, tr) &
!$OMP          SHARED (matrix_b,&
!$OMP                  part_row, part_f_row, part_l_row, f_row_f, l_row_l, &
!$OMP                  part_col, part_f_col, part_l_col, f_col_f, l_col_l,&
!$OMP                  block_row_bounds, block_col_bounds)
         CALL dbcsr_data_init(data_block)
         CALL dbcsr_data_new(data_block, dbcsr_type_1d_to_2d(dbcsr_get_data_type(matrix_b)))
         CALL dbcsr_iterator_start(iter, matrix_b, &
                                   dynamic=.TRUE., dynamic_byrows=.TRUE.)
         DO WHILE (dbcsr_iterator_blocks_left(iter))
            CALL dbcsr_iterator_next_block(iter, row, col, data_block, tr)
            IF (part_row) THEN
               IF (row .LT. block_row_bounds(1)) CYCLE
               IF (row .GT. block_row_bounds(2)) CYCLE
            END IF
            IF (part_col) THEN
               IF (col .LT. block_col_bounds(1)) CYCLE
               IF (col .GT. block_col_bounds(2)) CYCLE
            END IF
            IF (part_row) THEN
               IF (part_f_row .AND. row .EQ. block_row_bounds(1)) THEN
                  CALL dbcsr_data_clear(data_block, ub=f_row_f, tr=tr)
               END IF
               IF (part_l_row .AND. row .EQ. block_row_bounds(2)) THEN
                  CALL dbcsr_data_clear(data_block, lb=l_row_l, tr=tr)
               END IF
            END IF
            IF (part_col) THEN
               IF (part_f_col .AND. col .EQ. block_col_bounds(1)) THEN
                  CALL dbcsr_data_clear(data_block, ub2=f_col_f, tr=tr)
               END IF
               IF (part_l_col .AND. col .EQ. block_col_bounds(2)) THEN
                  CALL dbcsr_data_clear(data_block, lb2=l_col_l, tr=tr)
               END IF
            END IF
         END DO
         CALL dbcsr_iterator_stop(iter)
         CALL dbcsr_data_clear_pointer(data_block)
         CALL dbcsr_data_release(data_block)
         CALL dbcsr_finalize(matrix_b)
!$OMP END PARALLEL
      END IF
      !
      CALL timestop(handle)
   END SUBROUTINE dbcsr_crop_matrix