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