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