Crop and copies a matrix.
Type | Intent | Optional | 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 |
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