1!--------------------------------------------------------------------------------------------------! 2! Copyright (C) by the DBCSR developers group - All rights reserved ! 3! This file is part of the DBCSR library. ! 4! ! 5! For information on the license, see the LICENSE file. ! 6! For further information please visit https://dbcsr.cp2k.org ! 7! SPDX-License-Identifier: GPL-2.0+ ! 8!--------------------------------------------------------------------------------------------------! 9 10MODULE dbcsr_csr_conversions 11 !! DBCSR to CSR matrix format conversion 12 USE dbcsr_block_access, ONLY: dbcsr_put_block 13 USE dbcsr_data_methods, ONLY: dbcsr_data_clear_pointer, & 14 dbcsr_data_init, & 15 dbcsr_data_new, & 16 dbcsr_data_release 17 USE dbcsr_data_types, ONLY: dbcsr_type_complex_4, & 18 dbcsr_type_complex_8, & 19 dbcsr_type_real_4, & 20 dbcsr_type_real_8, & 21 dbcsr_type_real_default 22 USE dbcsr_dist_methods, ONLY: dbcsr_distribution_col_dist, & 23 dbcsr_distribution_mp, & 24 dbcsr_distribution_new, & 25 dbcsr_distribution_release 26 USE dbcsr_iterator_operations, ONLY: dbcsr_iterator_blocks_left, & 27 dbcsr_iterator_next_block, & 28 dbcsr_iterator_start, & 29 dbcsr_iterator_stop 30 USE dbcsr_kinds, ONLY: default_string_length, & 31 dp, & 32 int_8, & 33 real_4, & 34 real_8, & 35 sp 36 USE dbcsr_methods, ONLY: & 37 dbcsr_col_block_sizes, dbcsr_distribution, dbcsr_get_data_type, dbcsr_get_num_blocks, & 38 dbcsr_get_nze, dbcsr_has_symmetry, dbcsr_name, dbcsr_nblkcols_total, dbcsr_nblkrows_total, & 39 dbcsr_nfullrows_local, dbcsr_release, dbcsr_row_block_sizes, dbcsr_valid_index 40 USE dbcsr_mp_methods, ONLY: dbcsr_mp_group, & 41 dbcsr_mp_mynode, & 42 dbcsr_mp_new, & 43 dbcsr_mp_numnodes, & 44 dbcsr_mp_release 45 USE dbcsr_mpiwrap, ONLY: mp_environ, & 46 mp_gather, & 47 mp_recv, & 48 mp_send, & 49 mp_sum 50 USE dbcsr_operations, ONLY: dbcsr_copy, & 51 dbcsr_get_info, & 52 dbcsr_set 53 USE dbcsr_transformations, ONLY: dbcsr_complete_redistribute, & 54 dbcsr_desymmetrize_deep 55 USE dbcsr_types, ONLY: dbcsr_data_obj, & 56 dbcsr_distribution_obj, & 57 dbcsr_iterator, & 58 dbcsr_mp_obj, & 59 dbcsr_type 60 USE dbcsr_work_operations, ONLY: dbcsr_create 61#include "base/dbcsr_base_uses.f90" 62 63 IMPLICIT NONE 64 PRIVATE 65 66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbcsr_csr_conversions' 67 68 LOGICAL, PARAMETER, PRIVATE :: careful_mod = .FALSE. 69 70 INTEGER, PARAMETER, PUBLIC :: csr_dbcsr_blkrow_dist = 1, csr_eqrow_ceil_dist = 2, csr_eqrow_floor_dist = 3 71 72 TYPE csr_mapping_data 73 !! Mapping data relating local CSR indices to local indices of a block-row 74 !! distributed (BRD) DBCSR matrix, and containing the block structure 75 !! of the original DBCSR matrix from which the CSR matrix was created. 76 77 PRIVATE 78 INTEGER, DIMENSION(:), POINTER :: csr_to_brd_ind => NULL(), & 79 brd_to_csr_ind => NULL() 80 !! csr_to_brd_ind(csr_ind) gives the location of a matrix element with CSR index csr_ind (location in nzval_local) inside 81 !! the data_area of the corresponding BRD matrix. If an element of the DBCSR matrix is treated as 0 in the CSR format, the 82 !! index of this value is not in csr_to_brd_ind. 83 !! same as csr_to_brd_ind but inverse mapping. If a given DBCSR index dbcsr_ind points to a zero element, then 84 !! brd_to_csr_ind(dbcsr_ind) is -1. 85 TYPE(dbcsr_type) :: brd_mat 86 !! DBCSR BRD matrix acting as an intermediate step in any conversion from and to DBCSR format. 87 88 LOGICAL :: has_dbcsr_block_data = .FALSE. 89 !! whether dbcsr_* fields are defined 90 INTEGER :: dbcsr_nblkcols_total, & 91 dbcsr_nblkrows_total, & 92 dbcsr_nblks_local 93 !! data from original DBCSR matrix (not block-row distributed), 94 !! representing the original block structure. 95 INTEGER, DIMENSION(:), POINTER :: dbcsr_row_p, dbcsr_col_i, & 96 dbcsr_row_blk_size, dbcsr_col_blk_size 97 !! data from original DBCSR matrix (not block-row distributed), 98 !! representing the original block structure. 99 ENDTYPE 100 101 TYPE csr_data_area_type 102 !! Data type of CSR matrices 103 104 REAL(KIND=real_4), DIMENSION(:), POINTER :: r_sp => Null() 105 !! real, single precision data array 106 REAL(KIND=real_8), DIMENSION(:), POINTER :: r_dp => Null() 107 !! real, double precision data array 108 COMPLEX(KIND=real_4), DIMENSION(:), POINTER :: c_sp => Null() 109 !! complex, double precision data array 110 COMPLEX(KIND=real_8), DIMENSION(:), POINTER :: c_dp => Null() 111 INTEGER :: data_type = -1 112 !! data type of CSR matrix 113 ENDTYPE 114 115 TYPE csr_type 116 !! Type for CSR matrices 117 118 INTEGER :: nrows_total, ncols_total, & 119 nze_local, nrows_local, mp_group 120 !! total number of rows 121 !! total number of columns 122 !! local number of nonzero elements 123 !! local number of rows 124 !! message-passing group ID 125 INTEGER(KIND=int_8) :: nze_total 126 !! total number of nonzero elements 127 INTEGER, DIMENSION(:), POINTER :: rowptr_local => NULL(), & 128 colind_local => NULL(), & 129 nzerow_local => NULL() 130 !! indices of elements inside nzval_local starting a row 131 !! column indices of elements inside nzval_local 132 TYPE(csr_data_area_type) :: nzval_local 133 !! values of local non-zero elements, row-wise ordering. 134 TYPE(csr_mapping_data) :: dbcsr_mapping 135 !! mapping data relating indices of nzval_local to indices of a block-row distributed DBCSR matrix 136 LOGICAL :: has_mapping = .FALSE. 137 !! whether dbcsr_mapping is defined 138 LOGICAL :: valid = .FALSE. 139 !! whether essential data (excluding dbcsr_mapping) is completely allocated 140 LOGICAL :: has_indices = .FALSE. 141 !! whether rowptr_local and colind_local are defined 142 ENDTYPE csr_type 143 144 TYPE csr_p_type 145 TYPE(csr_type), POINTER :: csr_mat 146 END TYPE csr_p_type 147 148 PUBLIC :: csr_type, csr_p_type, convert_csr_to_dbcsr, & 149 csr_create_from_dbcsr, & 150 csr_destroy, & 151 convert_dbcsr_to_csr, & 152 csr_create_new, csr_create_template, & 153 csr_print_sparsity, dbcsr_to_csr_filter, & 154 csr_write 155 156 INTERFACE csr_create 157 MODULE PROCEDURE csr_create_new, csr_create_template 158 END INTERFACE 159 160CONTAINS 161 162 SUBROUTINE csr_create_new(csr_mat, nrows_total, ncols_total, nze_total, & 163 nze_local, nrows_local, mp_group, data_type) 164 !! Create a new CSR matrix and allocate all internal data (excluding dbcsr_mapping) 165 166 TYPE(csr_type), INTENT(OUT) :: csr_mat 167 !! CSR matrix to return 168 INTEGER, INTENT(IN) :: nrows_total, ncols_total 169 !! total number of rows 170 !! total number of columns 171 INTEGER(KIND=int_8) :: nze_total 172 !! total number of non-zero elements 173 INTEGER, INTENT(IN) :: nze_local, nrows_local, mp_group 174 !! local number of non-zero elements 175 !! local number of rows 176 INTEGER, INTENT(IN), OPTIONAL :: data_type 177 !! data type of the CSR matrix (default real double prec.) 178 179 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_create_new' 180 INTEGER :: handle 181 182 CALL timeset(routineN, handle) 183 184 IF (nrows_total .LT. nrows_local) & 185 DBCSR_ABORT("local number of rows must not exceed total number of rows") 186 187 IF (nze_total .LT. nze_local) CALL dbcsr_abort(__LOCATION__, "local number of non-zero "// & 188 "elements must not exceed total number of non-zero elements") 189 190 IF (INT(nrows_total, kind=int_8)*INT(ncols_total, kind=int_8) .LT. nze_total) & 191 DBCSR_ABORT("Total number of non-zero elements must not exceed total matrix size") 192 193 IF (INT(nrows_local, kind=int_8)*INT(ncols_total, kind=int_8) .LT. nze_local) & 194 DBCSR_ABORT("Local number of non-zero elements must not exceed local matrix size") 195 196 csr_mat%ncols_total = ncols_total 197 csr_mat%nrows_total = nrows_total 198 csr_mat%nze_total = nze_total 199 csr_mat%nze_local = nze_local 200 ALLOCATE (csr_mat%colind_local(nze_local)) 201 csr_mat%nrows_local = nrows_local 202 ALLOCATE (csr_mat%rowptr_local(nrows_local + 1)) 203 ALLOCATE (csr_mat%nzerow_local(nrows_local)) 204 205 IF (PRESENT(data_type)) THEN 206 csr_mat%nzval_local%data_type = data_type 207 ELSE 208 csr_mat%nzval_local%data_type = dbcsr_type_real_default 209 ENDIF 210 211 SELECT CASE (csr_mat%nzval_local%data_type) 212 CASE (dbcsr_type_real_4) 213 ALLOCATE (csr_mat%nzval_local%r_sp(nze_local)) 214 CASE (dbcsr_type_real_8) 215 ALLOCATE (csr_mat%nzval_local%r_dp(nze_local)) 216 CASE (dbcsr_type_complex_4) 217 ALLOCATE (csr_mat%nzval_local%c_sp(nze_local)) 218 CASE (dbcsr_type_complex_8) 219 ALLOCATE (csr_mat%nzval_local%c_dp(nze_local)) 220 CASE DEFAULT 221 DBCSR_ABORT("Invalid matrix type") 222 END SELECT 223 224 csr_mat%mp_group = mp_group 225 226 csr_mat%valid = .TRUE. 227 csr_mat%has_mapping = .FALSE. 228 csr_mat%has_indices = .FALSE. 229 230 CALL timestop(handle) 231 232 END SUBROUTINE csr_create_new 233 234 SUBROUTINE csr_create_template(matrix_b, matrix_a) 235 !! Create a new CSR matrix and allocate all internal data using 236 !! an existing CSR matrix. Copies the indices but no actual matrix data. 237 238 TYPE(csr_type), INTENT(OUT) :: matrix_b 239 !! Target CSR matrix 240 TYPE(csr_type), INTENT(IN) :: matrix_a 241 !! Source CSR matrix 242 243 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_create_template' 244 245 INTEGER :: handle 246 TYPE(csr_mapping_data) :: map 247 248 CALL timeset(routineN, handle) 249 250 IF (.NOT. matrix_a%valid) & 251 DBCSR_ABORT("Source CSR matrix must be created.") 252 253 CALL csr_create_new(matrix_b, matrix_a%nrows_total, matrix_a%ncols_total, & 254 matrix_a%nze_total, matrix_a%nze_local, matrix_a%nrows_local, & 255 matrix_a%mp_group, matrix_a%nzval_local%data_type) 256 257 matrix_b%mp_group = matrix_a%mp_group 258 matrix_b%has_mapping = matrix_a%has_mapping 259 matrix_b%has_indices = matrix_a%has_indices 260 261 IF (matrix_a%has_indices) THEN 262 matrix_b%rowptr_local = matrix_a%rowptr_local 263 matrix_b%nzerow_local = matrix_a%nzerow_local 264 matrix_b%colind_local = matrix_a%colind_local 265 ENDIF 266 267 IF (matrix_a%has_mapping) THEN 268 map = matrix_a%dbcsr_mapping 269 ALLOCATE (matrix_b%dbcsr_mapping%csr_to_brd_ind(SIZE(map%csr_to_brd_ind))) 270 ALLOCATE (matrix_b%dbcsr_mapping%brd_to_csr_ind(SIZE(map%brd_to_csr_ind))) 271 matrix_b%dbcsr_mapping%csr_to_brd_ind = map%csr_to_brd_ind 272 matrix_b%dbcsr_mapping%brd_to_csr_ind = map%brd_to_csr_ind 273 matrix_b%dbcsr_mapping%has_dbcsr_block_data = map%has_dbcsr_block_data 274 IF (matrix_b%dbcsr_mapping%has_dbcsr_block_data) THEN 275 matrix_b%dbcsr_mapping%dbcsr_nblkcols_total = map%dbcsr_nblkcols_total 276 matrix_b%dbcsr_mapping%dbcsr_nblkrows_total = map%dbcsr_nblkrows_total 277 ALLOCATE (matrix_b%dbcsr_mapping%dbcsr_row_blk_size(map%dbcsr_nblkrows_total)) 278 ALLOCATE (matrix_b%dbcsr_mapping%dbcsr_col_blk_size(map%dbcsr_nblkcols_total)) 279 ALLOCATE (matrix_b%dbcsr_mapping%dbcsr_row_p(map%dbcsr_nblkrows_total + 1)) 280 ALLOCATE (matrix_b%dbcsr_mapping%dbcsr_col_i(map%dbcsr_nblks_local)) 281 matrix_b%dbcsr_mapping%dbcsr_nblks_local = map%dbcsr_nblks_local 282 matrix_b%dbcsr_mapping%dbcsr_row_p = map%dbcsr_row_p 283 matrix_b%dbcsr_mapping%dbcsr_col_i = map%dbcsr_col_i 284 matrix_b%dbcsr_mapping%dbcsr_row_blk_size = map%dbcsr_row_blk_size 285 matrix_b%dbcsr_mapping%dbcsr_col_blk_size = map%dbcsr_col_blk_size 286 ENDIF 287 288 CALL dbcsr_copy(matrix_b%dbcsr_mapping%brd_mat, map%brd_mat) 289 290 matrix_b%valid = .TRUE. 291 292 ENDIF 293 294 CALL timestop(handle) 295 END SUBROUTINE csr_create_template 296 297 SUBROUTINE csr_create_nzerow(csr_mat, nzerow) 298 !! create a vector containing the number of non-zero elements in each 299 !! row of a CSR matrix 300 301 TYPE(csr_type), INTENT(IN) :: csr_mat 302 !! CSR matrix 303 INTEGER, DIMENSION(:), INTENT(INOUT), POINTER :: nzerow 304 !! number of non-zero elements in each row 305 306 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_create_nzerow' 307 308 INTEGER :: handle, k 309 310 CALL timeset(routineN, handle) 311 312 IF (.NOT. csr_mat%valid) & 313 DBCSR_ABORT("CSR matrix must be created.") 314 315 DO k = 1, csr_mat%nrows_local 316 nzerow(k) = csr_mat%rowptr_local(k + 1) - csr_mat%rowptr_local(k) 317 ENDDO 318 319 CALL timestop(handle) 320 END SUBROUTINE csr_create_nzerow 321 322 SUBROUTINE csr_destroy(csr_mat) 323 !! destroy a CSR matrix 324 TYPE(csr_type), INTENT(INOUT) :: csr_mat 325 326 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_destroy' 327 INTEGER :: handle 328 TYPE(csr_mapping_data) :: map 329 330 CALL timeset(routineN, handle) 331 332 IF (.NOT. csr_mat%valid) & 333 DBCSR_ABORT("CSR matrix must be created before destroying it.") 334 335 IF (ASSOCIATED(csr_mat%rowptr_local)) DEALLOCATE (csr_mat%rowptr_local) 336 IF (ASSOCIATED(csr_mat%nzerow_local)) DEALLOCATE (csr_mat%nzerow_local) 337 IF (ASSOCIATED(csr_mat%colind_local)) DEALLOCATE (csr_mat%colind_local) 338 339 IF (csr_mat%has_mapping) THEN 340 map = csr_mat%dbcsr_mapping 341 IF (ASSOCIATED(map%csr_to_brd_ind)) & 342 DEALLOCATE (map%csr_to_brd_ind) 343 IF (ASSOCIATED(map%brd_to_csr_ind)) & 344 DEALLOCATE (map%brd_to_csr_ind) 345 IF (ASSOCIATED(map%dbcsr_row_blk_size)) & 346 DEALLOCATE (map%dbcsr_row_blk_size) 347 IF (ASSOCIATED(map%dbcsr_col_blk_size)) & 348 DEALLOCATE (map%dbcsr_col_blk_size) 349 IF (ASSOCIATED(map%dbcsr_row_p)) & 350 DEALLOCATE (map%dbcsr_row_p) 351 IF (ASSOCIATED(map%dbcsr_col_i)) & 352 DEALLOCATE (map%dbcsr_col_i) 353 354 CALL dbcsr_release(map%brd_mat) 355 ENDIF 356 357 IF (ASSOCIATED(csr_mat%nzval_local%r_dp)) & 358 DEALLOCATE (csr_mat%nzval_local%r_dp) 359 IF (ASSOCIATED(csr_mat%nzval_local%r_sp)) & 360 DEALLOCATE (csr_mat%nzval_local%r_sp) 361 IF (ASSOCIATED(csr_mat%nzval_local%c_sp)) & 362 DEALLOCATE (csr_mat%nzval_local%c_sp) 363 IF (ASSOCIATED(csr_mat%nzval_local%c_dp)) & 364 DEALLOCATE (csr_mat%nzval_local%c_dp) 365 366 csr_mat%has_mapping = .FALSE. 367 csr_mat%valid = .FALSE. 368 csr_mat%dbcsr_mapping%has_dbcsr_block_data = .FALSE. 369 csr_mat%has_indices = .FALSE. 370 csr_mat%nzval_local%data_type = -1 371 372 CALL timestop(handle) 373 END SUBROUTINE csr_destroy 374 375 SUBROUTINE csr_create_from_brd(brd_mat, csr_mat, csr_sparsity_brd) 376 !! Allocate the internals of a CSR matrix using data from a block-row 377 !! distributed DBCSR matrix 378 379 TYPE(dbcsr_type), INTENT(IN) :: brd_mat 380 !! block-row-distributed DBCSR matrix 381 TYPE(csr_type), INTENT(OUT) :: csr_mat 382 !! CSR matrix 383 TYPE(dbcsr_type), INTENT(IN) :: csr_sparsity_brd 384 !! BRD matrix representing sparsity pattern of CSR matrix 385 386 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_create_from_brd' 387 388 INTEGER :: data_type, handle, mp_group, & 389 nfullcols_total, nfullrows, & 390 nfullrows_total, nze_local 391 INTEGER(KIND=int_8) :: nze_total 392 INTEGER, DIMENSION(:), POINTER :: cdist, csr_index, dbcsr_index 393 TYPE(dbcsr_distribution_obj) :: dist_current 394 395 CALL timeset(routineN, handle) 396 NULLIFY (dbcsr_index, csr_index, cdist) 397 398 dist_current = dbcsr_distribution(brd_mat) 399 400 mp_group = dbcsr_mp_group(dbcsr_distribution_mp(dist_current)) 401 cdist => dbcsr_distribution_col_dist(dist_current) 402 403 IF (ANY(cdist .NE. 0)) & 404 DBCSR_ABORT("DBCSR matrix not block-row distributed.") 405 406 ! Calculate mapping between BRD and CSR indices 407 CALL csr_get_dbcsr_mapping(brd_mat, dbcsr_index, csr_index, nze_local, & 408 csr_sparsity_brd) 409 410 CALL dbcsr_get_info(brd_mat, nfullrows_total=nfullrows_total, & 411 nfullcols_total=nfullcols_total) 412 413 ! Sum up local number of non-zero elements to get total number 414 nze_total = nze_local 415 CALL mp_sum(nze_total, mp_group) 416 417 nfullrows = dbcsr_nfullrows_local(brd_mat) 418 data_type = dbcsr_get_data_type(brd_mat) 419 420 ! Allocate CSR matrix 421 CALL csr_create_new(csr_mat, nfullrows_total, nfullcols_total, nze_total, & 422 nze_local, nfullrows, mp_group, data_type) 423 424 csr_mat%dbcsr_mapping%brd_to_csr_ind => csr_index 425 csr_mat%dbcsr_mapping%csr_to_brd_ind => dbcsr_index 426 427 csr_mat%has_mapping = .TRUE. 428 csr_mat%dbcsr_mapping%brd_mat = brd_mat 429 430 CALL timestop(handle) 431 END SUBROUTINE csr_create_from_brd 432 433 SUBROUTINE csr_get_dbcsr_mapping(brd_mat, dbcsr_index, csr_index, csr_nze_local, & 434 csr_sparsity_brd) 435 !! create the mapping information between a block-row distributed DBCSR 436 !! matrix and the corresponding CSR matrix 437 438 TYPE(dbcsr_type), INTENT(IN) :: brd_mat 439 !! the block-row distributed DBCSR matrix 440 INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: dbcsr_index, csr_index 441 !! csr to dbcsr index mapping 442 !! dbcsr to csr index mapping 443 INTEGER, INTENT(OUT) :: csr_nze_local 444 !! number of local non-zero elements 445 TYPE(dbcsr_type), INTENT(IN) :: csr_sparsity_brd 446 !! sparsity of CSR matrix represented in BRD format 447 448 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_get_dbcsr_mapping' 449 450 INTEGER :: blk, blkcol, blkrow, col_blk_size, csr_ind, data_type, dbcsr_ind, el_sum, & 451 fullcol_sum_blkrow, handle, l, m, n, nblkrows_total, nze, prev_blk, prev_blkrow, & 452 prev_row_blk_size, row_blk_offset, row_blk_size 453 INTEGER, ALLOCATABLE, DIMENSION(:) :: csr_nze, nfullcol_blkrow 454 INTEGER, DIMENSION(:), POINTER :: dbcsr_index_nozeroes 455 LOGICAL :: tr 456 TYPE(dbcsr_iterator) :: iter 457 458 CALL timeset(routineN, handle) 459 460 m = 0 461 dbcsr_ind = 0 462 fullcol_sum_blkrow = 0 463 NULLIFY (dbcsr_index, csr_index) 464 465 CALL dbcsr_get_info(brd_mat, nblkrows_total=nblkrows_total) 466 nze = dbcsr_get_nze(brd_mat) 467 468 ALLOCATE (nfullcol_blkrow(nblkrows_total)) 469 ALLOCATE (dbcsr_index(nze)) 470 ALLOCATE (csr_index(nze)) 471 472 CALL dbcsr_iterator_start(iter, brd_mat, read_only=.TRUE.) 473 nfullcol_blkrow = 0 ! number of non-zero full columns in each block row 474 prev_blk = 0 475 476 DO WHILE (dbcsr_iterator_blocks_left(iter)) 477 CALL dbcsr_iterator_next_block(iter, blkrow, blkcol, blk, transposed=tr, & 478 col_size=col_blk_size) 479 480 IF (blk /= prev_blk + 1) & 481 DBCSR_ABORT("iterator is required to traverse the blocks in a row-wise fashion") 482 483 prev_blk = blk 484 485 nfullcol_blkrow(blkrow) = nfullcol_blkrow(blkrow) + col_blk_size 486 IF (tr) & 487 DBCSR_ABORT("DBCSR block data must not be transposed") 488 ENDDO 489 CALL dbcsr_iterator_stop(iter) 490 491 el_sum = 0 ! number of elements above current block row 492 493 prev_blkrow = 0 ! store number and size of previous block row 494 prev_row_blk_size = 0 495 496 CALL dbcsr_iterator_start(iter, brd_mat, read_only=.TRUE.) 497 498 DO WHILE (dbcsr_iterator_blocks_left(iter)) 499 500 CALL dbcsr_iterator_next_block(iter, blkrow, blkcol, blk, transposed=tr, & 501 row_size=row_blk_size, col_size=col_blk_size, row_offset=row_blk_offset) 502 503 IF (blkrow .GT. prev_blkrow) THEN ! new block row 504 IF (prev_blkrow .GT. 0) THEN 505 el_sum = el_sum + nfullcol_blkrow(prev_blkrow)*prev_row_blk_size 506 ENDIF 507 508 ! number of non-zero full columns on the left of current block: 509 fullcol_sum_blkrow = 0 510 511 dbcsr_ind = el_sum 512 ENDIF 513 DO n = 1, col_blk_size !nr of columns 514 DO m = 1, row_blk_size !nr of rows 515 dbcsr_ind = dbcsr_ind + 1 516 csr_ind = (m - 1)*nfullcol_blkrow(blkrow) + fullcol_sum_blkrow + n + el_sum 517 dbcsr_index(csr_ind) = dbcsr_ind 518 csr_index(dbcsr_ind) = csr_ind 519 ENDDO 520 ENDDO 521 fullcol_sum_blkrow = fullcol_sum_blkrow + col_blk_size 522 prev_blkrow = blkrow 523 prev_row_blk_size = row_blk_size 524 ENDDO 525 CALL dbcsr_iterator_stop(iter) 526 527 ! remove BRD zero elements from CSR format 528 data_type = dbcsr_get_data_type(csr_sparsity_brd) 529 ALLOCATE (csr_nze(nze)) 530 531 SELECT CASE (data_type) 532 CASE (dbcsr_type_real_4) 533 csr_nze(:) = INT(csr_sparsity_brd%data_area%d%r_sp(1:nze)) 534 CASE (dbcsr_type_real_8) 535 csr_nze(:) = INT(csr_sparsity_brd%data_area%d%r_dp(1:nze)) 536 CASE DEFAULT 537 DBCSR_ABORT("CSR sparsity matrix must have a real datatype") 538 END SELECT 539 540 IF (ANY(csr_nze .EQ. 0)) THEN 541 ALLOCATE (dbcsr_index_nozeroes(SUM(csr_nze))) 542 m = 0 ! csr index if zeroes are excluded from CSR data 543 DO l = 1, nze ! csr index if zeroes are included in CSR data 544 IF (csr_nze(dbcsr_index(l)) .EQ. 0) THEN 545 csr_index(dbcsr_index(l)) = -1 546 ELSE 547 m = m + 1 548 dbcsr_index_nozeroes(m) = dbcsr_index(l) 549 csr_index(dbcsr_index(l)) = m 550 ENDIF 551 ENDDO 552 DEALLOCATE (dbcsr_index) 553 dbcsr_index => dbcsr_index_nozeroes 554 ENDIF 555 556 IF (ANY(csr_nze .EQ. 0)) THEN 557 csr_nze_local = m 558 ELSE 559 csr_nze_local = nze 560 ENDIF 561 562 CALL timestop(handle) 563 END SUBROUTINE csr_get_dbcsr_mapping 564 565 SUBROUTINE convert_csr_to_brd(brd_mat, csr_mat) 566 !! Copies data from a CSR matrix to a block-row distributed DBCSR matrix. 567 !! The DBCSR matrix must have a block structure consistent with the CSR matrix. 568 569 TYPE(dbcsr_type), INTENT(INOUT) :: brd_mat 570 !! block-row distributed DBCSR matrix 571 TYPE(csr_type), INTENT(IN) :: csr_mat 572 !! CSR matrix 573 574 CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_csr_to_brd' 575 576 INTEGER :: data_type, handle, ind, k, nze 577 578 CALL timeset(routineN, handle) 579 580 data_type = dbcsr_get_data_type(brd_mat) 581 nze = dbcsr_get_nze(brd_mat) 582 CALL dbcsr_data_release(brd_mat%data_area) 583 CALL dbcsr_data_new(brd_mat%data_area, data_type, nze) 584 585 SELECT CASE (data_type) 586 CASE (dbcsr_type_real_4) 587 brd_mat%data_area%d%r_sp(1:nze) = 0.0_sp 588 CASE (dbcsr_type_real_8) 589 brd_mat%data_area%d%r_dp(1:nze) = 0.0_dp 590 CASE (dbcsr_type_complex_4) 591 brd_mat%data_area%d%c_sp(1:nze) = 0.0_sp 592 CASE (dbcsr_type_complex_8) 593 brd_mat%data_area%d%c_dp(1:nze) = 0.0_dp 594 END SELECT 595 596 DO k = 1, csr_mat%nze_local 597 ind = csr_mat%dbcsr_mapping%csr_to_brd_ind(k) 598 SELECT CASE (data_type) 599 CASE (dbcsr_type_real_4) 600 brd_mat%data_area%d%r_sp(ind) = csr_mat%nzval_local%r_sp(k) 601 CASE (dbcsr_type_real_8) 602 brd_mat%data_area%d%r_dp(ind) = csr_mat%nzval_local%r_dp(k) 603 CASE (dbcsr_type_complex_4) 604 brd_mat%data_area%d%c_sp(ind) = csr_mat%nzval_local%c_sp(k) 605 CASE (dbcsr_type_complex_8) 606 brd_mat%data_area%d%c_dp(ind) = csr_mat%nzval_local%c_dp(k) 607 END SELECT 608 ENDDO 609 610 CALL timestop(handle) 611 END SUBROUTINE convert_csr_to_brd 612 613 SUBROUTINE convert_brd_to_csr(brd_mat, csr_mat) 614 !! Convert a block-row distributed DBCSR matrix to a CSR matrix 615 !! The DBCSR matrix must have a block structure consistent with the CSR matrix. 616 617 TYPE(dbcsr_type), INTENT(IN) :: brd_mat 618 !! block-row distributed DBCSR matrix 619 TYPE(csr_type), INTENT(INOUT) :: csr_mat 620 !! CSR matrix 621 622 CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_brd_to_csr' 623 624 INTEGER :: blk, blkcol, blkrow, col_blk_offset, col_blk_size, csr_ind, data_type, dbcsr_ind, & 625 el_sum, handle, ind_blk_data, k, local_row_ind, m, n, nblkrows_total, node_row_offset, & 626 prev_blkrow, prev_row_blk_size, row_blk_offset, row_blk_size 627 INTEGER, ALLOCATABLE, DIMENSION(:) :: nfullcol_blkrow 628 INTEGER, DIMENSION(:), POINTER :: colind, csr_index, dbcsr_index, nzerow, & 629 rowptr 630 LOGICAL :: new_ind, tr 631 TYPE(dbcsr_data_obj) :: block 632 TYPE(dbcsr_iterator) :: iter 633 634 CALL timeset(routineN, handle) 635 local_row_ind = 0 636 dbcsr_ind = 0 637 node_row_offset = 0 638 NULLIFY (rowptr, colind, dbcsr_index, csr_index) 639 640 dbcsr_index => csr_mat%dbcsr_mapping%csr_to_brd_ind 641 csr_index => csr_mat%dbcsr_mapping%brd_to_csr_ind 642 643 ! CSR indices are not recalculated if indices are already defined 644 new_ind = .NOT. (csr_mat%has_indices) 645 646 IF (.NOT. csr_mat%has_mapping) & 647 DBCSR_ABORT("DBCSR mapping of CSR matrix must be defined") 648 649 ! Calculate mapping between CSR matrix and DBCSR matrix if not yet defined 650 !IF (.NOT. csr_mat%has_mapping ) THEN 651 ! CALL csr_get_dbcsr_mapping (brd_mat, dbcsr_index, csr_index, nze) 652 !ENDIF 653 654 CALL dbcsr_get_info(brd_mat, nblkrows_total=nblkrows_total) 655 ALLOCATE (nfullcol_blkrow(nblkrows_total)) 656 657 ! iteration over blocks without touching data, 658 ! in order to get number of non-zero full columns in each block row 659 CALL dbcsr_iterator_start(iter, brd_mat, read_only=.TRUE.) 660 blkrow = 0 661 nfullcol_blkrow = 0 ! number of non-zero full columns in each block row 662 data_type = dbcsr_get_data_type(brd_mat) 663 664 DO WHILE (dbcsr_iterator_blocks_left(iter)) 665 CALL dbcsr_iterator_next_block(iter, blkrow, blkcol, blk, col_size=col_blk_size, & 666 row_offset=row_blk_offset) 667 nfullcol_blkrow(blkrow) = nfullcol_blkrow(blkrow) + col_blk_size 668 IF (blk .EQ. 1) THEN 669 node_row_offset = row_blk_offset 670 ENDIF 671 ENDDO 672 673 CALL dbcsr_iterator_stop(iter) 674 675 ! Copy data from BRD matrix to CSR matrix and calculate CSR indices 676 prev_blkrow = 0 677 prev_row_blk_size = 0 678 el_sum = 0 ! number of elements above current block row 679 colind => csr_mat%colind_local 680 rowptr => csr_mat%rowptr_local 681 nzerow => csr_mat%nzerow_local 682 CALL dbcsr_data_init(block) 683 CALL dbcsr_data_new(block, data_type) 684 685 CALL dbcsr_iterator_start(iter, brd_mat, read_only=.TRUE.) 686 687 IF (new_ind) rowptr(:) = 0 ! initialize to 0 in order to check which rows are 0 at a later time 688 DO WHILE (dbcsr_iterator_blocks_left(iter)) 689 CALL dbcsr_iterator_next_block(iter, blkrow, blkcol, block, tr, & 690 col_size=col_blk_size, row_size=row_blk_size, row_offset=row_blk_offset, & 691 col_offset=col_blk_offset) 692 693 IF (tr) & 694 DBCSR_ABORT("DBCSR block data must not be transposed") 695 696 IF (blkrow > prev_blkrow) THEN ! new block row 697 local_row_ind = row_blk_offset - node_row_offset ! actually: local row index - 1 698 IF (prev_blkrow .GT. 0) THEN 699 el_sum = el_sum + nfullcol_blkrow(prev_blkrow)*prev_row_blk_size 700 ENDIF 701 dbcsr_ind = el_sum 702 ENDIF 703 DO n = 1, col_blk_size !nr of columns 704 DO m = 1, row_blk_size !nr of rows 705 dbcsr_ind = dbcsr_ind + 1 706 csr_ind = csr_index(dbcsr_ind) ! get CSR index for current element 707 IF (csr_ind .GT. 0) THEN ! is non-zero element if csr_ind > 0 708 IF (new_ind) THEN 709 ! Calculate CSR column index 710 colind(csr_ind) = col_blk_offset + n - 1 711 ! Calculate CSR row pointer 712 ! (not accounting for zero elements inside non-zero blocks) 713 IF (rowptr(local_row_ind + m) .LE. 0) rowptr(local_row_ind + m) = & 714 rowptr(local_row_ind + m) + el_sum + 1 + nfullcol_blkrow(blkrow)*(m - 1) 715 ENDIF 716 ind_blk_data = (m + row_blk_size*(n - 1)) ! index of data inside DBCSR blocks 717 SELECT CASE (csr_mat%nzval_local%data_type) 718 CASE (dbcsr_type_real_4) 719 csr_mat%nzval_local%r_sp(csr_ind) = block%d%r_sp(ind_blk_data) 720 CASE (dbcsr_type_real_8) 721 csr_mat%nzval_local%r_dp(csr_ind) = block%d%r_dp(ind_blk_data) 722 CASE (dbcsr_type_complex_4) 723 csr_mat%nzval_local%c_sp(csr_ind) = block%d%c_sp(ind_blk_data) 724 CASE (dbcsr_type_complex_8) 725 csr_mat%nzval_local%c_dp(csr_ind) = block%d%c_dp(ind_blk_data) 726 END SELECT 727 ELSE ! is zero element if ind = -1 728 ! CSR row pointer has to be corrected if element is zero 729 ! (subtract 1 from all subsequent row pointers) 730 IF (new_ind) rowptr(local_row_ind + m + 1:) = rowptr(local_row_ind + m + 1:) - 1 731 ENDIF 732 ENDDO 733 ENDDO 734 prev_blkrow = blkrow 735 prev_row_blk_size = row_blk_size 736 ENDDO 737 738 IF (new_ind) THEN 739 ! Repeat previous row pointer for row pointers to zero rows 740 IF (csr_mat%nrows_local .GT. 0) rowptr(1) = 1 741 DO k = 1, csr_mat%nrows_local 742 IF (rowptr(k) .LE. 0) rowptr(k) = rowptr(k - 1) 743 ENDDO 744 745 rowptr(csr_mat%nrows_local + 1) = csr_mat%nze_local + 1 746 ENDIF 747 748 CALL csr_create_nzerow(csr_mat, nzerow) 749 750 CALL dbcsr_iterator_stop(iter) 751 CALL dbcsr_data_clear_pointer(block) 752 CALL dbcsr_data_release(block) 753 754 IF (new_ind) csr_mat%has_indices = .TRUE. 755 756 CALL timestop(handle) 757 END SUBROUTINE convert_brd_to_csr 758 759 SUBROUTINE csr_create_from_dbcsr(dbcsr_mat, csr_mat, dist_format, csr_sparsity, numnodes) 760 !! create CSR matrix including dbcsr_mapping from arbitrary DBCSR matrix 761 !! in order to prepare conversion. 762 763 TYPE(dbcsr_type), INTENT(IN) :: dbcsr_mat 764 TYPE(csr_type), INTENT(OUT) :: csr_mat 765 INTEGER, INTENT(IN) :: dist_format 766 !! how to distribute CSR rows over processes: csr_dbcsr_blkrow_dist: the number of rows per process is adapted to the row 767 !! block sizes in the DBCSR format such that blocks are not split over different processes. csr_eqrow_ceil_dist: each 768 !! process holds ceiling(N/P) CSR rows. csr_eqrow_floor_dist: each process holds floor(N/P) CSR rows. 769 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: csr_sparsity 770 !! DBCSR matrix containing 0 and 1, representing CSR sparsity pattern 1: non-zero element 0: zero element (not present in 771 !! CSR format) Note: matrix must be of data_type dbcsr_type_real_4 or dbcsr_type_real_8 (integer types not supported) 772 INTEGER, INTENT(IN), OPTIONAL :: numnodes 773 !! number of nodes to use for distributing CSR matrix (optional, default is number of nodes used for DBCSR matrix) 774 775 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_create_from_dbcsr' 776 777 INTEGER :: dbcsr_numnodes, handle, nblkcols_total, & 778 nblkrows_total, nblks_local, num_p 779 LOGICAL :: equal_dist, floor_dist 780 TYPE(dbcsr_type) :: brd_mat, csr_sparsity_brd, & 781 csr_sparsity_nosym, dbcsr_mat_nosym 782 783 CALL timeset(routineN, handle) 784 785 IF (.NOT. dbcsr_valid_index(dbcsr_mat)) & 786 DBCSR_ABORT("Invalid DBCSR matrix") 787 788 SELECT CASE (dist_format) 789 CASE (csr_dbcsr_blkrow_dist) 790 equal_dist = .FALSE. 791 floor_dist = .FALSE. 792 CASE (csr_eqrow_ceil_dist) 793 equal_dist = .TRUE. 794 floor_dist = .FALSE. 795 CASE (csr_eqrow_floor_dist) 796 equal_dist = .TRUE. 797 floor_dist = .TRUE. 798 END SELECT 799 800 ! Conversion does not support matrices in symmetric format, therefore desymmetrize 801 IF (dbcsr_has_symmetry(dbcsr_mat)) THEN 802 CALL dbcsr_desymmetrize_deep(dbcsr_mat, dbcsr_mat_nosym, untransposed_data=.TRUE.) 803 ELSE 804 CALL dbcsr_copy(dbcsr_mat_nosym, dbcsr_mat) 805 ENDIF 806 807 IF (PRESENT(csr_sparsity)) THEN 808 IF (dbcsr_has_symmetry(csr_sparsity)) THEN 809 CALL dbcsr_desymmetrize_deep(csr_sparsity, csr_sparsity_nosym, & 810 untransposed_data=.TRUE.) 811 ELSE 812 CALL dbcsr_copy(csr_sparsity_nosym, csr_sparsity) 813 ENDIF 814 ELSE 815 CALL dbcsr_create(csr_sparsity_nosym, & 816 template=dbcsr_mat_nosym, & 817 name="CSR sparsity matrix", & 818 data_type=dbcsr_type_real_8) 819 CALL dbcsr_copy(csr_sparsity_nosym, dbcsr_mat_nosym) 820 CALL dbcsr_set(csr_sparsity_nosym, 1.0_dp) 821 ENDIF 822 823 IF (.NOT. dbcsr_has_same_block_structure(dbcsr_mat_nosym, csr_sparsity_nosym)) & 824 DBCSR_ABORT("csr_sparsity and dbcsr_mat have different sparsity pattern") 825 826 dbcsr_numnodes = dbcsr_mp_numnodes(dbcsr_distribution_mp(dbcsr_distribution(dbcsr_mat))) 827 IF (PRESENT(numnodes)) THEN 828 IF (numnodes .GT. dbcsr_numnodes) & 829 CALL dbcsr_abort(__LOCATION__, "Number of nodes used for CSR matrix "// & 830 "must not exceed total number of nodes") 831 832 num_p = numnodes 833 ELSE 834 num_p = dbcsr_numnodes 835 ENDIF 836 837 CALL dbcsr_create_brd(dbcsr_mat_nosym, brd_mat, equal_dist, floor_dist, & 838 num_p) 839 CALL dbcsr_create_brd(csr_sparsity_nosym, csr_sparsity_brd, equal_dist, floor_dist, & 840 num_p) 841 842 ! Create CSR matrix from BRD matrix 843 CALL csr_create_from_brd(brd_mat, csr_mat, csr_sparsity_brd) 844 845 ! Store DBCSR block data inside CSR matrix 846 ! (otherwise, this data is lost when converting from DBCSR to CSR) 847 nblks_local = dbcsr_get_num_blocks(dbcsr_mat_nosym) 848 nblkrows_total = dbcsr_nblkrows_total(dbcsr_mat_nosym) 849 nblkcols_total = dbcsr_nblkcols_total(dbcsr_mat_nosym) 850 851 csr_mat%dbcsr_mapping%dbcsr_nblkcols_total = nblkcols_total 852 csr_mat%dbcsr_mapping%dbcsr_nblkrows_total = nblkrows_total 853 csr_mat%dbcsr_mapping%dbcsr_nblks_local = nblks_local 854 ALLOCATE (csr_mat%dbcsr_mapping%dbcsr_row_p(nblkrows_total + 1)) 855 csr_mat%dbcsr_mapping%dbcsr_row_p = dbcsr_mat_nosym%row_p 856 ALLOCATE (csr_mat%dbcsr_mapping%dbcsr_col_i(nblks_local)) 857 csr_mat%dbcsr_mapping%dbcsr_col_i = dbcsr_mat_nosym%col_i 858 859 ALLOCATE (csr_mat%dbcsr_mapping%dbcsr_row_blk_size(nblkrows_total)) 860 ALLOCATE (csr_mat%dbcsr_mapping%dbcsr_col_blk_size(nblkcols_total)) 861 862 csr_mat%dbcsr_mapping%dbcsr_row_blk_size = dbcsr_row_block_sizes(dbcsr_mat_nosym) 863 csr_mat%dbcsr_mapping%dbcsr_col_blk_size = dbcsr_col_block_sizes(dbcsr_mat_nosym) 864 865 csr_mat%dbcsr_mapping%has_dbcsr_block_data = .TRUE. 866 867 CALL dbcsr_release(dbcsr_mat_nosym) 868 CALL dbcsr_release(csr_sparsity_nosym) 869 CALL dbcsr_release(csr_sparsity_brd) 870 871 CALL timestop(handle) 872 END SUBROUTINE csr_create_from_dbcsr 873 874 FUNCTION dbcsr_has_same_block_structure(matrix_a, matrix_b) RESULT(is_equal) 875 !! Helper function to assert that two DBCSR matrices have the same block 876 !! structure and same sparsity pattern 877 878 TYPE(dbcsr_type), INTENT(IN) :: matrix_a, matrix_b 879 LOGICAL :: is_equal 880 !! whether matrix_a and matrix_b have the same block structure 881 882 is_equal = .TRUE. 883 884 IF (dbcsr_nblkcols_total(matrix_a) .NE. dbcsr_nblkcols_total(matrix_b)) is_equal = .FALSE. 885 IF (dbcsr_nblkrows_total(matrix_a) .NE. dbcsr_nblkrows_total(matrix_b)) is_equal = .FALSE. 886 IF ((matrix_a%nblks) .NE. (matrix_b%nblks)) is_equal = .FALSE. 887 IF (ANY(matrix_a%row_p .NE. matrix_b%row_p)) is_equal = .FALSE. 888 IF (ANY(matrix_a%col_i .NE. matrix_b%col_i)) is_equal = .FALSE. 889 IF (ANY(dbcsr_row_block_sizes(matrix_a) .NE. & 890 dbcsr_row_block_sizes(matrix_b))) is_equal = .FALSE. 891 IF (ANY(dbcsr_row_block_sizes(matrix_a) .NE. & 892 dbcsr_row_block_sizes(matrix_b))) is_equal = .FALSE. 893 894 END FUNCTION dbcsr_has_same_block_structure 895 896 SUBROUTINE csr_assert_consistency_with_dbcsr(csr_mat, dbcsr_mat) 897 !! Helper function to assert that a given CSR matrix and a given DBCSR 898 !! matrix are consistent before doing the conversion 899 900 TYPE(csr_type), INTENT(IN) :: csr_mat 901 TYPE(dbcsr_type), INTENT(IN) :: dbcsr_mat 902 903 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_assert_consistency_with_dbcsr' 904 905 INTEGER :: handle 906 TYPE(csr_mapping_data) :: map 907 908 CALL timeset(routineN, handle) 909 map = csr_mat%dbcsr_mapping 910 IF (map%has_dbcsr_block_data) THEN 911 IF (map%dbcsr_nblkcols_total .NE. dbcsr_nblkcols_total(dbcsr_mat)) & 912 CALL dbcsr_abort(__LOCATION__, & 913 "field nblkcols_total of DBCSR matrix not consistent with CSR matrix") 914 IF (map%dbcsr_nblkrows_total .NE. dbcsr_nblkrows_total(dbcsr_mat)) & 915 CALL dbcsr_abort(__LOCATION__, & 916 "field nblkrows_total of DBCSR matrix not consistent with CSR matrix") 917 IF (map%dbcsr_nblks_local .NE. dbcsr_mat%nblks) & 918 CALL dbcsr_abort(__LOCATION__, & 919 "field nblks of DBCSR matrix not consistent with CSR matrix") 920 IF (ANY(map%dbcsr_row_p .NE. dbcsr_mat%row_p)) & 921 CALL dbcsr_abort(__LOCATION__, & 922 "field row_p of DBCSR matrix not consistent with CSR matrix") 923 IF (ANY(map%dbcsr_col_i .NE. dbcsr_mat%col_i)) & 924 CALL dbcsr_abort(__LOCATION__, & 925 "field dbcsr_col_i of DBCSR matrix not consistent with CSR matrix") 926 IF (ANY(map%dbcsr_row_blk_size .NE. dbcsr_row_block_sizes(dbcsr_mat))) & 927 CALL dbcsr_abort(__LOCATION__, & 928 "field row_blk_size of DBCSR matrix not consistent with CSR matrix") 929 IF (ANY(map%dbcsr_col_blk_size .NE. dbcsr_col_block_sizes(dbcsr_mat))) & 930 CALL dbcsr_abort(__LOCATION__, & 931 "field col_blk_size of DBCSR matrix not consistent with CSR matrix") 932 ELSE 933 CALL dbcsr_warn(__LOCATION__, "Can not assert consistency of the matrices "// & 934 "as no block data stored in CSR matrix.") 935 ENDIF 936 CALL timestop(handle) 937 END SUBROUTINE csr_assert_consistency_with_dbcsr 938 939 SUBROUTINE convert_dbcsr_to_csr(dbcsr_mat, csr_mat) 940 !! Convert a DBCSR matrix to a CSR matrix. 941 942 TYPE(dbcsr_type), INTENT(IN) :: dbcsr_mat 943 !! DBCSR matrix to convert 944 TYPE(csr_type), INTENT(INOUT) :: csr_mat 945 !! correctly allocated CSR matrix 946 947 CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_dbcsr_to_csr' 948 949 INTEGER :: handle 950 TYPE(dbcsr_type) :: dbcsr_mat_nosym 951 952 CALL timeset(routineN, handle) 953 954 IF (.NOT. dbcsr_valid_index(dbcsr_mat)) & 955 DBCSR_ABORT("Invalid DBCSR matrix") 956 IF (dbcsr_get_data_type(dbcsr_mat) /= csr_mat%nzval_local%data_type) & 957 DBCSR_ABORT("DBCSR and CSR matrix must have same type") 958 959 IF (.NOT. csr_mat%has_mapping) & 960 DBCSR_ABORT("CSR_mat must contain mapping to DBCSR format") 961 962 IF (dbcsr_has_symmetry(dbcsr_mat)) THEN 963 CALL dbcsr_desymmetrize_deep(dbcsr_mat, dbcsr_mat_nosym, untransposed_data=.TRUE.) 964 ELSE 965 dbcsr_mat_nosym = dbcsr_mat 966 ENDIF 967 968 CALL csr_assert_consistency_with_dbcsr(csr_mat, dbcsr_mat_nosym) 969 970 ! 1) DBCSR -> BRD 971 CALL dbcsr_complete_redistribute(dbcsr_mat_nosym, csr_mat%dbcsr_mapping%brd_mat) 972 ! 2) BRD -> CSR 973 CALL convert_brd_to_csr(csr_mat%dbcsr_mapping%brd_mat, csr_mat) 974 975 IF (dbcsr_has_symmetry(dbcsr_mat)) CALL dbcsr_release(dbcsr_mat_nosym) 976 977 CALL timestop(handle) 978 END SUBROUTINE convert_dbcsr_to_csr 979 980 SUBROUTINE convert_csr_to_dbcsr(dbcsr_mat, csr_mat) 981 !! convert a CSR matrix to a DBCSR matrix 982 983 TYPE(dbcsr_type), INTENT(INOUT) :: dbcsr_mat 984 !! correctly allocated DBCSR matrix 985 TYPE(csr_type), INTENT(INOUT) :: csr_mat 986 !! CSR matrix to convert 987 988 CHARACTER(LEN=*), PARAMETER :: routineN = 'convert_csr_to_dbcsr' 989 990 INTEGER :: handle 991 TYPE(dbcsr_type) :: dbcsr_mat_nosym 992 993 CALL timeset(routineN, handle) 994 995 IF (.NOT. dbcsr_valid_index(dbcsr_mat)) & 996 DBCSR_ABORT("Invalid DBCSR matrix") 997 998 IF (dbcsr_get_data_type(dbcsr_mat) /= csr_mat%nzval_local%data_type) & 999 DBCSR_ABORT("DBCSR and CSR matrix must have same type") 1000 1001 IF (.NOT. csr_mat%has_mapping) & 1002 DBCSR_ABORT("CSR_mat must contain mapping to DBCSR format") 1003 1004 ! Desymmetrize to assert that DBCSR matrix has sparsity pattern consistent with CSR matrix 1005 IF (dbcsr_has_symmetry(dbcsr_mat)) THEN 1006 CALL dbcsr_desymmetrize_deep(dbcsr_mat, dbcsr_mat_nosym, untransposed_data=.TRUE.) 1007 ELSE 1008 dbcsr_mat_nosym = dbcsr_mat 1009 ENDIF 1010 1011 CALL csr_assert_consistency_with_dbcsr(csr_mat, dbcsr_mat_nosym) 1012 1013 IF (dbcsr_has_symmetry(dbcsr_mat)) CALL dbcsr_release(dbcsr_mat_nosym) 1014 1015 ! 1) CSR -> BRD 1016 CALL convert_csr_to_brd(csr_mat%dbcsr_mapping%brd_mat, csr_mat) 1017 1018 ! 2) BRD -> DBCSR 1019 CALL dbcsr_complete_redistribute(csr_mat%dbcsr_mapping%brd_mat, dbcsr_mat) 1020 1021 CALL timestop(handle) 1022 END SUBROUTINE convert_csr_to_dbcsr 1023 1024 SUBROUTINE dbcsr_to_csr_filter(dbcsr_mat, csr_sparsity, eps) 1025 !! Apply filtering threshold eps to DBCSR blocks in order to improve 1026 !! CSR sparsity (currently only used for testing purposes) 1027 1028 TYPE(dbcsr_type), INTENT(IN) :: dbcsr_mat 1029 TYPE(dbcsr_type), INTENT(OUT) :: csr_sparsity 1030 REAL(kind=real_8), INTENT(IN) :: eps 1031 1032 INTEGER :: blkcol, blkrow, col_blk_size, data_type, & 1033 row_blk_size 1034 LOGICAL :: tr 1035 REAL(kind=real_8), ALLOCATABLE, DIMENSION(:) :: block_abs, csr_sparsity_blk 1036 TYPE(dbcsr_data_obj) :: block 1037 TYPE(dbcsr_iterator) :: iter 1038 1039!REAL(kind=real_8), DIMENSION(:), POINTER :: block 1040 1041 CALL dbcsr_create(csr_sparsity, & 1042 template=dbcsr_mat, & 1043 name="CSR sparsity", & 1044 data_type=dbcsr_type_real_8) 1045 CALL dbcsr_copy(csr_sparsity, dbcsr_mat) 1046 CALL dbcsr_set(csr_sparsity, 1.0_dp) 1047 1048 IF (eps .GT. 0.0_dp) THEN 1049 data_type = dbcsr_get_data_type(dbcsr_mat) 1050 CALL dbcsr_data_init(block) 1051 CALL dbcsr_data_new(block, data_type) 1052 CALL dbcsr_iterator_start(iter, dbcsr_mat, read_only=.TRUE.) 1053 DO WHILE (dbcsr_iterator_blocks_left(iter)) 1054 CALL dbcsr_iterator_next_block(iter, blkrow, blkcol, block, transposed=tr, & 1055 row_size=row_blk_size, col_size=col_blk_size) 1056 1057 ALLOCATE (block_abs(row_blk_size*col_blk_size)) 1058 ALLOCATE (csr_sparsity_blk(row_blk_size*col_blk_size)) 1059 SELECT CASE (data_type) 1060 CASE (dbcsr_type_real_4) 1061 block_abs(:) = REAL(ABS(block%d%r_sp(:)), KIND=real_8) 1062 CASE (dbcsr_type_real_8) 1063 block_abs(:) = REAL(ABS(block%d%r_dp(:)), KIND=real_8) 1064 CASE (dbcsr_type_complex_4) 1065 block_abs(:) = REAL(ABS(block%d%c_sp(:)), KIND=real_8) 1066 CASE (dbcsr_type_complex_8) 1067 block_abs(:) = REAL(ABS(block%d%c_dp(:)), KIND=real_8) 1068 END SELECT 1069 1070 csr_sparsity_blk = 1.0_dp 1071 WHERE (block_abs .LT. eps) csr_sparsity_blk = 0.0_dp 1072 CALL dbcsr_put_block(csr_sparsity, blkrow, blkcol, csr_sparsity_blk, transposed=tr) 1073 DEALLOCATE (csr_sparsity_blk, block_abs) 1074 ENDDO 1075 CALL dbcsr_iterator_stop(iter) 1076 CALL dbcsr_data_clear_pointer(block) 1077 CALL dbcsr_data_release(block) 1078 ENDIF 1079 1080 END SUBROUTINE dbcsr_to_csr_filter 1081 1082 SUBROUTINE csr_write(csr_mat, unit_nr, upper_triangle, threshold, binary) 1083 !! Write a CSR matrix to file 1084 1085 TYPE(csr_type), INTENT(IN) :: csr_mat 1086 INTEGER, INTENT(IN) :: unit_nr 1087 !! unit number to which output is written 1088 LOGICAL, INTENT(IN), OPTIONAL :: upper_triangle 1089 !! If true (default: false), write only upper triangular part of matrix 1090 REAL(KIND=real_8), INTENT(IN), OPTIONAL :: threshold 1091 !! threshold on the absolute value of the elements to be printed 1092 LOGICAL, INTENT(IN), OPTIONAL :: binary 1093 1094 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_write' 1095 CHARACTER(LEN=default_string_length) :: data_format 1096 COMPLEX(KIND=real_4), ALLOCATABLE, DIMENSION(:) :: nzval_to_master_c_sp 1097 COMPLEX(KIND=real_8), ALLOCATABLE, DIMENSION(:) :: nzval_to_master_c_dp 1098 INTEGER :: handle, i, ii, k, l, m, mynode, & 1099 numnodes, rowind, tag1, tag2, tag3 1100 INTEGER, ALLOCATABLE, DIMENSION(:) :: colind_to_master, nzerow_to_master, & 1101 sizes_numrowlocal, sizes_nzelocal 1102 LOGICAL :: bin, ut 1103 REAL(KIND=real_4), ALLOCATABLE, DIMENSION(:) :: nzval_to_master_r_sp 1104 REAL(KIND=real_8) :: thld 1105 REAL(KIND=real_8), ALLOCATABLE, DIMENSION(:) :: nzval_to_master_r_dp 1106 1107 CALL timeset(routineN, handle) 1108 1109 IF (PRESENT(upper_triangle)) THEN 1110 ut = upper_triangle 1111 ELSE 1112 ut = .FALSE. 1113 ENDIF 1114 1115 IF (PRESENT(threshold)) THEN 1116 thld = threshold 1117 ELSE 1118 thld = 0.0_dp 1119 ENDIF 1120 1121 IF (PRESENT(binary)) THEN 1122 bin = binary 1123 ELSE 1124 bin = .FALSE. 1125 ENDIF 1126 1127 IF (.NOT. csr_mat%valid) & 1128 DBCSR_ABORT("can not write invalid CSR matrix") 1129 1130 tag1 = 0 1131 tag2 = 1 1132 tag3 = 2 1133 1134 CALL mp_environ(numnodes, mynode, csr_mat%mp_group) 1135 1136 ! gather sizes (number of local non-zero elements and number of local rows) 1137 ALLOCATE (sizes_nzelocal(numnodes)) 1138 ALLOCATE (sizes_numrowlocal(numnodes)) 1139 1140 CALL mp_gather(csr_mat%nze_local, sizes_nzelocal, 0, csr_mat%mp_group) 1141 CALL mp_gather(csr_mat%nrows_local, sizes_numrowlocal, 0, csr_mat%mp_group) 1142 1143 ! for each node, send matrix data to node 0 (master) and write data 1144 DO i = 0, numnodes - 1 1145 ii = i 1146 IF (mynode .EQ. 0) THEN ! allocations for receiving data from node i 1147 ALLOCATE (colind_to_master(sizes_nzelocal(ii + 1))) 1148 ALLOCATE (nzerow_to_master(sizes_numrowlocal(ii + 1))) 1149 1150 SELECT CASE (csr_mat%nzval_local%data_type) 1151 CASE (dbcsr_type_real_4) 1152 data_format = "(2(I8),E23.6E2)" 1153 ALLOCATE (nzval_to_master_r_sp(sizes_nzelocal(ii + 1))) 1154 CASE (dbcsr_type_real_8) 1155 data_format = "(2(I8),E23.14E3)" 1156 ALLOCATE (nzval_to_master_r_dp(sizes_nzelocal(ii + 1))) 1157 CASE (dbcsr_type_complex_4) 1158 data_format = "(2(I8),2(E23.6E2))" 1159 ALLOCATE (nzval_to_master_c_sp(sizes_nzelocal(ii + 1))) 1160 CASE (dbcsr_type_complex_8) 1161 data_format = "(2(I8),2(E23.14E3))" 1162 ALLOCATE (nzval_to_master_c_dp(sizes_nzelocal(ii + 1))) 1163 END SELECT 1164 ENDIF 1165 1166 IF (mynode .EQ. 0) THEN ! receive at node 0 1167 IF (ii .EQ. 0) THEN ! data from node 0, no need for mpi routines 1168 colind_to_master(:) = csr_mat%colind_local(:) 1169 nzerow_to_master(:) = csr_mat%nzerow_local(:) 1170 SELECT CASE (csr_mat%nzval_local%data_type) 1171 CASE (dbcsr_type_real_4) 1172 nzval_to_master_r_sp(:) = csr_mat%nzval_local%r_sp(:) 1173 CASE (dbcsr_type_real_8) 1174 nzval_to_master_r_dp(:) = csr_mat%nzval_local%r_dp(:) 1175 CASE (dbcsr_type_complex_4) 1176 nzval_to_master_c_sp(:) = csr_mat%nzval_local%c_sp(:) 1177 CASE (dbcsr_type_complex_8) 1178 nzval_to_master_c_dp(:) = csr_mat%nzval_local%c_dp(:) 1179 END SELECT 1180 ELSE ! receive data from nodes with rank > 0 1181 CALL mp_recv(colind_to_master, ii, tag1, csr_mat%mp_group) 1182 CALL mp_recv(nzerow_to_master, ii, tag2, csr_mat%mp_group) 1183 SELECT CASE (csr_mat%nzval_local%data_type) 1184 CASE (dbcsr_type_real_4) 1185 CALL mp_recv(nzval_to_master_r_sp, ii, tag3, csr_mat%mp_group) 1186 CASE (dbcsr_type_real_8) 1187 CALL mp_recv(nzval_to_master_r_dp, ii, tag3, csr_mat%mp_group) 1188 CASE (dbcsr_type_complex_4) 1189 CALL mp_recv(nzval_to_master_c_sp, ii, tag3, csr_mat%mp_group) 1190 CASE (dbcsr_type_complex_8) 1191 CALL mp_recv(nzval_to_master_c_dp, ii, tag3, csr_mat%mp_group) 1192 END SELECT 1193 ENDIF 1194 ENDIF 1195 1196 IF ((mynode .EQ. ii) .AND. (ii .NE. 0)) THEN ! send from nodes with rank > 0 1197 CALL mp_send(csr_mat%colind_local, 0, tag1, csr_mat%mp_group) 1198 CALL mp_send(csr_mat%nzerow_local, 0, tag2, csr_mat%mp_group) 1199 SELECT CASE (csr_mat%nzval_local%data_type) 1200 CASE (dbcsr_type_real_4) 1201 CALL mp_send(csr_mat%nzval_local%r_sp, 0, tag3, csr_mat%mp_group) 1202 CASE (dbcsr_type_real_8) 1203 CALL mp_send(csr_mat%nzval_local%r_dp, 0, tag3, csr_mat%mp_group) 1204 CASE (dbcsr_type_complex_4) 1205 CALL mp_send(csr_mat%nzval_local%c_sp, 0, tag3, csr_mat%mp_group) 1206 CASE (dbcsr_type_complex_8) 1207 CALL mp_send(csr_mat%nzval_local%c_dp, 0, tag3, csr_mat%mp_group) 1208 END SELECT 1209 ENDIF 1210 1211 IF (mynode .EQ. 0) THEN ! write data received at node 0 1212 !WRITE(unit_nr,"(A27)") "#row ind, col ind, value" 1213 m = 0 1214 DO k = 1, sizes_numrowlocal(ii + 1) 1215 rowind = k + SUM(sizes_numrowlocal(1:ii)) ! row index: local to global 1216 DO l = 1, nzerow_to_master(k) 1217 m = m + 1 1218 IF ((.NOT. ut) .OR. (rowind .LE. colind_to_master(m))) THEN 1219 SELECT CASE (csr_mat%nzval_local%data_type) 1220 CASE (dbcsr_type_real_4) 1221 IF (ABS(nzval_to_master_r_sp(m)) .GE. thld) THEN 1222 IF (bin) THEN 1223 WRITE (unit_nr) rowind, colind_to_master(m), nzval_to_master_r_sp(m) 1224 ELSE 1225 WRITE (unit_nr, data_format) rowind, colind_to_master(m), & 1226 nzval_to_master_r_sp(m) 1227 END IF 1228 END IF 1229 CASE (dbcsr_type_real_8) 1230 IF (ABS(nzval_to_master_r_dp(m)) .GE. thld) THEN 1231 IF (bin) THEN 1232 WRITE (unit_nr) rowind, colind_to_master(m), nzval_to_master_r_dp(m) 1233 ELSE 1234 WRITE (unit_nr, data_format) rowind, colind_to_master(m), & 1235 nzval_to_master_r_dp(m) 1236 END IF 1237 END IF 1238 CASE (dbcsr_type_complex_4) 1239 IF (ABS(nzval_to_master_c_sp(m)) .GE. thld) THEN 1240 IF (bin) THEN 1241 WRITE (unit_nr) rowind, colind_to_master(m), nzval_to_master_c_sp(m) 1242 ELSE 1243 WRITE (unit_nr, data_format) rowind, colind_to_master(m), & 1244 nzval_to_master_c_sp(m) 1245 END IF 1246 END IF 1247 CASE (dbcsr_type_complex_8) 1248 IF (ABS(nzval_to_master_c_dp(m)) .GE. thld) THEN 1249 IF (bin) THEN 1250 WRITE (unit_nr) rowind, colind_to_master(m), nzval_to_master_c_dp(m) 1251 ELSE 1252 WRITE (unit_nr, data_format) rowind, colind_to_master(m), & 1253 nzval_to_master_c_dp(m) 1254 END IF 1255 END IF 1256 END SELECT 1257 ENDIF 1258 ENDDO 1259 ENDDO 1260 1261 DEALLOCATE (colind_to_master) 1262 DEALLOCATE (nzerow_to_master) 1263 1264 SELECT CASE (csr_mat%nzval_local%data_type) 1265 CASE (dbcsr_type_real_4) 1266 DEALLOCATE (nzval_to_master_r_sp) 1267 CASE (dbcsr_type_real_8) 1268 DEALLOCATE (nzval_to_master_r_dp) 1269 CASE (dbcsr_type_complex_4) 1270 DEALLOCATE (nzval_to_master_c_sp) 1271 CASE (dbcsr_type_complex_8) 1272 DEALLOCATE (nzval_to_master_c_dp) 1273 END SELECT 1274 ENDIF 1275 ENDDO 1276 1277 CALL timestop(handle) 1278 1279 END SUBROUTINE csr_write 1280 1281 SUBROUTINE csr_print_sparsity(csr_mat, unit_nr) 1282 !! Print CSR sparsity 1283 TYPE(csr_type), INTENT(IN) :: csr_mat 1284 INTEGER, INTENT(IN) :: unit_nr 1285 1286 CHARACTER(LEN=*), PARAMETER :: routineN = 'csr_print_sparsity' 1287 1288 INTEGER :: handle, mynode, numnodes 1289 INTEGER(KIND=int_8) :: dbcsr_nze_total 1290 REAL(KIND=real_8) :: dbcsr_nze_percentage, nze_percentage 1291 1292 CALL timeset(routineN, handle) 1293 1294 IF (.NOT. csr_mat%valid) & 1295 DBCSR_ABORT("CSR matrix must be created first") 1296 1297 nze_percentage = 100.0_dp*(REAL(csr_mat%nze_total, KIND=real_8) & 1298 /REAL(csr_mat%nrows_total, KIND=real_8)) & 1299 /REAL(csr_mat%ncols_total, KIND=real_8) 1300 1301 IF (csr_mat%has_mapping) THEN 1302 dbcsr_nze_total = dbcsr_get_nze(csr_mat%dbcsr_mapping%brd_mat) 1303 CALL mp_sum(dbcsr_nze_total, csr_mat%mp_group) 1304 dbcsr_nze_percentage = 100.0_dp*(REAL(dbcsr_nze_total, KIND=real_8) & 1305 /REAL(csr_mat%nrows_total, KIND=real_8)) & 1306 /REAL(csr_mat%ncols_total, KIND=real_8) 1307 ENDIF 1308 1309 CALL mp_environ(numnodes, mynode, csr_mat%mp_group) 1310 1311 IF (mynode .EQ. 0) THEN 1312 WRITE (unit_nr, "(T15,A,T68,I13)") "Number of CSR non-zero elements:", csr_mat%nze_total 1313 WRITE (unit_nr, "(T15,A,T75,F6.2)") "Percentage CSR non-zero elements:", nze_percentage 1314 !IF(csr_mat%has_mapping) THEN 1315 ! WRITE(unit_nr,"(T15,A,T75,F6.2/)") "Percentage DBCSR non-zero elements:", dbcsr_nze_percentage 1316 !ENDIF 1317 ENDIF 1318 1319 CALL timestop(handle) 1320 END SUBROUTINE csr_print_sparsity 1321 1322 SUBROUTINE dbcsr_create_brd(dbcsr_mat, brd_mat, equal_dist, floor_dist, numnodes) 1323 !! Converts a DBCSR matrix to a block row distributed matrix. 1324 1325 TYPE(dbcsr_type), INTENT(IN) :: dbcsr_mat 1326 !! DBCSR matrix to be converted 1327 TYPE(dbcsr_type), INTENT(OUT) :: brd_mat 1328 !! converted matrix 1329 LOGICAL, INTENT(IN) :: equal_dist, floor_dist 1330 !! see documentation of csr_create_from_dbcsr 1331 !! see documentation of csr_create_from_dbcsr 1332 INTEGER, INTENT(IN) :: numnodes 1333 !! number of nodes to use for block row distribution 1334 1335 CHARACTER(LEN=*), PARAMETER :: routineN = 'dbcsr_create_brd' 1336 1337 CHARACTER :: matrix_type 1338 CHARACTER(LEN=default_string_length) :: matrix_name 1339 INTEGER :: cs, data_type, end_ind, handle, i, k, l, m, mp_group, mynode, nblkcols_total, & 1340 nblkrows_total, nfullrows_local, nfullrows_total, node_size, numnodes_total, row_index, & 1341 split_row, start_ind 1342 INTEGER, ALLOCATABLE, DIMENSION(:) :: rdist_tmp, row_blk_size_new_tmp 1343 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: pgrid 1344 INTEGER, DIMENSION(:), POINTER, CONTIGUOUS :: cdist, col_blk_size, rdist, & 1345 row_blk_size, row_blk_size_new 1346 REAL(KIND=real_8) :: chunk_size 1347 TYPE(dbcsr_distribution_obj) :: dist_current, dist_new 1348 TYPE(dbcsr_mp_obj) :: mp_obj_current, mp_obj_new 1349 1350 CALL timeset(routineN, handle) 1351 1352 NULLIFY (row_blk_size, rdist, row_blk_size_new) 1353 CALL dbcsr_get_info(dbcsr_mat, & 1354 nblkrows_total=nblkrows_total, & 1355 nblkcols_total=nblkcols_total, & 1356 nfullrows_total=nfullrows_total, & 1357 row_blk_size=row_blk_size, & 1358 col_blk_size=col_blk_size, & 1359 matrix_type=matrix_type, & 1360 data_type=data_type) 1361 1362 matrix_name = dbcsr_name(dbcsr_mat) 1363 1364 ALLOCATE (cdist(nblkcols_total)) 1365 cdist = 0 1366 1367 dist_current = dbcsr_distribution(dbcsr_mat) 1368 mp_obj_current = dbcsr_distribution_mp(dist_current) 1369 mp_group = dbcsr_mp_group(mp_obj_current) 1370 mynode = dbcsr_mp_mynode(mp_obj_current) 1371 numnodes_total = dbcsr_mp_numnodes(mp_obj_current) 1372 1373 ALLOCATE (pgrid(numnodes_total, 1)) 1374 1375 IF (equal_dist) THEN ! Equally distribute rows over processors -> cut blocks 1376 1377 ! Calculate the number of rows a processor can hold 1378 IF (floor_dist) THEN 1379 nfullrows_local = FLOOR(REAL(nfullrows_total, KIND=dp)/numnodes) 1380 ELSE 1381 nfullrows_local = CEILING(REAL(nfullrows_total, KIND=dp)/numnodes) 1382 ENDIF 1383 1384 ! allocate maximum amount of memory possibly needed 1385 ALLOCATE (rdist_tmp(nblkrows_total + numnodes - 1)) ! row distribution 1386 ALLOCATE (row_blk_size_new_tmp(nblkrows_total + numnodes - 1)) ! new sizes of block rows 1387 1388 k = 0 ! counter for block rows 1389 m = 0 ! node counter 1390 node_size = nfullrows_local ! space available on current node in number of rows 1391 IF (node_size .GT. 0) THEN 1392 DO l = 1, nblkrows_total 1393 split_row = row_blk_size(l) ! size of current block row (number of rows) 1394 DO WHILE (split_row .GE. node_size) ! cut block row and send it to two nodes 1395 k = k + 1 1396 m = m + 1 1397 row_blk_size_new_tmp(k) = node_size ! size of first part of block row 1398 rdist_tmp(k) = m - 1 ! send first part to node m 1399 split_row = split_row - node_size ! size of remaining part of block rows 1400 node_size = nfullrows_local ! space available on next node 1401 IF (floor_dist .AND. (m .EQ. numnodes - 1)) THEN ! send all remaining rows to last node 1402 node_size = nfullrows_total - (numnodes - 1)*node_size 1403 ENDIF 1404 ENDDO 1405 IF (split_row .GT. 0) THEN ! enough space left on next node for remaining rows 1406 k = k + 1 1407 row_blk_size_new_tmp(k) = split_row ! size of remaining part of block row 1408 rdist_tmp(k) = m ! send to next node 1409 node_size = node_size - split_row ! remaining space on next node 1410 ENDIF 1411 ENDDO 1412 ELSE ! send everything to last node if node_size = 0 1413 rdist_tmp(1:nblkrows_total) = numnodes - 1 1414 row_blk_size_new_tmp(1:nblkrows_total) = row_blk_size ! row blocks unchanged 1415 k = nblkrows_total 1416 ENDIF 1417 1418 ! Copy data to correctly allocated variables 1419 ALLOCATE (row_blk_size_new(k)) 1420 row_blk_size_new = row_blk_size_new_tmp(1:k) 1421 ALLOCATE (rdist(k)) 1422 rdist = rdist_tmp(1:k) 1423 1424 ELSE ! Leave block rows intact (do not cut) 1425 ALLOCATE (rdist(nblkrows_total)) 1426 rdist = 0 1427 IF (numnodes .GT. nblkrows_total) THEN 1428 rdist = (/(i, i=0, nblkrows_total - 1)/) 1429 ELSE 1430 chunk_size = REAL(nblkrows_total, KIND=dp)/numnodes 1431 row_index = 0 1432 start_ind = 1 1433 DO i = 0, numnodes - 1 1434 cs = NINT(i*chunk_size) - NINT((i - 1)*chunk_size) 1435 end_ind = MIN(start_ind - 1 + cs, nblkrows_total) 1436 rdist(start_ind:end_ind) = row_index 1437 start_ind = end_ind + 1 1438 row_index = row_index + 1 1439 END DO 1440 END IF 1441 row_blk_size_new => row_blk_size 1442 ENDIF 1443 1444 pgrid(:, :) = RESHAPE((/(i, i=0, numnodes_total - 1)/), (/numnodes_total, 1/)) 1445 CALL dbcsr_mp_new(mp_obj_new, mp_group, pgrid, mynode, numnodes=numnodes_total) 1446 CALL dbcsr_distribution_new(dist_new, mp_obj_new, rdist, cdist, reuse_arrays=.TRUE.) 1447 1448 CALL dbcsr_create(brd_mat, TRIM(matrix_name)//" row-block distributed", & 1449 dist_new, matrix_type, row_blk_size_new, col_blk_size, data_type=data_type) 1450 CALL dbcsr_complete_redistribute(dbcsr_mat, brd_mat) 1451 1452 DEALLOCATE (pgrid) 1453 1454 IF (equal_dist) DEALLOCATE (row_blk_size_new) 1455 1456 CALL dbcsr_distribution_release(dist_new) 1457 CALL dbcsr_mp_release(mp_obj_new) 1458 1459 CALL timestop(handle) 1460 1461 END SUBROUTINE dbcsr_create_brd 1462 1463END MODULE dbcsr_csr_conversions 1464