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 various cholesky decomposition related routines
      10             : !> \par History
      11             : !>      12.2002 Moved routines from cp_cfm_basic_linalg to this new module [Rocco Meli]
      12             : ! **************************************************************************************************
      13             : MODULE cp_cfm_cholesky
      14             :    USE cp_cfm_types, ONLY: cp_cfm_type
      15             :    USE kinds, ONLY: dp
      16             : #if defined(__DLAF)
      17             :    USE cp_cfm_dlaf_api, ONLY: cp_cfm_pzpotrf_dlaf
      18             :    USE cp_dlaf_utils_api, ONLY: cp_dlaf_initialize, cp_dlaf_create_grid
      19             :    USE cp_fm_cholesky, ONLY: cholesky_type, dlaf_cholesky_n_min, FM_CHOLESKY_TYPE_DLAF
      20             : #endif
      21             : 
      22             : #include "../base/base_uses.f90"
      23             : 
      24             :    IMPLICIT NONE
      25             :    PRIVATE
      26             : 
      27             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      28             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_cfm_cholesky'
      29             : 
      30             :    PUBLIC :: cp_cfm_cholesky_decompose, &
      31             :              cp_cfm_cholesky_invert
      32             : 
      33             : ! **************************************************************************************************
      34             : 
      35             : CONTAINS
      36             : 
      37             : ! **************************************************************************************************
      38             : !> \brief Used to replace a symmetric positive definite matrix M with its Cholesky
      39             : !>      decomposition U: M = U^T * U, with U upper triangular.
      40             : !> \param matrix   the matrix to replace with its Cholesky decomposition
      41             : !> \param n        the number of row (and columns) of the matrix &
      42             : !>                 (defaults to the min(size(matrix)))
      43             : !> \param info_out if present, outputs info from (p)zpotrf
      44             : !> \par History
      45             : !>      05.2002 created [JVdV]
      46             : !>      12.2002 updated, added n optional parm [fawzi]
      47             : !>      09.2021 removed CPASSERT(info == 0) since there is already check of info [Jan Wilhelm]
      48             : !>      12.2024 Added DLA-Future support [Rocco Meli]
      49             : !>      03.2025 Moved DLA-Future support to cp_cfm_cholesky_decompose_dlaf [Rocco Meli]
      50             : !> \author Joost
      51             : ! **************************************************************************************************
      52       20994 :    SUBROUTINE cp_cfm_cholesky_decompose(matrix, n, info_out)
      53             :       TYPE(cp_cfm_type), INTENT(IN)                      :: matrix
      54             :       INTEGER, INTENT(in), OPTIONAL                      :: n
      55             :       INTEGER, INTENT(out), OPTIONAL                     :: info_out
      56             : 
      57             :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cholesky_decompose'
      58             : 
      59       20994 :       COMPLEX(kind=dp), DIMENSION(:, :), POINTER         :: a
      60             :       INTEGER                                            :: handle, info, my_n
      61             : #if defined(__parallel)
      62             :       INTEGER, DIMENSION(9)                              :: desca
      63             : #else
      64             :       INTEGER                                            :: lda
      65             : #endif
      66             : 
      67       20994 :       CALL timeset(routineN, handle)
      68             : 
      69             :       my_n = MIN(matrix%matrix_struct%nrow_global, &
      70       20994 :                  matrix%matrix_struct%ncol_global)
      71       20994 :       IF (PRESENT(n)) THEN
      72        4506 :          CPASSERT(n <= my_n)
      73        4506 :          my_n = n
      74             :       END IF
      75             : 
      76       20994 :       a => matrix%local_data
      77             : 
      78             : #if defined(__parallel)
      79      209940 :       desca(:) = matrix%matrix_struct%descriptor(:)
      80             : #if defined(__DLAF)
      81             :       IF (cholesky_type == FM_CHOLESKY_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global .GE. dlaf_cholesky_n_min) THEN
      82             :          ! Initialize DLA-Future on-demand; if already initialized, does nothing
      83             :          CALL cp_dlaf_initialize()
      84             : 
      85             :          ! Create DLAF grid from BLACS context; if already present, does nothing
      86             :          CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
      87             : 
      88             :          CALL cp_cfm_pzpotrf_dlaf('U', my_n, a(:, :), 1, 1, desca, info)
      89             :       ELSE
      90             : #endif
      91       20994 :          CALL pzpotrf('U', my_n, a(1, 1), 1, 1, desca, info)
      92             : #if defined(__DLAF)
      93             :       END IF
      94             : #endif
      95             : 
      96             : #else
      97             :       lda = SIZE(a, 1)
      98             :       CALL zpotrf('U', my_n, a(1, 1), lda, info)
      99             : #endif
     100             : 
     101       20994 :       IF (PRESENT(info_out)) THEN
     102        4506 :          info_out = info
     103             :       ELSE
     104       16488 :          IF (info /= 0) &
     105             :             CALL cp_abort(__LOCATION__, &
     106           0 :                           "Cholesky decompose failed: matrix is not positive definite or ill-conditioned")
     107             :       END IF
     108             : 
     109       20994 :       CALL timestop(handle)
     110             : 
     111       20994 :    END SUBROUTINE cp_cfm_cholesky_decompose
     112             : 
     113             : ! **************************************************************************************************
     114             : !> \brief Used to replace Cholesky decomposition by the inverse.
     115             : !> \param matrix : the matrix to invert (must be an upper triangular matrix),
     116             : !>                 and is the output of Cholesky decomposition
     117             : !> \param n : size of the matrix to invert (defaults to the min(size(matrix)))
     118             : !> \param info_out : if present, outputs info of (p)zpotri
     119             : !> \par History
     120             : !>      05.2002 created Lianheng Tong, based on cp_fm_cholesky_invert
     121             : !> \author Lianheng Tong
     122             : ! **************************************************************************************************
     123        4602 :    SUBROUTINE cp_cfm_cholesky_invert(matrix, n, info_out)
     124             :       TYPE(cp_cfm_type), INTENT(IN)              :: matrix
     125             :       INTEGER, INTENT(in), OPTIONAL              :: n
     126             :       INTEGER, INTENT(out), OPTIONAL             :: info_out
     127             : 
     128             :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_cfm_cholesky_invert'
     129        4602 :       COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: aa
     130             :       INTEGER                                    :: info, handle
     131             :       INTEGER                                    :: my_n
     132             : #if defined(__parallel)
     133             :       INTEGER, DIMENSION(9)                      :: desca
     134             : #endif
     135             : 
     136        4602 :       CALL timeset(routineN, handle)
     137             : 
     138             :       my_n = MIN(matrix%matrix_struct%nrow_global, &
     139        4602 :                  matrix%matrix_struct%ncol_global)
     140        4602 :       IF (PRESENT(n)) THEN
     141           0 :          CPASSERT(n <= my_n)
     142           0 :          my_n = n
     143             :       END IF
     144             : 
     145        4602 :       aa => matrix%local_data
     146             : 
     147             : #if defined(__parallel)
     148       46020 :       desca = matrix%matrix_struct%descriptor
     149        4602 :       CALL pzpotri('U', my_n, aa(1, 1), 1, 1, desca, info)
     150             : #else
     151             :       CALL zpotri('U', my_n, aa(1, 1), SIZE(aa, 1), info)
     152             : #endif
     153             : 
     154        4602 :       IF (PRESENT(info_out)) THEN
     155           0 :          info_out = info
     156             :       ELSE
     157        4602 :          IF (info /= 0) &
     158             :             CALL cp_abort(__LOCATION__, &
     159           0 :                           "Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
     160             :       END IF
     161             : 
     162        4602 :       CALL timestop(handle)
     163             : 
     164        4602 :    END SUBROUTINE cp_cfm_cholesky_invert
     165             : 
     166             : END MODULE cp_cfm_cholesky