(git:ccc2433)
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 ! **************************************************************************************************
12 MODULE dbt_block
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: &
23  allocate_any
24  USE dbt_tas_types, ONLY: &
25  dbt_tas_iterator
26  USE dbt_tas_base, ONLY: &
27  dbt_tas_iterator_next_block, dbt_tas_iterator_blocks_left, dbt_tas_iterator_start, &
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, &
49  create_block, &
50  dbt_get_block, &
56  dbt_iterator_type, &
57  dbt_put_block, &
58  dbt_reserve_blocks, &
59  destroy_block, &
60  checker_tr, &
62 
63  TYPE dbt_iterator_type
64  TYPE(dbt_tas_iterator) :: iter
65  TYPE(dbt_type), POINTER :: tensor => null()
66  END TYPE dbt_iterator_type
67 
68  TYPE block_nd
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 
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  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 
925 END 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...
Definition: dbt_tas_base.F:13
integer function, public dbt_tas_iterator_num_blocks(iter)
As dbm_iterator_num_blocks.
Definition: dbt_tas_base.F:685
subroutine, public dbt_tas_iterator_start(iter, matrix_in)
As dbm_iterator_start.
Definition: dbt_tas_base.F:670
logical function, public dbt_tas_iterator_blocks_left(iter)
As dbm_iterator_blocks_left.
Definition: dbt_tas_base.F:698
subroutine, public dbt_tas_iterator_stop(iter)
As dbm_iterator_stop.
Definition: dbt_tas_base.F:710
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.
Definition: dbt_tas_types.F:13
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