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             : !>      09.2002 created [fawzi]
      12             : !> \author Fawzi Mohamed
      13             : ! **************************************************************************************************
      14             : MODULE cp_fm_cholesky
      15             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      16             :    USE cp_dlaf_utils_api,               ONLY: cp_dlaf_create_grid,&
      17             :                                               cp_dlaf_initialize
      18             :    USE cp_fm_dlaf_api,                  ONLY: cp_pdpotrf_dlaf,&
      19             :                                               cp_pspotrf_dlaf
      20             :    USE cp_fm_types,                     ONLY: cp_fm_type
      21             :    USE cp_log_handling,                 ONLY: cp_to_string
      22             :    USE kinds,                           ONLY: dp,&
      23             :                                               sp
      24             : #include "../base/base_uses.f90"
      25             : 
      26             :    IMPLICIT NONE
      27             :    PRIVATE
      28             : 
      29             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      30             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_cholesky'
      31             : 
      32             :    PUBLIC :: cp_fm_cholesky_decompose, cp_fm_cholesky_invert, &
      33             :              cp_fm_cholesky_reduce, cp_fm_cholesky_restore
      34             : 
      35             :    ! The following saved variables are Cholesky decomposition global
      36             :    ! Stores the default library for Cholesky decomposition
      37             :    INTEGER, SAVE, PUBLIC   :: cholesky_type = 0
      38             :    ! Minimum matrix size for the use of the DLAF Cholesky decomposition.
      39             :    ! ScaLAPACK is used as fallback for all smaller cases.
      40             :    INTEGER, SAVE, PUBLIC    :: dlaf_cholesky_n_min = 0
      41             :    ! Constants for the diag_type above
      42             :    INTEGER, PARAMETER, PUBLIC  :: FM_CHOLESKY_TYPE_SCALAPACK = 101, &
      43             :                                   FM_CHOLESKY_TYPE_DLAF = 104
      44             :    INTEGER, PARAMETER, PUBLIC :: FM_CHOLESKY_TYPE_DEFAULT = FM_CHOLESKY_TYPE_SCALAPACK
      45             : 
      46             : !***
      47             : CONTAINS
      48             : 
      49             : ! **************************************************************************************************
      50             : !> \brief used to replace a symmetric positive def. matrix M with its cholesky
      51             : !>      decomposition U: M = U^T * U, with U upper triangular
      52             : !> \param matrix the matrix to replace with its cholesky decomposition
      53             : !> \param n the number of row (and columns) of the matrix &
      54             : !>        (defaults to the min(size(matrix)))
      55             : !> \param info_out ...
      56             : !> \par History
      57             : !>      05.2002 created [JVdV]
      58             : !>      12.2002 updated, added n optional parm [fawzi]
      59             : !> \author Joost
      60             : ! **************************************************************************************************
      61       55923 :    SUBROUTINE cp_fm_cholesky_decompose(matrix, n, info_out)
      62             :       TYPE(cp_fm_type), INTENT(IN)             :: matrix
      63             :       INTEGER, INTENT(in), OPTIONAL            :: n
      64             :       INTEGER, INTENT(out), OPTIONAL           :: info_out
      65             : 
      66             :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_decompose'
      67             : 
      68             :       INTEGER                                  :: handle, info, my_n
      69       55923 :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a
      70       55923 :       REAL(KIND=sp), DIMENSION(:, :), POINTER  :: a_sp
      71             : #if defined(__parallel)
      72             :       INTEGER, DIMENSION(9)                    :: desca
      73             : #endif
      74             : 
      75       55923 :       CALL timeset(routineN, handle)
      76             : 
      77             :       my_n = MIN(matrix%matrix_struct%nrow_global, &
      78       55923 :                  matrix%matrix_struct%ncol_global)
      79       55923 :       IF (PRESENT(n)) THEN
      80       15438 :          CPASSERT(n <= my_n)
      81       15438 :          my_n = n
      82             :       END IF
      83             : 
      84       55923 :       a => matrix%local_data
      85       55923 :       a_sp => matrix%local_data_sp
      86             : 
      87             : #if defined(__parallel)
      88      559230 :       desca(:) = matrix%matrix_struct%descriptor(:)
      89             : #if defined(__DLAF)
      90             :       IF (cholesky_type == FM_CHOLESKY_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global >= dlaf_cholesky_n_min) THEN
      91             :          ! Initialize DLA-Future on-demand; if already initialized, does nothing
      92             :          CALL cp_dlaf_initialize()
      93             : 
      94             :          ! Create DLAF grid from BLACS context; if already present, does nothing
      95             :          CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
      96             : 
      97             :          IF (matrix%use_sp) THEN
      98             :             CALL cp_pspotrf_dlaf('U', my_n, a_sp(:, :), 1, 1, desca, info)
      99             :          ELSE
     100             :             CALL cp_pdpotrf_dlaf('U', my_n, a(:, :), 1, 1, desca, info)
     101             :          END IF
     102             :       ELSE
     103             : #endif
     104       55923 :          IF (matrix%use_sp) THEN
     105           0 :             CALL pspotrf('U', my_n, a_sp(1, 1), 1, 1, desca, info)
     106             :          ELSE
     107       55923 :             CALL pdpotrf('U', my_n, a(1, 1), 1, 1, desca, info)
     108             :          END IF
     109             : #if defined(__DLAF)
     110             :       END IF
     111             : #endif
     112             : #else
     113             :       IF (matrix%use_sp) THEN
     114             :          CALL spotrf('U', my_n, a_sp(1, 1), SIZE(a_sp, 1), info)
     115             :       ELSE
     116             :          CALL dpotrf('U', my_n, a(1, 1), SIZE(a, 1), info)
     117             :       END IF
     118             : #endif
     119             : 
     120       55923 :       IF (PRESENT(info_out)) THEN
     121       20442 :          info_out = info
     122             :       ELSE
     123       35481 :          IF (info /= 0) &
     124             :             CALL cp_abort(__LOCATION__, &
     125           0 :                           "Cholesky decompose failed: the matrix is not positive definite or ill-conditioned.")
     126             :       END IF
     127             : 
     128       55923 :       CALL timestop(handle)
     129             : 
     130       55923 :    END SUBROUTINE cp_fm_cholesky_decompose
     131             : 
     132             : ! **************************************************************************************************
     133             : !> \brief used to replace the cholesky decomposition by the inverse
     134             : !> \param matrix the matrix to invert (must be an upper triangular matrix)
     135             : !> \param n size of the matrix to invert (defaults to the min(size(matrix)))
     136             : !> \param info_out ...
     137             : !> \par History
     138             : !>      05.2002 created [JVdV]
     139             : !> \author Joost VandeVondele
     140             : ! **************************************************************************************************
     141       17497 :    SUBROUTINE cp_fm_cholesky_invert(matrix, n, info_out)
     142             :       TYPE(cp_fm_type), INTENT(IN)           :: matrix
     143             :       INTEGER, INTENT(in), OPTIONAL             :: n
     144             :       INTEGER, INTENT(OUT), OPTIONAL            :: info_out
     145             : 
     146             :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_invert'
     147       17497 :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a
     148       17497 :       REAL(KIND=sp), DIMENSION(:, :), POINTER  :: a_sp
     149             :       INTEGER                                   :: info, handle
     150             :       INTEGER                                   :: my_n
     151             : #if defined(__parallel)
     152             :       INTEGER, DIMENSION(9)                     :: desca
     153             : #endif
     154             : 
     155       17497 :       CALL timeset(routineN, handle)
     156             : 
     157             :       my_n = MIN(matrix%matrix_struct%nrow_global, &
     158       17497 :                  matrix%matrix_struct%ncol_global)
     159       17497 :       IF (PRESENT(n)) THEN
     160           0 :          CPASSERT(n <= my_n)
     161           0 :          my_n = n
     162             :       END IF
     163             : 
     164       17497 :       a => matrix%local_data
     165       17497 :       a_sp => matrix%local_data_sp
     166             : 
     167             : #if defined(__parallel)
     168             : 
     169      174970 :       desca(:) = matrix%matrix_struct%descriptor(:)
     170             : 
     171       17497 :       IF (matrix%use_sp) THEN
     172           0 :          CALL pspotri('U', my_n, a_sp(1, 1), 1, 1, desca, info)
     173             :       ELSE
     174       17497 :          CALL pdpotri('U', my_n, a(1, 1), 1, 1, desca, info)
     175             :       END IF
     176             : 
     177             : #else
     178             : 
     179             :       IF (matrix%use_sp) THEN
     180             :          CALL spotri('U', my_n, a_sp(1, 1), SIZE(a_sp, 1), info)
     181             :       ELSE
     182             :          CALL dpotri('U', my_n, a(1, 1), SIZE(a, 1), info)
     183             :       END IF
     184             : 
     185             : #endif
     186             : 
     187       17497 :       IF (PRESENT(info_out)) THEN
     188           0 :          info_out = info
     189             :       ELSE
     190       17497 :          IF (info /= 0) &
     191           0 :             CPABORT("Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
     192             :       END IF
     193             : 
     194       17497 :       CALL timestop(handle)
     195             : 
     196       17497 :    END SUBROUTINE cp_fm_cholesky_invert
     197             : 
     198             : ! **************************************************************************************************
     199             : !> \brief reduce a matrix pencil A,B to normal form
     200             : !>      B has to be cholesky decomposed with  cp_fm_cholesky_decompose
     201             : !>      before calling this routine
     202             : !>      A,B -> inv(U^T)*A*inv(U),1
     203             : !>      (AX=BX -> inv(U^T)*A*inv(U)*U*X=U*X hence evecs U*X)
     204             : !> \param matrix the symmetric matrix A
     205             : !> \param matrixb the cholesky decomposition of matrix B
     206             : !> \param itype ...
     207             : !> \par History
     208             : !>      05.2002 created [JVdV]
     209             : !> \author Joost VandeVondele
     210             : ! **************************************************************************************************
     211        4002 :    SUBROUTINE cp_fm_cholesky_reduce(matrix, matrixb, itype)
     212             :       TYPE(cp_fm_type), INTENT(IN)     :: matrix, matrixb
     213             :       INTEGER, OPTIONAL                   :: itype
     214             : 
     215             :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_reduce'
     216        4002 :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: a, b
     217             :       INTEGER                                   :: info, handle
     218             :       INTEGER                                   :: n, my_itype
     219             : #if defined(__parallel)
     220             :       REAL(KIND=dp)                             :: scale
     221             :       INTEGER, DIMENSION(9)                     :: desca, descb
     222             : #endif
     223             : 
     224        4002 :       CALL timeset(routineN, handle)
     225             : 
     226        4002 :       n = matrix%matrix_struct%nrow_global
     227             : 
     228        4002 :       my_itype = 1
     229        4002 :       IF (PRESENT(itype)) my_itype = itype
     230             : 
     231        4002 :       a => matrix%local_data
     232        4002 :       b => matrixb%local_data
     233             : 
     234             : #if defined(__parallel)
     235             : 
     236       40020 :       desca(:) = matrix%matrix_struct%descriptor(:)
     237       40020 :       descb(:) = matrixb%matrix_struct%descriptor(:)
     238             : 
     239        4002 :       CALL pdsygst(my_itype, 'U', n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, scale, info)
     240             : 
     241             :       ! this is supposed to be one in current version of lapack
     242             :       ! if not, eigenvalues have to be scaled by this number
     243        4002 :       IF (scale /= 1.0_dp) &
     244           0 :          CPABORT("scale not equal 1 (scale="//cp_to_string(scale)//")")
     245             : #else
     246             : 
     247             :       CALL dsygst(my_itype, 'U', n, a(1, 1), n, b(1, 1), n, info)
     248             : 
     249             : #endif
     250             : 
     251        4002 :       CPASSERT(info == 0)
     252             : 
     253        4002 :       CALL timestop(handle)
     254             : 
     255        4002 :    END SUBROUTINE cp_fm_cholesky_reduce
     256             : 
     257             : !
     258             : ! op can be "SOLVE" (out = U^-1 * in ) or "MULTIPLY"   (out = U * in )
     259             : ! pos can be "LEFT" or "RIGHT" (U at the left or at the right)
     260             : !
     261             : ! DEPRECATED, see cp_fm_basic_linalg:cp_fm_triangular_multiply
     262             : !
     263             : ! **************************************************************************************************
     264             : !> \brief ...
     265             : !> \param matrix ...
     266             : !> \param neig ...
     267             : !> \param matrixb ...
     268             : !> \param matrixout ...
     269             : !> \param op ...
     270             : !> \param pos ...
     271             : !> \param transa ...
     272             : ! **************************************************************************************************
     273      231649 :    SUBROUTINE cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
     274             :       TYPE(cp_fm_type), INTENT(IN)         :: matrix, matrixb, matrixout
     275             :       INTEGER, INTENT(IN)                     :: neig
     276             :       CHARACTER(LEN=*), INTENT(IN)            :: op
     277             :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL  :: pos
     278             :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL  :: transa
     279             : 
     280             :       CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_restore'
     281      231649 :       REAL(KIND=dp), DIMENSION(:, :), POINTER         :: a, b, out
     282      231649 :       REAL(KIND=sp), DIMENSION(:, :), POINTER         :: a_sp, b_sp, out_sp
     283             :       INTEGER                                   :: itype, handle
     284             :       INTEGER                                   :: n
     285             :       REAL(KIND=dp)                           :: alpha
     286             :       INTEGER                                   :: myprow, mypcol
     287             :       TYPE(cp_blacs_env_type), POINTER          :: context
     288             :       CHARACTER                                 :: chol_pos, chol_transa
     289             : #if defined(__parallel)
     290             :       INTEGER                                   :: i
     291             :       INTEGER, DIMENSION(9)                     :: desca, descb, descout
     292             : #endif
     293             : 
     294      231649 :       CALL timeset(routineN, handle)
     295             : 
     296      231649 :       context => matrix%matrix_struct%context
     297      231649 :       myprow = context%mepos(1)
     298      231649 :       mypcol = context%mepos(2)
     299      231649 :       n = matrix%matrix_struct%nrow_global
     300      231649 :       itype = 1
     301      231649 :       IF (op /= "SOLVE" .AND. op /= "MULTIPLY") &
     302           0 :          CPABORT("wrong argument op")
     303             : 
     304      231649 :       IF (PRESENT(pos)) THEN
     305       75881 :          SELECT CASE (pos)
     306             :          CASE ("LEFT")
     307       75881 :             chol_pos = 'L'
     308             :          CASE ("RIGHT")
     309       74913 :             chol_pos = 'R'
     310             :          CASE DEFAULT
     311      150794 :             CPABORT("wrong argument pos")
     312             :          END SELECT
     313             :       ELSE
     314       80855 :          chol_pos = 'L'
     315             :       END IF
     316             : 
     317      231649 :       chol_transa = 'N'
     318      231649 :       IF (PRESENT(transa)) chol_transa = transa
     319             : 
     320      231649 :       IF ((matrix%use_sp .NEQV. matrixb%use_sp) .OR. (matrix%use_sp .NEQV. matrixout%use_sp)) &
     321           0 :          CPABORT("not the same precision")
     322             : 
     323             :       ! notice b is the cholesky guy
     324      231649 :       a => matrix%local_data
     325      231649 :       b => matrixb%local_data
     326      231649 :       out => matrixout%local_data
     327      231649 :       a_sp => matrix%local_data_sp
     328      231649 :       b_sp => matrixb%local_data_sp
     329      231649 :       out_sp => matrixout%local_data_sp
     330             : 
     331             : #if defined(__parallel)
     332             : 
     333     2316490 :       desca(:) = matrix%matrix_struct%descriptor(:)
     334     2316490 :       descb(:) = matrixb%matrix_struct%descriptor(:)
     335     2316490 :       descout(:) = matrixout%matrix_struct%descriptor(:)
     336      231649 :       alpha = 1.0_dp
     337     4452097 :       DO i = 1, neig
     338     4452097 :          IF (matrix%use_sp) THEN
     339           0 :             CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, out_sp(1, 1), 1, i, descout, 1)
     340             :          ELSE
     341     4220448 :             CALL pdcopy(n, a(1, 1), 1, i, desca, 1, out(1, 1), 1, i, descout, 1)
     342             :          END IF
     343             :       END DO
     344      231649 :       IF (op .EQ. "SOLVE") THEN
     345      230219 :          IF (matrix%use_sp) THEN
     346             :             CALL pstrsm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, &
     347           0 :                         out_sp(1, 1), 1, 1, descout)
     348             :          ELSE
     349      230219 :             CALL pdtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
     350             :          END IF
     351             :       ELSE
     352        1430 :          IF (matrix%use_sp) THEN
     353             :             CALL pstrmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, &
     354           0 :                         out_sp(1, 1), 1, 1, descout)
     355             :          ELSE
     356        1430 :             CALL pdtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
     357             :          END IF
     358             :       END IF
     359             : #else
     360             : 
     361             :       alpha = 1.0_dp
     362             :       IF (matrix%use_sp) THEN
     363             :          CALL scopy(neig*n, a_sp(1, 1), 1, out_sp(1, 1), 1)
     364             :       ELSE
     365             :          CALL dcopy(neig*n, a(1, 1), 1, out(1, 1), 1)
     366             :       END IF
     367             :       IF (op .EQ. "SOLVE") THEN
     368             :          IF (matrix%use_sp) THEN
     369             :             CALL strsm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), SIZE(b_sp, 1), out_sp(1, 1), n)
     370             :          ELSE
     371             :             CALL dtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, 1), out(1, 1), n)
     372             :          END IF
     373             :       ELSE
     374             :          IF (matrix%use_sp) THEN
     375             :             CALL strmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), n, out_sp(1, 1), n)
     376             :          ELSE
     377             :             CALL dtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), n, out(1, 1), n)
     378             :          END IF
     379             :       END IF
     380             : 
     381             : #endif
     382             : 
     383      231649 :       CALL timestop(handle)
     384             : 
     385      231649 :    END SUBROUTINE cp_fm_cholesky_restore
     386             : 
     387             : END MODULE cp_fm_cholesky