Line data    Source code

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