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_index_operations 11 !! Operations on the DBCSR index 12 13 USE dbcsr_array_types, ONLY: array_data, & 14 array_size 15 USE dbcsr_config, ONLY: default_resize_factor 16 USE dbcsr_dist_methods, ONLY: dbcsr_distribution_local_cols, & 17 dbcsr_distribution_local_rows, & 18 dbcsr_distribution_ncols, & 19 dbcsr_distribution_nlocal_cols, & 20 dbcsr_distribution_nlocal_rows, & 21 dbcsr_distribution_nrows, & 22 dbcsr_distribution_thread_dist 23 USE dbcsr_dist_operations, ONLY: get_stored_canonical 24 USE dbcsr_kinds, ONLY: int_4, & 25 int_8 26 USE dbcsr_methods, ONLY: dbcsr_distribution, & 27 dbcsr_get_index_memory_type 28 USE dbcsr_ptr_util, ONLY: ensure_array_size, & 29 memory_allocate, & 30 memory_deallocate 31 USE dbcsr_toollib, ONLY: joaat_hash, & 32 sort, & 33 swap 34 USE dbcsr_types, ONLY: & 35 dbcsr_distribution_obj, dbcsr_meta_size, dbcsr_num_slots, dbcsr_slot_blk_p, & 36 dbcsr_slot_col_i, dbcsr_slot_coo_l, dbcsr_slot_home_prow, dbcsr_slot_home_vprow, & 37 dbcsr_slot_nblkcols_local, dbcsr_slot_nblkcols_total, dbcsr_slot_nblkrows_local, & 38 dbcsr_slot_nblkrows_total, dbcsr_slot_nblks, dbcsr_slot_nfullcols_local, dbcsr_slot_nze, & 39 dbcsr_slot_row_p, dbcsr_slot_size, dbcsr_slot_thr_c, dbcsr_type 40#include "base/dbcsr_base_uses.f90" 41 42!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads 43 44 IMPLICIT NONE 45 46 PRIVATE 47 48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbcsr_index_operations' 49 50 LOGICAL, PARAMETER :: careful_mod = .FALSE. 51 LOGICAL, PARAMETER :: debug_mod = .FALSE. 52 53 ! Index transformations 54 PUBLIC :: dbcsr_make_index_canonical, & 55 transpose_index_local 56 PUBLIC :: dbcsr_make_index_local_row, dbcsr_has_local_row_index 57 PUBLIC :: dbcsr_make_index_list 58 ! Dense/Sparse 59 PUBLIC :: make_dense_index, make_undense_index 60 ! Working with DBCSR and linear indices 61 PUBLIC :: dbcsr_make_dbcsr_index, dbcsr_sort_indices, & 62 merge_index_arrays, & 63 dbcsr_expand_row_index, & 64 dbcsr_count_row_index, dbcsr_build_row_index, & 65 dbcsr_index_prune_deleted, & 66 dbcsr_index_compact 67 PUBLIC :: dbcsr_index_checksum 68 ! Index array manipulation 69 PUBLIC :: dbcsr_addto_index_array, dbcsr_clearfrom_index_array, & 70 dbcsr_repoint_index, dbcsr_make_index_exist 71 72 INTERFACE dbcsr_count_row_index 73 MODULE PROCEDURE dbcsr_count_row_index_copy, & 74 dbcsr_count_row_index_inplace 75 END INTERFACE 76 77 INTERFACE dbcsr_build_row_index 78 MODULE PROCEDURE dbcsr_build_row_index_copy, & 79 dbcsr_build_row_index_inplace 80 END INTERFACE 81 82CONTAINS 83 84 SUBROUTINE make_index_canonical(new_row_p, new_col_i, new_blk_p, & 85 old_row_p, old_col_i, old_blk_p, matrix) 86 !! Makes a canonical index given the index arrays 87 !! 88 !! Description of canonical ordering 89 !! A non-(anti)symmetric matrix is left as is. Otherwise, the row and column 90 !! are stored in the position prescribed by the distribution. 91 !! @note 92 !! This routine uses hard-coded logic as to what constitutes a 93 !! canonical ordering 94 95 INTEGER, DIMENSION(:), INTENT(OUT) :: new_row_p, new_col_i, new_blk_p 96 INTEGER, DIMENSION(:), INTENT(IN) :: old_row_p, old_col_i, old_blk_p 97 TYPE(dbcsr_type), INTENT(IN) :: matrix 98 99 CHARACTER(len=*), PARAMETER :: routineN = 'make_index_canonical' 100 101 INTEGER :: blk, col, nblks, row, stored_col, & 102 stored_row 103 INTEGER, ALLOCATABLE, DIMENSION(:) :: row_i 104 LOGICAL :: tr 105 106! --------------------------------------------------------------------------- 107 108 nblks = SIZE(old_blk_p) 109 ALLOCATE (row_i(nblks)) 110 IF (debug_mod) THEN 111 WRITE (*, *) "old row_p", old_row_p 112 WRITE (*, *) "old col_i", old_col_i 113 WRITE (*, *) "old blk_p", old_blk_p 114 ENDIF 115 DO row = 1, SIZE(old_row_p) - 1 116 DO blk = old_row_p(row) + 1, old_row_p(row + 1) 117 col = old_col_i(blk) 118 stored_row = row 119 stored_col = col 120 tr = .FALSE. 121 CALL get_stored_canonical(matrix, stored_row, stored_col, tr) 122 IF (debug_mod) & 123 WRITE (*, '(A,2(1X,I5),A,2(1X,I5),";",I7,1X,L1)') & 124 routineN//" X->", row, col, "->", & 125 stored_row, stored_col, blk, tr 126 row_i(blk) = stored_row 127 new_col_i(blk) = stored_col 128 IF (.NOT. tr) THEN 129 new_blk_p(blk) = old_blk_p(blk) 130 ELSE 131 new_blk_p(blk) = -old_blk_p(blk) 132 ENDIF 133 ENDDO 134 ENDDO 135 CALL dbcsr_sort_indices(nblks, row_i, new_col_i, blk_p=new_blk_p) 136 ! Re-create the index 137 CALL dbcsr_make_dbcsr_index(new_row_p, row_i, SIZE(new_row_p) - 1, nblks) 138 IF (debug_mod) THEN 139 WRITE (*, *) "new row_p", new_row_p 140 WRITE (*, *) "new row_i", row_i 141 WRITE (*, *) "new col_i", new_col_i 142 WRITE (*, *) "new blk_p", new_blk_p 143 ENDIF 144 END SUBROUTINE make_index_canonical 145 146 SUBROUTINE make_index_triangular(new_row_p, new_col_i, new_blk_p, & 147 old_row_p, old_col_i, old_blk_p, matrix) 148 !! Makes a CP2K triangular index given the index arrays 149 !! 150 !! Description of canonical ordering 151 !! A non-(anti)symmetric matrix is left as is. Otherwise, the row and column 152 !! are stored in the position prescribed by the distribution. 153 !! @note 154 !! This routine uses hard-coded logic as to what constitutes a 155 !! canonical ordering 156 157 INTEGER, DIMENSION(:), INTENT(OUT) :: new_row_p, new_col_i, new_blk_p 158 INTEGER, DIMENSION(:), INTENT(IN) :: old_row_p, old_col_i, old_blk_p 159 TYPE(dbcsr_type), INTENT(IN) :: matrix 160 161 CHARACTER(len=*), PARAMETER :: routineN = 'make_index_triangular' 162 163 INTEGER :: blk, col, nblks, row, stored_col, & 164 stored_row 165 INTEGER, ALLOCATABLE, DIMENSION(:) :: row_i 166 LOGICAL :: tr 167 168! --------------------------------------------------------------------------- 169 170 nblks = SIZE(old_blk_p) 171 ALLOCATE (row_i(nblks)) 172 IF (debug_mod) THEN 173 WRITE (*, *) "old row_p", old_row_p 174 WRITE (*, *) "old col_i", old_col_i 175 WRITE (*, *) "old blk_p", old_blk_p 176 ENDIF 177 DO row = 1, SIZE(old_row_p) - 1 178 DO blk = old_row_p(row) + 1, old_row_p(row + 1) 179 col = old_col_i(blk) 180 stored_row = row 181 stored_col = col 182 tr = .FALSE. 183 CALL get_stored_canonical(matrix, stored_row, stored_col, tr) 184 IF (stored_row .GT. stored_col) THEN 185 CALL swap(stored_row, stored_col) 186 tr = .NOT. tr 187 ENDIF 188 IF (debug_mod) & 189 WRITE (*, '(A,2(1X,I5),A,2(1X,I5),";",I7,1X,L1)') & 190 routineN//" X->", row, col, "->", & 191 stored_row, stored_col, blk, tr 192 row_i(blk) = stored_row 193 new_col_i(blk) = stored_col 194 IF (.NOT. tr) THEN 195 new_blk_p(blk) = old_blk_p(blk) 196 ELSE 197 new_blk_p(blk) = -old_blk_p(blk) 198 ENDIF 199 ENDDO 200 ENDDO 201 CALL dbcsr_sort_indices(nblks, row_i, new_col_i, blk_p=new_blk_p) 202 ! Re-create the index 203 CALL dbcsr_make_dbcsr_index(new_row_p, row_i, SIZE(new_row_p) - 1, nblks) 204 IF (debug_mod) THEN 205 WRITE (*, *) "new row_p", new_row_p 206 WRITE (*, *) "new row_i", row_i 207 WRITE (*, *) "new col_i", new_col_i 208 WRITE (*, *) "new blk_p", new_blk_p 209 ENDIF 210 END SUBROUTINE make_index_triangular 211 212 SUBROUTINE dbcsr_make_dbcsr_index(row_p, row_i, nrows, nblks) 213 !! Collapses a row_p index 214 INTEGER, INTENT(in) :: nrows, nblks 215 INTEGER, DIMENSION(1:nrows + 1), INTENT(out) :: row_p 216 INTEGER, DIMENSION(1:nblks), INTENT(in) :: row_i 217 218 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_make_dbcsr_index' 219 220 INTEGER :: blk, error_handle, row 221 222 CALL timeset(routineN, error_handle) 223 224 row_p(1) = 0 225 row_p(nrows + 1) = nblks 226 row = 1 227 blk = 1 228 DO WHILE (row .LE. nrows) 229 IF (blk .LE. nblks) THEN 230 DO WHILE (row_i(blk) .EQ. row) 231 blk = blk + 1 232 IF (blk .GT. nblks) THEN 233 row_p(row + 1) = nblks - 1 234 EXIT 235 ENDIF 236 ENDDO 237 ENDIF 238 row_p(row + 1) = blk - 1 239 row = row + 1 240 ENDDO 241 242 CALL timestop(error_handle) 243 244 END SUBROUTINE dbcsr_make_dbcsr_index 245 246 PURE SUBROUTINE dbcsr_expand_row_index(row_p, row_i, nrows, nblks) 247 !! Expands a row_p index 248 INTEGER, INTENT(IN) :: nrows, nblks 249 INTEGER, DIMENSION(1:nrows + 1), INTENT(IN) :: row_p 250 INTEGER, DIMENSION(1:nblks), INTENT(OUT) :: row_i 251 252 INTEGER :: row 253 254 DO row = 1, nrows 255 row_i(row_p(row) + 1:row_p(row + 1)) = row 256 ENDDO 257 END SUBROUTINE dbcsr_expand_row_index 258 259 PURE SUBROUTINE dbcsr_expand_row_index_2d(row_p, row_i, nrows, dst_i) 260 !! Expands a row_p index 261 INTEGER, INTENT(IN) :: nrows, dst_i 262 INTEGER, DIMENSION(1:nrows + 1), INTENT(IN) :: row_p 263 INTEGER, DIMENSION(:, :), INTENT(OUT) :: row_i 264 265 INTEGER :: row 266 267 DO row = 1, nrows 268 row_i(dst_i, row_p(row) + 1:row_p(row + 1)) = row 269 ENDDO 270 END SUBROUTINE dbcsr_expand_row_index_2d 271 272 PURE SUBROUTINE dbcsr_count_row_index_inplace(rows, nrows) 273 !! Counts columns-per-row count from row index array, in-place. 274 275 INTEGER, INTENT(IN) :: nrows 276 !! number of rows 277 INTEGER, DIMENSION(1:nrows + 1), INTENT(INOUT) :: rows 278 !! the row_p index (input); the count of the number of columns per row (output) 279 280 INTEGER :: row 281 282 DO row = 1, nrows 283 rows(row) = rows(row + 1) - rows(row) 284 ENDDO 285 rows(nrows + 1) = 0 286 END SUBROUTINE dbcsr_count_row_index_inplace 287 288 PURE SUBROUTINE dbcsr_count_row_index_copy(rows, counts, nrows) 289 !! Counts columns-per-row count from row index array. 290 291 INTEGER, INTENT(IN) :: nrows 292 !! number of rows 293 INTEGER, DIMENSION(1:nrows), INTENT(OUT) :: counts 294 !! the count of the number of columns per row 295 INTEGER, DIMENSION(1:nrows + 1), INTENT(IN) :: rows 296 !! the row_p index (input) 297 298 INTEGER :: row 299 300 DO row = 1, nrows 301 counts(row) = rows(row + 1) - rows(row) 302 END DO 303 END SUBROUTINE dbcsr_count_row_index_copy 304 305 PURE SUBROUTINE dbcsr_build_row_index_inplace(rows, nrows) 306 !! Builds row index array from a columns-per-row count, in-place. 307 308 INTEGER, INTENT(IN) :: nrows 309 !! number of rows 310 INTEGER, DIMENSION(1:nrows + 1), INTENT(INOUT) :: rows 311 !! count of the number of columns per row (input); the row_p index (output) 312 313 INTEGER :: o, old_count, row 314 315 old_count = rows(1) 316 rows(1) = 0 317 IF (nrows .GE. 1) THEN 318 DO row = 2, nrows + 1 319 o = rows(row) 320 rows(row) = rows(row - 1) + old_count 321 old_count = o 322 ENDDO 323 ENDIF 324 END SUBROUTINE dbcsr_build_row_index_inplace 325 326 PURE SUBROUTINE dbcsr_build_row_index_copy(counts, rows, nrows) 327 !! Builds row index array from a columns-per-row count. 328 329 INTEGER, INTENT(IN) :: nrows 330 !! number of rows 331 INTEGER, DIMENSION(1:nrows + 1), INTENT(OUT) :: rows 332 !! count of the number of columns per row (input); the row_p index (output) 333 INTEGER, DIMENSION(1:nrows), INTENT(IN) :: counts 334 !! count of the number of columns per row 335 336!WTF?!rows(1) = 0 337!WTF?!IF (nrows .GE. 1) THEN 338!WTF?! DO row = 2, nrows+1 339!WTF?! rows(row) = rows(row-1) + counts(rows-1) 340!WTF?! ENDDO 341!WTF?!ENDIF 342 343 rows(1:nrows) = counts(1:nrows) 344 CALL dbcsr_build_row_index_inplace(rows, nrows) 345 END SUBROUTINE dbcsr_build_row_index_copy 346 347 SUBROUTINE dbcsr_addto_index_array(matrix, slot, DATA, reservation, extra) 348 !! Adds data to the index. Increases the index size when necessary. 349 350 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 351 !! bcsr matrix 352 INTEGER, INTENT(IN) :: slot 353 !! which index array to add (e.g., dbcsr_slot_row_blk_sizes) 354 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: DATA 355 !! array holding the index data to add to the index array 356 INTEGER, INTENT(IN), OPTIONAL :: reservation, extra 357 !! only reserve space for subsequent array 358 !! reserve extra space for later additions 359 360 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_addto_index_array', & 361 routineP = moduleN//':'//routineN 362 363 INTEGER :: deplus, space, ub, ub_new 364 365! --------------------------------------------------------------------------- 366 367 IF (debug_mod) THEN 368 IF (.NOT. ASSOCIATED(matrix%index)) & 369 DBCSR_ABORT("Index must be preallocated.") 370 IF (UBOUND(matrix%index, 1) < dbcsr_num_slots) & 371 DBCSR_ABORT("Actual index size less than declared size") 372 IF (.NOT. PRESENT(DATA) .AND. .NOT. PRESENT(reservation)) & 373 DBCSR_ABORT('Either an array or its size must be specified.') 374 WRITE (*, *) routineP//' index', matrix%index(:dbcsr_num_slots) 375 ENDIF 376 IF (PRESENT(reservation)) THEN 377 space = reservation 378 ELSE 379 space = SIZE(DATA) 380 ENDIF 381 IF (PRESENT(extra)) THEN 382 deplus = extra 383 ELSE 384 deplus = 0 385 ENDIF 386 ub = UBOUND(matrix%index, 1) 387 ! The data area was not defined or the new area is greater than the old: 388 IF (matrix%index(slot) .EQ. 0 .OR. & 389 space .GT. matrix%index(slot + 1) - matrix%index(slot) + 1) THEN 390 IF (debug_mod) WRITE (*, *) routineP//' Slot', slot, 'not filled, adding at', & 391 matrix%index(dbcsr_slot_size) + 1, 'sized', space 392 matrix%index(slot) = matrix%index(dbcsr_slot_size) + 1 393 matrix%index(slot + 1) = matrix%index(slot) + space - 1 394 matrix%index(dbcsr_slot_size) = matrix%index(slot + 1) 395 ENDIF 396 ! Shorten an index entry. 397 IF (space .LT. matrix%index(slot + 1) - matrix%index(slot) + 1) THEN 398 IF (debug_mod) WRITE (*, *) routineP//' Shortening index' 399 matrix%index(slot + 1) = matrix%index(slot) + space - 1 400 CALL dbcsr_repoint_index(matrix, slot) 401 ENDIF 402 ub_new = matrix%index(slot + 1) + deplus 403 IF (debug_mod) WRITE (*, *) routineP//' need', space, 'at', matrix%index(slot), & 404 'to', matrix%index(slot + 1), '(', ub_new, ')', 'have', ub 405 IF (ub_new .GT. ub) THEN 406 IF (debug_mod) WRITE (*, *) routineP//' Reallocating index to ubound', ub_new 407 !CALL reallocate(matrix%index, 1, ub_new) 408 CALL ensure_array_size(matrix%index, lb=1, ub=ub_new, & 409 factor=default_resize_factor, & 410 nocopy=.FALSE., & 411 memory_type=matrix%index_memory_type) 412 CALL dbcsr_repoint_index(matrix) 413 ENDIF 414 IF (debug_mod) WRITE (*, *) routineP//' Adding slot', slot, 'at', & 415 matrix%index(slot), 'sized', space 416 CALL dbcsr_repoint_index(matrix, slot) 417 IF (PRESENT(DATA)) & 418 matrix%index(matrix%index(slot):matrix%index(slot + 1)) = DATA(:) 419 END SUBROUTINE dbcsr_addto_index_array 420 421 SUBROUTINE dbcsr_clearfrom_index_array(matrix, slot) 422 !! Removes data from the index. 423 424 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 425 !! bcsr matrix 426 INTEGER, INTENT(IN) :: slot 427 !! which index array to remove (e.g., dbcsr_slot_row_blk_sizes) 428 429 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_clearfrom_index_array', & 430 routineP = moduleN//':'//routineN 431 432 INTEGER :: space 433 INTEGER, DIMENSION(5) :: max_extents 434 435! --------------------------------------------------------------------------- 436 437 IF (.NOT. ASSOCIATED(matrix%index)) & 438 DBCSR_ABORT("Index must be preallocated.") 439 IF (UBOUND(matrix%index, 1) < dbcsr_num_slots) & 440 DBCSR_ABORT("Actual index size less than declared size") 441 IF (debug_mod) WRITE (*, *) routineP//' index', & 442 matrix%index(:dbcsr_num_slots) 443 ! Clear index entry pointer 444 matrix%index(slot) = 1 445 matrix%index(slot + 1) = 0 446 CALL dbcsr_repoint_index(matrix, slot) 447 ! Update the declared index size 448 max_extents = (/ & 449 matrix%index(dbcsr_slot_row_p + 1), & 450 matrix%index(dbcsr_slot_col_i + 1), & 451 matrix%index(dbcsr_slot_blk_p + 1), & 452 matrix%index(dbcsr_slot_thr_c + 1), & 453 matrix%index(dbcsr_slot_coo_l + 1)/) 454 space = MAX(MAXVAL(max_extents), dbcsr_num_slots) 455 matrix%index(dbcsr_slot_size) = space 456 END SUBROUTINE dbcsr_clearfrom_index_array 457 458 SUBROUTINE dbcsr_repoint_index(m, slot) 459 !! Updates the index pointers of a bcsr matrix 460 461 TYPE(dbcsr_type), INTENT(INOUT) :: m 462 !! matrix for which index pointers are updated 463 INTEGER, INTENT(IN), OPTIONAL :: slot 464 !! only repoint this index 465 466 INTEGER :: s 467 LOGICAL :: all 468 469! --------------------------------------------------------------------------- 470 471 IF (m%nblks .NE. m%index(dbcsr_slot_nblks)) THEN 472 m%nblks = m%index(dbcsr_slot_nblks) 473 m%nze = m%index(dbcsr_slot_nze) 474 ENDIF 475 all = .TRUE. 476 IF (PRESENT(slot)) THEN 477 all = .FALSE. 478 s = slot 479 ELSE 480 s = 0 481 ENDIF 482 483 IF (m%index(dbcsr_slot_row_p) .GT. 0 .AND. all .OR. & 484 s .EQ. dbcsr_slot_row_p) THEN 485 IF (m%index(dbcsr_slot_row_p) .GT. 0) THEN 486 m%row_p => m%index(m%index(dbcsr_slot_row_p): & 487 m%index(dbcsr_slot_row_p + 1)) 488 ELSE 489 NULLIFY (m%row_p) 490 ENDIF 491 ENDIF 492 IF (m%index(dbcsr_slot_col_i) .GT. 0 .AND. all .OR. & 493 s .EQ. dbcsr_slot_col_i) THEN 494 IF (m%index(dbcsr_slot_col_i) .GT. 0) THEN 495 m%col_i => m%index(m%index(dbcsr_slot_col_i): & 496 m%index(dbcsr_slot_col_i + 1)) 497 ELSE 498 NULLIFY (m%col_i) 499 ENDIF 500 ENDIF 501 IF (m%index(dbcsr_slot_blk_p) .GT. 0 .AND. all .OR. & 502 s .EQ. dbcsr_slot_blk_p) THEN 503 IF (m%index(dbcsr_slot_blk_p) .GT. 0) THEN 504 m%blk_p => m%index(m%index(dbcsr_slot_blk_p): & 505 m%index(dbcsr_slot_blk_p + 1)) 506 ELSE 507 NULLIFY (m%blk_p) 508 ENDIF 509 ENDIF 510 IF (m%index(dbcsr_slot_thr_c) .GT. 0 .AND. all .OR. & 511 s .EQ. dbcsr_slot_thr_c) THEN 512 IF (m%index(dbcsr_slot_thr_c) .GT. 0) THEN 513 m%thr_c => m%index(m%index(dbcsr_slot_thr_c): & 514 m%index(dbcsr_slot_thr_c + 1)) 515 ELSE 516 NULLIFY (m%thr_c) 517 ENDIF 518 ENDIF 519 IF (m%index(dbcsr_slot_coo_l) .GT. 0 .AND. all .OR. & 520 s .EQ. dbcsr_slot_coo_l) THEN 521 IF (m%index(dbcsr_slot_coo_l) .GT. 0) THEN 522 m%coo_l => m%index(m%index(dbcsr_slot_coo_l): & 523 m%index(dbcsr_slot_coo_l + 1)) 524 ELSE 525 NULLIFY (m%coo_l) 526 ENDIF 527 ENDIF 528 IF (all) THEN 529 m%index(dbcsr_slot_nblks) = m%nblks 530 m%index(dbcsr_slot_nze) = m%nze 531 ENDIF 532 END SUBROUTINE dbcsr_repoint_index 533 534 SUBROUTINE dbcsr_make_index_exist(m) 535 536 TYPE(dbcsr_type), INTENT(INOUT) :: m 537 !! Create index for this matrix 538 539 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_make_index_exist' 540 541 INTEGER :: error_handle 542 543! --------------------------------------------------------------------------- 544 545 CALL timeset(routineN, error_handle) 546 IF (.NOT. ASSOCIATED(m%index)) & 547 DBCSR_ABORT("Index array does not yet exist.") 548 IF (.NOT. ASSOCIATED(m%row_p)) THEN 549 CALL dbcsr_addto_index_array(m, dbcsr_slot_row_p, & 550 reservation=m%nblkrows_total + 1) 551 m%row_p(:) = 0 552 ENDIF 553 IF (.NOT. ASSOCIATED(m%col_i)) THEN 554 CALL dbcsr_addto_index_array(m, dbcsr_slot_col_i, & 555 reservation=0) 556 ENDIF 557 IF (.NOT. ASSOCIATED(m%blk_p)) THEN 558 CALL dbcsr_addto_index_array(m, dbcsr_slot_blk_p, & 559 reservation=0) 560 ENDIF 561 CALL dbcsr_repoint_index(m) 562 CALL timestop(error_handle) 563 END SUBROUTINE dbcsr_make_index_exist 564 565 SUBROUTINE dbcsr_sort_indices(n, row_i, col_i, blk_p, blk_d) 566 !! Sorts the rows & columns of a work matrix 567 !! 568 !! Description 569 !! Sorts the row and column indices so that the rows monotonically 570 !! increase and the columns monotonically increase within each row. 571 !! Passing the blk_p array rearranges the block pointers accordingly. 572 !! This must be done if they are pointing to valid data, otherwise 573 !! they become invalid. 574 575 INTEGER, INTENT(IN) :: n 576 !! number of blocks (elements) to sort 577 INTEGER, DIMENSION(1:), INTENT(INOUT) :: row_i, col_i 578 !! row indices 579 !! column indices 580 INTEGER, DIMENSION(1:), INTENT(INOUT), OPTIONAL :: blk_p, blk_d 581 !! block pointers 582 !! data storage 583 584 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_sort_indices', & 585 routineP = moduleN//':'//routineN 586 INTEGER(KIND=int_8), PARAMETER :: lmask8 = 4294967295_int_8 587 588 INTEGER :: error_handle, i 589 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: sort_keys 590 INTEGER, ALLOCATABLE, DIMENSION(:) :: buf, buf_d 591 592! --------------------------------------------------------------------------- 593 594 IF (SIZE(row_i) .EQ. 0) RETURN 595 596 CALL timeset(routineN, error_handle) 597 598 IF (SIZE(row_i) < n) DBCSR_ABORT('row_i too small') 599 IF (SIZE(col_i) < n) DBCSR_ABORT('col_i too small') 600 IF (PRESENT(blk_p)) THEN 601 IF (SIZE(blk_p) < n) DBCSR_ABORT('blk_p too small') 602 ALLOCATE (buf(n)) 603 buf(1:n) = blk_p(1:n) 604 ENDIF 605 IF (PRESENT(blk_d)) THEN 606 ALLOCATE (buf_d(n)) 607 buf_d(1:n) = blk_d(1:n) 608 ENDIF 609 ! Create an ordering for both rows and columns. If the blk_p must 610 ! be rearranged, then the col_i array will be used as a 611 ! permutation vector. 612 ALLOCATE (sort_keys(n)) 613 sort_keys(:) = IOR(ISHFT(INT(row_i(1:n), int_8), 32), INT(col_i(1:n), int_8)) 614 IF (PRESENT(blk_p)) col_i(1:n) = (/(i, i=1, n)/) 615 ! Now do a nice quicksort. 616 CALL sort(sort_keys, n, col_i) 617 ! Since blk_d is usually not present we can have two loops that 618 ! are essentially the same. 619 IF (PRESENT(blk_p)) THEN 620 DO i = 1, n 621 blk_p(i) = buf(col_i(i)) 622 ENDDO 623 DEALLOCATE (buf) 624 END IF 625 IF (PRESENT(blk_d)) THEN 626 DO i = 1, n 627 blk_d(i) = buf_d(col_i(i)) 628 ENDDO 629 DEALLOCATE (buf_d) 630 ENDIF 631 DO i = 1, n 632 col_i(i) = INT(IAND(sort_keys(i), lmask8), int_4) 633 row_i(i) = INT(ISHFT(sort_keys(i), -32), int_4) 634 ENDDO 635 DEALLOCATE (sort_keys) 636 IF (debug_mod .AND. PRESENT(blk_p)) & 637 WRITE (*, *) routineP//' sort, blk_p =', blk_p 638 CALL timestop(error_handle) 639 640 END SUBROUTINE dbcsr_sort_indices 641 642 SUBROUTINE dbcsr_index_prune_deleted(matrix) 643 !! Removes the deleted blocks from the index. 644 !! 645 !! Description 646 647 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 648 !! Prune the index of this matrix. 649 650 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_index_prune_deleted' 651 652 INTEGER :: error_handle, nblks_max, new_blk, nrows, & 653 old_blk, row 654 INTEGER, ALLOCATABLE, DIMENSION(:) :: new_blk_p, new_col_i, new_row_count 655 INTEGER, DIMENSION(:), POINTER :: old_blk_p, old_col_i, old_row_p 656 657! --------------------------------------------------------------------------- 658 659 CALL timeset(routineN, error_handle) 660 ! 661 old_row_p => matrix%row_p 662 old_col_i => matrix%col_i 663 old_blk_p => matrix%blk_p 664 ! 665 nrows = matrix%nblkrows_total 666 nblks_max = old_row_p(nrows + 1) 667 ALLOCATE (new_row_count(nrows)) 668 ALLOCATE (new_col_i(nblks_max)) 669 ALLOCATE (new_blk_p(nblks_max)) 670 ! 671 ! Build up the new index from all non-deleted blocks in the 672 ! existing index. 673 new_blk = 0 674 DO row = 1, nrows 675 new_row_count(row) = 0 676 DO old_blk = old_row_p(row) + 1, old_row_p(row + 1) 677 IF (old_blk_p(old_blk) .GT. 0) THEN 678 new_blk = new_blk + 1 679 new_row_count(row) = new_row_count(row) + 1 680 new_col_i(new_blk) = old_col_i(old_blk) 681 new_blk_p(new_blk) = old_blk_p(old_blk) 682 ENDIF 683 ENDDO 684 ENDDO 685 ! 686 CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_row_p) 687 CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_col_i) 688 CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_blk_p) 689 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_row_p, & 690 reservation=nrows + 1, extra=2*new_blk) 691 old_row_p => matrix%row_p 692 CALL dbcsr_build_row_index(counts=new_row_count, rows=old_row_p, & 693 nrows=nrows) 694 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_col_i, DATA=new_col_i(1:new_blk)) 695 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_blk_p, DATA=new_blk_p(1:new_blk)) 696 matrix%nblks = new_blk 697 matrix%index(dbcsr_slot_nblks) = new_blk 698 ! 699 DEALLOCATE (new_col_i, new_blk_p, new_row_count) 700 ! 701 CALL timestop(error_handle) 702 END SUBROUTINE dbcsr_index_prune_deleted 703 704 SUBROUTINE transpose_index_local(new_col_p, new_row_i, old_row_p, & 705 old_col_i, new_blk_p, old_blk_p) 706 !! Re-indexes row_p and blk_i according to columns. 707 !! 708 !! The re-indexing is equivalent to a local-only transpose. 709 710 INTEGER, DIMENSION(:), INTENT(OUT) :: new_col_p, new_row_i 711 !! new column pointer 712 !! new row index 713 INTEGER, DIMENSION(:), INTENT(IN) :: old_row_p, old_col_i 714 !! old row pointer 715 !! old column index 716 INTEGER, DIMENSION(:), INTENT(OUT), OPTIONAL :: new_blk_p 717 !! new block pointer 718 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: old_blk_p 719 !! old block pointer 720 721 CHARACTER(len=*), PARAMETER :: routineN = 'transpose_index_local' 722 723 INTEGER :: error_handle, nblks, ncols_new, nrows_old 724 INTEGER, ALLOCATABLE, DIMENSION(:) :: new_col_i 725 LOGICAL :: blks 726 727! --------------------------------------------------------------------------- 728 729 CALL timeset(routineN, error_handle) 730 blks = PRESENT(new_blk_p) .AND. PRESENT(old_blk_p) 731 nblks = SIZE(old_col_i) 732 nrows_old = SIZE(old_row_p) - 1 733 ncols_new = SIZE(new_col_p) - 1 734 IF (blks) new_blk_p(:) = old_blk_p(:) 735 ALLOCATE (new_col_i(nblks)) 736 CALL dbcsr_expand_row_index(old_row_p, new_row_i, nrows_old, nblks) 737 new_col_i(:) = old_col_i(:) 738 CALL dbcsr_sort_indices(nblks, new_col_i, new_row_i, new_blk_p) 739 CALL dbcsr_make_dbcsr_index(new_col_p, new_col_i, ncols_new, nblks) 740 DEALLOCATE (new_col_i) 741 CALL timestop(error_handle) 742 END SUBROUTINE transpose_index_local 743 744 SUBROUTINE dbcsr_make_index_canonical(matrix, cp2k) 745 !! Makes a canonical index to the distribution. 746 !! 747 !! Canonical means that it respects the distribution. 748 749 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 750 !! matrix for which to make canonical index 751 LOGICAL, INTENT(IN), OPTIONAL :: cp2k 752 !! make CP2K triangular index from canonical; default is false 753 754 INTEGER :: nb, nc, nr 755 INTEGER, ALLOCATABLE, DIMENSION(:) :: new_blk_p, new_col_i, new_row_p 756 LOGICAL :: rev 757 758! --------------------------------------------------------------------------- 759 760 rev = .FALSE. 761 IF (PRESENT(cp2k)) rev = cp2k 762 nr = SIZE(matrix%row_p) 763 ALLOCATE (new_row_p(nr)) 764 nc = SIZE(matrix%col_i) 765 ALLOCATE (new_col_i(nc)) 766 nb = SIZE(matrix%blk_p) 767 ALLOCATE (new_blk_p(nb)) 768 IF (rev) THEN 769 CALL make_index_triangular(new_row_p, new_col_i, new_blk_p, & 770 matrix%row_p, matrix%col_i, matrix%blk_p, matrix) 771 ELSE 772 CALL make_index_canonical(new_row_p, new_col_i, new_blk_p, & 773 matrix%row_p, matrix%col_i, matrix%blk_p, matrix) 774 ENDIF 775 matrix%row_p(:) = new_row_p 776 matrix%col_i(:) = new_col_i 777 matrix%blk_p(:) = new_blk_p 778 END SUBROUTINE dbcsr_make_index_canonical 779 780 SUBROUTINE make_dense_index(row_p, col_i, blk_p, & 781 nblkrows_total, nblkcols_total, myblkrows, myblkcols, & 782 row_blk_offsets, col_blk_offsets, meta, make_tr) 783 !! Makes the index for a dense matrix 784 !! @note Used for making matrices dense/undense 785 786 !INTEGER, DIMENSION(:), INTENT(OUT) :: row_p, col_i, blk_p 787 INTEGER, INTENT(IN) :: nblkrows_total 788 !! Total blocked rows 789 INTEGER, DIMENSION(:), INTENT(OUT) :: blk_p, col_i 790 !! Storage for new index 791 !! Storage for new index 792 INTEGER, DIMENSION(1:nblkrows_total + 1), & 793 INTENT(OUT) :: row_p 794 !! Storage for new index 795 INTEGER, INTENT(IN) :: nblkcols_total 796 !! Total blocked columns 797 INTEGER, DIMENSION(:), INTENT(IN) :: myblkrows, myblkcols, row_blk_offsets, & 798 col_blk_offsets 799 !! List of blocked rows in my process row 800 !! List of blocked columns in my process column 801 INTEGER, DIMENSION(dbcsr_meta_size), INTENT(INOUT) :: meta 802 !! Metadata updates for new index 803 LOGICAL, INTENT(IN), OPTIONAL :: make_tr 804 !! Dense blocks are transposed 805 806 CHARACTER(len=*), PARAMETER :: routineN = 'make_dense_index' 807 808 INTEGER :: blk, c, col_l, mynblkcols, mynblkrows, & 809 nblks, nze, prev_row, row, row_l, & 810 sign_carrier, sz 811 812! --------------------------------------------------------------------------- 813 814 sign_carrier = 1 815 IF (PRESENT(make_tr)) THEN 816 IF (make_tr) sign_carrier = -1 817 ENDIF 818 mynblkrows = SIZE(myblkrows) 819 mynblkcols = SIZE(myblkcols) 820 meta(dbcsr_slot_nblkrows_local) = mynblkrows 821 meta(dbcsr_slot_nblkcols_local) = mynblkcols 822 nblks = mynblkrows*mynblkcols 823 nze = 1 824 IF (nblks .EQ. 0) THEN 825 row_p(1:) = 0 826 ELSE 827 row_p(1) = 0 828 !row_p(nrows+1) = nblks 829 prev_row = 1 830 blk = 0 831 DO row_l = 1, mynblkrows 832 row = myblkrows(row_l) 833 row_p(prev_row + 1:row) = blk 834 DO col_l = 1, mynblkcols 835 c = myblkcols(col_l) 836 col_i(blk + col_l) = c 837 sz = (row_blk_offsets(row + 1) - row_blk_offsets(row))* & 838 (col_blk_offsets(c + 1) - col_blk_offsets(c)) 839 IF (sz .GT. 0) THEN 840 blk_p(blk + col_l) = SIGN(nze, sign_carrier) 841 nze = nze + sz 842 ELSE 843 blk_p(blk + col_l) = 0 844 ENDIF 845 ENDDO 846 prev_row = row 847 blk = blk + mynblkcols 848 END DO 849 IF (blk /= nblks) DBCSR_ABORT("Block mismatch") 850 row_p(prev_row + 1:nblkrows_total + 1) = nblks 851 ENDIF 852 IF (debug_mod) THEN 853 WRITE (*, *) routineN//" new index" 854 WRITE (*, *) "row_p=", row_p 855 WRITE (*, *) "col_i=", col_i 856 WRITE (*, *) "blk_p=", blk_p 857 ENDIF 858 meta(dbcsr_slot_nblkrows_total) = nblkrows_total 859 meta(dbcsr_slot_nblkcols_total) = nblkcols_total 860 END SUBROUTINE make_dense_index 861 862 SUBROUTINE make_undense_index( & 863 row_p, col_i, blk_p, & 864 distribution, local_row_offsets, local_col_offsets, & 865 meta) 866 !! Makes a blocked index from a dense matrix 867 !! @note Used for making matrices dense/undense 868 869 INTEGER, DIMENSION(:), INTENT(OUT) :: row_p, col_i, blk_p 870 !! Storage for new index 871 !! Storage for new index 872 !! Storage for new index 873 TYPE(dbcsr_distribution_obj) :: distribution 874 !! Blocked distribution 875 INTEGER, DIMENSION(:), INTENT(IN) :: local_row_offsets, local_col_offsets 876 INTEGER, DIMENSION(dbcsr_meta_size), INTENT(INOUT) :: meta 877 !! Metadata updates for new index 878 879 INTEGER :: col, lr, lrow, nblkcols_local, & 880 nblkrows_local, nblkrows_total, & 881 nfullcols_local, prev_row, row 882 INTEGER, DIMENSION(:), POINTER :: local_cols, local_rows 883 884! --------------------------------------------------------------------------- 885 886 local_cols => dbcsr_distribution_local_cols(distribution) 887 local_rows => dbcsr_distribution_local_rows(distribution) 888 meta(dbcsr_slot_nblkrows_total) = dbcsr_distribution_nrows(distribution) 889 meta(dbcsr_slot_nblkcols_total) = dbcsr_distribution_ncols(distribution) 890 meta(dbcsr_slot_nblkrows_local) = dbcsr_distribution_nlocal_rows(distribution) 891 meta(dbcsr_slot_nblkcols_local) = dbcsr_distribution_nlocal_cols(distribution) 892 nblkrows_total = meta(dbcsr_slot_nblkrows_total) 893 nblkcols_local = meta(dbcsr_slot_nblkcols_local) 894 nblkrows_local = meta(dbcsr_slot_nblkrows_local) 895 nfullcols_local = meta(dbcsr_slot_nfullcols_local) 896 ! Fill the row_p array. 897 lr = 0 898 row_p(1) = 0 899 prev_row = 1 900 DO lrow = 1, nblkrows_local 901 row = local_rows(lrow) 902 row_p(prev_row + 1:row) = lr 903 lr = lr + nblkcols_local 904 row_p(row + 1) = lr 905 prev_row = row 906 ENDDO 907 row_p(prev_row + 1:nblkrows_total + 1) = lr 908 ! 909 DO row = 1, nblkrows_local 910 DO col = 1, nblkcols_local 911 col_i(nblkcols_local*(row - 1) + col) = local_cols(col) 912 blk_p(nblkcols_local*(row - 1) + col) = 1 + & 913 (local_row_offsets(row) - 1)*nfullcols_local & 914 + (local_col_offsets(col) - 1)* & 915 (local_row_offsets(row + 1) - local_row_offsets(row)) 916 END DO 917 END DO 918 END SUBROUTINE make_undense_index 919 920 SUBROUTINE merge_index_arrays(new_row_i, new_col_i, new_blk_p, new_size, & 921 old_row_i, old_col_i, old_blk_p, old_size, & 922 add_ip, add_size, new_blk_d, old_blk_d, & 923 added_size_offset, added_sizes, added_size, added_nblks) 924 !! Merges two indices 925 !! 926 !! Added sizes 927 !! added_size_offset and added_sizes can be optionally 928 !! specified. This is meant for cases where the added blocks may 929 !! be duplicates of existing blocks. In this way it is possible 930 !! to recalculate new block pointers to avoid wasted space. 931 !! @note Used in local multiply 932 !! Assumes they are both pre-sorted 933 934 INTEGER, INTENT(IN) :: new_size 935 !! size of merged index 936 INTEGER, DIMENSION(new_size), INTENT(OUT) :: new_blk_p, new_col_i, new_row_i 937 !! merged result 938 !! merged result 939 !! merged result 940 INTEGER, INTENT(IN) :: old_size 941 !! size of current index 942 INTEGER, DIMENSION(old_size), INTENT(IN) :: old_blk_p, old_col_i, old_row_i 943 !! current index 944 !! current index 945 !! current index 946 INTEGER, INTENT(IN) :: add_size 947 !! size of index to add into the current index 948 INTEGER, DIMENSION(3, add_size), INTENT(IN) :: add_ip 949 !! index to add into the current index 950 INTEGER, DIMENSION(new_size), INTENT(OUT), & 951 OPTIONAL :: new_blk_d 952 INTEGER, DIMENSION(old_size), INTENT(IN), OPTIONAL :: old_blk_d 953 INTEGER, INTENT(IN), OPTIONAL :: added_size_offset 954 !! specify base of added sizes 955 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: added_sizes 956 !! specify sizes of added blocks 957 INTEGER, INTENT(OUT), OPTIONAL :: added_size, added_nblks 958 !! counts number of sizes of added blocks 959 !! actual number of new elements 960 961 INTEGER :: add_blk, bp, i, merge_from_whom, & 962 new_blk, old_blk 963 LOGICAL :: multidata 964 965! --------------------------------------------------------------------------- 966 967 bp = 0 968 multidata = PRESENT(old_blk_d) .AND. PRESENT(new_blk_d) 969 IF (old_size + add_size .NE. new_size) & 970 DBCSR_WARN("Mismatch of new and old size") 971 IF (PRESENT(added_size_offset) .NEQV. PRESENT(added_sizes)) & 972 DBCSR_ABORT("Must specify a set of arguments") 973 IF (PRESENT(added_sizes) .NEQV. PRESENT(added_size)) & 974 DBCSR_ABORT("Must specify a set of arguments") 975 IF (debug_mod) THEN 976 WRITE (*, *) " Old array", old_size 977 DO i = 1, old_size 978 WRITE (*, '(I7,2X,I7,2X,I7)') old_row_i(i), old_col_i(i), old_blk_p(i) 979 ENDDO 980 WRITE (*, *) " Add array", add_size 981 DO i = 1, add_size 982 WRITE (*, '(I7,2X,I7,2X,I7)') add_ip(1:3, i) 983 ENDDO 984 ENDIF 985 IF (PRESENT(added_nblks)) added_nblks = 0 986 IF (PRESENT(added_size)) THEN 987 added_size = 0 988 bp = added_size_offset 989 ENDIF 990 IF (add_size .GT. 0) THEN 991 old_blk = 1 992 add_blk = 1 993 new_blk = 1 994 IF (old_size .EQ. 0) THEN 995 new_row_i(1:add_size) = add_ip(1, 1:add_size) 996 new_col_i(1:add_size) = add_ip(2, 1:add_size) 997 new_blk_p(1:add_size) = add_ip(3, 1:add_size) 998 !IF (multidata) new_blk_d(1:add_size) = add_ip(4, 1:add_size) 999 IF (PRESENT(added_nblks)) added_nblks = add_size 1000 IF (PRESENT(added_size)) added_size = SUM(added_sizes) 1001 ELSE 1002 DO WHILE (new_blk .LE. new_size) 1003 merge_from_whom = 0 1004 IF (old_blk .LE. old_size .AND. add_blk .LE. add_size) THEN 1005 IF (add_ip(1, add_blk) .EQ. old_row_i(old_blk) & 1006 .AND. add_ip(2, add_blk) .EQ. old_col_i(old_blk)) THEN 1007 IF (debug_mod) THEN 1008 WRITE (*, *) "Duplicate block! addblk", & 1009 add_blk, "oldblk", old_blk 1010 ENDIF 1011 ENDIF 1012 ! Rows come first 1013 IF (add_ip(1, add_blk) .LT. old_row_i(old_blk)) THEN 1014 merge_from_whom = 2 1015 ELSEIF (add_ip(1, add_blk) .GT. old_row_i(old_blk)) THEN 1016 merge_from_whom = 1 1017 ELSE ! Same rows, so now come the columns 1018 IF (add_ip(2, add_blk) .LT. old_col_i(old_blk)) THEN 1019 ! Merges from the add array 1020 merge_from_whom = 2 1021 ELSEIF (add_ip(2, add_blk) .GT. old_col_i(old_blk)) THEN 1022 ! Merges from the old array 1023 merge_from_whom = 1 1024 ELSE 1025 ! Merge from old array and skip one in the new array 1026 IF (debug_mod) THEN 1027 WRITE (*, *) "Duplicate, keeping old", & 1028 add_ip(1, add_blk), add_ip(2, add_blk) 1029 ENDIF 1030 merge_from_whom = 1 1031 add_blk = add_blk + 1 1032 ENDIF 1033 ENDIF 1034 ELSE 1035 IF (add_blk .LE. add_size) THEN 1036 ! Merges from the add array 1037 merge_from_whom = 2 1038 ELSEIF (old_blk .LE. old_size) THEN 1039 ! Merges from the old array 1040 merge_from_whom = 1 1041 ELSE 1042 ! Hmmm, nothing to merge... 1043 merge_from_whom = 0 1044 !WRITE(*,*)"Ran out of data to merge" 1045 ENDIF 1046 ENDIF 1047 SELECT CASE (merge_from_whom) 1048 CASE (2) 1049 ! Merges from the add array 1050 new_row_i(new_blk) = add_ip(1, add_blk) 1051 new_col_i(new_blk) = add_ip(2, add_blk) 1052 new_blk_p(new_blk) = add_ip(3, add_blk) 1053 !IF (multidata) new_blk_d(new_blk) = add_ip(4, add_blk) 1054 IF (PRESENT(added_nblks)) added_nblks = added_nblks + 1 1055 IF (PRESENT(added_sizes)) THEN 1056 new_blk_p(new_blk) = bp 1057 bp = bp + added_sizes(add_blk) 1058 added_size = added_size + added_sizes(add_blk) 1059 ENDIF 1060 add_blk = add_blk + 1 1061 CASE (1) 1062 ! Merges from the old array 1063 new_row_i(new_blk) = old_row_i(old_blk) 1064 new_col_i(new_blk) = old_col_i(old_blk) 1065 new_blk_p(new_blk) = old_blk_p(old_blk) 1066 IF (multidata) new_blk_p(new_blk) = old_blk_d(old_blk) 1067 old_blk = old_blk + 1 1068 CASE DEFAULT 1069 !WRITE(*,*)"Nothing to merge" 1070 END SELECT 1071 new_blk = new_blk + 1 1072 ENDDO 1073 ENDIF 1074 ELSE 1075 new_row_i(1:old_size) = old_row_i(1:old_size) 1076 new_col_i(1:old_size) = old_col_i(1:old_size) 1077 new_blk_p(1:old_size) = old_blk_p(1:old_size) 1078 IF (multidata) new_blk_d(1:old_size) = old_blk_d(1:old_size) 1079 ENDIF 1080 IF (debug_mod) THEN 1081 WRITE (*, *) " New array" 1082 DO i = 1, new_size 1083 WRITE (*, '(4(2X,I7))') new_row_i(i), new_col_i(i), new_blk_p(i) 1084 ENDDO 1085 ENDIF 1086 END SUBROUTINE merge_index_arrays 1087 1088 SUBROUTINE dbcsr_make_index_local_row(matrix) 1089 !! Converts BCSR global row index to local row index. 1090 1091 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 1092 !! matrix for which to make canonical index 1093 1094 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_make_index_local_row' 1095 1096 INTEGER :: error_handle, lrow, nlocal_rows, & 1097 ntotal_rows, prow 1098 INTEGER, ALLOCATABLE, DIMENSION(:) :: local_row_p 1099 INTEGER, DIMENSION(:), POINTER :: local_rows 1100 1101! --------------------------------------------------------------------------- 1102 1103 CALL timeset(routineN, error_handle) 1104 IF (.NOT. ASSOCIATED(matrix%row_p)) & 1105 DBCSR_ABORT("The row index must be initialized.") 1106 IF (matrix%bcsc) & 1107 DBCSR_ABORT("Not support for BCSC yet.") 1108 ! 1109 prow = matrix%index(dbcsr_slot_home_vprow) 1110 IF (prow .LT. 0) THEN 1111 prow = matrix%index(dbcsr_slot_home_prow) 1112 ENDIF 1113 nlocal_rows = matrix%nblkrows_local 1114 ALLOCATE (local_row_p(nlocal_rows + 1)) 1115 ! The existing row_p is converted from an indexing array into a 1116 ! counting array. Because it is later discarded, the counting is 1117 ! done in-place. 1118 ntotal_rows = matrix%nblkrows_total 1119 CALL dbcsr_count_row_index(matrix%row_p, ntotal_rows) 1120 ! We first have to find the local rows for the given prow. 1121 local_rows => array_data(matrix%local_rows) 1122 IF (SIZE(local_rows) /= nlocal_rows) & 1123 DBCSR_ABORT("Mismatch in the number of local rows.") 1124 ! The counts are mapped to local rows, 1125 DO lrow = 1, nlocal_rows 1126 local_row_p(lrow) = matrix%row_p(local_rows(lrow)) 1127 ENDDO 1128 IF (SUM(matrix%row_p(1:ntotal_rows)) /= SUM(local_row_p(1:nlocal_rows))) & 1129 DBCSR_ABORT("Inconsistent row counts. Perhaps non-local rows contain data?.") 1130 ! then converted into an index. 1131 CALL dbcsr_build_row_index(local_row_p, nlocal_rows) 1132 ! The local row index replaces the global one. 1133 CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_row_p) 1134 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_row_p, DATA=local_row_p) 1135 ! Finally the matrix is designated as having a local-based index. 1136 matrix%local_indexing = .TRUE. 1137 IF (careful_mod) THEN 1138 IF (array_size(matrix%local_rows) /= nlocal_rows) & 1139 DBCSR_ABORT("Inconsistent local row counts.") 1140 IF (array_size(matrix%global_rows) /= ntotal_rows) & 1141 DBCSR_ABORT("Inconsistent global row counts.") 1142 IF (array_size(matrix%global_rows) .EQ. 0) THEN 1143 IF (nlocal_rows /= 0) & 1144 DBCSR_ABORT("Invalid number of local or global rows.") 1145 ENDIF 1146 ENDIF 1147 DEALLOCATE (local_row_p) 1148 ! 1149 CALL timestop(error_handle) 1150 END SUBROUTINE dbcsr_make_index_local_row 1151 1152 SUBROUTINE dbcsr_make_index_list(matrix, thread_redist) 1153 !! Converts BCSR index into list indices (similar to work matrices) 1154 1155 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 1156 !! matrix for which to make canonical index 1157 LOGICAL, INTENT(IN) :: thread_redist 1158 !! make a thread subdistribution 1159 1160 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_make_index_list' 1161 LOGICAL, PARAMETER :: dbg = .FALSE. 1162 1163 INTEGER :: blk, error_handle, ithread, my_cnt, & 1164 nblks, nrows, nthreads 1165 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: blki 1166 INTEGER, DIMENSION(0) :: zero_len_array 1167 INTEGER, DIMENSION(:), POINTER :: global_cols, local_rows, td, thr_c 1168 1169! --------------------------------------------------------------------------- 1170 1171 CALL timeset(routineN, error_handle) 1172 IF (.NOT. ASSOCIATED(matrix%row_p)) & 1173 DBCSR_ABORT("The row index must be initialized.") 1174 IF (matrix%list_indexing) & 1175 DBCSR_ABORT("List index already exists?") 1176 IF (matrix%bcsc) & 1177 DBCSR_ABORT("Not support for BCSC yet.") 1178 IF (matrix%nblks /= SIZE(matrix%col_i)) & 1179 DBCSR_ABORT("Block count mismatch.") 1180 IF (matrix%nblks /= SIZE(matrix%blk_p)) & 1181 DBCSR_ABORT("Block count mismatch") 1182 ! 1183 IF (matrix%local_indexing) THEN 1184 IF (SIZE(matrix%row_p) - 1 /= matrix%nblkrows_local) & 1185 DBCSR_ABORT("Local row index incorrectly sized.") 1186 ELSE 1187 IF (SIZE(matrix%row_p) - 1 /= matrix%nblkrows_total) & 1188 DBCSR_ABORT("Global row index incorrectly sized") 1189 ENDIF 1190 ! 1191 matrix%list_indexing = .TRUE. 1192 ! 1193 IF (matrix%local_indexing) THEN 1194 nrows = matrix%nblkrows_local 1195 ELSE 1196 nrows = matrix%nblkrows_total 1197 ENDIF 1198 ! 1199 nblks = matrix%nblks 1200 ALLOCATE (blki(3, nblks)) 1201 CALL dbcsr_expand_row_index_2d(matrix%row_p, blki, nrows, 1) 1202 IF (matrix%local_indexing) THEN 1203 global_cols => array_data(matrix%global_cols) 1204 ! If local indexing is enabled, then the rows but not the 1205 ! columns are already localized 1206 IF (dbg) THEN 1207 WRITE (*, *) routineN//" Making local columns" 1208 WRITE (*, '(10(1X,i7))') global_cols 1209 WRITE (*, *) 'local' 1210 WRITE (*, '(10(1X,i7))') array_data(matrix%local_cols) 1211 ENDIF 1212 DO blk = 1, nblks 1213 blki(2, blk) = global_cols(matrix%col_i(blk)) 1214 blki(3, blk) = matrix%blk_p(blk) 1215 ENDDO 1216 ELSE 1217 IF (dbg) WRITE (*, *) routineN//" Not making local columns" 1218 DO blk = 1, nblks 1219 blki(2, blk) = matrix%col_i(blk) 1220 blki(3, blk) = matrix%blk_p(blk) 1221 END DO 1222 ENDIF 1223 CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_row_p) 1224 CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_col_i) 1225 CALL dbcsr_clearfrom_index_array(matrix, dbcsr_slot_blk_p) 1226 nthreads = 0 1227!$ nthreads = OMP_GET_MAX_THREADS() 1228 ! 1229 ! Reshuffle according to threads 1230 IF (nthreads .GT. 0 .AND. thread_redist) THEN 1231 td => array_data(dbcsr_distribution_thread_dist(dbcsr_distribution(matrix))) 1232 IF (matrix%local_indexing) THEN 1233 local_rows => array_data(matrix%local_rows) 1234 ENDIF 1235!$OMP PARALLEL DEFAULT (NONE) & 1236!$OMP PRIVATE (my_cnt, ithread, blk) & 1237!$OMP SHARED (td, blki, nthreads, thr_c, nblks, matrix,local_rows) 1238 ! 1239 ithread = 0 1240!$ ithread = OMP_GET_THREAD_NUM() 1241!$OMP MASTER 1242!$ nthreads = OMP_GET_NUM_THREADS() 1243 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_thr_c, & 1244 reservation=nthreads + 1, extra=3*nblks) 1245 thr_c => matrix%thr_c 1246 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_coo_l, & 1247 reservation=3*nblks) 1248!$OMP END MASTER 1249!$OMP BARRIER 1250 my_cnt = 0 1251 IF (matrix%local_indexing) THEN 1252 my_cnt = COUNT(td(local_rows(blki(1, :))) .EQ. ithread) 1253 ELSE 1254 my_cnt = COUNT(td(blki(1, :)) .EQ. ithread) 1255 ENDIF 1256 !DO blk = 1, nblks 1257 ! IF (td(blki(1, blk)) .EQ. ithread) my_cnt = my_cnt+1 1258 !ENDDO 1259 thr_c(ithread + 1) = my_cnt 1260!$OMP BARRIER 1261!$OMP MASTER 1262 CALL dbcsr_build_row_index_inplace(thr_c, nthreads) 1263!$OMP END MASTER 1264!$OMP BARRIER 1265 my_cnt = (thr_c(ithread + 1) + 1)*3 - 2 1266 IF (matrix%local_indexing) THEN 1267 DO blk = 1, nblks 1268 IF (td(local_rows(blki(1, blk))) .EQ. ithread) THEN 1269 matrix%coo_l(my_cnt:my_cnt + 2) = blki(1:3, blk) 1270 my_cnt = my_cnt + 3 1271 ENDIF 1272 ENDDO 1273 ELSE 1274 DO blk = 1, nblks 1275 IF (td(blki(1, blk)) .EQ. ithread) THEN 1276 matrix%coo_l(my_cnt:my_cnt + 2) = blki(1:3, blk) 1277 my_cnt = my_cnt + 3 1278 ENDIF 1279 ENDDO 1280 ENDIF 1281!$OMP END PARALLEL 1282 ELSE 1283 ! Small price to pay for avoiding infinite recursions. 1284 DO blk = 2, nblks 1285 IF (blki(1, blk) .EQ. blki(1, blk - 1) .AND. blki(2, blk) .EQ. blki(2, blk - 1)) THEN 1286 ! Weird assertion to give some idea of the two blocks. 1287 IF (-blki(1, blk) /= blki(2, blk)) & 1288 CALL dbcsr_abort(__LOCATION__, & 1289 "Should not have duplicate matrix blocks. (-row, col) is duplicated.") 1290 ENDIF 1291 ENDDO 1292 ! 1293 IF (nblks > 0) THEN 1294 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_coo_l, & 1295 DATA=RESHAPE(blki, (/3*nblks/))) 1296 ELSE 1297 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_coo_l, & 1298 DATA=zero_len_array) 1299 ENDIF 1300 ENDIF 1301 1302 DEALLOCATE (blki) 1303 ! 1304 CALL timestop(error_handle) 1305 END SUBROUTINE dbcsr_make_index_list 1306 1307 SUBROUTINE dbcsr_index_compact(matrix) 1308 !! Compacts an index. 1309 1310 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 1311 !! matrix for which to make canonical index 1312 1313 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_index_compact' 1314 1315 INTEGER :: error_handle, new_size, size_blk_p, & 1316 size_col_i, size_coo_l, size_row_p, & 1317 size_thr_c 1318 INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_p, col_i, coo_l, meta, row_p, thr_c 1319 LOGICAL :: compact, has_blk_p, has_col_i, & 1320 has_coo_l, has_row_p, has_thr_c 1321 1322! --------------------------------------------------------------------------- 1323 1324 CALL timeset(routineN, error_handle) 1325 ! Ensures the index pointers are set. 1326 CALL dbcsr_repoint_index(matrix) 1327 ! Check that compaction is even needed. 1328 has_row_p = ASSOCIATED(matrix%row_p) 1329 IF (has_row_p) THEN 1330 size_row_p = SIZE(matrix%row_p) 1331 ELSE 1332 size_row_p = 0 1333 ENDIF 1334 has_col_i = ASSOCIATED(matrix%col_i) 1335 IF (has_col_i) THEN 1336 size_col_i = SIZE(matrix%col_i) 1337 ELSE 1338 size_col_i = 0 1339 ENDIF 1340 has_blk_p = ASSOCIATED(matrix%blk_p) 1341 IF (has_blk_p) THEN 1342 size_blk_p = SIZE(matrix%blk_p) 1343 ELSE 1344 size_blk_p = 0 1345 ENDIF 1346 has_thr_c = ASSOCIATED(matrix%thr_c) 1347 IF (has_thr_c) THEN 1348 size_thr_c = SIZE(matrix%thr_c) 1349 ELSE 1350 size_thr_c = 0 1351 ENDIF 1352 has_coo_l = ASSOCIATED(matrix%coo_l) 1353 IF (has_coo_l) THEN 1354 size_coo_l = SIZE(matrix%coo_l) 1355 ELSE 1356 size_coo_l = 0 1357 ENDIF 1358 ! 1359 new_size = dbcsr_num_slots + & 1360 size_row_p + size_col_i + size_blk_p + size_thr_c + size_coo_l 1361 compact = new_size .LT. SIZE(matrix%index) 1362 IF (compact) THEN 1363 ! Store old index arrays. 1364 IF (has_row_p) THEN 1365 ALLOCATE (row_p(size_row_p)) 1366 row_p(:) = matrix%row_p(:) 1367 ENDIF 1368 IF (has_col_i) THEN 1369 ALLOCATE (col_i(size_col_i)) 1370 col_i(:) = matrix%col_i(:) 1371 ENDIF 1372 IF (has_blk_p) THEN 1373 ALLOCATE (blk_p(size_blk_p)) 1374 blk_p(:) = matrix%blk_p(:) 1375 ENDIF 1376 IF (has_thr_c) THEN 1377 ALLOCATE (thr_c(size_thr_c)) 1378 thr_c(:) = matrix%thr_c(:) 1379 ENDIF 1380 IF (has_coo_l) THEN 1381 ALLOCATE (coo_l(size_coo_l)) 1382 coo_l(:) = matrix%coo_l(:) 1383 ENDIF 1384 ALLOCATE (meta(dbcsr_num_slots)) 1385 meta(:) = matrix%index(1:dbcsr_num_slots) 1386 ! Clear the index. 1387 CALL memory_deallocate(matrix%index, & 1388 dbcsr_get_index_memory_type(matrix)) 1389 NULLIFY (matrix%index) 1390 CALL memory_allocate(matrix%index, new_size, & 1391 dbcsr_get_index_memory_type(matrix)) 1392 ! 1393 ! Now copy the old index arrays into the index. We must not 1394 ! copy the positions of the old pointers. 1395 matrix%index(1:dbcsr_meta_size) = meta(1:dbcsr_meta_size) 1396 matrix%index(dbcsr_meta_size + 1:) = 0 1397 matrix%index(dbcsr_slot_size) = dbcsr_num_slots 1398 IF (has_thr_c) THEN 1399 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_thr_c, thr_c) 1400 DEALLOCATE (thr_c) 1401 ENDIF 1402 IF (has_row_p) THEN 1403 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_row_p, row_p) 1404 DEALLOCATE (row_p) 1405 ENDIF 1406 IF (has_col_i) THEN 1407 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_col_i, col_i) 1408 DEALLOCATE (col_i) 1409 ENDIF 1410 IF (has_blk_p) THEN 1411 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_blk_p, blk_p) 1412 DEALLOCATE (blk_p) 1413 ENDIF 1414 IF (has_coo_l) THEN 1415 CALL dbcsr_addto_index_array(matrix, dbcsr_slot_coo_l, coo_l) 1416 DEALLOCATE (coo_l) 1417 ENDIF 1418 DEALLOCATE (meta) 1419 IF (careful_mod) THEN 1420 ! This is pretty strong but it should be true. 1421 IF (matrix%index(dbcsr_slot_size) /= new_size) & 1422 DBCSR_ABORT("Unexpected index size.") 1423 IF (SIZE(matrix%index) /= new_size) & 1424 DBCSR_ABORT("Unexpected index size.") 1425 ENDIF 1426 CALL dbcsr_repoint_index(matrix) 1427 ENDIF 1428 CALL timestop(error_handle) 1429 END SUBROUTINE dbcsr_index_compact 1430 1431 SUBROUTINE dbcsr_index_checksum(matrix, checksum) 1432 !! Calculates the checksum of an index. 1433 1434 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 1435 !! matrix for which to make canonical index 1436 INTEGER, INTENT(OUT) :: checksum 1437 1438 CHARACTER(len=*), PARAMETER :: routineN = 'dbcsr_index_checksum' 1439 1440 INTEGER :: error_handle 1441 1442! --------------------------------------------------------------------------- 1443 1444 CALL timeset(routineN, error_handle) 1445 ! 1446 checksum = joaat_hash((/joaat_hash(matrix%row_p), & 1447 joaat_hash(matrix%col_i), & 1448 joaat_hash(matrix%blk_p)/)) 1449 ! 1450 CALL timestop(error_handle) 1451 END SUBROUTINE dbcsr_index_checksum 1452 1453 FUNCTION dbcsr_has_local_row_index(matrix) RESULT(local_indexing) 1454 TYPE(dbcsr_type), INTENT(INOUT) :: matrix 1455 LOGICAL :: local_indexing 1456 1457 local_indexing = matrix%local_indexing 1458 END FUNCTION dbcsr_has_local_row_index 1459 1460END MODULE dbcsr_index_operations 1461