Line data Source code
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 : ! **************************************************************************************************
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 : 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 : !***
48 : CONTAINS
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 56727 : 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 56727 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
71 56727 : REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp
72 : #if defined(__parallel)
73 : INTEGER, DIMENSION(9) :: desca
74 : #endif
75 :
76 56727 : CALL timeset(routineN, handle)
77 :
78 : my_n = MIN(matrix%matrix_struct%nrow_global, &
79 56727 : matrix%matrix_struct%ncol_global)
80 56727 : IF (PRESENT(n)) THEN
81 15496 : CPASSERT(n <= my_n)
82 15496 : my_n = n
83 : END IF
84 :
85 56727 : a => matrix%local_data
86 56727 : a_sp => matrix%local_data_sp
87 :
88 : #if defined(__parallel)
89 567270 : 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 56727 : IF (matrix%use_sp) THEN
106 0 : CALL pspotrf('U', my_n, a_sp(1, 1), 1, 1, desca, info)
107 : ELSE
108 56727 : 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 56727 : IF (PRESENT(info_out)) THEN
122 20602 : info_out = info
123 36125 : ELSE IF (info /= 0) THEN
124 : CALL cp_abort(__LOCATION__, &
125 0 : "Cholesky decompose failed: the matrix is not positive definite or ill-conditioned.")
126 : END IF
127 :
128 56727 : CALL timestop(handle)
129 :
130 56727 : 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 17635 : 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 17635 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
148 17635 : 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 17635 : CALL timeset(routineN, handle)
155 :
156 : my_n = MIN(matrix%matrix_struct%nrow_global, &
157 17635 : matrix%matrix_struct%ncol_global)
158 17635 : IF (PRESENT(n)) THEN
159 0 : CPASSERT(n <= my_n)
160 0 : my_n = n
161 : END IF
162 :
163 17635 : a => matrix%local_data
164 17635 : a_sp => matrix%local_data_sp
165 :
166 : #if defined(__parallel)
167 :
168 176350 : 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 17635 : IF (matrix%use_sp) THEN
186 0 : CALL pspotri('U', my_n, a_sp(1, 1), 1, 1, desca, info)
187 : ELSE
188 17635 : 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 17635 : IF (PRESENT(info_out)) THEN
205 0 : info_out = info
206 : ELSE
207 17635 : IF (info /= 0) &
208 0 : CPABORT("Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
209 : END IF
210 :
211 17635 : CALL timestop(handle)
212 :
213 17635 : 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 4152 : 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 4152 : 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 4152 : CALL timeset(routineN, handle)
242 :
243 4152 : n = matrix%matrix_struct%nrow_global
244 :
245 4152 : my_itype = 1
246 4152 : IF (PRESENT(itype)) my_itype = itype
247 :
248 4152 : a => matrix%local_data
249 4152 : b => matrixb%local_data
250 :
251 : #if defined(__parallel)
252 :
253 41520 : desca(:) = matrix%matrix_struct%descriptor(:)
254 41520 : descb(:) = matrixb%matrix_struct%descriptor(:)
255 :
256 4152 : 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 4152 : IF (scale /= 1.0_dp) &
261 0 : 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 4152 : CPASSERT(info == 0)
269 :
270 4152 : CALL timestop(handle)
271 :
272 4152 : 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 235137 : 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 235137 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, outm
295 235137 : 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 235137 : CALL timeset(routineN, handle)
307 :
308 : ! wrong argument op
309 235137 : CPASSERT((op == "SOLVE") .OR. (op == "MULTIPLY"))
310 : ! not the same precision
311 235137 : CPASSERT(fm_matrix%use_sp .EQV. fm_matrixb%use_sp)
312 235137 : CPASSERT(fm_matrix%use_sp .EQV. fm_matrixout%use_sp)
313 :
314 235137 : chol_pos = 'L'
315 235137 : IF (PRESENT(pos)) chol_pos = pos(1:1)
316 235137 : CPASSERT((chol_pos == 'L') .OR. (chol_pos == 'R'))
317 :
318 235137 : chol_transa = 'N'
319 235137 : IF (PRESENT(transa)) chol_transa = transa
320 :
321 : ! notice b is the cholesky guy
322 235137 : a => fm_matrix%local_data
323 235137 : b => fm_matrixb%local_data
324 235137 : outm => fm_matrixout%local_data
325 235137 : a_sp => fm_matrix%local_data_sp
326 235137 : b_sp => fm_matrixb%local_data_sp
327 235137 : outm_sp => fm_matrixout%local_data_sp
328 235137 : n = fm_matrix%matrix_struct%nrow_global
329 : mdim = MERGE(1, 2, 'L' == chol_pos)
330 235137 : alpha = 1.0_dp
331 :
332 : #if defined(__parallel)
333 2351370 : desca(:) = fm_matrix%matrix_struct%descriptor(:)
334 2351370 : descb(:) = fm_matrixb%matrix_struct%descriptor(:)
335 2351370 : descout(:) = fm_matrixout%matrix_struct%descriptor(:)
336 4497331 : DO i = 1, neig
337 4497331 : IF (fm_matrix%use_sp) THEN
338 0 : CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, outm_sp(1, 1), 1, i, descout, 1)
339 : ELSE
340 4262194 : CALL pdcopy(n, a(1, 1), 1, i, desca, 1, outm(1, 1), 1, i, descout, 1)
341 : END IF
342 : END DO
343 235137 : IF (op == "SOLVE") THEN
344 233707 : 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 0 : outm_sp(1, 1), 1, 1, descout)
347 : ELSE
348 233707 : 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 1430 : 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 0 : outm_sp(1, 1), 1, 1, descout)
354 : ELSE
355 1430 : 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 235137 : CALL timestop(handle)
395 :
396 235137 : END SUBROUTINE cp_fm_cholesky_restore
397 :
398 : END MODULE cp_fm_cholesky
|
