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