LCOV - code coverage report
Current view: top level - src/dbt - dbt_block.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 92.8 % 152 141
Test Date: 2025-07-25 12:55:17 Functions: 83.9 % 31 26

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Methods to operate on n-dimensional tensor blocks.
      10              : !> \author Patrick Seewald
      11              : ! **************************************************************************************************
      12              : MODULE dbt_block
      13              : 
      14              :    #:include "dbt_macros.fypp"
      15              :    #:set maxdim = maxrank
      16              :    #:set ndims = range(2,maxdim+1)
      17              : 
      18              :    USE OMP_LIB, ONLY: omp_get_thread_num, omp_get_num_threads
      19              :    USE cp_dbcsr_api, ONLY: &
      20              :       dbcsr_type, dbcsr_release, &
      21              :       dbcsr_iterator_type, dbcsr_iterator_start, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      22              :       dbcsr_has_symmetry, dbcsr_desymmetrize, dbcsr_get_num_blocks, dbcsr_iterator_stop, &
      23              :       dbcsr_reserve_blocks, dbcsr_finalize
      24              :    USE dbt_allocate_wrap, ONLY: &
      25              :       allocate_any
      26              :    USE dbt_tas_types, ONLY: &
      27              :       dbt_tas_iterator
      28              :    USE dbt_tas_base, ONLY: &
      29              :       dbt_tas_iterator_next_block, dbt_tas_iterator_blocks_left, dbt_tas_iterator_start, &
      30              :       dbt_tas_iterator_stop, dbt_tas_get_block_p, dbt_tas_put_block, dbt_tas_reserve_blocks, &
      31              :       dbt_tas_iterator_num_blocks
      32              :    USE kinds, ONLY: dp, int_8, dp
      33              :    USE dbt_index, ONLY: &
      34              :       nd_to_2d_mapping, ndims_mapping, get_nd_indices_tensor, destroy_nd_to_2d_mapping, get_2d_indices_tensor, &
      35              :       create_nd_to_2d_mapping
      36              :    USE dbt_array_list_methods, ONLY: &
      37              :       array_list, get_array_elements, destroy_array_list, sizes_of_arrays, create_array_list, &
      38              :       get_arrays
      39              :    USE dbt_types, ONLY: &
      40              :       dbt_type, ndims_tensor, dbt_blk_sizes, dbt_get_num_blocks, &
      41              :       dbt_finalize, ndims_matrix_row, ndims_matrix_column
      42              : 
      43              : #include "../base/base_uses.f90"
      44              : 
      45              :    IMPLICIT NONE
      46              :    PRIVATE
      47              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dbt_block'
      48              : 
      49              :    PUBLIC :: &
      50              :       block_nd, &
      51              :       create_block, &
      52              :       dbt_get_block, &
      53              :       dbt_iterator_num_blocks, &
      54              :       dbt_iterator_blocks_left, &
      55              :       dbt_iterator_next_block, &
      56              :       dbt_iterator_start, &
      57              :       dbt_iterator_stop, &
      58              :       dbt_iterator_type, &
      59              :       dbt_put_block, &
      60              :       dbt_reserve_blocks, &
      61              :       destroy_block, &
      62              :       checker_tr, &
      63              :       ndims_iterator
      64              : 
      65              :    TYPE dbt_iterator_type
      66              :       TYPE(dbt_tas_iterator)      :: iter
      67              :       TYPE(dbt_type), POINTER     :: tensor => NULL()
      68              :    END TYPE dbt_iterator_type
      69              : 
      70              :    TYPE block_nd
      71              :       INTEGER, DIMENSION(:), ALLOCATABLE   :: sizes
      72              :       REAL(dp), DIMENSION(:), ALLOCATABLE :: blk
      73              :    END TYPE
      74              : 
      75              :    INTERFACE create_block
      76              :       MODULE PROCEDURE create_block_data
      77              :       MODULE PROCEDURE create_block_nodata
      78              :    END INTERFACE
      79              : 
      80              :    INTERFACE dbt_put_block
      81              :       #:for ndim in ndims
      82              :          MODULE PROCEDURE dbt_put_${ndim}$d_block
      83              :       #:endfor
      84              :       MODULE PROCEDURE dbt_put_anyd_block
      85              :    END INTERFACE
      86              : 
      87              :    INTERFACE dbt_get_block
      88              :       #:for ndim in ndims
      89              :          MODULE PROCEDURE dbt_get_${ndim}$d_block
      90              :          MODULE PROCEDURE dbt_allocate_and_get_${ndim}$d_block
      91              :       #:endfor
      92              :       MODULE PROCEDURE dbt_get_anyd_block
      93              :    END INTERFACE
      94              : 
      95              :    INTERFACE dbt_reserve_blocks
      96              :       MODULE PROCEDURE dbt_reserve_blocks_index
      97              :       MODULE PROCEDURE dbt_reserve_blocks_index_array
      98              :       MODULE PROCEDURE dbt_reserve_blocks_template
      99              :       MODULE PROCEDURE dbt_reserve_blocks_tensor_to_matrix
     100              :       MODULE PROCEDURE dbt_reserve_blocks_matrix_to_tensor
     101              :    END INTERFACE
     102              : 
     103              : CONTAINS
     104              : 
     105              : ! **************************************************************************************************
     106              : !> \brief block size
     107              : !> \author Patrick Seewald
     108              : ! **************************************************************************************************
     109            0 :    FUNCTION block_size(block)
     110              :       TYPE(block_nd), INTENT(IN)         :: block
     111              :       INTEGER, ALLOCATABLE, DIMENSION(:) :: block_size
     112              : 
     113            0 :       ALLOCATE (block_size, source=block%sizes)
     114              :    END FUNCTION
     115              : 
     116              : ! **************************************************************************************************
     117              : !> \brief Generalization of block_iterator_start for tensors.
     118              : !> \author Patrick Seewald
     119              : ! **************************************************************************************************
     120      1793122 :    SUBROUTINE dbt_iterator_start(iterator, tensor)
     121              :       TYPE(dbt_iterator_type), INTENT(OUT)           :: iterator
     122              :       TYPE(dbt_type), INTENT(IN), TARGET             :: tensor
     123              : 
     124      1793122 :       CPASSERT(tensor%valid)
     125      1793122 :       CALL dbt_tas_iterator_start(iterator%iter, tensor%matrix_rep)
     126      1793122 :       iterator%tensor => tensor
     127      1793122 :    END SUBROUTINE
     128              : 
     129              : ! **************************************************************************************************
     130              : !> \brief Generalization of block_iterator_stop for tensors.
     131              : !> \author Patrick Seewald
     132              : ! **************************************************************************************************
     133      1793122 :    SUBROUTINE dbt_iterator_stop(iterator)
     134              :       TYPE(dbt_iterator_type), INTENT(INOUT) :: iterator
     135              : 
     136      1793122 :       CALL dbt_tas_iterator_stop(iterator%iter)
     137      1793122 :    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     36684840 :    PURE FUNCTION ndims_iterator(iterator)
     146              :       TYPE(dbt_iterator_type), INTENT(IN) :: iterator
     147              :       INTEGER                                 :: ndims_iterator
     148              : 
     149     36684840 :       ndims_iterator = iterator%tensor%nd_index%ndim_nd
     150     36684840 :    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     35642512 :    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     36684840 :       CALL dbt_tas_iterator_next_block(iterator%iter, ind_2d(1), ind_2d(2))
     171              : 
     172     36684840 :       ind_nd(:) = get_nd_indices_tensor(iterator%tensor%nd_index_blk, ind_2d)
     173     36684840 :       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     36684840 :       IF (PRESENT(blk_offset)) blk_offset(:) = get_array_elements(iterator%tensor%blk_offsets, ind_nd)
     177              : 
     178     59948259 :    END SUBROUTINE
     179              : 
     180              : ! **************************************************************************************************
     181              : !> \brief Generalization of block_iterator_num_blocks for tensors.
     182              : !> \author Ole Schuett
     183              : ! **************************************************************************************************
     184       256131 :    FUNCTION dbt_iterator_num_blocks(iterator)
     185              :       TYPE(dbt_iterator_type), INTENT(IN) :: iterator
     186              :       INTEGER                             :: dbt_iterator_num_blocks
     187              : 
     188       256131 :       dbt_iterator_num_blocks = dbt_tas_iterator_num_blocks(iterator%iter)
     189              : 
     190       256131 :    END FUNCTION
     191              : 
     192              : ! **************************************************************************************************
     193              : !> \brief Generalization of block_iterator_blocks_left for tensors.
     194              : !> \author Patrick Seewald
     195              : ! **************************************************************************************************
     196     37887163 :    FUNCTION dbt_iterator_blocks_left(iterator)
     197              :       TYPE(dbt_iterator_type), INTENT(IN) :: iterator
     198              :       LOGICAL                                 :: dbt_iterator_blocks_left
     199              : 
     200     37887163 :       dbt_iterator_blocks_left = dbt_tas_iterator_blocks_left(iterator%iter)
     201              : 
     202     37887163 :    END FUNCTION
     203              : 
     204              : ! **************************************************************************************************
     205              : !> \brief reserve blocks from indices as array object
     206              : !> \author Patrick Seewald
     207              : ! **************************************************************************************************
     208       592893 :    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       592893 :       CALL timeset(routineN, handle)
     215              :       #:for ndim in ndims
     216      1185494 :          IF (ndims_tensor(tensor) == ${ndim}$) THEN
     217       592893 :             CALL dbt_reserve_blocks(tensor, ${arrlist("blk_ind", nmax=ndim, ndim_pre=1)}$)
     218              :          END IF
     219              :       #:endfor
     220       592893 :       CALL timestop(handle)
     221              : 
     222       592893 :    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       679225 :    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       679225 :       TYPE(array_list)                            :: blks
     236      1358450 :       INTEGER, DIMENSION(ndims_tensor(tensor))   :: iblk_nd, ind_nd, nblk_tmp
     237              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'dbt_reserve_blocks_index'
     238              : 
     239       679225 :       CALL timeset(routineN, handle)
     240       679225 :       CPASSERT(tensor%valid)
     241              : 
     242              :       CALL create_array_list(blks, ndims_tensor(tensor), &
     243      1637840 :                              ${varlist("blk_ind")}$)
     244              : 
     245      2437510 :       nblk_tmp(:) = sizes_of_arrays(blks)
     246       679225 :       nblk = nblk_tmp(1)
     247      2471582 :       ALLOCATE (cols(nblk), rows(nblk))
     248     14546476 :       DO iblk = 1, nblk
     249     52271761 :          iblk_nd(:) = iblk
     250     13867251 :          ind_nd(:) = get_array_elements(blks, iblk_nd)
     251     13867251 :          ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind_nd)
     252     14546476 :          rows(iblk) = ind_2d(1); cols(iblk) = ind_2d(2)
     253              :       END DO
     254              : 
     255       679225 :       CALL dbt_tas_reserve_blocks(tensor%matrix_rep, rows=rows, columns=cols)
     256       679225 :       CALL dbt_finalize(tensor)
     257       679225 :       CALL timestop(handle)
     258      1358450 :    END SUBROUTINE
     259              : 
     260              : ! **************************************************************************************************
     261              : !> \brief reserve tensor blocks using template
     262              : !> \param tensor_in template tensor
     263              : !> \author Patrick Seewald
     264              : ! **************************************************************************************************
     265        11182 :    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        11182 :       INTEGER, DIMENSION(:, :), ALLOCATABLE :: blk_ind
     274              : 
     275        11182 :       CALL timeset(routineN, handle)
     276              : 
     277              : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,tensor_out) &
     278        11182 : !$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        11182 :       CALL timestop(handle)
     292        22364 :    END SUBROUTINE
     293              : 
     294              : ! **************************************************************************************************
     295              : !> \brief reserve tensor blocks using matrix template
     296              : !> \author Patrick Seewald
     297              : ! **************************************************************************************************
     298        76478 :    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                            :: iblk, nblk, nblk_per_thread, a, b
     304        76478 :       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        76478 :       CALL timeset(routineN, handle)
     311              : 
     312        76478 :       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        76478 :       nblk = dbcsr_get_num_blocks(matrix_in_desym)
     320       363970 :       ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
     321        76478 :       CALL dbcsr_iterator_start(iter, matrix_in_desym)
     322       330975 :       DO iblk = 1, nblk
     323       254497 :          CALL dbcsr_iterator_next_block(iter, ind_2d(1), ind_2d(2))
     324       330975 :          blk_ind_1(iblk) = ind_2d(1); blk_ind_2(iblk) = ind_2d(2)
     325              :       END DO
     326        76478 :       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        76478 : !$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        76478 :       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        76478 :       CALL timestop(handle)
     343       229434 :    END SUBROUTINE
     344              : 
     345              : ! **************************************************************************************************
     346              : !> \brief reserve matrix blocks using tensor template
     347              : !> \author Patrick Seewald
     348              : ! **************************************************************************************************
     349        48090 :    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        48090 :       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        48090 :       CALL timeset(routineN, handle)
     360              : 
     361        48090 :       nblk = dbt_get_num_blocks(tensor_in)
     362       170206 :       ALLOCATE (blk_ind_1(nblk), blk_ind_2(nblk))
     363        48090 :       iblk = 0
     364              : 
     365              : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor_in,matrix_out,iblk,blk_ind_1,blk_ind_2) &
     366        48090 : !$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        48090 :       CALL dbcsr_reserve_blocks(matrix_out, blk_ind_1(:iblk), blk_ind_2(:iblk))
     384        48090 :       CALL dbcsr_finalize(matrix_out)
     385              : 
     386        48090 :       CALL timestop(handle)
     387        96180 :    END SUBROUTINE
     388              : 
     389              : ! **************************************************************************************************
     390              : !> \brief Swaps two integers
     391              : !> \author Patrick Seewald
     392              : ! **************************************************************************************************
     393         4659 :    ELEMENTAL SUBROUTINE swap(a, b)
     394              :       INTEGER, INTENT(INOUT)                             :: a, b
     395              :       INTEGER                                            :: tmp
     396              : 
     397         4659 :       tmp = a
     398         4659 :       a = b
     399         4659 :       b = tmp
     400         4659 :    END SUBROUTINE swap
     401              : 
     402              : ! **************************************************************************************************
     403              : !> \brief Create block from array, array can be n-dimensional.
     404              : !> \author Patrick Seewald
     405              : ! **************************************************************************************************
     406     28719792 :    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     44001792 :       ALLOCATE (block%sizes, source=sizes)
     412   2635746587 :       ALLOCATE (block%blk, source=array)
     413      7641000 :    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      7638072 :    SUBROUTINE destroy_block(block)
     431              :       TYPE(block_nd), INTENT(INOUT) :: block
     432      7638072 :       DEALLOCATE (block%blk)
     433      7638072 :       DEALLOCATE (block%sizes)
     434      7638072 :    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        39774 :    ELEMENTAL FUNCTION checker_tr(row, column) RESULT(transpose)
     449              :       INTEGER, INTENT(IN)                                :: row, column
     450              :       LOGICAL                                            :: transpose
     451              : 
     452        39774 :       transpose = BTEST(column + row, 0) .EQV. column .GE. row
     453              : 
     454        39774 :    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      3891407 :    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       924840 :       SELECT CASE (ndims_tensor(tensor))
     476              :          #:for ndim in ndims
     477              :             CASE (${ndim}$)
     478      3891407 :             CALL dbt_put_${ndim}$d_block(tensor, ind, block%sizes, block%blk, summation=summation)
     479              :          #:endfor
     480              :       END SELECT
     481      3891407 :    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      3894335 :    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      7788670 :       INTEGER, DIMENSION(ndims_tensor(tensor))    :: blk_size
     502      3894335 :       REAL(dp), DIMENSION(:), ALLOCATABLE         :: block_arr
     503              : 
     504      3894335 :       CALL dbt_blk_sizes(tensor, ind, blk_size)
     505     22443721 :       ALLOCATE (block_arr(PRODUCT(blk_size)))
     506              : 
     507       924888 :       SELECT CASE (ndims_tensor(tensor))
     508              :          #:for ndim in ndims
     509              :             CASE (${ndim}$)
     510      3894335 :             CALL dbt_get_${ndim}$d_block(tensor, ind, blk_size, block_arr, found)
     511              :          #:endfor
     512              :       END SELECT
     513      3894335 :       CALL create_block(block, blk_size, block_arr)
     514      3894335 :    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     18651629 :       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     18651629 :          REAL(dp), POINTER, DIMENSION(:, :)                   :: block_2d
     539              :          INTEGER, DIMENSION(${ndim}$)                          :: shape_nd
     540              :          LOGICAL :: found, new_block
     541     18651629 :          REAL(dp), DIMENSION(${arrlist("sizes", nmax=ndim)}$) :: block_check
     542              : 
     543              :          LOGICAL, PARAMETER :: debug = .FALSE.
     544              :          INTEGER :: i
     545              : 
     546     18651629 :          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    108482812 :             shape_2d = [PRODUCT(sizes(map1_2d)), PRODUCT(sizes(map2_2d))]
     558              : 
     559    160841270 :             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     18578827 :                block_2d(1:shape_2d(1), 1:shape_2d(2)) => block(${shape_colon(ndim)}$)
     562              :             ELSE
     563              :                ! need reshape due to rank reordering
     564       291208 :                ALLOCATE (block_2d(shape_2d(1), shape_2d(2)))
     565        72802 :                new_block = .TRUE.
     566       286789 :                shape_nd(map_nd) = sizes
     567       646380 :                block_2d(:, :) = RESHAPE(RESHAPE(block, SHAPE=shape_nd, order=map_nd), SHAPE=shape_2d)
     568              :             END IF
     569              : 
     570     37303258 :             ind_2d(:) = get_2d_indices_tensor(tensor%nd_index_blk, ind)
     571              : 
     572              :          END ASSOCIATE
     573              : 
     574     18651629 :          CALL dbt_tas_put_block(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d, summation=summation)
     575              : 
     576     18651629 :          IF (new_block) DEALLOCATE (block_2d)
     577              : 
     578     18651629 :       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     17450718 :       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     17450718 :          CALL dbt_blk_sizes(tensor, ind, blk_size)
     603     17450718 :          CALL allocate_any(block, shape_spec=blk_size)
     604     17450718 :          CALL dbt_get_${ndim}$d_block(tensor, ind, blk_size, block, found)
     605              : 
     606     17450718 :       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     21345351 :       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     21345351 :          REAL(dp), DIMENSION(:, :), POINTER, CONTIGUOUS       :: block_2d_ptr
     628              :          INTEGER                                               :: i
     629     21345351 :          REAL(dp), DIMENSION(${shape_colon(ndim)}$), POINTER  :: block_ptr
     630              : 
     631     21345351 :          NULLIFY (block_2d_ptr)
     632              : 
     633     21345351 :          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     21345351 :             CALL dbt_tas_get_block_p(tensor%matrix_rep, ind_2d(1), ind_2d(2), block_2d_ptr)
     639     21345351 :             found = ASSOCIATED(block_2d_ptr)
     640              : 
     641     21345351 :             IF (found) THEN
     642    165725282 :                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     49726210 :                   block_ptr(${shape_explicit('block', ndim)}$) => block_2d_ptr(:, :)
     645   9125721172 :                   block(${shape_colon(ndim)}$) = block_ptr(${shape_colon(ndim)}$)
     646              :                ELSE
     647              :                   ! need reshape due to rank reordering
     648     25965606 :                   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     21345351 :       END SUBROUTINE
     655              :    #:endfor
     656              : 
     657            0 : END MODULE
        

Generated by: LCOV version 2.0-1