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