gw_utils_dbcsr.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Common DBCSR matrix operations used by GW modules.
!> \par History
!>      09.2026 created Jan Wilhelm
! **************************************************************************************************
MODULE gw_utils_dbcsr
   USE cp_dbcsr_api,                    ONLY: &
        dbcsr_copy, dbcsr_create, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
        dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
        dbcsr_multiply, dbcsr_release, dbcsr_type
   USE kinds,                           ONLY: dp
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   PUBLIC :: dbcsr_contract_aba, hadamard_product, hadamard_product_inplace
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils_dbcsr'
 
CONTAINS
 
! **************************************************************************************************
!> \brief Computes the scaled element-wise product C = factor (A ◦ B) while preserving the
!>        block structure of A. Blocks absent from B are retained in C with zero values.
!> \param matrix_A First factor and source of the block structure
!> \param matrix_B Second factor
!> \param matrix_C Scaled element-wise product
!> \param factor Scaling factor
! **************************************************************************************************
   SUBROUTINE hadamard_product(matrix_A, matrix_B, matrix_C, factor)
      TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_a, matrix_b, matrix_c
      REAL(kind=dp), INTENT(IN)                          :: factor
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'hadamard_product'
 
      INTEGER                                            :: handle, icol, irow
      LOGICAL                                            :: found
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: block_b, block_c
      TYPE(dbcsr_iterator_type)                          :: iterator
 
      CALL timeset(routinen, handle)
 
      CALL dbcsr_copy(matrix_c, matrix_a)
      CALL dbcsr_iterator_start(iterator, matrix_c)
      DO WHILE (dbcsr_iterator_blocks_left(iterator))
         CALL dbcsr_iterator_next_block(iterator, irow, icol, block_c)
         CALL dbcsr_get_block_p(matrix_b, irow, icol, block_b, found)
         IF (found) THEN
            block_c(:, :) = factor*block_c(:, :)*block_b(:, :)
         ELSE
            block_c(:, :) = 0.0_dp
         END IF
      END DO
      CALL dbcsr_iterator_stop(iterator)
 
      CALL timestop(handle)
 
   END SUBROUTINE hadamard_product
 
! **************************************************************************************************
!> \brief Form A = factor * (A element-wise B) without changing A's block structure.
!> \param matrix_A First factor, overwritten by the product.
!> \param matrix_B Second factor; a missing block represents zero.
!> \param factor Product scale factor.
! **************************************************************************************************
   SUBROUTINE hadamard_product_inplace(matrix_A, matrix_B, factor)
      TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_a, matrix_b
      REAL(kind=dp), INTENT(IN)                          :: factor
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'hadamard_product_inplace'
 
      INTEGER                                            :: handle, icol, irow
      LOGICAL                                            :: found
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: block_a, block_b
      TYPE(dbcsr_iterator_type)                          :: iterator
 
      CALL timeset(routinen, handle)
 
      CALL dbcsr_iterator_start(iterator, matrix_a)
      DO WHILE (dbcsr_iterator_blocks_left(iterator))
         CALL dbcsr_iterator_next_block(iterator, irow, icol, block_a)
         CALL dbcsr_get_block_p(matrix_b, irow, icol, block_b, found)
         IF (found) THEN
            block_a(:, :) = factor*block_a(:, :)*block_b(:, :)
         ELSE
            block_a(:, :) = 0.0_dp
         END IF
      END DO
      CALL dbcsr_iterator_stop(iterator)
 
      CALL timestop(handle)
   END SUBROUTINE hadamard_product_inplace
 
! **************************************************************************************************
!> \brief Computes C=A B A^T or C=A^T B A for DBCSR matrices.
!> \param trans_A_left transposition applied to the left occurrence of A
!> \param trans_A_right transposition applied to the right occurrence of A
!> \param matrix_A left and right matrix A
!> \param matrix_B input matrix B
!> \param matrix_C output matrix C
!> \param eps_filter filtering threshold for both matrix multiplications
!> \param retain_sparsity if true, only existing blocks of C are filled
! **************************************************************************************************
   SUBROUTINE dbcsr_contract_aba(trans_A_left, trans_A_right, matrix_A, matrix_B, matrix_C, &
                                 eps_filter, retain_sparsity)
      CHARACTER(LEN=1), INTENT(IN)                       :: trans_a_left, trans_a_right
      TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_a, matrix_b, matrix_c
      REAL(kind=dp), INTENT(IN)                          :: eps_filter
      LOGICAL, INTENT(IN), OPTIONAL                      :: retain_sparsity
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'dbcsr_contract_ABA'
 
      INTEGER                                            :: handle
      LOGICAL                                            :: my_retain_sparsity
      TYPE(dbcsr_type)                                   :: work
 
      CALL timeset(routinen, handle)
 
      my_retain_sparsity = .false.
      IF (PRESENT(retain_sparsity)) my_retain_sparsity = retain_sparsity
 
      CALL dbcsr_create(work, template=matrix_a)
 
      IF (trans_a_left == "N" .AND. trans_a_right == "T") THEN
         ! C = A B A^T
         CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_a, matrix_b, &
                             0.0_dp, work, filter_eps=eps_filter)
         CALL dbcsr_multiply("N", "T", 1.0_dp, work, matrix_a, &
                             0.0_dp, matrix_c, filter_eps=eps_filter, &
                             retain_sparsity=my_retain_sparsity)
      ELSE IF (trans_a_left == "T" .AND. trans_a_right == "N") THEN
         ! C = A^T B A
         CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_b, matrix_a, &
                             0.0_dp, work, filter_eps=eps_filter)
         CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_a, work, &
                             0.0_dp, matrix_c, filter_eps=eps_filter, &
                             retain_sparsity=my_retain_sparsity)
      ELSE
         cpabort("Unsupported transposition pair in dbcsr_contract_ABA")
      END IF
      CALL dbcsr_release(work)
 
      CALL timestop(handle)
 
   END SUBROUTINE dbcsr_contract_aba
 
END MODULE gw_utils_dbcsr