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