Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods to operate on n-dimensional tensor blocks.
10 : !> \author Patrick Seewald
11 : ! **************************************************************************************************
12 : MODULE dbt_block
13 :
14 : #:include "dbt_macros.fypp"
15 : #:set maxdim = maxrank
16 : #:set ndims = range(2,maxdim+1)
17 :
18 : USE OMP_LIB, ONLY: omp_get_thread_num, omp_get_num_threads, omp_get_level
19 : USE util, ONLY: get_limit
20 : USE cp_dbcsr_api, ONLY: &
21 : dbcsr_type, dbcsr_release, &
22 : dbcsr_iterator_type, dbcsr_iterator_start, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
23 : dbcsr_has_symmetry, dbcsr_desymmetrize, dbcsr_get_num_blocks, dbcsr_iterator_stop, &
24 : dbcsr_reserve_blocks, dbcsr_finalize
25 : USE dbt_allocate_wrap, ONLY: &
26 : allocate_any
27 : USE dbt_tas_types, ONLY: &
28 : dbt_tas_iterator
29 : USE dbt_tas_base, ONLY: &
30 : dbt_tas_iterator_next_block, dbt_tas_iterator_blocks_left, dbt_tas_iterator_start, &
31 : dbt_tas_iterator_stop, dbt_tas_get_block_p, dbt_tas_put_block, dbt_tas_reserve_blocks, &
32 : dbt_tas_iterator_num_blocks
33 : USE kinds, ONLY: dp, int_8, dp
34 : USE dbt_index, ONLY: &
35 : nd_to_2d_mapping, ndims_mapping, get_nd_indices_tensor, destroy_nd_to_2d_mapping, get_2d_indices_tensor, &
36 : create_nd_to_2d_mapping
37 : USE dbt_array_list_methods, ONLY: &
38 : array_list, get_array_elements, destroy_array_list, sizes_of_arrays, create_array_list, &
39 : get_arrays
40 : USE dbt_types, ONLY: &
41 : dbt_type, ndims_tensor, dbt_blk_sizes, dbt_get_num_blocks, dbt_get_info, &
42 : dbt_finalize, ndims_matrix_row, ndims_matrix_column
43 :
44 : #include "../base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 : PRIVATE
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbt_block'
49 :
50 : PUBLIC :: &
51 : dbt_get_block_pattern, &
52 : block_nd, &
53 : create_block, &
54 : dbt_get_block, &
55 : dbt_iterator_num_blocks, &
56 : dbt_iterator_blocks_left, &
57 : dbt_iterator_next_block, &
58 : dbt_iterator_start, &
59 : dbt_iterator_stop, &
60 : dbt_iterator_type, &
61 : dbt_put_block, &
62 : dbt_reserve_blocks, &
63 : destroy_block, &
64 : checker_tr, &
65 : dbt_copy_block, &
66 : dbt_get_block_data, &
67 : dbt_put_block_data, &
68 : ndims_iterator
69 :
70 : TYPE dbt_iterator_type
71 : TYPE(dbt_tas_iterator) :: iter
72 : TYPE(dbt_type), POINTER :: tensor => NULL()
73 : END TYPE dbt_iterator_type
74 :
75 : TYPE block_nd
76 : INTEGER, DIMENSION(:), ALLOCATABLE :: sizes
77 : REAL(dp), DIMENSION(:), ALLOCATABLE :: blk
78 : END TYPE block_nd
79 :
80 : INTERFACE create_block
81 : MODULE PROCEDURE create_block_data
82 : MODULE PROCEDURE create_block_nodata
83 : END INTERFACE create_block
84 :
85 : INTERFACE dbt_put_block
86 : #:for ndim in ndims
87 : MODULE PROCEDURE dbt_put_${ndim}$d_block
88 : #:endfor
89 : MODULE PROCEDURE dbt_put_anyd_block
90 : END INTERFACE dbt_put_block
91 :
92 : INTERFACE dbt_get_block
93 : #:for ndim in ndims
94 : MODULE PROCEDURE dbt_get_${ndim}$d_block
95 : MODULE PROCEDURE dbt_allocate_and_get_${ndim}$d_block
96 : #:endfor
97 : MODULE PROCEDURE dbt_get_anyd_block
98 : END INTERFACE dbt_get_block
99 :
100 : INTERFACE dbt_reserve_blocks
101 : MODULE PROCEDURE dbt_reserve_blocks_index
102 : MODULE PROCEDURE dbt_reserve_blocks_index_array
103 : MODULE PROCEDURE dbt_reserve_blocks_template
104 : MODULE PROCEDURE dbt_reserve_blocks_tensor_to_matrix
105 : MODULE PROCEDURE dbt_reserve_blocks_matrix_to_tensor
106 : END INTERFACE dbt_reserve_blocks
107 :
108 : CONTAINS
109 :
110 : ! **************************************************************************************************
111 : !> \brief block size
112 : !> \author Patrick Seewald
113 : ! **************************************************************************************************
114 0 : FUNCTION block_size(block)
115 : TYPE(block_nd), INTENT(IN) :: block
116 : INTEGER, ALLOCATABLE, DIMENSION(:) :: block_size
117 :
118 0 : ALLOCATE (block_size, source=block%sizes)
119 : END FUNCTION block_size
120 :
121 : ! **************************************************************************************************
122 : !> \brief Generalization of block_iterator_start for tensors.
123 : !> \author Patrick Seewald
124 : ! **************************************************************************************************
125 2418760 : SUBROUTINE dbt_iterator_start(iterator, tensor)
126 : TYPE(dbt_iterator_type), INTENT(OUT) :: iterator
127 : TYPE(dbt_type), INTENT(IN), TARGET :: tensor
128 :
129 2418760 : CPASSERT(tensor%valid)
130 2418760 : CALL dbt_tas_iterator_start(iterator%iter, tensor%matrix_rep)
131 2418760 : iterator%tensor => tensor
132 2418760 : END SUBROUTINE dbt_iterator_start
133 :
134 : ! **************************************************************************************************
135 : !> \brief Generalization of block_iterator_stop for tensors.
136 : !> \author Patrick Seewald
137 : ! **************************************************************************************************
138 2418760 : SUBROUTINE dbt_iterator_stop(iterator)
139 : TYPE(dbt_iterator_type), INTENT(INOUT) :: iterator
140 :
141 2418760 : CALL dbt_tas_iterator_stop(iterator%iter)
142 2418760 : END SUBROUTINE dbt_iterator_stop
143 :
144 : ! **************************************************************************************************
145 : !> \brief Number of dimensions.
146 : !> \note specification function below must be defined before it is used in
147 : !> the source due to a bug in the IBM XL Fortran compiler (compilation fails)
148 : !> \author Patrick Seewald
149 : ! **************************************************************************************************
150 57153330 : PURE FUNCTION ndims_iterator(iterator)
151 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
152 : INTEGER :: ndims_iterator
153 :
154 57153330 : ndims_iterator = iterator%tensor%nd_index%ndim_nd
155 57153330 : END FUNCTION ndims_iterator
156 :
157 : ! **************************************************************************************************
158 : !> \brief iterate over nd blocks of an nd rank tensor, index only (blocks must be retrieved by
159 : !> calling dbt_get_block on tensor).
160 : !> \param ind_nd nd index of block
161 : !> \param blk_size blk size in each dimension
162 : !> \param blk_offset blk offset in each dimension
163 : !> \author Patrick Seewald
164 : ! **************************************************************************************************
165 80896072 : SUBROUTINE dbt_iterator_next_block(iterator, ind_nd, blk_size, blk_offset)
166 : !!
167 : TYPE(dbt_iterator_type), INTENT(INOUT) :: iterator
168 : INTEGER, DIMENSION(ndims_iterator(iterator)), &
169 : INTENT(OUT) :: ind_nd
170 : INTEGER, DIMENSION(ndims_iterator(iterator)), &
171 : INTENT(OUT), OPTIONAL :: blk_size, blk_offset
172 :
173 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
174 :
175 57153330 : CALL dbt_tas_iterator_next_block(iterator%iter, ind_2d(1), ind_2d(2))
176 :
177 57153330 : ind_nd(:) = get_nd_indices_tensor(iterator%tensor%nd_index_blk, ind_2d)
178 57153330 : IF (PRESENT(blk_size)) blk_size(:) = get_array_elements(iterator%tensor%blk_sizes, ind_nd)
179 : ! note: blk_offset needs to be determined by tensor metadata, can not be derived from 2d row/col
180 : ! offset since block index mapping is not consistent with element index mapping
181 57153330 : IF (PRESENT(blk_offset)) blk_offset(:) = get_array_elements(iterator%tensor%blk_offsets, ind_nd)
182 :
183 99537147 : END SUBROUTINE dbt_iterator_next_block
184 :
185 : ! **************************************************************************************************
186 : !> \brief Generalization of block_iterator_num_blocks for tensors.
187 : !> \author Ole Schuett
188 : ! **************************************************************************************************
189 208793 : FUNCTION dbt_iterator_num_blocks(iterator)
190 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
191 : INTEGER :: dbt_iterator_num_blocks
192 :
193 208793 : dbt_iterator_num_blocks = dbt_tas_iterator_num_blocks(iterator%iter)
194 :
195 208793 : END FUNCTION dbt_iterator_num_blocks
196 :
197 : ! **************************************************************************************************
198 : !> \brief Generalization of block_iterator_blocks_left for tensors.
199 : !> \author Patrick Seewald
200 : ! **************************************************************************************************
201 58944714 : FUNCTION dbt_iterator_blocks_left(iterator)
202 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
203 : LOGICAL :: dbt_iterator_blocks_left
204 :
205 58944714 : dbt_iterator_blocks_left = dbt_tas_iterator_blocks_left(iterator%iter)
206 :
207 58944714 : END FUNCTION dbt_iterator_blocks_left
208 :
209 : ! **************************************************************************************************
210 : !> \brief reserve blocks from indices as array object
211 : !> A serial caller starts an OpenMP team internally. An existing team calls collectively,
212 : !> with each thread providing its own part of the block list.
213 : !> \author Patrick Seewald
214 : ! **************************************************************************************************
215 682592 : RECURSIVE SUBROUTINE dbt_reserve_blocks_index_array(tensor, blk_ind)
216 : TYPE(dbt_type), INTENT(INOUT) :: tensor
217 : INTEGER, DIMENSION(:, :), INTENT(IN) :: blk_ind
218 : INTEGER :: handle, block_range(2)
219 : LOGICAL :: inside_team
220 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_index_array'
221 :
222 682592 : CALL timeset(routineN, handle)
223 : ! The array interface accepts a serial caller and keeps collective team calls intact.
224 682592 : inside_team = .TRUE.
225 682592 : !$ inside_team = omp_get_level() > 0
226 682592 : IF (.NOT. inside_team) THEN
227 136 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,blk_ind) PRIVATE(block_range)
228 : block_range = get_limit(SIZE(blk_ind, 1), omp_get_num_threads(), omp_get_thread_num())
229 : CALL dbt_reserve_blocks_index_array(tensor, blk_ind(block_range(1):block_range(2), :))
230 : !$OMP END PARALLEL
231 : ELSE
232 : #:for ndim in ndims
233 1364620 : IF (ndims_tensor(tensor) == ${ndim}$) THEN
234 682456 : CALL dbt_reserve_blocks(tensor, ${arrlist("blk_ind", nmax=ndim, ndim_pre=1)}$)
235 : END IF
236 : #:endfor
237 : END IF
238 682592 : CALL timestop(handle)
239 :
240 682592 : END SUBROUTINE dbt_reserve_blocks_index_array
241 :
242 : ! **************************************************************************************************
243 : !> \brief Collect the global stored-block pattern projected onto two tensor dimensions.
244 : !> Called outside an OpenMP team; iteration and reduction are internal.
245 : !> \param tensor source tensor
246 : !> \param dimensions tensor dimensions retained in the projection
247 : !> \param pattern output pattern, one for present entries and zero otherwise
248 : ! **************************************************************************************************
249 136 : SUBROUTINE dbt_get_block_pattern(tensor, dimensions, pattern)
250 : TYPE(dbt_type), INTENT(IN) :: tensor
251 : INTEGER, INTENT(IN) :: dimensions(2)
252 : INTEGER, INTENT(OUT) :: pattern(:, :)
253 :
254 272 : INTEGER :: ind(ndims_tensor(tensor)), nblks(ndims_tensor(tensor))
255 : TYPE(dbt_iterator_type) :: iter
256 :
257 136 : CPASSERT(omp_get_level() == 0)
258 680 : CPASSERT(ALL(dimensions >= 1) .AND. ALL(dimensions <= SIZE(ind)))
259 136 : CALL dbt_get_info(tensor, nblks_total=nblks)
260 408 : CPASSERT(ALL(SHAPE(pattern) == nblks(dimensions)))
261 38280 : pattern = 0
262 136 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,dimensions,pattern) PRIVATE(iter,ind)
263 : CALL dbt_iterator_start(iter, tensor)
264 : DO WHILE (dbt_iterator_blocks_left(iter))
265 : CALL dbt_iterator_next_block(iter, ind)
266 : !$OMP ATOMIC WRITE
267 : pattern(ind(dimensions(1)), ind(dimensions(2))) = 1
268 : END DO
269 : CALL dbt_iterator_stop(iter)
270 : !$OMP END PARALLEL
271 76424 : CALL tensor%pgrid%mp_comm_2d%max(pattern)
272 136 : END SUBROUTINE dbt_get_block_pattern
273 :
274 : ! **************************************************************************************************
275 : !> \brief reserve tensor blocks using block indices
276 : !> \param blk_ind index of blocks to reserve in each dimension
277 : !> \author Patrick Seewald
278 : ! **************************************************************************************************
279 798919 : SUBROUTINE dbt_reserve_blocks_index(tensor, ${varlist("blk_ind")}$)
280 : TYPE(dbt_type), INTENT(INOUT) :: tensor
281 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: ${varlist("blk_ind")}$
282 : INTEGER :: iblk, nblk, handle
283 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: cols, rows
284 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
285 798919 : TYPE(array_list) :: blks
286 1597838 : INTEGER, DIMENSION(ndims_tensor(tensor)) :: iblk_nd, ind_nd, nblk_tmp
287 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_index'
288 :
289 798919 : CALL timeset(routineN, handle)
290 798919 : CPASSERT(tensor%valid)
291 :
292 : CALL create_array_list(blks, ndims_tensor(tensor), &
293 1916251 : ${varlist("blk_ind")}$)
294 :
295 2877263 : nblk_tmp(:) = sizes_of_arrays(blks)
296 798919 : nblk = nblk_tmp(1)
297 2929916 : ALLOCATE (cols(nblk), rows(nblk))
298 15783896 : DO iblk = 1, nblk
299 56718504 : iblk_nd(:) = iblk
300 14984977 : ind_nd(:) = get_array_elements(blks, iblk_nd)
301 14984977 : ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind_nd)
302 15783896 : rows(iblk) = ind_2d(1); cols(iblk) = ind_2d(2)
303 : END DO
304 :
305 798919 : CALL dbt_tas_reserve_blocks(tensor%matrix_rep, rows=rows, columns=cols)
306 798919 : CALL dbt_finalize(tensor)
307 798919 : CALL timestop(handle)
308 1597838 : END SUBROUTINE dbt_reserve_blocks_index
309 :
310 : ! **************************************************************************************************
311 : !> \brief reserve tensor blocks using template
312 : !> \param tensor_in template tensor
313 : !> \author Patrick Seewald
314 : ! **************************************************************************************************
315 24382 : SUBROUTINE dbt_reserve_blocks_template(tensor_in, tensor_out)
316 : TYPE(dbt_type), INTENT(IN) :: tensor_in
317 : TYPE(dbt_type), INTENT(INOUT) :: tensor_out
318 :
319 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_template'
320 :
321 : TYPE(dbt_iterator_type) :: iter
322 : INTEGER :: handle, nblk, iblk
323 24382 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: blk_ind
324 :
325 24382 : CALL timeset(routineN, handle)
326 :
327 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,tensor_out) &
328 24382 : !$OMP PRIVATE(iter,nblk,iblk,blk_ind)
329 : CALL dbt_iterator_start(iter, tensor_in)
330 : nblk = dbt_iterator_num_blocks(iter)
331 : ALLOCATE (blk_ind(nblk, ndims_tensor(tensor_in)))
332 : DO iblk = 1, nblk
333 : CALL dbt_iterator_next_block(iter, ind_nd=blk_ind(iblk, :))
334 : END DO
335 : CPASSERT(.NOT. dbt_iterator_blocks_left(iter))
336 : CALL dbt_iterator_stop(iter)
337 :
338 : CALL dbt_reserve_blocks(tensor_out, blk_ind)
339 : !$OMP END PARALLEL
340 :
341 24382 : CALL timestop(handle)
342 48764 : END SUBROUTINE dbt_reserve_blocks_template
343 :
344 : ! **************************************************************************************************
345 : !> \brief reserve tensor blocks using matrix template
346 : !> \author Patrick Seewald
347 : ! **************************************************************************************************
348 97109 : SUBROUTINE dbt_reserve_blocks_matrix_to_tensor(matrix_in, tensor_out)
349 : TYPE(dbcsr_type), TARGET, INTENT(IN) :: matrix_in
350 : TYPE(dbt_type), INTENT(INOUT) :: tensor_out
351 : TYPE(dbcsr_type), POINTER :: matrix_in_desym
352 :
353 : INTEGER :: iblk, nblk, nblk_per_thread, a, b
354 97109 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_ind_1, blk_ind_2
355 : INTEGER, DIMENSION(2) :: ind_2d
356 : TYPE(dbcsr_iterator_type) :: iter
357 : INTEGER :: handle
358 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_matrix_to_tensor'
359 :
360 97109 : CALL timeset(routineN, handle)
361 :
362 97109 : IF (dbcsr_has_symmetry(matrix_in)) THEN
363 0 : ALLOCATE (matrix_in_desym)
364 0 : CALL dbcsr_desymmetrize(matrix_in, matrix_in_desym)
365 : ELSE
366 : matrix_in_desym => matrix_in
367 : END IF
368 :
369 97109 : nblk = dbcsr_get_num_blocks(matrix_in_desym)
370 378704 : ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
371 97109 : CALL dbcsr_iterator_start(iter, matrix_in_desym)
372 392654 : DO iblk = 1, nblk
373 295545 : CALL dbcsr_iterator_next_block(iter, ind_2d(1), ind_2d(2))
374 392654 : blk_ind_1(iblk) = ind_2d(1); blk_ind_2(iblk) = ind_2d(2)
375 : END DO
376 97109 : CALL dbcsr_iterator_stop(iter)
377 :
378 : !TODO: Parallelize creation of block list.
379 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_out,nblk,blk_ind_1,blk_ind_2) &
380 97109 : !$OMP PRIVATE(nblk_per_thread,a,b)
381 : nblk_per_thread = nblk/omp_get_num_threads() + 1
382 : a = omp_get_thread_num()*nblk_per_thread + 1
383 : b = MIN(a + nblk_per_thread, nblk)
384 : CALL dbt_reserve_blocks(tensor_out, blk_ind_1(a:b), blk_ind_2(a:b))
385 : !$OMP END PARALLEL
386 :
387 97109 : IF (dbcsr_has_symmetry(matrix_in)) THEN
388 0 : CALL dbcsr_release(matrix_in_desym)
389 0 : DEALLOCATE (matrix_in_desym)
390 : END IF
391 :
392 97109 : CALL timestop(handle)
393 291327 : END SUBROUTINE dbt_reserve_blocks_matrix_to_tensor
394 :
395 : ! **************************************************************************************************
396 : !> \brief reserve matrix blocks using tensor template
397 : !> \author Patrick Seewald
398 : ! **************************************************************************************************
399 65292 : SUBROUTINE dbt_reserve_blocks_tensor_to_matrix(tensor_in, matrix_out)
400 : TYPE(dbt_type), INTENT(IN) :: tensor_in
401 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_out
402 : TYPE(dbt_iterator_type) :: iter
403 65292 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_ind_1, blk_ind_2
404 :
405 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_tensor_to_matrix'
406 : INTEGER :: handle, iblk, nblk
407 : INTEGER, DIMENSION(2) :: ind_2d
408 :
409 65292 : CALL timeset(routineN, handle)
410 :
411 65292 : nblk = dbt_get_num_blocks(tensor_in)
412 224420 : ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
413 65292 : iblk = 0
414 :
415 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,matrix_out,iblk,blk_ind_1,blk_ind_2) &
416 65292 : !$OMP PRIVATE(iter,ind_2d)
417 : CALL dbt_iterator_start(iter, tensor_in)
418 : DO WHILE (dbt_iterator_blocks_left(iter))
419 : CALL dbt_iterator_next_block(iter, ind_2d)
420 : IF (dbcsr_has_symmetry(matrix_out)) THEN
421 : IF (checker_tr(ind_2d(1), ind_2d(2))) CYCLE
422 : IF (ind_2d(1) > ind_2d(2)) CALL swap(ind_2d(1), ind_2d(2))
423 : END IF
424 : !$OMP CRITICAL(omp_dbt_reserve_blocks_tensor_to_matrix)
425 : iblk = iblk + 1
426 : blk_ind_1(iblk) = ind_2d(1)
427 : blk_ind_2(iblk) = ind_2d(2)
428 : !$OMP END CRITICAL(omp_dbt_reserve_blocks_tensor_to_matrix)
429 : END DO
430 : CALL dbt_iterator_stop(iter)
431 : !$OMP END PARALLEL
432 :
433 65292 : CALL dbcsr_reserve_blocks(matrix_out, blk_ind_1(:iblk), blk_ind_2(:iblk))
434 65292 : CALL dbcsr_finalize(matrix_out)
435 :
436 65292 : CALL timestop(handle)
437 130584 : END SUBROUTINE dbt_reserve_blocks_tensor_to_matrix
438 :
439 : ! **************************************************************************************************
440 : !> \brief Swaps two integers
441 : !> \author Patrick Seewald
442 : ! **************************************************************************************************
443 5008 : ELEMENTAL SUBROUTINE swap(a, b)
444 : INTEGER, INTENT(INOUT) :: a, b
445 : INTEGER :: tmp
446 :
447 5008 : tmp = a
448 5008 : a = b
449 5008 : b = tmp
450 5008 : END SUBROUTINE swap
451 :
452 : ! **************************************************************************************************
453 : !> \brief Create block from array, array can be n-dimensional.
454 : !> \author Patrick Seewald
455 : ! **************************************************************************************************
456 0 : SUBROUTINE create_block_data(block, sizes, array)
457 : TYPE(block_nd), INTENT(OUT) :: block
458 : INTEGER, DIMENSION(:), INTENT(IN) :: sizes
459 : REAL(dp), DIMENSION(PRODUCT(sizes)), INTENT(IN) :: array
460 :
461 0 : ALLOCATE (block%sizes, source=sizes)
462 0 : ALLOCATE (block%blk, source=array)
463 0 : END SUBROUTINE create_block_data
464 :
465 : ! **************************************************************************************************
466 : !> \brief Create and allocate block, but no data.
467 : !> \author Patrick Seewald
468 : ! **************************************************************************************************
469 206995 : SUBROUTINE create_block_nodata(block, sizes)
470 : INTEGER, INTENT(IN), DIMENSION(:) :: sizes
471 : TYPE(block_nd), INTENT(OUT) :: block
472 1244696 : ALLOCATE (block%sizes, source=sizes)
473 1244098 : ALLOCATE (block%blk(PRODUCT(sizes)))
474 206995 : END SUBROUTINE create_block_nodata
475 :
476 : ! **************************************************************************************************
477 : !> \brief
478 : !> \author Patrick Seewald
479 : ! **************************************************************************************************
480 204067 : SUBROUTINE destroy_block(block)
481 : TYPE(block_nd), INTENT(INOUT) :: block
482 204067 : DEALLOCATE (block%blk)
483 204067 : DEALLOCATE (block%sizes)
484 204067 : END SUBROUTINE destroy_block
485 :
486 : ! **************************************************************************************************
487 : !> \brief Determines whether a transpose must be applied
488 : !> \param row The absolute matrix row.
489 : !> \param column The absolute matrix column
490 : !> \param
491 : !> \param
492 : !> \param
493 : !> \param
494 : !> \param
495 : !> \param
496 : !> \author Patrick Seewald
497 : ! **************************************************************************************************
498 43292 : ELEMENTAL FUNCTION checker_tr(row, column) RESULT(transpose)
499 : INTEGER, INTENT(IN) :: row, column
500 : LOGICAL :: transpose
501 :
502 43292 : transpose = BTEST(column + row, 0) .EQV. column >= row
503 :
504 43292 : END FUNCTION checker_tr
505 :
506 : ! **************************************************************************************************
507 : !> \brief Generic implementation of dbt_put_block, template for datatype
508 : !> \param block block to put
509 : !> \param summation whether block should be summed to existing block
510 : !> \param ind block index
511 : !> \param
512 : !> \param
513 : !> \param
514 : !> \param
515 : !> \param
516 : !> \author Patrick Seewald
517 : ! **************************************************************************************************
518 204067 : SUBROUTINE dbt_put_anyd_block(tensor, ind, block, summation)
519 : TYPE(block_nd), INTENT(IN) :: block
520 : TYPE(dbt_type), INTENT(INOUT) :: tensor
521 : LOGICAL, INTENT(IN), OPTIONAL :: summation
522 : INTEGER, DIMENSION(ndims_tensor(tensor)), &
523 : INTENT(IN) :: ind
524 :
525 204067 : CALL dbt_put_block_data(tensor, ind, block%sizes, block%blk, summation=summation)
526 204067 : END SUBROUTINE dbt_put_anyd_block
527 :
528 : ! **************************************************************************************************
529 : !> \brief Generic implementation of dbt_get_block (arbitrary tensor rank)
530 : !> \param block block to get
531 : !> \param found whether block was found
532 : !> \param ind block index
533 : !> \param
534 : !> \param
535 : !> \param
536 : !> \param
537 : !> \param
538 : !> \author Patrick Seewald
539 : ! **************************************************************************************************
540 206995 : SUBROUTINE dbt_get_anyd_block(tensor, ind, block, found)
541 : TYPE(block_nd), INTENT(OUT) :: block
542 : LOGICAL, INTENT(OUT) :: found
543 : TYPE(dbt_type), INTENT(INOUT) :: tensor
544 : INTEGER, DIMENSION(ndims_tensor(tensor)), &
545 : INTENT(IN) :: ind
546 206995 : INTEGER, DIMENSION(ndims_tensor(tensor)) :: blk_size
547 :
548 206995 : CALL dbt_blk_sizes(tensor, ind, blk_size)
549 206995 : CALL create_block(block, blk_size)
550 206995 : CALL dbt_get_block_data(tensor, ind, blk_size, block%blk, found)
551 206995 : END SUBROUTINE dbt_get_anyd_block
552 :
553 : ! **************************************************************************************************
554 : !> \brief Get a block into the given contiguous data (arbitrary tensor rank), without allocating
555 : !> any intermediate block.
556 : !> \param ind block index
557 : !> \param sizes block size
558 : !> \param data block data
559 : !> \param found whether block was found
560 : !> \author Hans Pabst
561 : ! **************************************************************************************************
562 19535032 : SUBROUTINE dbt_get_block_data(tensor, ind, sizes, data, found)
563 : TYPE(dbt_type), INTENT(INOUT) :: tensor
564 : INTEGER, DIMENSION(ndims_tensor(tensor)), INTENT(IN) :: ind, sizes
565 : REAL(dp), DIMENSION(PRODUCT(sizes)), INTENT(OUT) :: data
566 : LOGICAL, INTENT(OUT) :: found
567 :
568 5150372 : found = .FALSE.
569 1070603 : SELECT CASE (ndims_tensor(tensor))
570 : #:for ndim in ndims
571 : CASE (${ndim}$)
572 5150372 : CALL dbt_get_${ndim}$d_block(tensor, ind, sizes, data, found)
573 : #:endfor
574 : END SELECT
575 5150372 : END SUBROUTINE dbt_get_block_data
576 :
577 : ! **************************************************************************************************
578 : !> \brief Put a block from the given contiguous data (arbitrary tensor rank), without allocating
579 : !> any intermediate block.
580 : !> \param ind block index
581 : !> \param sizes block size
582 : !> \param data block data
583 : !> \param summation whether block should be summed to existing block
584 : !> \author Hans Pabst
585 : ! **************************************************************************************************
586 19922532 : SUBROUTINE dbt_put_block_data(tensor, ind, sizes, data, summation)
587 : TYPE(dbt_type), INTENT(INOUT) :: tensor
588 : INTEGER, DIMENSION(ndims_tensor(tensor)), INTENT(IN) :: ind, sizes
589 : REAL(dp), DIMENSION(PRODUCT(sizes)), INTENT(IN) :: data
590 : LOGICAL, INTENT(IN), OPTIONAL :: summation
591 :
592 1070555 : SELECT CASE (ndims_tensor(tensor))
593 : #:for ndim in ndims
594 : CASE (${ndim}$)
595 5247919 : CALL dbt_put_${ndim}$d_block(tensor, ind, sizes, data, summation=summation)
596 : #:endfor
597 : END SELECT
598 5247919 : END SUBROUTINE dbt_put_block_data
599 :
600 : ! **************************************************************************************************
601 : !> \brief Copy a block between tensors with the same index mapping and block sizes (e.g., one
602 : !> created from the other), which needs no reshape since the matrix blocks are equal.
603 : !> \param ind block index
604 : !> \param summation whether block should be summed to existing block
605 : !> \author Hans Pabst
606 : ! **************************************************************************************************
607 1185885 : SUBROUTINE dbt_copy_block(tensor_in, tensor_out, ind, summation)
608 : TYPE(dbt_type), INTENT(INOUT) :: tensor_in, tensor_out
609 : INTEGER, DIMENSION(ndims_tensor(tensor_in)), INTENT(IN) :: ind
610 : LOGICAL, INTENT(IN), OPTIONAL :: summation
611 :
612 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
613 1185885 : REAL(dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: block_2d
614 :
615 1185885 : NULLIFY (block_2d)
616 1185885 : ind_2d(:) = get_2d_indices_tensor(tensor_in%nd_index_blk, ind)
617 1185885 : CALL dbt_tas_get_block_p(tensor_in%matrix_rep, ind_2d(1), ind_2d(2), block_2d)
618 1185885 : IF (ASSOCIATED(block_2d)) THEN
619 1185885 : CALL dbt_tas_put_block(tensor_out%matrix_rep, ind_2d(1), ind_2d(2), block_2d, summation=summation)
620 : END IF
621 1185885 : END SUBROUTINE dbt_copy_block
622 :
623 : #:for ndim in ndims
624 : ! **************************************************************************************************
625 : !> \brief Template for dbt_put_block.
626 : !> \param ind block index
627 : !> \param sizes block size
628 : !> \param block block to put
629 : !> \param summation whether block should be summed to existing block
630 : !> \param
631 : !> \param
632 : !> \param
633 : !> \param
634 : !> \author Patrick Seewald
635 : ! **************************************************************************************************
636 17592085 : SUBROUTINE dbt_put_${ndim}$d_block(tensor, ind, sizes, block, summation)
637 : TYPE(dbt_type), INTENT(INOUT) :: tensor
638 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: ind
639 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: sizes
640 : REAL(dp), DIMENSION(${arrlist("sizes", nmax=ndim)}$), &
641 : INTENT(IN), TARGET :: block
642 : LOGICAL, INTENT(IN), OPTIONAL :: summation
643 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
644 : INTEGER, DIMENSION(2) :: shape_2d
645 17592085 : REAL(dp), POINTER, DIMENSION(:, :) :: block_2d
646 : INTEGER, DIMENSION(${ndim}$) :: shape_nd
647 : LOGICAL :: found, new_block
648 17592085 : REAL(dp), DIMENSION(${arrlist("sizes", nmax=ndim)}$) :: block_check
649 :
650 : LOGICAL, PARAMETER :: debug = .FALSE.
651 : INTEGER :: i
652 :
653 17592085 : new_block = .FALSE.
654 :
655 : IF (debug) THEN
656 : CALL dbt_get_block(tensor, ind, sizes, block_check, found=found)
657 : CPASSERT(found)
658 : END IF
659 :
660 : ASSOCIATE (map_nd => tensor%nd_index_blk%map_nd, &
661 : map1_2d => tensor%nd_index_blk%map1_2d, &
662 : map2_2d => tensor%nd_index_blk%map2_2d)
663 :
664 102267627 : shape_2d = [PRODUCT(sizes(map1_2d)), PRODUCT(sizes(map2_2d))]
665 :
666 150660896 : IF (ALL([map1_2d, map2_2d] == [(i, i=1, ${ndim}$)])) THEN
667 : ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
668 17134430 : block_2d(1:shape_2d(1), 1:shape_2d(2)) => block(${shape_colon(ndim)}$)
669 : ELSE
670 : ! need reshape due to rank reordering
671 1830620 : ALLOCATE (block_2d(shape_2d(1), shape_2d(2)))
672 457655 : new_block = .TRUE.
673 1819705 : shape_nd(map_nd) = sizes
674 4097065 : block_2d(:, :) = RESHAPE(RESHAPE(block, SHAPE=shape_nd, order=map_nd), SHAPE=shape_2d)
675 : END IF
676 :
677 35184170 : ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
678 :
679 : END ASSOCIATE
680 :
681 17592085 : CALL dbt_tas_put_block(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d, summation=summation)
682 :
683 17592085 : IF (new_block) DEALLOCATE (block_2d)
684 :
685 17592085 : END SUBROUTINE dbt_put_${ndim}$d_block
686 : #:endfor
687 :
688 : #:for ndim in ndims
689 : ! **************************************************************************************************
690 : !> \brief allocate and get block
691 : !> \param ind block index
692 : !> \param block block to get
693 : !> \param found whether block was found
694 : !> \param
695 : !> \param
696 : !> \param
697 : !> \param
698 : !> \param
699 : !> \author Patrick Seewald
700 : ! **************************************************************************************************
701 26031678 : SUBROUTINE dbt_allocate_and_get_${ndim}$d_block(tensor, ind, block, found)
702 : TYPE(dbt_type), INTENT(INOUT) :: tensor
703 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: ind
704 : REAL(dp), DIMENSION(${shape_colon(ndim)}$), &
705 : ALLOCATABLE, INTENT(OUT) :: block
706 : LOGICAL, INTENT(OUT) :: found
707 : INTEGER, DIMENSION(${ndim}$) :: blk_size
708 :
709 13015839 : CALL dbt_blk_sizes(tensor, ind, blk_size)
710 13015839 : CALL allocate_any(block, shape_spec=blk_size)
711 13015839 : CALL dbt_get_${ndim}$d_block(tensor, ind, blk_size, block, found)
712 :
713 13015839 : END SUBROUTINE dbt_allocate_and_get_${ndim}$d_block
714 : #:endfor
715 :
716 : #:for ndim in ndims
717 : ! **************************************************************************************************
718 : !> \brief Template for dbt_get_block.
719 : !> \param ind block index
720 : !> \param sizes block size
721 : !> \param block block to get
722 : !> \param found whether block was found
723 : !> \author Patrick Seewald
724 : ! **************************************************************************************************
725 18166723 : SUBROUTINE dbt_get_${ndim}$d_block(tensor, ind, sizes, block, found)
726 : TYPE(dbt_type), INTENT(INOUT) :: tensor
727 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: ind
728 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: sizes
729 : REAL(dp), DIMENSION(${arrlist("sizes", nmax=ndim)}$), &
730 : INTENT(OUT) :: block
731 : LOGICAL, INTENT(OUT) :: found
732 :
733 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
734 18166723 : REAL(dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: block_2d_ptr
735 : INTEGER :: i
736 18166723 : REAL(dp), DIMENSION(${shape_colon(ndim)}$), POINTER :: block_ptr
737 :
738 18166723 : NULLIFY (block_2d_ptr)
739 :
740 18166723 : ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
741 :
742 : ASSOCIATE (map1_2d => tensor%nd_index_blk%map1_2d, &
743 : map2_2d => tensor%nd_index_blk%map2_2d)
744 :
745 18166723 : CALL dbt_tas_get_block_p(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d_ptr)
746 18166723 : found = ASSOCIATED(block_2d_ptr)
747 :
748 18166723 : IF (found) THEN
749 137202733 : IF (ALL([map1_2d, map2_2d] == [(i, i=1, ${ndim}$)])) THEN
750 : ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
751 40348680 : block_ptr(${shape_explicit('block', ndim)}$) => block_2d_ptr(:, :)
752 6755277801 : block(${shape_colon(ndim)}$) = block_ptr(${shape_colon(ndim)}$)
753 : ELSE
754 : ! need reshape due to rank reordering
755 24982536 : block(${shape_colon(ndim)}$) = RESHAPE(block_2d_ptr, SHAPE=SHAPE(block), ORDER=[map1_2d, map2_2d])
756 : END IF
757 : END IF
758 :
759 : END ASSOCIATE
760 :
761 18166723 : END SUBROUTINE dbt_get_${ndim}$d_block
762 : #:endfor
763 :
764 0 : END MODULE dbt_block
|
