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