(git:374b731)
Loading...
Searching...
No Matches
dbt_block.F
Go to the documentation of this file.
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! **************************************************************************************************
13
14
15
16 USE omp_lib, ONLY: omp_get_thread_num, omp_get_num_threads
17 USE dbcsr_api, ONLY: &
18 dbcsr_type, dbcsr_release, &
19 dbcsr_iterator_type, dbcsr_iterator_start, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
20 dbcsr_has_symmetry, dbcsr_desymmetrize, dbcsr_get_num_blocks, dbcsr_iterator_stop, &
21 dbcsr_reserve_blocks, dbcsr_finalize
22 USE dbt_allocate_wrap, ONLY: &
24 USE dbt_tas_types, ONLY: &
26 USE dbt_tas_base, ONLY: &
30 USE kinds, ONLY: dp, int_8, dp
31 USE dbt_index, ONLY: &
34 USE dbt_array_list_methods, ONLY: &
37 USE dbt_types, ONLY: &
40
41#include "../base/base_uses.f90"
42
43 IMPLICIT NONE
44 PRIVATE
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbt_block'
46
47 PUBLIC :: &
48 block_nd, &
60 checker_tr, &
62
64 TYPE(dbt_tas_iterator) :: iter
65 TYPE(dbt_type), POINTER :: tensor => null()
66 END TYPE dbt_iterator_type
67
69 INTEGER, DIMENSION(:), ALLOCATABLE :: sizes
70 REAL(dp), DIMENSION(:), ALLOCATABLE :: blk
71 END TYPE
72
73 INTERFACE create_block
74 MODULE PROCEDURE create_block_data
75 MODULE PROCEDURE create_block_nodata
76 END INTERFACE
77
78 INTERFACE dbt_put_block
79 MODULE PROCEDURE dbt_put_2d_block
80 MODULE PROCEDURE dbt_put_3d_block
81 MODULE PROCEDURE dbt_put_4d_block
82 MODULE PROCEDURE dbt_put_anyd_block
83 END INTERFACE
84
85 INTERFACE dbt_get_block
86 MODULE PROCEDURE dbt_get_2d_block
87 MODULE PROCEDURE dbt_allocate_and_get_2d_block
88 MODULE PROCEDURE dbt_get_3d_block
89 MODULE PROCEDURE dbt_allocate_and_get_3d_block
90 MODULE PROCEDURE dbt_get_4d_block
91 MODULE PROCEDURE dbt_allocate_and_get_4d_block
92 MODULE PROCEDURE dbt_get_anyd_block
93 END INTERFACE
94
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
103CONTAINS
104
105! **************************************************************************************************
106!> \brief block size
107!> \author Patrick Seewald
108! **************************************************************************************************
109 FUNCTION block_size(block)
110 TYPE(block_nd), INTENT(IN) :: block
111 INTEGER, ALLOCATABLE, DIMENSION(:) :: block_size
112
113 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 SUBROUTINE dbt_iterator_start(iterator, tensor)
121 TYPE(dbt_iterator_type), INTENT(OUT) :: iterator
122 TYPE(dbt_type), INTENT(IN), TARGET :: tensor
123
124 cpassert(tensor%valid)
125 CALL dbt_tas_iterator_start(iterator%iter, tensor%matrix_rep)
126 iterator%tensor => tensor
127 END SUBROUTINE
128
129! **************************************************************************************************
130!> \brief Generalization of block_iterator_stop for tensors.
131!> \author Patrick Seewald
132! **************************************************************************************************
133 SUBROUTINE dbt_iterator_stop(iterator)
134 TYPE(dbt_iterator_type), INTENT(INOUT) :: iterator
135
136 CALL dbt_tas_iterator_stop(iterator%iter)
137 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 PURE FUNCTION ndims_iterator(iterator)
146 TYPE(dbt_iterator_type), INTENT(IN) :: iterator
147 INTEGER :: ndims_iterator
148
149 ndims_iterator = iterator%tensor%nd_index%ndim_nd
150 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 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 CALL dbt_tas_iterator_next_block(iterator%iter, ind_2d(1), ind_2d(2))
171
172 ind_nd(:) = get_nd_indices_tensor(iterator%tensor%nd_index_blk, ind_2d)
173 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 IF (PRESENT(blk_offset)) blk_offset(:) = get_array_elements(iterator%tensor%blk_offsets, ind_nd)
177
178 END SUBROUTINE
179
180! **************************************************************************************************
181!> \brief Generalization of block_iterator_num_blocks for tensors.
182!> \author Ole Schuett
183! **************************************************************************************************
184 FUNCTION dbt_iterator_num_blocks(iterator)
185 TYPE(dbt_iterator_type), INTENT(IN) :: iterator
186 INTEGER :: dbt_iterator_num_blocks
187
189
190 END FUNCTION
191
192! **************************************************************************************************
193!> \brief Generalization of block_iterator_blocks_left for tensors.
194!> \author Patrick Seewald
195! **************************************************************************************************
196 FUNCTION dbt_iterator_blocks_left(iterator)
197 TYPE(dbt_iterator_type), INTENT(IN) :: iterator
198 LOGICAL :: dbt_iterator_blocks_left
199
201
202 END FUNCTION
203
204! **************************************************************************************************
205!> \brief reserve blocks from indices as array object
206!> \author Patrick Seewald
207! **************************************************************************************************
208 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 CALL timeset(routinen, handle)
215 IF (ndims_tensor(tensor) == 2) THEN
216 CALL dbt_reserve_blocks(tensor, blk_ind(:,1), blk_ind(:,2))
217 END IF
218 IF (ndims_tensor(tensor) == 3) THEN
219 CALL dbt_reserve_blocks(tensor, blk_ind(:,1), blk_ind(:,2), blk_ind(:,3))
220 END IF
221 IF (ndims_tensor(tensor) == 4) THEN
222 CALL dbt_reserve_blocks(tensor, blk_ind(:,1), blk_ind(:,2), blk_ind(:,3), blk_ind(:,4))
223 END IF
224 CALL timestop(handle)
225
226 END SUBROUTINE
227
228! **************************************************************************************************
229!> \brief reserve tensor blocks using block indices
230!> \param blk_ind index of blocks to reserve in each dimension
231!> \author Patrick Seewald
232! **************************************************************************************************
233 SUBROUTINE dbt_reserve_blocks_index(tensor, blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4)
234 TYPE(dbt_type), INTENT(INOUT) :: tensor
235 INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4
236 INTEGER :: iblk, nblk, handle
237 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: cols, rows
238 INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
239 TYPE(array_list) :: blks
240 INTEGER, DIMENSION(ndims_tensor(tensor)) :: iblk_nd, ind_nd, nblk_tmp
241 CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_index'
242
243 CALL timeset(routinen, handle)
244 cpassert(tensor%valid)
245
246 CALL create_array_list(blks, ndims_tensor(tensor), &
247 blk_ind_1, blk_ind_2, blk_ind_3, blk_ind_4)
248
249 nblk_tmp(:) = sizes_of_arrays(blks)
250 nblk = nblk_tmp(1)
251 ALLOCATE (cols(nblk), rows(nblk))
252 DO iblk = 1, nblk
253 iblk_nd(:) = iblk
254 ind_nd(:) = get_array_elements(blks, iblk_nd)
255 ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind_nd)
256 rows(iblk) = ind_2d(1); cols(iblk) = ind_2d(2)
257 END DO
258
259 CALL dbt_tas_reserve_blocks(tensor%matrix_rep, rows=rows, columns=cols)
260 CALL dbt_finalize(tensor)
261 CALL timestop(handle)
262 END SUBROUTINE
263
264! **************************************************************************************************
265!> \brief reserve tensor blocks using template
266!> \param tensor_in template tensor
267!> \author Patrick Seewald
268! **************************************************************************************************
269 SUBROUTINE dbt_reserve_blocks_template(tensor_in, tensor_out)
270 TYPE(dbt_type), INTENT(IN) :: tensor_in
271 TYPE(dbt_type), INTENT(INOUT) :: tensor_out
272
273 CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_template'
274
275 TYPE(dbt_iterator_type) :: iter
276 INTEGER :: handle, nblk, iblk
277 INTEGER, DIMENSION(:, :), ALLOCATABLE :: blk_ind
278
279 CALL timeset(routinen, handle)
280
281!$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,tensor_out) &
282!$OMP PRIVATE(iter,nblk,iblk,blk_ind)
283 CALL dbt_iterator_start(iter, tensor_in)
284 nblk = dbt_iterator_num_blocks(iter)
285 ALLOCATE (blk_ind(nblk, ndims_tensor(tensor_in)))
286 DO iblk = 1, nblk
287 CALL dbt_iterator_next_block(iter, ind_nd=blk_ind(iblk, :))
288 END DO
289 cpassert(.NOT. dbt_iterator_blocks_left(iter))
290 CALL dbt_iterator_stop(iter)
291
292 CALL dbt_reserve_blocks(tensor_out, blk_ind)
293!$OMP END PARALLEL
294
295 CALL timestop(handle)
296 END SUBROUTINE
297
298! **************************************************************************************************
299!> \brief reserve tensor blocks using matrix template
300!> \author Patrick Seewald
301! **************************************************************************************************
302 SUBROUTINE dbt_reserve_blocks_matrix_to_tensor(matrix_in, tensor_out)
303 TYPE(dbcsr_type), TARGET, INTENT(IN) :: matrix_in
304 TYPE(dbt_type), INTENT(INOUT) :: tensor_out
305 TYPE(dbcsr_type), POINTER :: matrix_in_desym
306
307 INTEGER :: blk, iblk, nblk, nblk_per_thread, a, b
308 INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_ind_1, blk_ind_2
309 INTEGER, DIMENSION(2) :: ind_2d
310 TYPE(dbcsr_iterator_type) :: iter
311 INTEGER :: handle
312 CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_matrix_to_tensor'
313
314 CALL timeset(routinen, handle)
315
316 IF (dbcsr_has_symmetry(matrix_in)) THEN
317 ALLOCATE (matrix_in_desym)
318 CALL dbcsr_desymmetrize(matrix_in, matrix_in_desym)
319 ELSE
320 matrix_in_desym => matrix_in
321 END IF
322
323 nblk = dbcsr_get_num_blocks(matrix_in_desym)
324 ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
325 CALL dbcsr_iterator_start(iter, matrix_in_desym)
326 DO iblk = 1, nblk
327 CALL dbcsr_iterator_next_block(iter, ind_2d(1), ind_2d(2), blk)
328 blk_ind_1(iblk) = ind_2d(1); blk_ind_2(iblk) = ind_2d(2)
329 END DO
330 CALL dbcsr_iterator_stop(iter)
331
332!TODO: Parallelize creation of block list.
333!$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_out,nblk,blk_ind_1,blk_ind_2) &
334!$OMP PRIVATE(nblk_per_thread,a,b)
335 nblk_per_thread = nblk/omp_get_num_threads() + 1
336 a = omp_get_thread_num()*nblk_per_thread + 1
337 b = min(a + nblk_per_thread, nblk)
338 CALL dbt_reserve_blocks(tensor_out, blk_ind_1(a:b), blk_ind_2(a:b))
339!$OMP END PARALLEL
340
341 IF (dbcsr_has_symmetry(matrix_in)) THEN
342 CALL dbcsr_release(matrix_in_desym)
343 DEALLOCATE (matrix_in_desym)
344 END IF
345
346 CALL timestop(handle)
347 END SUBROUTINE
348
349! **************************************************************************************************
350!> \brief reserve matrix blocks using tensor template
351!> \author Patrick Seewald
352! **************************************************************************************************
353 SUBROUTINE dbt_reserve_blocks_tensor_to_matrix(tensor_in, matrix_out)
354 TYPE(dbt_type), INTENT(IN) :: tensor_in
355 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_out
356 TYPE(dbt_iterator_type) :: iter
357 INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_ind_1, blk_ind_2
358
359 CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_tensor_to_matrix'
360 INTEGER :: handle, iblk, nblk
361 INTEGER, DIMENSION(2) :: ind_2d
362
363 CALL timeset(routinen, handle)
364
365 nblk = dbt_get_num_blocks(tensor_in)
366 ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
367 iblk = 0
368
369!$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,matrix_out,iblk,blk_ind_1,blk_ind_2) &
370!$OMP PRIVATE(iter,ind_2d)
371 CALL dbt_iterator_start(iter, tensor_in)
372 DO WHILE (dbt_iterator_blocks_left(iter))
373 CALL dbt_iterator_next_block(iter, ind_2d)
374 IF (dbcsr_has_symmetry(matrix_out)) THEN
375 IF (checker_tr(ind_2d(1), ind_2d(2))) cycle
376 IF (ind_2d(1) > ind_2d(2)) CALL swap(ind_2d(1), ind_2d(2))
377 END IF
378!$OMP CRITICAL
379 iblk = iblk + 1
380 blk_ind_1(iblk) = ind_2d(1)
381 blk_ind_2(iblk) = ind_2d(2)
382!$OMP END CRITICAL
383 END DO
384 CALL dbt_iterator_stop(iter)
385!$OMP END PARALLEL
386
387 CALL dbcsr_reserve_blocks(matrix_out, blk_ind_1(:iblk), blk_ind_2(:iblk))
388 CALL dbcsr_finalize(matrix_out)
389
390 CALL timestop(handle)
391 END SUBROUTINE
392
393! **************************************************************************************************
394!> \brief Swaps two integers
395!> \author Patrick Seewald
396! **************************************************************************************************
397 ELEMENTAL SUBROUTINE swap(a, b)
398 INTEGER, INTENT(INOUT) :: a, b
399 INTEGER :: tmp
400
401 tmp = a
402 a = b
403 b = tmp
404 END SUBROUTINE swap
405
406! **************************************************************************************************
407!> \brief Create block from array, array can be n-dimensional.
408!> \author Patrick Seewald
409! **************************************************************************************************
410 SUBROUTINE create_block_data(block, sizes, array)
411 TYPE(block_nd), INTENT(OUT) :: block
412 INTEGER, DIMENSION(:), INTENT(IN) :: sizes
413 REAL(dp), DIMENSION(PRODUCT(sizes)), INTENT(IN) :: array
414
415 ALLOCATE (block%sizes, source=sizes)
416 ALLOCATE (block%blk, source=array)
417 END SUBROUTINE
418
419! **************************************************************************************************
420!> \brief Create and allocate block, but no data.
421!> \author Patrick Seewald
422! **************************************************************************************************
423 SUBROUTINE create_block_nodata(block, sizes)
424 INTEGER, INTENT(IN), DIMENSION(:) :: sizes
425 TYPE(block_nd), INTENT(OUT) :: block
426 ALLOCATE (block%sizes, source=sizes)
427 ALLOCATE (block%blk(product(sizes)))
428 END SUBROUTINE
429
430! **************************************************************************************************
431!> \brief
432!> \author Patrick Seewald
433! **************************************************************************************************
434 SUBROUTINE destroy_block(block)
435 TYPE(block_nd), INTENT(INOUT) :: block
436 DEALLOCATE (block%blk)
437 DEALLOCATE (block%sizes)
438 END SUBROUTINE
439
440! **************************************************************************************************
441!> \brief Determines whether a transpose must be applied
442!> \param row The absolute matrix row.
443!> \param column The absolute matrix column
444!> \param
445!> \param
446!> \param
447!> \param
448!> \param
449!> \param
450!> \author Patrick Seewald
451! **************************************************************************************************
452 ELEMENTAL FUNCTION checker_tr(row, column) RESULT(transpose)
453 INTEGER, INTENT(IN) :: row, column
454 LOGICAL :: transpose
455
456 transpose = btest(column + row, 0) .EQV. column .GE. row
457
458 END FUNCTION checker_tr
459
460! **************************************************************************************************
461!> \brief Generic implementation of dbt_put_block, template for datatype
462!> \param block block to put
463!> \param summation whether block should be summed to existing block
464!> \param ind block index
465!> \param
466!> \param
467!> \param
468!> \param
469!> \param
470!> \author Patrick Seewald
471! **************************************************************************************************
472 SUBROUTINE dbt_put_anyd_block(tensor, ind, block, summation)
473 TYPE(block_nd), INTENT(IN) :: block
474 TYPE(dbt_type), INTENT(INOUT) :: tensor
475 LOGICAL, INTENT(IN), OPTIONAL :: summation
476 INTEGER, DIMENSION(ndims_tensor(tensor)), &
477 INTENT(IN) :: ind
478
479 SELECT CASE (ndims_tensor(tensor))
480 CASE (2)
481 CALL dbt_put_2d_block(tensor, ind, block%sizes, block%blk, summation=summation)
482 CASE (3)
483 CALL dbt_put_3d_block(tensor, ind, block%sizes, block%blk, summation=summation)
484 CASE (4)
485 CALL dbt_put_4d_block(tensor, ind, block%sizes, block%blk, summation=summation)
486 END SELECT
487 END SUBROUTINE
488
489! **************************************************************************************************
490!> \brief Generic implementation of dbt_get_block (arbitrary tensor rank)
491!> \param block block to get
492!> \param found whether block was found
493!> \param ind block index
494!> \param
495!> \param
496!> \param
497!> \param
498!> \param
499!> \author Patrick Seewald
500! **************************************************************************************************
501 SUBROUTINE dbt_get_anyd_block(tensor, ind, block, found)
502 TYPE(block_nd), INTENT(OUT) :: block
503 LOGICAL, INTENT(OUT) :: found
504 TYPE(dbt_type), INTENT(INOUT) :: tensor
505 INTEGER, DIMENSION(ndims_tensor(tensor)), &
506 INTENT(IN) :: ind
507 INTEGER, DIMENSION(ndims_tensor(tensor)) :: blk_size
508 REAL(dp), DIMENSION(:), ALLOCATABLE :: block_arr
509
510 CALL dbt_blk_sizes(tensor, ind, blk_size)
511 ALLOCATE (block_arr(product(blk_size)))
512
513 SELECT CASE (ndims_tensor(tensor))
514 CASE (2)
515 CALL dbt_get_2d_block(tensor, ind, blk_size, block_arr, found)
516 CASE (3)
517 CALL dbt_get_3d_block(tensor, ind, blk_size, block_arr, found)
518 CASE (4)
519 CALL dbt_get_4d_block(tensor, ind, blk_size, block_arr, found)
520 END SELECT
521 CALL create_block(block, blk_size, block_arr)
522 END SUBROUTINE
523
524! **************************************************************************************************
525!> \brief Template for dbt_put_block.
526!> \param ind block index
527!> \param sizes block size
528!> \param block block to put
529!> \param summation whether block should be summed to existing block
530!> \param
531!> \param
532!> \param
533!> \param
534!> \author Patrick Seewald
535! **************************************************************************************************
536 SUBROUTINE dbt_put_2d_block(tensor, ind, sizes, block, summation)
537 TYPE(dbt_type), INTENT(INOUT) :: tensor
538 INTEGER, DIMENSION(2), INTENT(IN) :: ind
539 INTEGER, DIMENSION(2), INTENT(IN) :: sizes
540 REAL(dp), DIMENSION(sizes(1), sizes(2)), &
541 INTENT(IN), TARGET :: block
542 LOGICAL, INTENT(IN), OPTIONAL :: summation
543 INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
544 INTEGER, DIMENSION(2) :: shape_2d
545 REAL(dp), POINTER, DIMENSION(:, :) :: block_2d
546 INTEGER, DIMENSION(2) :: shape_nd
547 LOGICAL :: found, new_block
548 REAL(dp), DIMENSION(sizes(1), sizes(2)) :: block_check
549
550 LOGICAL, PARAMETER :: debug = .false.
551 INTEGER :: i
552
553 new_block = .false.
554
555 IF (debug) THEN
556 CALL dbt_get_block(tensor, ind, sizes, block_check, found=found)
557 cpassert(found)
558 END IF
559
560 associate(map_nd => tensor%nd_index_blk%map_nd, &
561 map1_2d => tensor%nd_index_blk%map1_2d, &
562 map2_2d => tensor%nd_index_blk%map2_2d)
563
564 shape_2d = [product(sizes(map1_2d)), product(sizes(map2_2d))]
565
566 IF (all([map1_2d, map2_2d] == (/(i, i=1, 2)/))) THEN
567 ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
568 block_2d(1:shape_2d(1), 1:shape_2d(2)) => block(:,:)
569 ELSE
570 ! need reshape due to rank reordering
571 ALLOCATE (block_2d(shape_2d(1), shape_2d(2)))
572 new_block = .true.
573 shape_nd(map_nd) = sizes
574 block_2d(:, :) = reshape(reshape(block, shape=shape_nd, order=map_nd), shape=shape_2d)
575 END IF
576
577 ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
578
579 END associate
580
581 CALL dbt_tas_put_block(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d, summation=summation)
582
583 IF (new_block) DEALLOCATE (block_2d)
584
585 END SUBROUTINE
586! **************************************************************************************************
587!> \brief Template for dbt_put_block.
588!> \param ind block index
589!> \param sizes block size
590!> \param block block to put
591!> \param summation whether block should be summed to existing block
592!> \param
593!> \param
594!> \param
595!> \param
596!> \author Patrick Seewald
597! **************************************************************************************************
598 SUBROUTINE dbt_put_3d_block(tensor, ind, sizes, block, summation)
599 TYPE(dbt_type), INTENT(INOUT) :: tensor
600 INTEGER, DIMENSION(3), INTENT(IN) :: ind
601 INTEGER, DIMENSION(3), INTENT(IN) :: sizes
602 REAL(dp), DIMENSION(sizes(1), sizes(2), sizes(3)), &
603 INTENT(IN), TARGET :: block
604 LOGICAL, INTENT(IN), OPTIONAL :: summation
605 INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
606 INTEGER, DIMENSION(2) :: shape_2d
607 REAL(dp), POINTER, DIMENSION(:, :) :: block_2d
608 INTEGER, DIMENSION(3) :: shape_nd
609 LOGICAL :: found, new_block
610 REAL(dp), DIMENSION(sizes(1), sizes(2), sizes(3)) :: block_check
611
612 LOGICAL, PARAMETER :: debug = .false.
613 INTEGER :: i
614
615 new_block = .false.
616
617 IF (debug) THEN
618 CALL dbt_get_block(tensor, ind, sizes, block_check, found=found)
619 cpassert(found)
620 END IF
621
622 associate(map_nd => tensor%nd_index_blk%map_nd, &
623 map1_2d => tensor%nd_index_blk%map1_2d, &
624 map2_2d => tensor%nd_index_blk%map2_2d)
625
626 shape_2d = [product(sizes(map1_2d)), product(sizes(map2_2d))]
627
628 IF (all([map1_2d, map2_2d] == (/(i, i=1, 3)/))) THEN
629 ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
630 block_2d(1:shape_2d(1), 1:shape_2d(2)) => block(:,:,:)
631 ELSE
632 ! need reshape due to rank reordering
633 ALLOCATE (block_2d(shape_2d(1), shape_2d(2)))
634 new_block = .true.
635 shape_nd(map_nd) = sizes
636 block_2d(:, :) = reshape(reshape(block, shape=shape_nd, order=map_nd), shape=shape_2d)
637 END IF
638
639 ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
640
641 END associate
642
643 CALL dbt_tas_put_block(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d, summation=summation)
644
645 IF (new_block) DEALLOCATE (block_2d)
646
647 END SUBROUTINE
648! **************************************************************************************************
649!> \brief Template for dbt_put_block.
650!> \param ind block index
651!> \param sizes block size
652!> \param block block to put
653!> \param summation whether block should be summed to existing block
654!> \param
655!> \param
656!> \param
657!> \param
658!> \author Patrick Seewald
659! **************************************************************************************************
660 SUBROUTINE dbt_put_4d_block(tensor, ind, sizes, block, summation)
661 TYPE(dbt_type), INTENT(INOUT) :: tensor
662 INTEGER, DIMENSION(4), INTENT(IN) :: ind
663 INTEGER, DIMENSION(4), INTENT(IN) :: sizes
664 REAL(dp), DIMENSION(sizes(1), sizes(2), sizes(3), sizes(4)), &
665 INTENT(IN), TARGET :: block
666 LOGICAL, INTENT(IN), OPTIONAL :: summation
667 INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
668 INTEGER, DIMENSION(2) :: shape_2d
669 REAL(dp), POINTER, DIMENSION(:, :) :: block_2d
670 INTEGER, DIMENSION(4) :: shape_nd
671 LOGICAL :: found, new_block
672 REAL(dp), DIMENSION(sizes(1), sizes(2), sizes(3), sizes(4)) :: block_check
673
674 LOGICAL, PARAMETER :: debug = .false.
675 INTEGER :: i
676
677 new_block = .false.
678
679 IF (debug) THEN
680 CALL dbt_get_block(tensor, ind, sizes, block_check, found=found)
681 cpassert(found)
682 END IF
683
684 associate(map_nd => tensor%nd_index_blk%map_nd, &
685 map1_2d => tensor%nd_index_blk%map1_2d, &
686 map2_2d => tensor%nd_index_blk%map2_2d)
687
688 shape_2d = [product(sizes(map1_2d)), product(sizes(map2_2d))]
689
690 IF (all([map1_2d, map2_2d] == (/(i, i=1, 4)/))) THEN
691 ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
692 block_2d(1:shape_2d(1), 1:shape_2d(2)) => block(:,:,:,:)
693 ELSE
694 ! need reshape due to rank reordering
695 ALLOCATE (block_2d(shape_2d(1), shape_2d(2)))
696 new_block = .true.
697 shape_nd(map_nd) = sizes
698 block_2d(:, :) = reshape(reshape(block, shape=shape_nd, order=map_nd), shape=shape_2d)
699 END IF
700
701 ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
702
703 END associate
704
705 CALL dbt_tas_put_block(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d, summation=summation)
706
707 IF (new_block) DEALLOCATE (block_2d)
708
709 END SUBROUTINE
710
711! **************************************************************************************************
712!> \brief allocate and get block
713!> \param ind block index
714!> \param block block to get
715!> \param found whether block was found
716!> \param
717!> \param
718!> \param
719!> \param
720!> \param
721!> \author Patrick Seewald
722! **************************************************************************************************
723 SUBROUTINE dbt_allocate_and_get_2d_block(tensor, ind, block, found)
724 TYPE(dbt_type), INTENT(INOUT) :: tensor
725 INTEGER, DIMENSION(2), INTENT(IN) :: ind
726 REAL(dp), DIMENSION(:,:), &
727 ALLOCATABLE, INTENT(OUT) :: block
728 LOGICAL, INTENT(OUT) :: found
729 INTEGER, DIMENSION(2) :: blk_size
730
731 CALL dbt_blk_sizes(tensor, ind, blk_size)
732 CALL allocate_any(block, shape_spec=blk_size)
733 CALL dbt_get_2d_block(tensor, ind, blk_size, block, found)
734
735 END SUBROUTINE
736! **************************************************************************************************
737!> \brief allocate and get block
738!> \param ind block index
739!> \param block block to get
740!> \param found whether block was found
741!> \param
742!> \param
743!> \param
744!> \param
745!> \param
746!> \author Patrick Seewald
747! **************************************************************************************************
748 SUBROUTINE dbt_allocate_and_get_3d_block(tensor, ind, block, found)
749 TYPE(dbt_type), INTENT(INOUT) :: tensor
750 INTEGER, DIMENSION(3), INTENT(IN) :: ind
751 REAL(dp), DIMENSION(:,:,:), &
752 ALLOCATABLE, INTENT(OUT) :: block
753 LOGICAL, INTENT(OUT) :: found
754 INTEGER, DIMENSION(3) :: blk_size
755
756 CALL dbt_blk_sizes(tensor, ind, blk_size)
757 CALL allocate_any(block, shape_spec=blk_size)
758 CALL dbt_get_3d_block(tensor, ind, blk_size, block, found)
759
760 END SUBROUTINE
761! **************************************************************************************************
762!> \brief allocate and get block
763!> \param ind block index
764!> \param block block to get
765!> \param found whether block was found
766!> \param
767!> \param
768!> \param
769!> \param
770!> \param
771!> \author Patrick Seewald
772! **************************************************************************************************
773 SUBROUTINE dbt_allocate_and_get_4d_block(tensor, ind, block, found)
774 TYPE(dbt_type), INTENT(INOUT) :: tensor
775 INTEGER, DIMENSION(4), INTENT(IN) :: ind
776 REAL(dp), DIMENSION(:,:,:,:), &
777 ALLOCATABLE, INTENT(OUT) :: block
778 LOGICAL, INTENT(OUT) :: found
779 INTEGER, DIMENSION(4) :: blk_size
780
781 CALL dbt_blk_sizes(tensor, ind, blk_size)
782 CALL allocate_any(block, shape_spec=blk_size)
783 CALL dbt_get_4d_block(tensor, ind, blk_size, block, found)
784
785 END SUBROUTINE
786
787! **************************************************************************************************
788!> \brief Template for dbt_get_block.
789!> \param ind block index
790!> \param sizes block size
791!> \param block block to get
792!> \param found whether block was found
793!> \author Patrick Seewald
794! **************************************************************************************************
795 SUBROUTINE dbt_get_2d_block(tensor, ind, sizes, block, found)
796 TYPE(dbt_type), INTENT(INOUT) :: tensor
797 INTEGER, DIMENSION(2), INTENT(IN) :: ind
798 INTEGER, DIMENSION(2), INTENT(IN) :: sizes
799 REAL(dp), DIMENSION(sizes(1), sizes(2)), &
800 INTENT(OUT) :: block
801 LOGICAL, INTENT(OUT) :: found
802
803 INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
804 REAL(dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: block_2d_ptr
805 INTEGER :: i
806 REAL(dp), DIMENSION(:,:), POINTER :: block_ptr
807
808 NULLIFY (block_2d_ptr)
809
810 ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
811
812 associate(map1_2d => tensor%nd_index_blk%map1_2d, &
813 map2_2d => tensor%nd_index_blk%map2_2d)
814
815 CALL dbt_tas_get_block_p(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d_ptr)
816 found = ASSOCIATED(block_2d_ptr)
817
818 IF (found) THEN
819 IF (all([map1_2d, map2_2d] == (/(i, i=1, 2)/))) THEN
820 ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
821 block_ptr(lbound(block, 1):ubound(block, 1), lbound(block, 2):ubound(block, 2)) => block_2d_ptr(:, :)
822 block(:,:) = block_ptr(:,:)
823 ELSE
824 ! need reshape due to rank reordering
825 block(:,:) = reshape(block_2d_ptr, shape=shape(block), order=[map1_2d, map2_2d])
826 END IF
827 END IF
828
829 END associate
830
831 END SUBROUTINE
832! **************************************************************************************************
833!> \brief Template for dbt_get_block.
834!> \param ind block index
835!> \param sizes block size
836!> \param block block to get
837!> \param found whether block was found
838!> \author Patrick Seewald
839! **************************************************************************************************
840 SUBROUTINE dbt_get_3d_block(tensor, ind, sizes, block, found)
841 TYPE(dbt_type), INTENT(INOUT) :: tensor
842 INTEGER, DIMENSION(3), INTENT(IN) :: ind
843 INTEGER, DIMENSION(3), INTENT(IN) :: sizes
844 REAL(dp), DIMENSION(sizes(1), sizes(2), sizes(3)), &
845 INTENT(OUT) :: block
846 LOGICAL, INTENT(OUT) :: found
847
848 INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
849 REAL(dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: block_2d_ptr
850 INTEGER :: i
851 REAL(dp), DIMENSION(:,:,:), POINTER :: block_ptr
852
853 NULLIFY (block_2d_ptr)
854
855 ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
856
857 associate(map1_2d => tensor%nd_index_blk%map1_2d, &
858 map2_2d => tensor%nd_index_blk%map2_2d)
859
860 CALL dbt_tas_get_block_p(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d_ptr)
861 found = ASSOCIATED(block_2d_ptr)
862
863 IF (found) THEN
864 IF (all([map1_2d, map2_2d] == (/(i, i=1, 3)/))) THEN
865 ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
866 block_ptr(lbound(block, 1):ubound(block, 1), lbound(block, 2):ubound(block, 2), lbound(block, 3):ubound(block,&
867 & 3)) => block_2d_ptr(:, :)
868 block(:,:,:) = block_ptr(:,:,:)
869 ELSE
870 ! need reshape due to rank reordering
871 block(:,:,:) = reshape(block_2d_ptr, shape=shape(block), order=[map1_2d, map2_2d])
872 END IF
873 END IF
874
875 END associate
876
877 END SUBROUTINE
878! **************************************************************************************************
879!> \brief Template for dbt_get_block.
880!> \param ind block index
881!> \param sizes block size
882!> \param block block to get
883!> \param found whether block was found
884!> \author Patrick Seewald
885! **************************************************************************************************
886 SUBROUTINE dbt_get_4d_block(tensor, ind, sizes, block, found)
887 TYPE(dbt_type), INTENT(INOUT) :: tensor
888 INTEGER, DIMENSION(4), INTENT(IN) :: ind
889 INTEGER, DIMENSION(4), INTENT(IN) :: sizes
890 REAL(dp), DIMENSION(sizes(1), sizes(2), sizes(3), sizes(4)), &
891 INTENT(OUT) :: block
892 LOGICAL, INTENT(OUT) :: found
893
894 INTEGER(KIND=int_8), DIMENSION(2) :: ind_2d
895 REAL(dp), DIMENSION(:, :), POINTER, CONTIGUOUS :: block_2d_ptr
896 INTEGER :: i
897 REAL(dp), DIMENSION(:,:,:,:), POINTER :: block_ptr
898
899 NULLIFY (block_2d_ptr)
900
901 ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
902
903 associate(map1_2d => tensor%nd_index_blk%map1_2d, &
904 map2_2d => tensor%nd_index_blk%map2_2d)
905
906 CALL dbt_tas_get_block_p(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d_ptr)
907 found = ASSOCIATED(block_2d_ptr)
908
909 IF (found) THEN
910 IF (all([map1_2d, map2_2d] == (/(i, i=1, 4)/))) THEN
911 ! to avoid costly reshape can do pointer bounds remapping as long as arrays are equivalent in memory
912 block_ptr(lbound(block, 1):ubound(block, 1), lbound(block, 2):ubound(block, 2), lbound(block, 3):ubound(block,&
913 & 3), lbound(block, 4):ubound(block, 4)) => block_2d_ptr(:, :)
914 block(:,:,:,:) = block_ptr(:,:,:,:)
915 ELSE
916 ! need reshape due to rank reordering
917 block(:,:,:,:) = reshape(block_2d_ptr, shape=shape(block), order=[map1_2d, map2_2d])
918 END IF
919 END IF
920
921 END associate
922
923 END SUBROUTINE
924
925END MODULE
struct tensor_ tensor
Wrapper for allocating, copying and reshaping arrays.
Representation of arbitrary number of 1d integer arrays with arbitrary sizes. This is needed for gene...
subroutine, public get_arrays(list, data_1, data_2, data_3, data_4, i_selected)
Get all arrays contained in list.
subroutine, public create_array_list(list, ndata, data_1, data_2, data_3, data_4)
collects any number of arrays of different sizes into a single array (listcol_data),...
subroutine, public destroy_array_list(list)
destroy array list.
integer function, dimension(:), allocatable, public sizes_of_arrays(list)
sizes of arrays stored in list
pure integer function, dimension(number_of_arrays(list)), public get_array_elements(list, indices)
Get an element for each array.
Methods to operate on n-dimensional tensor blocks.
Definition dbt_block.F:12
elemental logical function, public checker_tr(row, column)
Determines whether a transpose must be applied.
Definition dbt_block.F:453
integer function, public dbt_iterator_num_blocks(iterator)
Generalization of block_iterator_num_blocks for tensors.
Definition dbt_block.F:185
logical function, public dbt_iterator_blocks_left(iterator)
Generalization of block_iterator_blocks_left for tensors.
Definition dbt_block.F:197
subroutine, public destroy_block(block)
Definition dbt_block.F:435
pure integer function, public ndims_iterator(iterator)
Number of dimensions.
Definition dbt_block.F:146
subroutine, public dbt_iterator_stop(iterator)
Generalization of block_iterator_stop for tensors.
Definition dbt_block.F:134
subroutine, public dbt_iterator_start(iterator, tensor)
Generalization of block_iterator_start for tensors.
Definition dbt_block.F:121
subroutine, public dbt_iterator_next_block(iterator, ind_nd, blk_size, blk_offset)
iterate over nd blocks of an nd rank tensor, index only (blocks must be retrieved by calling dbt_get_...
Definition dbt_block.F:161
tensor index and mapping to DBM index
Definition dbt_index.F:12
subroutine, public create_nd_to_2d_mapping(map, dims, map1_2d, map2_2d, base, col_major)
Create all data needed to quickly map between nd index and 2d index.
Definition dbt_index.F:73
pure integer(kind=int_8) function, dimension(2), public get_2d_indices_tensor(map, ind_in)
transform nd index to 2d index, using info from index mapping.
Definition dbt_index.F:318
pure integer function, public ndims_mapping(map)
Definition dbt_index.F:130
pure integer function, dimension(map%ndim_nd), public get_nd_indices_tensor(map, ind_in)
transform 2d index to nd index, using info from index mapping.
Definition dbt_index.F:368
subroutine, public destroy_nd_to_2d_mapping(map)
Definition dbt_index.F:115
Tall-and-skinny matrices: base routines similar to DBM API, mostly wrappers around existing DBM routi...
integer function, public dbt_tas_iterator_num_blocks(iter)
As dbm_iterator_num_blocks.
subroutine, public dbt_tas_iterator_start(iter, matrix_in)
As dbm_iterator_start.
logical function, public dbt_tas_iterator_blocks_left(iter)
As dbm_iterator_blocks_left.
subroutine, public dbt_tas_iterator_stop(iter)
As dbm_iterator_stop.
subroutine, public dbt_tas_put_block(matrix, row, col, block, summation)
As dbm_put_block.
subroutine, public dbt_tas_get_block_p(matrix, row, col, block, row_size, col_size)
As dbm_get_block_p.
DBT tall-and-skinny base types. Mostly wrappers around existing DBM routines.
DBT tensor framework for block-sparse tensor contraction: Types and create/destroy routines.
Definition dbt_types.F:12
subroutine, public dbt_blk_sizes(tensor, ind, blk_size)
Size of tensor block.
Definition dbt_types.F:1479
pure integer function, public ndims_tensor(tensor)
tensor rank
Definition dbt_types.F:1227
pure integer function, public dbt_get_num_blocks(tensor)
As block_get_num_blocks: get number of local blocks.
Definition dbt_types.F:1759
subroutine, public dbt_finalize(tensor)
Finalize tensor, as block_finalize. This should be taken care of internally in DBT tensors,...
Definition dbt_types.F:1790
pure integer(int_8) function, public ndims_matrix_row(tensor)
how many tensor dimensions are mapped to matrix row
Definition dbt_types.F:1204
pure integer(int_8) function, public ndims_matrix_column(tensor)
how many tensor dimensions are mapped to matrix column
Definition dbt_types.F:1216
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34