Converts BCSR index into list indices (similar to work matrices)
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(dbcsr_type), | intent(inout) | :: | matrix |
matrix for which to make canonical index |
||
logical, | intent(in) | :: | thread_redist |
make a thread subdistribution |
SUBROUTINE dbcsr_make_index_list(matrix, thread_redist)
!! Converts BCSR index into list indices (similar to work matrices)
TYPE(dbcsr_type), INTENT(INOUT) :: matrix
!! matrix for which to make canonical index
LOGICAL, INTENT(IN) :: thread_redist
!! make a thread subdistribution
CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_make_index_list'
LOGICAL, PARAMETER :: dbg = .FALSE.
INTEGER :: blk, error_handle, ithread, my_cnt, &
nblks, nrows, nthreads
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: blki
INTEGER, DIMENSION(0) :: zero_len_array
INTEGER, DIMENSION(:), POINTER :: global_cols, local_rows, td, thr_c
! ---------------------------------------------------------------------------
CALL timeset(routineN, error_handle)
IF (.NOT. ASSOCIATED(matrix%row_p)) &
DBCSR_ABORT("The row index must be initialized.")
IF (matrix%list_indexing) &
DBCSR_ABORT("List index already exists?")
IF (matrix%bcsc) &
DBCSR_ABORT("Not support for BCSC yet.")
IF (matrix%nblks /= SIZE(matrix%col_i)) &
DBCSR_ABORT("Block count mismatch.")
IF (matrix%nblks /= SIZE(matrix%blk_p)) &
DBCSR_ABORT("Block count mismatch")
!
IF (matrix%local_indexing) THEN
IF (SIZE(matrix%row_p) - 1 /= matrix%nblkrows_local) &
DBCSR_ABORT("Local row index incorrectly sized.")
ELSE
IF (SIZE(matrix%row_p) - 1 /= matrix%nblkrows_total) &
DBCSR_ABORT("Global row index incorrectly sized")
END IF
!
matrix%list_indexing = .TRUE.
!
IF (matrix%local_indexing) THEN
nrows = matrix%nblkrows_local
ELSE
nrows = matrix%nblkrows_total
END IF
!
nblks = matrix%nblks
ALLOCATE (blki(3, nblks))
CALL dbcsr_expand_row_index_2d(matrix%row_p, blki, nrows, 1)
IF (matrix%local_indexing) THEN
global_cols => array_data(matrix%global_cols)
! If local indexing is enabled, then the rows but not the
! columns are already localized
IF (dbg) THEN
WRITE (*, *) routineN//" Making local columns"
WRITE (*, '(10(1X,i7))') global_cols
WRITE (*, *) 'local'
WRITE (*, '(10(1X,i7))') array_data(matrix%local_cols)
END IF
DO blk = 1, nblks
blki(2, blk) = global_cols(matrix%col_i(blk))
blki(3, blk) = matrix%blk_p(blk)
END DO
ELSE
IF (dbg) WRITE (*, *) routineN//" Not making local columns"
DO blk = 1, nblks
blki(2, blk) = matrix%col_i(blk)
blki(3, blk) = matrix%blk_p(blk)
END DO
END IF
CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_row_p)
CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_col_i)
CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_blk_p)
nthreads = 0
!$ nthreads = OMP_GET_MAX_THREADS()
!
! Reshuffle according to threads
IF (nthreads .GT. 0 .AND. thread_redist) THEN
td => array_data(dbcsr_distribution_thread_dist(dbcsr_distribution(matrix)))
IF (matrix%local_indexing) THEN
local_rows => array_data(matrix%local_rows)
END IF
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (my_cnt, ithread, blk) &
!$OMP SHARED (td, blki, nthreads, thr_c, nblks, matrix,local_rows)
!
ithread = 0
!$ ithread = OMP_GET_THREAD_NUM()
!$OMP MASTER
!$ nthreads = OMP_GET_NUM_THREADS()
CALL dbcsr_addto_index_array(matrix, dbcsr_slot_thr_c, &
reservation=nthreads + 1, extra=3*nblks)
thr_c => matrix%thr_c
CALL dbcsr_addto_index_array(matrix, dbcsr_slot_coo_l, &
reservation=3*nblks)
!$OMP END MASTER
!$OMP BARRIER
my_cnt = 0
IF (matrix%local_indexing) THEN
my_cnt = COUNT(td(local_rows(blki(1, :))) .EQ. ithread)
ELSE
my_cnt = COUNT(td(blki(1, :)) .EQ. ithread)
END IF
!DO blk = 1, nblks
! IF (td(blki(1, blk)) .EQ. ithread) my_cnt = my_cnt+1
!ENDDO
thr_c(ithread + 1) = my_cnt
!$OMP BARRIER
!$OMP MASTER
CALL dbcsr_build_row_index_inplace(thr_c, nthreads)
!$OMP END MASTER
!$OMP BARRIER
my_cnt = (thr_c(ithread + 1) + 1)*3 - 2
IF (matrix%local_indexing) THEN
DO blk = 1, nblks
IF (td(local_rows(blki(1, blk))) .EQ. ithread) THEN
matrix%coo_l(my_cnt:my_cnt + 2) = blki(1:3, blk)
my_cnt = my_cnt + 3
END IF
END DO
ELSE
DO blk = 1, nblks
IF (td(blki(1, blk)) .EQ. ithread) THEN
matrix%coo_l(my_cnt:my_cnt + 2) = blki(1:3, blk)
my_cnt = my_cnt + 3
END IF
END DO
END IF
!$OMP END PARALLEL
ELSE
! Small price to pay for avoiding infinite recursions.
DO blk = 2, nblks
IF (blki(1, blk) .EQ. blki(1, blk - 1) .AND. blki(2, blk) .EQ. blki(2, blk - 1)) THEN
! Weird assertion to give some idea of the two blocks.
IF (-blki(1, blk) /= blki(2, blk)) &
CALL dbcsr_abort(__LOCATION__, &
"Should not have duplicate matrix blocks. (-row, col) is duplicated.")
END IF
END DO
!
IF (nblks > 0) THEN
CALL dbcsr_addto_index_array(matrix, dbcsr_slot_coo_l, &
DATA=RESHAPE(blki, (/3*nblks/)))
ELSE
CALL dbcsr_addto_index_array(matrix, dbcsr_slot_coo_l, &
DATA=zero_len_array)
END IF
END IF
DEALLOCATE (blki)
!
CALL timestop(error_handle)
END SUBROUTINE dbcsr_make_index_list