Convert a block-row distributed DBCSR matrix to a CSR matrix The DBCSR matrix must have a block structure consistent with the CSR matrix.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(dbcsr_type), | intent(in) | :: | brd_mat |
block-row distributed DBCSR matrix |
||
type(csr_type), | intent(inout) | :: | csr_mat |
CSR matrix |
SUBROUTINE convert_brd_to_csr(brd_mat, csr_mat)
!! Convert a block-row distributed DBCSR matrix to a CSR matrix
!! The DBCSR matrix must have a block structure consistent with the CSR matrix.
TYPE(dbcsr_type), INTENT(IN) :: brd_mat
!! block-row distributed DBCSR matrix
TYPE(csr_type), INTENT(INOUT) :: csr_mat
!! CSR matrix
CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_brd_to_csr'
INTEGER :: blk, blkcol, blkrow, col_blk_offset, col_blk_size, csr_ind, data_type, dbcsr_ind, &
el_sum, handle, ind_blk_data, k, local_row_ind, m, n, nblkrows_total, node_row_offset, &
prev_blkrow, prev_row_blk_size, row_blk_offset, row_blk_size
INTEGER, ALLOCATABLE, DIMENSION(:) :: nfullcol_blkrow
INTEGER, DIMENSION(:), POINTER :: colind, csr_index, dbcsr_index, nzerow, &
rowptr
LOGICAL :: new_ind, tr
TYPE(dbcsr_data_obj) :: block
TYPE(dbcsr_iterator) :: iter
CALL timeset(routineN, handle)
local_row_ind = 0
dbcsr_ind = 0
node_row_offset = 0
NULLIFY (rowptr, colind, dbcsr_index, csr_index)
dbcsr_index => csr_mat%dbcsr_mapping%csr_to_brd_ind
csr_index => csr_mat%dbcsr_mapping%brd_to_csr_ind
! CSR indices are not recalculated if indices are already defined
new_ind = .NOT. (csr_mat%has_indices)
IF (.NOT. csr_mat%has_mapping) &
DBCSR_ABORT("DBCSR mapping of CSR matrix must be defined")
! Calculate mapping between CSR matrix and DBCSR matrix if not yet defined
!IF (.NOT. csr_mat%has_mapping ) THEN
! CALL csr_get_dbcsr_mapping (brd_mat, dbcsr_index, csr_index, nze)
!ENDIF
CALL dbcsr_get_info(brd_mat, nblkrows_total=nblkrows_total)
ALLOCATE (nfullcol_blkrow(nblkrows_total))
! iteration over blocks without touching data,
! in order to get number of non-zero full columns in each block row
CALL dbcsr_iterator_start(iter, brd_mat, read_only=.TRUE.)
blkrow = 0
nfullcol_blkrow = 0 ! number of non-zero full columns in each block row
data_type = dbcsr_get_data_type(brd_mat)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, blkrow, blkcol, blk, col_size=col_blk_size, &
row_offset=row_blk_offset)
nfullcol_blkrow(blkrow) = nfullcol_blkrow(blkrow) + col_blk_size
IF (blk .EQ. 1) THEN
node_row_offset = row_blk_offset
END IF
END DO
CALL dbcsr_iterator_stop(iter)
! Copy data from BRD matrix to CSR matrix and calculate CSR indices
prev_blkrow = 0
prev_row_blk_size = 0
el_sum = 0 ! number of elements above current block row
colind => csr_mat%colind_local
rowptr => csr_mat%rowptr_local
nzerow => csr_mat%nzerow_local
CALL dbcsr_data_init(block)
CALL dbcsr_data_new(block, data_type)
CALL dbcsr_iterator_start(iter, brd_mat, read_only=.TRUE.)
IF (new_ind) rowptr(:) = 0 ! initialize to 0 in order to check which rows are 0 at a later time
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, blkrow, blkcol, block, tr, &
col_size=col_blk_size, row_size=row_blk_size, row_offset=row_blk_offset, &
col_offset=col_blk_offset)
IF (tr) &
DBCSR_ABORT("DBCSR block data must not be transposed")
IF (blkrow > prev_blkrow) THEN ! new block row
local_row_ind = row_blk_offset - node_row_offset ! actually: local row index - 1
IF (prev_blkrow .GT. 0) THEN
el_sum = el_sum + nfullcol_blkrow(prev_blkrow)*prev_row_blk_size
END IF
dbcsr_ind = el_sum
END IF
DO n = 1, col_blk_size !nr of columns
DO m = 1, row_blk_size !nr of rows
dbcsr_ind = dbcsr_ind + 1
csr_ind = csr_index(dbcsr_ind) ! get CSR index for current element
IF (csr_ind .GT. 0) THEN ! is non-zero element if csr_ind > 0
IF (new_ind) THEN
! Calculate CSR column index
colind(csr_ind) = col_blk_offset + n - 1
! Calculate CSR row pointer
! (not accounting for zero elements inside non-zero blocks)
IF (rowptr(local_row_ind + m) .LE. 0) rowptr(local_row_ind + m) = &
rowptr(local_row_ind + m) + el_sum + 1 + nfullcol_blkrow(blkrow)*(m - 1)
END IF
ind_blk_data = (m + row_blk_size*(n - 1)) ! index of data inside DBCSR blocks
SELECT CASE (csr_mat%nzval_local%data_type)
CASE (dbcsr_type_real_4)
csr_mat%nzval_local%r_sp(csr_ind) = block%d%r_sp(ind_blk_data)
CASE (dbcsr_type_real_8)
csr_mat%nzval_local%r_dp(csr_ind) = block%d%r_dp(ind_blk_data)
CASE (dbcsr_type_complex_4)
csr_mat%nzval_local%c_sp(csr_ind) = block%d%c_sp(ind_blk_data)
CASE (dbcsr_type_complex_8)
csr_mat%nzval_local%c_dp(csr_ind) = block%d%c_dp(ind_blk_data)
END SELECT
ELSE ! is zero element if ind = -1
! CSR row pointer has to be corrected if element is zero
! (subtract 1 from all subsequent row pointers)
IF (new_ind) rowptr(local_row_ind + m + 1:) = rowptr(local_row_ind + m + 1:) - 1
END IF
END DO
END DO
prev_blkrow = blkrow
prev_row_blk_size = row_blk_size
END DO
IF (new_ind) THEN
! Repeat previous row pointer for row pointers to zero rows
IF (csr_mat%nrows_local .GT. 0) rowptr(1) = 1
DO k = 1, csr_mat%nrows_local
IF (rowptr(k) .LE. 0) rowptr(k) = rowptr(k - 1)
END DO
rowptr(csr_mat%nrows_local + 1) = csr_mat%nze_local + 1
END IF
CALL csr_create_nzerow(csr_mat, nzerow)
CALL dbcsr_iterator_stop(iter)
CALL dbcsr_data_clear_pointer(block)
CALL dbcsr_data_release(block)
IF (new_ind) csr_mat%has_indices = .TRUE.
CALL timestop(handle)
END SUBROUTINE convert_brd_to_csr