(git:44e3845)
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' == 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 == "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 == "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 == "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