Line data Source code
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 : ! **************************************************************************************************
14 : MODULE 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 : !***
47 : CONTAINS
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 68874 : 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 68874 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
70 : #if defined(__parallel)
71 : INTEGER, DIMENSION(9) :: desca
72 : #endif
73 :
74 68874 : CALL timeset(routineN, handle)
75 :
76 : my_n = MIN(matrix%matrix_struct%nrow_global, &
77 68874 : matrix%matrix_struct%ncol_global)
78 68874 : IF (PRESENT(n)) THEN
79 17255 : CPASSERT(n <= my_n)
80 17255 : my_n = n
81 : END IF
82 :
83 68874 : a => matrix%local_data
84 :
85 : #if defined(__parallel)
86 688740 : 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 68874 : 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 68874 : IF (PRESENT(info_out)) THEN
107 23549 : info_out = info
108 45325 : ELSE IF (info /= 0) THEN
109 : CALL cp_abort(__LOCATION__, &
110 0 : "Cholesky decompose failed: the matrix is not positive definite or ill-conditioned.")
111 : END IF
112 :
113 68874 : CALL timestop(handle)
114 :
115 68874 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
136 : #if defined(__parallel)
137 : INTEGER, DIMENSION(9) :: desca, descb
138 : #endif
139 :
140 0 : CALL timeset(routineN, handle)
141 :
142 : my_n = MIN(matrix%matrix_struct%nrow_global, &
143 0 : matrix%matrix_struct%ncol_global)
144 0 : IF (PRESENT(n)) THEN
145 0 : CPASSERT(n <= my_n)
146 0 : my_n = n
147 : END IF
148 0 : nrhs = matrixb%matrix_struct%ncol_global
149 :
150 0 : a => matrix%local_data
151 0 : b => matrixb%local_data
152 :
153 : #if defined(__parallel)
154 0 : desca(:) = matrix%matrix_struct%descriptor(:)
155 0 : descb(:) = matrixb%matrix_struct%descriptor(:)
156 0 : 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 0 : IF (PRESENT(info_out)) THEN
162 0 : info_out = info
163 0 : ELSE IF (info /= 0) THEN
164 : CALL cp_abort(__LOCATION__, &
165 : "Cholesky solve failed: "// &
166 0 : "the matrix is not positive definite or ill-conditioned, info="//cp_to_string(info))
167 : END IF
168 :
169 0 : CALL timestop(handle)
170 :
171 0 : 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 20575 : 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 20575 : 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 20575 : CALL timeset(routineN, handle)
198 :
199 : my_n = MIN(matrix%matrix_struct%nrow_global, &
200 20575 : matrix%matrix_struct%ncol_global)
201 20575 : IF (PRESENT(n)) THEN
202 12212 : CPASSERT(n <= my_n)
203 12212 : my_n = n
204 : END IF
205 :
206 20575 : a => matrix%local_data
207 :
208 : #if defined(__parallel)
209 :
210 205750 : 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 20575 : 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 20575 : IF (PRESENT(info_out)) THEN
235 0 : info_out = info
236 : ELSE
237 20575 : IF (info /= 0) &
238 0 : CPABORT("Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
239 : END IF
240 :
241 20575 : IF (PRESENT(work) .AND. info == 0) CALL cp_fm_uplo_to_full(matrix, work)
242 :
243 20575 : CALL timestop(handle)
244 :
245 20575 : 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 4802 : 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 4802 : 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 4802 : CALL timeset(routineN, handle)
274 :
275 4802 : n = matrix%matrix_struct%nrow_global
276 :
277 4802 : my_itype = 1
278 4802 : IF (PRESENT(itype)) my_itype = itype
279 :
280 4802 : a => matrix%local_data
281 4802 : b => matrixb%local_data
282 :
283 : #if defined(__parallel)
284 :
285 48020 : desca(:) = matrix%matrix_struct%descriptor(:)
286 48020 : descb(:) = matrixb%matrix_struct%descriptor(:)
287 :
288 4802 : 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 4802 : IF (scale /= 1.0_dp) &
293 0 : 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 4802 : CPASSERT(info == 0)
301 :
302 4802 : CALL timestop(handle)
303 :
304 4802 : 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 372299 : 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 372299 : 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 372299 : CALL timeset(routineN, handle)
338 :
339 : ! wrong argument op
340 372299 : CPASSERT((op == "SOLVE") .OR. (op == "MULTIPLY"))
341 :
342 372299 : chol_pos = 'L'
343 372299 : IF (PRESENT(pos)) chol_pos = pos(1:1)
344 372299 : CPASSERT((chol_pos == 'L') .OR. (chol_pos == 'R'))
345 :
346 372299 : chol_transa = 'N'
347 372299 : IF (PRESENT(transa)) chol_transa = transa
348 :
349 : ! notice b is the cholesky guy
350 372299 : a => fm_matrix%local_data
351 372299 : b => fm_matrixb%local_data
352 372299 : outm => fm_matrixout%local_data
353 372299 : n = fm_matrix%matrix_struct%nrow_global
354 : mdim = MERGE(1, 2, 'L' == chol_pos)
355 372299 : alpha = 1.0_dp
356 :
357 : #if defined(__parallel)
358 3722990 : desca(:) = fm_matrix%matrix_struct%descriptor(:)
359 3722990 : descb(:) = fm_matrixb%matrix_struct%descriptor(:)
360 3722990 : descout(:) = fm_matrixout%matrix_struct%descriptor(:)
361 7498485 : DO i = 1, neig
362 7498485 : CALL pdcopy(n, a(1, 1), 1, i, desca, 1, outm(1, 1), 1, i, descout, 1)
363 : END DO
364 372299 : IF (op == "SOLVE") THEN
365 370681 : 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 1618 : 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 372299 : CALL timestop(handle)
386 :
387 372299 : END SUBROUTINE cp_fm_cholesky_restore
388 :
389 : END MODULE cp_fm_cholesky
|
