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