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_pspotrf_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 : sp
24 : #include "../base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 : PRIVATE
28 :
29 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_fm_cholesky'
31 :
32 : PUBLIC :: cp_fm_cholesky_decompose, cp_fm_cholesky_invert, &
33 : cp_fm_cholesky_reduce, cp_fm_cholesky_restore
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 55929 : 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 55929 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
70 55929 : REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp
71 : #if defined(__parallel)
72 : INTEGER, DIMENSION(9) :: desca
73 : #endif
74 :
75 55929 : CALL timeset(routineN, handle)
76 :
77 : my_n = MIN(matrix%matrix_struct%nrow_global, &
78 55929 : matrix%matrix_struct%ncol_global)
79 55929 : IF (PRESENT(n)) THEN
80 15438 : CPASSERT(n <= my_n)
81 15438 : my_n = n
82 : END IF
83 :
84 55929 : a => matrix%local_data
85 55929 : a_sp => matrix%local_data_sp
86 :
87 : #if defined(__parallel)
88 559290 : desca(:) = matrix%matrix_struct%descriptor(:)
89 : #if defined(__DLAF)
90 : IF (cholesky_type == FM_CHOLESKY_TYPE_DLAF .AND. matrix%matrix_struct%nrow_global >= dlaf_cholesky_n_min) THEN
91 : ! Initialize DLA-Future on-demand; if already initialized, does nothing
92 : CALL cp_dlaf_initialize()
93 :
94 : ! Create DLAF grid from BLACS context; if already present, does nothing
95 : CALL cp_dlaf_create_grid(matrix%matrix_struct%context%get_handle())
96 :
97 : IF (matrix%use_sp) THEN
98 : CALL cp_pspotrf_dlaf('U', my_n, a_sp(:, :), 1, 1, desca, info)
99 : ELSE
100 : CALL cp_pdpotrf_dlaf('U', my_n, a(:, :), 1, 1, desca, info)
101 : END IF
102 : ELSE
103 : #endif
104 55929 : IF (matrix%use_sp) THEN
105 0 : CALL pspotrf('U', my_n, a_sp(1, 1), 1, 1, desca, info)
106 : ELSE
107 55929 : CALL pdpotrf('U', my_n, a(1, 1), 1, 1, desca, info)
108 : END IF
109 : #if defined(__DLAF)
110 : END IF
111 : #endif
112 : #else
113 : IF (matrix%use_sp) THEN
114 : CALL spotrf('U', my_n, a_sp(1, 1), SIZE(a_sp, 1), info)
115 : ELSE
116 : CALL dpotrf('U', my_n, a(1, 1), SIZE(a, 1), info)
117 : END IF
118 : #endif
119 :
120 55929 : IF (PRESENT(info_out)) THEN
121 20442 : info_out = info
122 35487 : ELSE IF (info /= 0) THEN
123 : CALL cp_abort(__LOCATION__, &
124 0 : "Cholesky decompose failed: the matrix is not positive definite or ill-conditioned.")
125 : END IF
126 :
127 55929 : CALL timestop(handle)
128 :
129 55929 : END SUBROUTINE cp_fm_cholesky_decompose
130 :
131 : ! **************************************************************************************************
132 : !> \brief used to replace the cholesky decomposition by the inverse
133 : !> \param matrix the matrix to invert (must be an upper triangular matrix)
134 : !> \param n size of the matrix to invert (defaults to the min(size(matrix)))
135 : !> \param info_out ...
136 : !> \par History
137 : !> 05.2002 created [JVdV]
138 : !> \author Joost VandeVondele
139 : ! **************************************************************************************************
140 17497 : SUBROUTINE cp_fm_cholesky_invert(matrix, n, info_out)
141 : TYPE(cp_fm_type), INTENT(IN) :: matrix
142 : INTEGER, INTENT(in), OPTIONAL :: n
143 : INTEGER, INTENT(OUT), OPTIONAL :: info_out
144 :
145 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_invert'
146 17497 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
147 17497 : REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp
148 : INTEGER :: info, handle
149 : INTEGER :: my_n
150 : #if defined(__parallel)
151 : INTEGER, DIMENSION(9) :: desca
152 : #endif
153 :
154 17497 : CALL timeset(routineN, handle)
155 :
156 : my_n = MIN(matrix%matrix_struct%nrow_global, &
157 17497 : matrix%matrix_struct%ncol_global)
158 17497 : IF (PRESENT(n)) THEN
159 0 : CPASSERT(n <= my_n)
160 0 : my_n = n
161 : END IF
162 :
163 17497 : a => matrix%local_data
164 17497 : a_sp => matrix%local_data_sp
165 :
166 : #if defined(__parallel)
167 :
168 174970 : desca(:) = matrix%matrix_struct%descriptor(:)
169 :
170 17497 : IF (matrix%use_sp) THEN
171 0 : CALL pspotri('U', my_n, a_sp(1, 1), 1, 1, desca, info)
172 : ELSE
173 17497 : CALL pdpotri('U', my_n, a(1, 1), 1, 1, desca, info)
174 : END IF
175 :
176 : #else
177 :
178 : IF (matrix%use_sp) THEN
179 : CALL spotri('U', my_n, a_sp(1, 1), SIZE(a_sp, 1), info)
180 : ELSE
181 : CALL dpotri('U', my_n, a(1, 1), SIZE(a, 1), info)
182 : END IF
183 :
184 : #endif
185 :
186 17497 : IF (PRESENT(info_out)) THEN
187 0 : info_out = info
188 : ELSE
189 17497 : IF (info /= 0) &
190 0 : CPABORT("Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
191 : END IF
192 :
193 17497 : CALL timestop(handle)
194 :
195 17497 : END SUBROUTINE cp_fm_cholesky_invert
196 :
197 : ! **************************************************************************************************
198 : !> \brief reduce a matrix pencil A,B to normal form
199 : !> B has to be cholesky decomposed with cp_fm_cholesky_decompose
200 : !> before calling this routine
201 : !> A,B -> inv(U^T)*A*inv(U),1
202 : !> (AX=BX -> inv(U^T)*A*inv(U)*U*X=U*X hence evecs U*X)
203 : !> \param matrix the symmetric matrix A
204 : !> \param matrixb the cholesky decomposition of matrix B
205 : !> \param itype ...
206 : !> \par History
207 : !> 05.2002 created [JVdV]
208 : !> \author Joost VandeVondele
209 : ! **************************************************************************************************
210 4002 : SUBROUTINE cp_fm_cholesky_reduce(matrix, matrixb, itype)
211 : TYPE(cp_fm_type), INTENT(IN) :: matrix, matrixb
212 : INTEGER, OPTIONAL :: itype
213 :
214 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_reduce'
215 4002 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
216 : INTEGER :: info, handle
217 : INTEGER :: n, my_itype
218 : #if defined(__parallel)
219 : REAL(KIND=dp) :: scale
220 : INTEGER, DIMENSION(9) :: desca, descb
221 : #endif
222 :
223 4002 : CALL timeset(routineN, handle)
224 :
225 4002 : n = matrix%matrix_struct%nrow_global
226 :
227 4002 : my_itype = 1
228 4002 : IF (PRESENT(itype)) my_itype = itype
229 :
230 4002 : a => matrix%local_data
231 4002 : b => matrixb%local_data
232 :
233 : #if defined(__parallel)
234 :
235 40020 : desca(:) = matrix%matrix_struct%descriptor(:)
236 40020 : descb(:) = matrixb%matrix_struct%descriptor(:)
237 :
238 4002 : CALL pdsygst(my_itype, 'U', n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, scale, info)
239 :
240 : ! this is supposed to be one in current version of lapack
241 : ! if not, eigenvalues have to be scaled by this number
242 4002 : IF (scale /= 1.0_dp) &
243 0 : CPABORT("scale not equal 1 (scale="//cp_to_string(scale)//")")
244 : #else
245 :
246 : CALL dsygst(my_itype, 'U', n, a(1, 1), n, b(1, 1), n, info)
247 :
248 : #endif
249 :
250 4002 : CPASSERT(info == 0)
251 :
252 4002 : CALL timestop(handle)
253 :
254 4002 : END SUBROUTINE cp_fm_cholesky_reduce
255 :
256 : !
257 : ! op can be "SOLVE" (out = U^-1 * in ) or "MULTIPLY" (out = U * in )
258 : ! pos can be "LEFT" or "RIGHT" (U at the left or at the right)
259 : !
260 : ! DEPRECATED, see cp_fm_basic_linalg:cp_fm_triangular_multiply
261 : !
262 : ! **************************************************************************************************
263 : !> \brief ...
264 : !> \param matrix ...
265 : !> \param neig ...
266 : !> \param matrixb ...
267 : !> \param matrixout ...
268 : !> \param op ...
269 : !> \param pos ...
270 : !> \param transa ...
271 : ! **************************************************************************************************
272 231979 : SUBROUTINE cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
273 : TYPE(cp_fm_type), INTENT(IN) :: matrix, matrixb, matrixout
274 : INTEGER, INTENT(IN) :: neig
275 : CHARACTER(LEN=*), INTENT(IN) :: op
276 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pos
277 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: transa
278 :
279 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_restore'
280 231979 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, out
281 231979 : REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp, out_sp
282 : INTEGER :: itype, handle
283 : INTEGER :: n
284 : REAL(KIND=dp) :: alpha
285 : INTEGER :: myprow, mypcol
286 : TYPE(cp_blacs_env_type), POINTER :: context
287 : CHARACTER :: chol_pos, chol_transa
288 : #if defined(__parallel)
289 : INTEGER :: i
290 : INTEGER, DIMENSION(9) :: desca, descb, descout
291 : #endif
292 :
293 231979 : CALL timeset(routineN, handle)
294 :
295 231979 : context => matrix%matrix_struct%context
296 231979 : myprow = context%mepos(1)
297 231979 : mypcol = context%mepos(2)
298 231979 : n = matrix%matrix_struct%nrow_global
299 231979 : itype = 1
300 231979 : IF (op /= "SOLVE" .AND. op /= "MULTIPLY") &
301 0 : CPABORT("wrong argument op")
302 :
303 231979 : IF (PRESENT(pos)) THEN
304 75991 : SELECT CASE (pos)
305 : CASE ("LEFT")
306 75991 : chol_pos = 'L'
307 : CASE ("RIGHT")
308 75023 : chol_pos = 'R'
309 : CASE DEFAULT
310 151014 : CPABORT("wrong argument pos")
311 : END SELECT
312 : ELSE
313 80965 : chol_pos = 'L'
314 : END IF
315 :
316 231979 : chol_transa = 'N'
317 231979 : IF (PRESENT(transa)) chol_transa = transa
318 :
319 231979 : IF ((matrix%use_sp .NEQV. matrixb%use_sp) .OR. (matrix%use_sp .NEQV. matrixout%use_sp)) &
320 0 : CPABORT("not the same precision")
321 :
322 : ! notice b is the cholesky guy
323 231979 : a => matrix%local_data
324 231979 : b => matrixb%local_data
325 231979 : out => matrixout%local_data
326 231979 : a_sp => matrix%local_data_sp
327 231979 : b_sp => matrixb%local_data_sp
328 231979 : out_sp => matrixout%local_data_sp
329 :
330 : #if defined(__parallel)
331 :
332 2319790 : desca(:) = matrix%matrix_struct%descriptor(:)
333 2319790 : descb(:) = matrixb%matrix_struct%descriptor(:)
334 2319790 : descout(:) = matrixout%matrix_struct%descriptor(:)
335 231979 : alpha = 1.0_dp
336 4460667 : DO i = 1, neig
337 4460667 : IF (matrix%use_sp) THEN
338 0 : CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, out_sp(1, 1), 1, i, descout, 1)
339 : ELSE
340 4228688 : CALL pdcopy(n, a(1, 1), 1, i, desca, 1, out(1, 1), 1, i, descout, 1)
341 : END IF
342 : END DO
343 231979 : IF (op .EQ. "SOLVE") THEN
344 230549 : IF (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 : out_sp(1, 1), 1, 1, descout)
347 : ELSE
348 230549 : CALL pdtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
349 : END IF
350 : ELSE
351 1430 : IF (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 : out_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, out(1, 1), 1, 1, descout)
356 : END IF
357 : END IF
358 : #else
359 :
360 : alpha = 1.0_dp
361 : IF (matrix%use_sp) THEN
362 : CALL scopy(neig*n, a_sp(1, 1), 1, out_sp(1, 1), 1)
363 : ELSE
364 : CALL dcopy(neig*n, a(1, 1), 1, out(1, 1), 1)
365 : END IF
366 : IF (op .EQ. "SOLVE") THEN
367 : IF (matrix%use_sp) THEN
368 : 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)
369 : ELSE
370 : CALL dtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, 1), out(1, 1), n)
371 : END IF
372 : ELSE
373 : IF (matrix%use_sp) THEN
374 : CALL strmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), n, out_sp(1, 1), n)
375 : ELSE
376 : CALL dtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), n, out(1, 1), n)
377 : END IF
378 : END IF
379 :
380 : #endif
381 :
382 231979 : CALL timestop(handle)
383 :
384 231979 : END SUBROUTINE cp_fm_cholesky_restore
385 :
386 : END MODULE cp_fm_cholesky
|
