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_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 : !***
46 : CONTAINS
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 63636 : 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 63636 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
69 : #if defined(__parallel)
70 : INTEGER, DIMENSION(9) :: desca
71 : #endif
72 :
73 63636 : CALL timeset(routineN, handle)
74 :
75 : my_n = MIN(matrix%matrix_struct%nrow_global, &
76 63636 : matrix%matrix_struct%ncol_global)
77 63636 : IF (PRESENT(n)) THEN
78 15899 : CPASSERT(n <= my_n)
79 15899 : my_n = n
80 : END IF
81 :
82 63636 : a => matrix%local_data
83 :
84 : #if defined(__parallel)
85 636360 : 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 63636 : 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 63636 : IF (PRESENT(info_out)) THEN
106 21343 : info_out = info
107 42293 : ELSE IF (info /= 0) THEN
108 : CALL cp_abort(__LOCATION__, &
109 0 : "Cholesky decompose failed: the matrix is not positive definite or ill-conditioned.")
110 : END IF
111 :
112 63636 : CALL timestop(handle)
113 :
114 63636 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
135 : #if defined(__parallel)
136 : INTEGER, DIMENSION(9) :: desca, descb
137 : #endif
138 :
139 0 : CALL timeset(routineN, handle)
140 :
141 : my_n = MIN(matrix%matrix_struct%nrow_global, &
142 0 : matrix%matrix_struct%ncol_global)
143 0 : IF (PRESENT(n)) THEN
144 0 : CPASSERT(n <= my_n)
145 0 : my_n = n
146 : END IF
147 0 : nrhs = matrixb%matrix_struct%ncol_global
148 :
149 0 : a => matrix%local_data
150 0 : b => matrixb%local_data
151 :
152 : #if defined(__parallel)
153 0 : desca(:) = matrix%matrix_struct%descriptor(:)
154 0 : descb(:) = matrixb%matrix_struct%descriptor(:)
155 0 : 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 0 : IF (PRESENT(info_out)) THEN
161 0 : info_out = info
162 0 : ELSE IF (info /= 0) THEN
163 : CALL cp_abort(__LOCATION__, &
164 : "Cholesky solve failed: "// &
165 0 : "the matrix is not positive definite or ill-conditioned, info="//cp_to_string(info))
166 : END IF
167 :
168 0 : CALL timestop(handle)
169 :
170 0 : 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 18477 : 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 18477 : 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 18477 : CALL timeset(routineN, handle)
194 :
195 : my_n = MIN(matrix%matrix_struct%nrow_global, &
196 18477 : matrix%matrix_struct%ncol_global)
197 18477 : IF (PRESENT(n)) THEN
198 0 : CPASSERT(n <= my_n)
199 0 : my_n = n
200 : END IF
201 :
202 18477 : a => matrix%local_data
203 :
204 : #if defined(__parallel)
205 :
206 184770 : 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 18477 : 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 18477 : IF (PRESENT(info_out)) THEN
231 0 : info_out = info
232 : ELSE
233 18477 : IF (info /= 0) &
234 0 : CPABORT("Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
235 : END IF
236 :
237 18477 : CALL timestop(handle)
238 :
239 18477 : 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 4592 : 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 4592 : 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 4592 : CALL timeset(routineN, handle)
268 :
269 4592 : n = matrix%matrix_struct%nrow_global
270 :
271 4592 : my_itype = 1
272 4592 : IF (PRESENT(itype)) my_itype = itype
273 :
274 4592 : a => matrix%local_data
275 4592 : b => matrixb%local_data
276 :
277 : #if defined(__parallel)
278 :
279 45920 : desca(:) = matrix%matrix_struct%descriptor(:)
280 45920 : descb(:) = matrixb%matrix_struct%descriptor(:)
281 :
282 4592 : 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 4592 : IF (scale /= 1.0_dp) &
287 0 : 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 4592 : CPASSERT(info == 0)
295 :
296 4592 : CALL timestop(handle)
297 :
298 4592 : 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 340855 : 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 340855 : 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 340855 : CALL timeset(routineN, handle)
332 :
333 : ! wrong argument op
334 340855 : CPASSERT((op == "SOLVE") .OR. (op == "MULTIPLY"))
335 :
336 340855 : chol_pos = 'L'
337 340855 : IF (PRESENT(pos)) chol_pos = pos(1:1)
338 340855 : CPASSERT((chol_pos == 'L') .OR. (chol_pos == 'R'))
339 :
340 340855 : chol_transa = 'N'
341 340855 : IF (PRESENT(transa)) chol_transa = transa
342 :
343 : ! notice b is the cholesky guy
344 340855 : a => fm_matrix%local_data
345 340855 : b => fm_matrixb%local_data
346 340855 : outm => fm_matrixout%local_data
347 340855 : n = fm_matrix%matrix_struct%nrow_global
348 : mdim = MERGE(1, 2, 'L' == chol_pos)
349 340855 : alpha = 1.0_dp
350 :
351 : #if defined(__parallel)
352 3408550 : desca(:) = fm_matrix%matrix_struct%descriptor(:)
353 3408550 : descb(:) = fm_matrixb%matrix_struct%descriptor(:)
354 3408550 : descout(:) = fm_matrixout%matrix_struct%descriptor(:)
355 7275737 : DO i = 1, neig
356 7275737 : CALL pdcopy(n, a(1, 1), 1, i, desca, 1, outm(1, 1), 1, i, descout, 1)
357 : END DO
358 340855 : IF (op == "SOLVE") THEN
359 339329 : 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 1526 : 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 340855 : CALL timestop(handle)
380 :
381 340855 : END SUBROUTINE cp_fm_cholesky_restore
382 :
383 : END MODULE cp_fm_cholesky
|
