Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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
19 : USE cp_dbcsr_api, ONLY: &
20 : dbcsr_type, dbcsr_release, &
21 : dbcsr_iterator_type, dbcsr_iterator_start, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
22 : dbcsr_has_symmetry, dbcsr_desymmetrize, dbcsr_get_num_blocks, dbcsr_iterator_stop, &
23 : dbcsr_reserve_blocks, dbcsr_finalize
24 : USE dbt_allocate_wrap, ONLY: &
25 : allocate_any
26 : USE dbt_tas_types, ONLY: &
27 : dbt_tas_iterator
28 : USE dbt_tas_base, ONLY: &
29 : dbt_tas_iterator_next_block, dbt_tas_iterator_blocks_left, dbt_tas_iterator_start, &
30 : dbt_tas_iterator_stop, dbt_tas_get_block_p, dbt_tas_put_block, dbt_tas_reserve_blocks, &
31 : dbt_tas_iterator_num_blocks
32 : USE kinds, ONLY: dp, int_8, dp
33 : USE dbt_index, ONLY: &
34 : nd_to_2d_mapping, ndims_mapping, get_nd_indices_tensor, destroy_nd_to_2d_mapping, get_2d_indices_tensor, &
35 : create_nd_to_2d_mapping
36 : USE dbt_array_list_methods, ONLY: &
37 : array_list, get_array_elements, destroy_array_list, sizes_of_arrays, create_array_list, &
38 : get_arrays
39 : USE dbt_types, ONLY: &
40 : dbt_type, ndims_tensor, dbt_blk_sizes, dbt_get_num_blocks, &
41 : dbt_finalize, ndims_matrix_row, ndims_matrix_column
42 :
43 : #include "../base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 : PRIVATE
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbt_block'
48 :
49 : PUBLIC :: &
50 : block_nd, &
51 : create_block, &
52 : dbt_get_block, &
53 : dbt_iterator_num_blocks, &
54 : dbt_iterator_blocks_left, &
55 : dbt_iterator_next_block, &
56 : dbt_iterator_start, &
57 : dbt_iterator_stop, &
58 : dbt_iterator_type, &
59 : dbt_put_block, &
60 : dbt_reserve_blocks, &
61 : destroy_block, &
62 : checker_tr, &
63 : ndims_iterator
64 :
65 : TYPE dbt_iterator_type
66 : TYPE(dbt_tas_iterator) :: iter
67 : TYPE(dbt_type), POINTER :: tensor => NULL()
68 : END TYPE dbt_iterator_type
69 :
70 : TYPE block_nd
71 : INTEGER, DIMENSION(:), ALLOCATABLE :: sizes
72 : REAL(dp), DIMENSION(:), ALLOCATABLE :: blk
73 : END TYPE
74 :
75 : INTERFACE create_block
76 : MODULE PROCEDURE create_block_data
77 : MODULE PROCEDURE create_block_nodata
78 : END INTERFACE
79 :
80 : INTERFACE dbt_put_block
81 : #:for ndim in ndims
82 : MODULE PROCEDURE dbt_put_${ndim}$d_block
83 : #:endfor
84 : MODULE PROCEDURE dbt_put_anyd_block
85 : END INTERFACE
86 :
87 : INTERFACE dbt_get_block
88 : #:for ndim in ndims
89 : MODULE PROCEDURE dbt_get_${ndim}$d_block
90 : MODULE PROCEDURE dbt_allocate_and_get_${ndim}$d_block
91 : #:endfor
92 : MODULE PROCEDURE dbt_get_anyd_block
93 : END INTERFACE
94 :
95 : INTERFACE dbt_reserve_blocks
96 : MODULE PROCEDURE dbt_reserve_blocks_index
97 : MODULE PROCEDURE dbt_reserve_blocks_index_array
98 : MODULE PROCEDURE dbt_reserve_blocks_template
99 : MODULE PROCEDURE dbt_reserve_blocks_tensor_to_matrix
100 : MODULE PROCEDURE dbt_reserve_blocks_matrix_to_tensor
101 : END INTERFACE
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief block size
107 : !> \author Patrick Seewald
108 : ! **************************************************************************************************
109 0 : FUNCTION block_size(block)
110 : TYPE(block_nd), INTENT(IN) :: block
111 : INTEGER, ALLOCATABLE, DIMENSION(:) :: block_size
112 :
113 0 : ALLOCATE (block_size, source=block%sizes)
114 : END FUNCTION
115 :
116 : ! **************************************************************************************************
117 : !> \brief Generalization of block_iterator_start for tensors.
118 : !> \author Patrick Seewald
119 : ! **************************************************************************************************
120 1656446 : SUBROUTINE dbt_iterator_start(iterator, tensor)
121 : TYPE(dbt_iterator_type), INTENT(OUT) :: iterator
122 : TYPE(dbt_type), INTENT(IN), TARGET :: tensor
123 :
124 1656446 : CPASSERT(tensor%valid)
125 1656446 : CALL dbt_tas_iterator_start(iterator%iter, tensor%matrix_rep)
126 1656446 : iterator%tensor => tensor
127 1656446 : END SUBROUTINE
128 :
129 : ! **************************************************************************************************
130 : !> \brief Generalization of block_iterator_stop for tensors.
131 : !> \author Patrick Seewald
132 : ! **************************************************************************************************
133 1656446 : SUBROUTINE dbt_iterator_stop(iterator)
134 : TYPE(dbt_iterator_type), INTENT(INOUT) :: iterator
135 :
136 1656446 : CALL dbt_tas_iterator_stop(iterator%iter)
137 1656446 : END SUBROUTINE
138 :
139 : ! **************************************************************************************************
140 : !> \brief Number of dimensions.
141 : !> \note specification function below must be defined before it is used in
142 : !> the source due to a bug in the IBM XL Fortran compiler (compilation fails)
143 : !> \author Patrick Seewald
144 : ! **************************************************************************************************
145 34394680 : PURE FUNCTION ndims_iterator(iterator)
146 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
147 : INTEGER :: ndims_iterator
148 :
149 34394680 : ndims_iterator = iterator%tensor%nd_index%ndim_nd
150 34394680 : END FUNCTION
151 :
152 : ! **************************************************************************************************
153 : !> \brief iterate over nd blocks of an nd rank tensor, index only (blocks must be retrieved by
154 : !> calling dbt_get_block on tensor).
155 : !> \param ind_nd nd index of block
156 : !> \param blk_size blk size in each dimension
157 : !> \param blk_offset blk offset in each dimension
158 : !> \author Patrick Seewald
159 : ! **************************************************************************************************
160 28966828 : SUBROUTINE dbt_iterator_next_block(iterator, ind_nd, blk_size, blk_offset)
161 : !!
162 : TYPE(dbt_iterator_type), INTENT(INOUT) :: iterator
163 : INTEGER, DIMENSION(ndims_iterator(iterator)), &
164 : INTENT(OUT) :: ind_nd
165 : INTEGER, DIMENSION(ndims_iterator(iterator)), &
166 : INTENT(OUT), OPTIONAL :: blk_size, blk_offset
167 :
168 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
169 :
170 34394680 : CALL dbt_tas_iterator_next_block(iterator%iter, ind_2d(1), ind_2d(2))
171 :
172 34394680 : ind_nd(:) = get_nd_indices_tensor(iterator%tensor%nd_index_blk, ind_2d)
173 34394680 : IF (PRESENT(blk_size)) blk_size(:) = get_array_elements(iterator%tensor%blk_sizes, ind_nd)
174 : ! note: blk_offset needs to be determined by tensor metadata, can not be derived from 2d row/col
175 : ! offset since block index mapping is not consistent with element index mapping
176 34394680 : IF (PRESENT(blk_offset)) blk_offset(:) = get_array_elements(iterator%tensor%blk_offsets, ind_nd)
177 :
178 53567770 : END SUBROUTINE
179 :
180 : ! **************************************************************************************************
181 : !> \brief Generalization of block_iterator_num_blocks for tensors.
182 : !> \author Ole Schuett
183 : ! **************************************************************************************************
184 249001 : FUNCTION dbt_iterator_num_blocks(iterator)
185 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
186 : INTEGER :: dbt_iterator_num_blocks
187 :
188 249001 : dbt_iterator_num_blocks = dbt_tas_iterator_num_blocks(iterator%iter)
189 :
190 249001 : END FUNCTION
191 :
192 : ! **************************************************************************************************
193 : !> \brief Generalization of block_iterator_blocks_left for tensors.
194 : !> \author Patrick Seewald
195 : ! **************************************************************************************************
196 35584313 : FUNCTION dbt_iterator_blocks_left(iterator)
197 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
198 : LOGICAL :: dbt_iterator_blocks_left
199 :
200 35584313 : dbt_iterator_blocks_left = dbt_tas_iterator_blocks_left(iterator%iter)
201 :
202 35584313 : END FUNCTION
203 :
204 : ! **************************************************************************************************
205 : !> \brief reserve blocks from indices as array object
206 : !> \author Patrick Seewald
207 : ! **************************************************************************************************
208 566081 : SUBROUTINE dbt_reserve_blocks_index_array(tensor, blk_ind)
209 : TYPE(dbt_type), INTENT(INOUT) :: tensor
210 : INTEGER, DIMENSION(:, :), INTENT(IN) :: blk_ind
211 : INTEGER :: handle
212 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_index_array'
213 :
214 566081 : CALL timeset(routineN, handle)
215 : #:for ndim in ndims
216 1131870 : IF (ndims_tensor(tensor) == ${ndim}$) THEN
217 566081 : CALL dbt_reserve_blocks(tensor, ${arrlist("blk_ind", nmax=ndim, ndim_pre=1)}$)
218 : END IF
219 : #:endfor
220 566081 : CALL timestop(handle)
221 :
222 566081 : END SUBROUTINE
223 :
224 : ! **************************************************************************************************
225 : !> \brief reserve tensor blocks using block indices
226 : !> \param blk_ind index of blocks to reserve in each dimension
227 : !> \author Patrick Seewald
228 : ! **************************************************************************************************
229 639741 : SUBROUTINE dbt_reserve_blocks_index(tensor, ${varlist("blk_ind")}$)
230 : TYPE(dbt_type), INTENT(INOUT) :: tensor
231 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: ${varlist("blk_ind")}$
232 : INTEGER :: iblk, nblk, handle
233 : INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: cols, rows
234 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
235 639741 : TYPE(array_list) :: blks
236 1279482 : INTEGER, DIMENSION(ndims_tensor(tensor)) :: iblk_nd, ind_nd, nblk_tmp
237 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_index'
238 :
239 639741 : CALL timeset(routineN, handle)
240 639741 : CPASSERT(tensor%valid)
241 :
242 : CALL create_array_list(blks, ndims_tensor(tensor), &
243 1524482 : ${varlist("blk_ind")}$)
244 :
245 2313964 : nblk_tmp(:) = sizes_of_arrays(blks)
246 639741 : nblk = nblk_tmp(1)
247 2316886 : ALLOCATE (cols(nblk), rows(nblk))
248 12383809 : DO iblk = 1, nblk
249 44272306 : iblk_nd(:) = iblk
250 11744068 : ind_nd(:) = get_array_elements(blks, iblk_nd)
251 11744068 : ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind_nd)
252 12383809 : rows(iblk) = ind_2d(1); cols(iblk) = ind_2d(2)
253 : END DO
254 :
255 639741 : CALL dbt_tas_reserve_blocks(tensor%matrix_rep, rows=rows, columns=cols)
256 639741 : CALL dbt_finalize(tensor)
257 639741 : CALL timestop(handle)
258 1279482 : END SUBROUTINE
259 :
260 : ! **************************************************************************************************
261 : !> \brief reserve tensor blocks using template
262 : !> \param tensor_in template tensor
263 : !> \author Patrick Seewald
264 : ! **************************************************************************************************
265 15050 : SUBROUTINE dbt_reserve_blocks_template(tensor_in, tensor_out)
266 : TYPE(dbt_type), INTENT(IN) :: tensor_in
267 : TYPE(dbt_type), INTENT(INOUT) :: tensor_out
268 :
269 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_template'
270 :
271 : TYPE(dbt_iterator_type) :: iter
272 : INTEGER :: handle, nblk, iblk
273 15050 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: blk_ind
274 :
275 15050 : CALL timeset(routineN, handle)
276 :
277 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,tensor_out) &
278 15050 : !$OMP PRIVATE(iter,nblk,iblk,blk_ind)
279 : CALL dbt_iterator_start(iter, tensor_in)
280 : nblk = dbt_iterator_num_blocks(iter)
281 : ALLOCATE (blk_ind(nblk, ndims_tensor(tensor_in)))
282 : DO iblk = 1, nblk
283 : CALL dbt_iterator_next_block(iter, ind_nd=blk_ind(iblk, :))
284 : END DO
285 : CPASSERT(.NOT. dbt_iterator_blocks_left(iter))
286 : CALL dbt_iterator_stop(iter)
287 :
288 : CALL dbt_reserve_blocks(tensor_out, blk_ind)
289 : !$OMP END PARALLEL
290 :
291 15050 : CALL timestop(handle)
292 30100 : END SUBROUTINE
293 :
294 : ! **************************************************************************************************
295 : !> \brief reserve tensor blocks using matrix template
296 : !> \author Patrick Seewald
297 : ! **************************************************************************************************
298 64766 : SUBROUTINE dbt_reserve_blocks_matrix_to_tensor(matrix_in, tensor_out)
299 : TYPE(dbcsr_type), TARGET, INTENT(IN) :: matrix_in
300 : TYPE(dbt_type), INTENT(INOUT) :: tensor_out
301 : TYPE(dbcsr_type), POINTER :: matrix_in_desym
302 :
303 : INTEGER :: blk, iblk, nblk, nblk_per_thread, a, b
304 64766 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_ind_1, blk_ind_2
305 : INTEGER, DIMENSION(2) :: ind_2d
306 : TYPE(dbcsr_iterator_type) :: iter
307 : INTEGER :: handle
308 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_matrix_to_tensor'
309 :
310 64766 : CALL timeset(routineN, handle)
311 :
312 64766 : IF (dbcsr_has_symmetry(matrix_in)) THEN
313 0 : ALLOCATE (matrix_in_desym)
314 0 : CALL dbcsr_desymmetrize(matrix_in, matrix_in_desym)
315 : ELSE
316 : matrix_in_desym => matrix_in
317 : END IF
318 :
319 64766 : nblk = dbcsr_get_num_blocks(matrix_in_desym)
320 248416 : ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
321 64766 : CALL dbcsr_iterator_start(iter, matrix_in_desym)
322 280485 : DO iblk = 1, nblk
323 215719 : CALL dbcsr_iterator_next_block(iter, ind_2d(1), ind_2d(2), blk)
324 280485 : blk_ind_1(iblk) = ind_2d(1); blk_ind_2(iblk) = ind_2d(2)
325 : END DO
326 64766 : CALL dbcsr_iterator_stop(iter)
327 :
328 : !TODO: Parallelize creation of block list.
329 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_out,nblk,blk_ind_1,blk_ind_2) &
330 64766 : !$OMP PRIVATE(nblk_per_thread,a,b)
331 : nblk_per_thread = nblk/omp_get_num_threads() + 1
332 : a = omp_get_thread_num()*nblk_per_thread + 1
333 : b = MIN(a + nblk_per_thread, nblk)
334 : CALL dbt_reserve_blocks(tensor_out, blk_ind_1(a:b), blk_ind_2(a:b))
335 : !$OMP END PARALLEL
336 :
337 64766 : IF (dbcsr_has_symmetry(matrix_in)) THEN
338 0 : CALL dbcsr_release(matrix_in_desym)
339 0 : DEALLOCATE (matrix_in_desym)
340 : END IF
341 :
342 64766 : CALL timestop(handle)
343 194298 : END SUBROUTINE
344 :
345 : ! **************************************************************************************************
346 : !> \brief reserve matrix blocks using tensor template
347 : !> \author Patrick Seewald
348 : ! **************************************************************************************************
349 42244 : SUBROUTINE dbt_reserve_blocks_tensor_to_matrix(tensor_in, matrix_out)
350 : TYPE(dbt_type), INTENT(IN) :: tensor_in
351 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_out
352 : TYPE(dbt_iterator_type) :: iter
353 42244 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_ind_1, blk_ind_2
354 :
355 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_tensor_to_matrix'
356 : INTEGER :: handle, iblk, nblk
357 : INTEGER, DIMENSION(2) :: ind_2d
358 :
359 42244 : CALL timeset(routineN, handle)
360 :
361 42244 : nblk = dbt_get_num_blocks(tensor_in)
362 151532 : ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
363 42244 : iblk = 0
364 :
365 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,matrix_out,iblk,blk_ind_1,blk_ind_2) &
366 42244 : !$OMP PRIVATE(iter,ind_2d)
367 : CALL dbt_iterator_start(iter, tensor_in)
368 : DO WHILE (dbt_iterator_blocks_left(iter))
369 : CALL dbt_iterator_next_block(iter, ind_2d)
370 : IF (dbcsr_has_symmetry(matrix_out)) THEN
371 : IF (checker_tr(ind_2d(1), ind_2d(2))) CYCLE
372 : IF (ind_2d(1) > ind_2d(2)) CALL swap(ind_2d(1), ind_2d(2))
373 : END IF
374 : !$OMP CRITICAL
375 : iblk = iblk + 1
376 : blk_ind_1(iblk) = ind_2d(1)
377 : blk_ind_2(iblk) = ind_2d(2)
378 : !$OMP END CRITICAL
379 : END DO
380 : CALL dbt_iterator_stop(iter)
381 : !$OMP END PARALLEL
382 :
383 42244 : CALL dbcsr_reserve_blocks(matrix_out, blk_ind_1(:iblk), blk_ind_2(:iblk))
384 42244 : CALL dbcsr_finalize(matrix_out)
385 :
386 42244 : CALL timestop(handle)
387 84488 : END SUBROUTINE
388 :
389 : ! **************************************************************************************************
390 : !> \brief Swaps two integers
391 : !> \author Patrick Seewald
392 : ! **************************************************************************************************
393 3421 : ELEMENTAL SUBROUTINE swap(a, b)
394 : INTEGER, INTENT(INOUT) :: a, b
395 : INTEGER :: tmp
396 :
397 3421 : tmp = a
398 3421 : a = b
399 3421 : b = tmp
400 3421 : END SUBROUTINE swap
401 :
402 : ! **************************************************************************************************
403 : !> \brief Create block from array, array can be n-dimensional.
404 : !> \author Patrick Seewald
405 : ! **************************************************************************************************
406 24336258 : SUBROUTINE create_block_data(block, sizes, array)
407 : TYPE(block_nd), INTENT(OUT) :: block
408 : INTEGER, DIMENSION(:), INTENT(IN) :: sizes
409 : REAL(dp), DIMENSION(PRODUCT(sizes)), INTENT(IN) :: array
410 :
411 37284448 : ALLOCATE (block%sizes, source=sizes)
412 2428217028 : ALLOCATE (block%blk, source=array)
413 6474095 : END SUBROUTINE
414 :
415 : ! **************************************************************************************************
416 : !> \brief Create and allocate block, but no data.
417 : !> \author Patrick Seewald
418 : ! **************************************************************************************************
419 0 : SUBROUTINE create_block_nodata(block, sizes)
420 : INTEGER, INTENT(IN), DIMENSION(:) :: sizes
421 : TYPE(block_nd), INTENT(OUT) :: block
422 0 : ALLOCATE (block%sizes, source=sizes)
423 0 : ALLOCATE (block%blk(PRODUCT(sizes)))
424 0 : END SUBROUTINE
425 :
426 : ! **************************************************************************************************
427 : !> \brief
428 : !> \author Patrick Seewald
429 : ! **************************************************************************************************
430 6471167 : SUBROUTINE destroy_block(block)
431 : TYPE(block_nd), INTENT(INOUT) :: block
432 6471167 : DEALLOCATE (block%blk)
433 6471167 : DEALLOCATE (block%sizes)
434 6471167 : END SUBROUTINE
435 :
436 : ! **************************************************************************************************
437 : !> \brief Determines whether a transpose must be applied
438 : !> \param row The absolute matrix row.
439 : !> \param column The absolute matrix column
440 : !> \param
441 : !> \param
442 : !> \param
443 : !> \param
444 : !> \param
445 : !> \param
446 : !> \author Patrick Seewald
447 : ! **************************************************************************************************
448 29432 : ELEMENTAL FUNCTION checker_tr(row, column) RESULT(transpose)
449 : INTEGER, INTENT(IN) :: row, column
450 : LOGICAL :: transpose
451 :
452 29432 : transpose = BTEST(column + row, 0) .EQV. column .GE. row
453 :
454 29432 : END FUNCTION checker_tr
455 :
456 : ! **************************************************************************************************
457 : !> \brief Generic implementation of dbt_put_block, template for datatype
458 : !> \param block block to put
459 : !> \param summation whether block should be summed to existing block
460 : !> \param ind block index
461 : !> \param
462 : !> \param
463 : !> \param
464 : !> \param
465 : !> \param
466 : !> \author Patrick Seewald
467 : ! **************************************************************************************************
468 3310695 : SUBROUTINE dbt_put_anyd_block(tensor, ind, block, summation)
469 : TYPE(block_nd), INTENT(IN) :: block
470 : TYPE(dbt_type), INTENT(INOUT) :: tensor
471 : LOGICAL, INTENT(IN), OPTIONAL :: summation
472 : INTEGER, DIMENSION(ndims_tensor(tensor)), &
473 : INTENT(IN) :: ind
474 :
475 782797 : SELECT CASE (ndims_tensor(tensor))
476 : #:for ndim in ndims
477 : CASE (${ndim}$)
478 3310695 : CALL dbt_put_${ndim}$d_block(tensor, ind, block%sizes, block%blk, summation=summation)
479 : #:endfor
480 : END SELECT
481 3310695 : END SUBROUTINE
482 :
483 : ! **************************************************************************************************
484 : !> \brief Generic implementation of dbt_get_block (arbitrary tensor rank)
485 : !> \param block block to get
486 : !> \param found whether block was found
487 : !> \param ind block index
488 : !> \param
489 : !> \param
490 : !> \param
491 : !> \param
492 : !> \param
493 : !> \author Patrick Seewald
494 : ! **************************************************************************************************
495 3313623 : SUBROUTINE dbt_get_anyd_block(tensor, ind, block, found)
496 : TYPE(block_nd), INTENT(OUT) :: block
497 : LOGICAL, INTENT(OUT) :: found
498 : TYPE(dbt_type), INTENT(INOUT) :: tensor
499 : INTEGER, DIMENSION(ndims_tensor(tensor)), &
500 : INTENT(IN) :: ind
501 6627246 : INTEGER, DIMENSION(ndims_tensor(tensor)) :: blk_size
502 3313623 : REAL(dp), DIMENSION(:), ALLOCATABLE :: block_arr
503 :
504 3313623 : CALL dbt_blk_sizes(tensor, ind, blk_size)
505 19101580 : ALLOCATE (block_arr(PRODUCT(blk_size)))
506 :
507 782845 : SELECT CASE (ndims_tensor(tensor))
508 : #:for ndim in ndims
509 : CASE (${ndim}$)
510 3313623 : CALL dbt_get_${ndim}$d_block(tensor, ind, blk_size, block_arr, found)
511 : #:endfor
512 : END SELECT
513 3313623 : CALL create_block(block, blk_size, block_arr)
514 3313623 : END SUBROUTINE
515 :
516 : #:for ndim in ndims
517 : ! **************************************************************************************************
518 : !> \brief Template for dbt_put_block.
519 : !> \param ind block index
520 : !> \param sizes block size
521 : !> \param block block to put
522 : !> \param summation whether block should be summed to existing block
523 : !> \param
524 : !> \param
525 : !> \param
526 : !> \param
527 : !> \author Patrick Seewald
528 : ! **************************************************************************************************
529 18684525 : SUBROUTINE dbt_put_${ndim}$d_block(tensor, ind, sizes, block, summation)
530 : TYPE(dbt_type), INTENT(INOUT) :: tensor
531 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: ind
532 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: sizes
533 : REAL(dp), DIMENSION(${arrlist("sizes", nmax=ndim)}$), &
534 : INTENT(IN), TARGET :: block
535 : LOGICAL, INTENT(IN), OPTIONAL :: summation
536 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
537 : INTEGER, DIMENSION(2) :: shape_2d
538 18684525 : REAL(dp), POINTER, DIMENSION(:, :) :: block_2d
539 : INTEGER, DIMENSION(${ndim}$) :: shape_nd
540 : LOGICAL :: found, new_block
541 18684525 : REAL(dp), DIMENSION(${arrlist("sizes", nmax=ndim)}$) :: block_check
542 :
543 : LOGICAL, PARAMETER :: debug = .FALSE.
544 : INTEGER :: i
545 :
546 18684525 : new_block = .FALSE.
547 :
548 : IF (debug) THEN
549 : CALL dbt_get_block(tensor, ind, sizes, block_check, found=found)
550 : CPASSERT(found)
551 : END IF
552 :
553 : ASSOCIATE (map_nd => tensor%nd_index_blk%map_nd, &
554 : map1_2d => tensor%nd_index_blk%map1_2d, &
555 : map2_2d => tensor%nd_index_blk%map2_2d)
556 :
557 109219124 : shape_2d = [PRODUCT(sizes(map1_2d)), PRODUCT(sizes(map2_2d))]
558 :
559 162192894 : IF (ALL([map1_2d, map2_2d] == (/(i, i=1, ${ndim}$)/))) THEN
560 : ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
561 18599139 : block_2d(1:shape_2d(1), 1:shape_2d(2)) => block(${shape_colon(ndim)}$)
562 : ELSE
563 : ! need reshape due to rank reordering
564 341544 : ALLOCATE (block_2d(shape_2d(1), shape_2d(2)))
565 85386 : new_block = .TRUE.
566 336757 : shape_nd(map_nd) = sizes
567 758900 : block_2d(:, :) = RESHAPE(RESHAPE(block, SHAPE=shape_nd, order=map_nd), SHAPE=shape_2d)
568 : END IF
569 :
570 37369050 : ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
571 :
572 : END ASSOCIATE
573 :
574 18684525 : CALL dbt_tas_put_block(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d, summation=summation)
575 :
576 18684525 : IF (new_block) DEALLOCATE (block_2d)
577 :
578 18684525 : END SUBROUTINE
579 : #:endfor
580 :
581 : #:for ndim in ndims
582 : ! **************************************************************************************************
583 : !> \brief allocate and get block
584 : !> \param ind block index
585 : !> \param block block to get
586 : !> \param found whether block was found
587 : !> \param
588 : !> \param
589 : !> \param
590 : !> \param
591 : !> \param
592 : !> \author Patrick Seewald
593 : ! **************************************************************************************************
594 19176017 : SUBROUTINE dbt_allocate_and_get_${ndim}$d_block(tensor, ind, block, found)
595 : TYPE(dbt_type), INTENT(INOUT) :: tensor
596 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: ind
597 : REAL(dp), DIMENSION(${shape_colon(ndim)}$), &
598 : ALLOCATABLE, INTENT(OUT) :: block
599 : LOGICAL, INTENT(OUT) :: found
600 : INTEGER, DIMENSION(${ndim}$) :: blk_size
601 :
602 19176017 : CALL dbt_blk_sizes(tensor, ind, blk_size)
603 19176017 : CALL allocate_any(block, shape_spec=blk_size)
604 19176017 : CALL dbt_get_${ndim}$d_block(tensor, ind, blk_size, block, found)
605 :
606 19176017 : END SUBROUTINE
607 : #:endfor
608 :
609 : #:for ndim in ndims
610 : ! **************************************************************************************************
611 : !> \brief Template for dbt_get_block.
612 : !> \param ind block index
613 : !> \param sizes block size
614 : !> \param block block to get
615 : !> \param found whether block was found
616 : !> \author Patrick Seewald
617 : ! **************************************************************************************************
618 22489946 : SUBROUTINE dbt_get_${ndim}$d_block(tensor, ind, sizes, block, found)
619 : TYPE(dbt_type), INTENT(INOUT) :: tensor
620 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: ind
621 : INTEGER, DIMENSION(${ndim}$), INTENT(IN) :: sizes
622 : REAL(dp), DIMENSION(${arrlist("sizes", nmax=ndim)}$), &
623 : INTENT(OUT) :: block
624 : LOGICAL, INTENT(OUT) :: found
625 :
626 : INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
627 22489946 : REAL(dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: block_2d_ptr
628 : INTEGER :: i
629 22489946 : REAL(dp), DIMENSION(${shape_colon(ndim)}$), POINTER :: block_ptr
630 :
631 22489946 : NULLIFY (block_2d_ptr)
632 :
633 22489946 : ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
634 :
635 : ASSOCIATE (map1_2d => tensor%nd_index_blk%map1_2d, &
636 : map2_2d => tensor%nd_index_blk%map2_2d)
637 :
638 22489946 : CALL dbt_tas_get_block_p(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d_ptr)
639 22489946 : found = ASSOCIATED(block_2d_ptr)
640 :
641 22489946 : IF (found) THEN
642 180981446 : IF (ALL([map1_2d, map2_2d] == (/(i, i=1, ${ndim}$)/))) THEN
643 : ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
644 55795225 : block_ptr(${shape_explicit('block', ndim)}$) => block_2d_ptr(:, :)
645 11842825381 : block(${shape_colon(ndim)}$) = block_ptr(${shape_colon(ndim)}$)
646 : ELSE
647 : ! need reshape due to rank reordering
648 20963106 : block(${shape_colon(ndim)}$) = RESHAPE(block_2d_ptr, SHAPE=SHAPE(block), ORDER=[map1_2d, map2_2d])
649 : END IF
650 : END IF
651 :
652 : END ASSOCIATE
653 :
654 22489946 : END SUBROUTINE
655 : #:endfor
656 :
657 0 : END MODULE
|
