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 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 1530957 : SUBROUTINE dbt_iterator_start(iterator, tensor)
121 : TYPE(dbt_iterator_type), INTENT(OUT) :: iterator
122 : TYPE(dbt_type), INTENT(IN), TARGET :: tensor
123 :
124 1530957 : CPASSERT(tensor%valid)
125 1530957 : CALL dbt_tas_iterator_start(iterator%iter, tensor%matrix_rep)
126 1530957 : iterator%tensor => tensor
127 1530957 : END SUBROUTINE
128 :
129 : ! **************************************************************************************************
130 : !> \brief Generalization of block_iterator_stop for tensors.
131 : !> \author Patrick Seewald
132 : ! **************************************************************************************************
133 1530957 : SUBROUTINE dbt_iterator_stop(iterator)
134 : TYPE(dbt_iterator_type), INTENT(INOUT) :: iterator
135 :
136 1530957 : CALL dbt_tas_iterator_stop(iterator%iter)
137 1530957 : 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 28562006 : PURE FUNCTION ndims_iterator(iterator)
146 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
147 : INTEGER :: ndims_iterator
148 :
149 28562006 : ndims_iterator = iterator%tensor%nd_index%ndim_nd
150 28562006 : 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 26549239 : 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 28562006 : CALL dbt_tas_iterator_next_block(iterator%iter, ind_2d(1), ind_2d(2))
171 :
172 28562006 : ind_nd(:) = get_nd_indices_tensor(iterator%tensor%nd_index_blk, ind_2d)
173 28562006 : 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 28562006 : IF (PRESENT(blk_offset)) blk_offset(:) = get_array_elements(iterator%tensor%blk_offsets, ind_nd)
177 :
178 46017290 : END SUBROUTINE
179 :
180 : ! **************************************************************************************************
181 : !> \brief Generalization of block_iterator_num_blocks for tensors.
182 : !> \author Ole Schuett
183 : ! **************************************************************************************************
184 237873 : FUNCTION dbt_iterator_num_blocks(iterator)
185 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
186 : INTEGER :: dbt_iterator_num_blocks
187 :
188 237873 : dbt_iterator_num_blocks = dbt_tas_iterator_num_blocks(iterator%iter)
189 :
190 237873 : END FUNCTION
191 :
192 : ! **************************************************************************************************
193 : !> \brief Generalization of block_iterator_blocks_left for tensors.
194 : !> \author Patrick Seewald
195 : ! **************************************************************************************************
196 29706715 : FUNCTION dbt_iterator_blocks_left(iterator)
197 : TYPE(dbt_iterator_type), INTENT(IN) :: iterator
198 : LOGICAL :: dbt_iterator_blocks_left
199 :
200 29706715 : dbt_iterator_blocks_left = dbt_tas_iterator_blocks_left(iterator%iter)
201 :
202 29706715 : END FUNCTION
203 :
204 : ! **************************************************************************************************
205 : !> \brief reserve blocks from indices as array object
206 : !> \author Patrick Seewald
207 : ! **************************************************************************************************
208 513922 : 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 513922 : CALL timeset(routineN, handle)
215 : #:for ndim in ndims
216 1027552 : IF (ndims_tensor(tensor) == ${ndim}$) THEN
217 513922 : CALL dbt_reserve_blocks(tensor, ${arrlist("blk_ind", nmax=ndim, ndim_pre=1)}$)
218 : END IF
219 : #:endfor
220 513922 : CALL timestop(handle)
221 :
222 513922 : 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 578730 : 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 578730 : TYPE(array_list) :: blks
236 1157460 : INTEGER, DIMENSION(ndims_tensor(tensor)) :: iblk_nd, ind_nd, nblk_tmp
237 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_index'
238 :
239 578730 : CALL timeset(routineN, handle)
240 578730 : CPASSERT(tensor%valid)
241 :
242 : CALL create_array_list(blks, ndims_tensor(tensor), &
243 1400052 : ${varlist("blk_ind")}$)
244 :
245 2072328 : nblk_tmp(:) = sizes_of_arrays(blks)
246 578730 : nblk = nblk_tmp(1)
247 2096374 : ALLOCATE (cols(nblk), rows(nblk))
248 11712817 : DO iblk = 1, nblk
249 41928753 : iblk_nd(:) = iblk
250 11134087 : ind_nd(:) = get_array_elements(blks, iblk_nd)
251 11134087 : ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind_nd)
252 11712817 : rows(iblk) = ind_2d(1); cols(iblk) = ind_2d(2)
253 : END DO
254 :
255 578730 : CALL dbt_tas_reserve_blocks(tensor%matrix_rep, rows=rows, columns=cols)
256 578730 : CALL dbt_finalize(tensor)
257 578730 : CALL timestop(handle)
258 1157460 : END SUBROUTINE
259 :
260 : ! **************************************************************************************************
261 : !> \brief reserve tensor blocks using template
262 : !> \param tensor_in template tensor
263 : !> \author Patrick Seewald
264 : ! **************************************************************************************************
265 3134 : 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 3134 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: blk_ind
274 :
275 3134 : CALL timeset(routineN, handle)
276 :
277 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,tensor_out) &
278 3134 : !$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 3134 : CALL timestop(handle)
292 6268 : END SUBROUTINE
293 :
294 : ! **************************************************************************************************
295 : !> \brief reserve tensor blocks using matrix template
296 : !> \author Patrick Seewald
297 : ! **************************************************************************************************
298 57004 : 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 57004 : 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 57004 : CALL timeset(routineN, handle)
311 :
312 57004 : 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 57004 : nblk = dbcsr_get_num_blocks(matrix_in_desym)
320 219874 : ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
321 57004 : CALL dbcsr_iterator_start(iter, matrix_in_desym)
322 263922 : DO iblk = 1, nblk
323 206918 : CALL dbcsr_iterator_next_block(iter, ind_2d(1), ind_2d(2), blk)
324 263922 : blk_ind_1(iblk) = ind_2d(1); blk_ind_2(iblk) = ind_2d(2)
325 : END DO
326 57004 : 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 57004 : !$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 57004 : 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 57004 : CALL timestop(handle)
343 171012 : END SUBROUTINE
344 :
345 : ! **************************************************************************************************
346 : !> \brief reserve matrix blocks using tensor template
347 : !> \author Patrick Seewald
348 : ! **************************************************************************************************
349 35298 : 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 35298 : 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 35298 : CALL timeset(routineN, handle)
360 :
361 35298 : nblk = dbt_get_num_blocks(tensor_in)
362 124798 : ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
363 35298 : iblk = 0
364 :
365 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,matrix_out,iblk,blk_ind_1,blk_ind_2) &
366 35298 : !$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 35298 : CALL dbcsr_reserve_blocks(matrix_out, blk_ind_1(:iblk), blk_ind_2(:iblk))
384 35298 : CALL dbcsr_finalize(matrix_out)
385 :
386 35298 : CALL timestop(handle)
387 70596 : END SUBROUTINE
388 :
389 : ! **************************************************************************************************
390 : !> \brief Swaps two integers
391 : !> \author Patrick Seewald
392 : ! **************************************************************************************************
393 3502 : ELEMENTAL SUBROUTINE swap(a, b)
394 : INTEGER, INTENT(INOUT) :: a, b
395 : INTEGER :: tmp
396 :
397 3502 : tmp = a
398 3502 : a = b
399 3502 : b = tmp
400 3502 : END SUBROUTINE swap
401 :
402 : ! **************************************************************************************************
403 : !> \brief Create block from array, array can be n-dimensional.
404 : !> \author Patrick Seewald
405 : ! **************************************************************************************************
406 22219808 : 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 34075724 : ALLOCATE (block%sizes, source=sizes)
412 2374636888 : ALLOCATE (block%blk, source=array)
413 5927958 : 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 5925030 : SUBROUTINE destroy_block(block)
431 : TYPE(block_nd), INTENT(INOUT) :: block
432 5925030 : DEALLOCATE (block%blk)
433 5925030 : DEALLOCATE (block%sizes)
434 5925030 : 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 30056 : ELEMENTAL FUNCTION checker_tr(row, column) RESULT(transpose)
449 : INTEGER, INTENT(IN) :: row, column
450 : LOGICAL :: transpose
451 :
452 30056 : transpose = BTEST(column + row, 0) .EQV. column .GE. row
453 :
454 30056 : 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 3008436 : 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 748748 : SELECT CASE (ndims_tensor(tensor))
476 : #:for ndim in ndims
477 : CASE (${ndim}$)
478 3008436 : CALL dbt_put_${ndim}$d_block(tensor, ind, block%sizes, block%blk, summation=summation)
479 : #:endfor
480 : END SELECT
481 3008436 : 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 3011364 : 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 6022728 : INTEGER, DIMENSION(ndims_tensor(tensor)) :: blk_size
502 3011364 : REAL(dp), DIMENSION(:), ALLOCATABLE :: block_arr
503 :
504 3011364 : CALL dbt_blk_sizes(tensor, ind, blk_size)
505 17322075 : ALLOCATE (block_arr(PRODUCT(blk_size)))
506 :
507 748796 : SELECT CASE (ndims_tensor(tensor))
508 : #:for ndim in ndims
509 : CASE (${ndim}$)
510 3011364 : CALL dbt_get_${ndim}$d_block(tensor, ind, blk_size, block_arr, found)
511 : #:endfor
512 : END SELECT
513 3011364 : CALL create_block(block, blk_size, block_arr)
514 3011364 : 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 15379071 : 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 15379071 : REAL(dp), POINTER, DIMENSION(:, :) :: block_2d
539 : INTEGER, DIMENSION(${ndim}$) :: shape_nd
540 : LOGICAL :: found, new_block
541 15379071 : REAL(dp), DIMENSION(${arrlist("sizes", nmax=ndim)}$) :: block_check
542 :
543 : LOGICAL, PARAMETER :: debug = .FALSE.
544 : INTEGER :: i
545 :
546 15379071 : 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 89465318 : shape_2d = [PRODUCT(sizes(map1_2d)), PRODUCT(sizes(map2_2d))]
558 :
559 132764014 : 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 15368719 : block_2d(1:shape_2d(1), 1:shape_2d(2)) => block(${shape_colon(ndim)}$)
562 : ELSE
563 : ! need reshape due to rank reordering
564 41408 : ALLOCATE (block_2d(shape_2d(1), shape_2d(2)))
565 10352 : new_block = .TRUE.
566 40697 : shape_nd(map_nd) = sizes
567 91746 : block_2d(:, :) = RESHAPE(RESHAPE(block, SHAPE=shape_nd, order=map_nd), SHAPE=shape_2d)
568 : END IF
569 :
570 30758142 : ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
571 :
572 : END ASSOCIATE
573 :
574 15379071 : CALL dbt_tas_put_block(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d, summation=summation)
575 :
576 15379071 : IF (new_block) DEALLOCATE (block_2d)
577 :
578 15379071 : 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 14471167 : 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 14471167 : CALL dbt_blk_sizes(tensor, ind, blk_size)
603 14471167 : CALL allocate_any(block, shape_spec=blk_size)
604 14471167 : CALL dbt_get_${ndim}$d_block(tensor, ind, blk_size, block, found)
605 :
606 14471167 : 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 17482713 : 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 17482713 : REAL(dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: block_2d_ptr
628 : INTEGER :: i
629 17482713 : REAL(dp), DIMENSION(${shape_colon(ndim)}$), POINTER :: block_ptr
630 :
631 17482713 : NULLIFY (block_2d_ptr)
632 :
633 17482713 : 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 17482713 : CALL dbt_tas_get_block_p(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d_ptr)
639 17482713 : found = ASSOCIATED(block_2d_ptr)
640 :
641 17482713 : IF (found) THEN
642 134896292 : 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 40792230 : block_ptr(${shape_explicit('block', ndim)}$) => block_2d_ptr(:, :)
645 8857247164 : block(${shape_colon(ndim)}$) = block_ptr(${shape_colon(ndim)}$)
646 : ELSE
647 : ! need reshape due to rank reordering
648 19126674 : 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 17482713 : END SUBROUTINE
655 : #:endfor
656 :
657 0 : END MODULE
|
