Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_fm_types, ONLY: cp_fm_type
17 : USE cp_log_handling, ONLY: cp_to_string
18 : USE kinds, ONLY: dp, &
19 : sp
20 : #ifdef __DLAF
21 : USE cp_fm_dlaf_api, ONLY: cp_pdpotrf_dlaf, cp_pspotrf_dlaf
22 : #endif
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 :
34 : !***
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief used to replace a symmetric positive def. matrix M with its cholesky
39 : !> decomposition U: M = U^T * U, with U upper triangular
40 : !> \param matrix the matrix to replace with its cholesky decomposition
41 : !> \param n the number of row (and columns) of the matrix &
42 : !> (defaults to the min(size(matrix)))
43 : !> \param info_out ...
44 : !> \par History
45 : !> 05.2002 created [JVdV]
46 : !> 12.2002 updated, added n optional parm [fawzi]
47 : !> \author Joost
48 : ! **************************************************************************************************
49 49677 : SUBROUTINE cp_fm_cholesky_decompose(matrix, n, info_out)
50 : TYPE(cp_fm_type), INTENT(IN) :: matrix
51 : INTEGER, INTENT(in), OPTIONAL :: n
52 : INTEGER, INTENT(out), OPTIONAL :: info_out
53 :
54 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_decompose'
55 :
56 : INTEGER :: handle, info, my_n
57 49677 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
58 49677 : REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp
59 : #if defined(__parallel)
60 : INTEGER, DIMENSION(9) :: desca
61 : #endif
62 :
63 49677 : CALL timeset(routineN, handle)
64 :
65 : my_n = MIN(matrix%matrix_struct%nrow_global, &
66 49677 : matrix%matrix_struct%ncol_global)
67 49677 : IF (PRESENT(n)) THEN
68 12436 : CPASSERT(n <= my_n)
69 12436 : my_n = n
70 : END IF
71 :
72 49677 : a => matrix%local_data
73 49677 : a_sp => matrix%local_data_sp
74 :
75 : #if defined(__parallel)
76 496770 : desca(:) = matrix%matrix_struct%descriptor(:)
77 : #if defined(__DLAF)
78 : IF (matrix%use_sp) THEN
79 : CALL cp_pspotrf_dlaf('U', my_n, a_sp(:, :), 1, 1, desca, info)
80 : ELSE
81 : CALL cp_pdpotrf_dlaf('U', my_n, a(:, :), 1, 1, desca, info)
82 : END IF
83 : #else
84 49677 : IF (matrix%use_sp) THEN
85 0 : CALL pspotrf('U', my_n, a_sp(1, 1), 1, 1, desca, info)
86 : ELSE
87 49677 : CALL pdpotrf('U', my_n, a(1, 1), 1, 1, desca, info)
88 : END IF
89 : #endif
90 : #else
91 :
92 : IF (matrix%use_sp) THEN
93 : CALL spotrf('U', my_n, a_sp(1, 1), SIZE(a_sp, 1), info)
94 : ELSE
95 : CALL dpotrf('U', my_n, a(1, 1), SIZE(a, 1), info)
96 : END IF
97 :
98 : #endif
99 :
100 49677 : IF (PRESENT(info_out)) THEN
101 17315 : info_out = info
102 : ELSE
103 32362 : IF (info /= 0) &
104 : CALL cp_abort(__LOCATION__, &
105 0 : "Cholesky decompose failed: the matrix is not positive definite or ill-conditioned.")
106 : END IF
107 :
108 49677 : CALL timestop(handle)
109 :
110 49677 : END SUBROUTINE cp_fm_cholesky_decompose
111 :
112 : ! **************************************************************************************************
113 : !> \brief used to replace the cholesky decomposition by the inverse
114 : !> \param matrix the matrix to invert (must be an upper triangular matrix)
115 : !> \param n size of the matrix to invert (defaults to the min(size(matrix)))
116 : !> \param info_out ...
117 : !> \par History
118 : !> 05.2002 created [JVdV]
119 : !> \author Joost VandeVondele
120 : ! **************************************************************************************************
121 14206 : SUBROUTINE cp_fm_cholesky_invert(matrix, n, info_out)
122 : TYPE(cp_fm_type), INTENT(IN) :: matrix
123 : INTEGER, INTENT(in), OPTIONAL :: n
124 : INTEGER, INTENT(OUT), OPTIONAL :: info_out
125 :
126 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_invert'
127 14206 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
128 14206 : REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp
129 : INTEGER :: info, handle
130 : INTEGER :: my_n
131 : #if defined(__parallel)
132 : INTEGER, DIMENSION(9) :: desca
133 : #endif
134 :
135 14206 : CALL timeset(routineN, handle)
136 :
137 : my_n = MIN(matrix%matrix_struct%nrow_global, &
138 14206 : matrix%matrix_struct%ncol_global)
139 14206 : IF (PRESENT(n)) THEN
140 0 : CPASSERT(n <= my_n)
141 0 : my_n = n
142 : END IF
143 :
144 14206 : a => matrix%local_data
145 14206 : a_sp => matrix%local_data_sp
146 :
147 : #if defined(__parallel)
148 :
149 142060 : desca(:) = matrix%matrix_struct%descriptor(:)
150 :
151 14206 : IF (matrix%use_sp) THEN
152 0 : CALL pspotri('U', my_n, a_sp(1, 1), 1, 1, desca, info)
153 : ELSE
154 14206 : CALL pdpotri('U', my_n, a(1, 1), 1, 1, desca, info)
155 : END IF
156 :
157 : #else
158 :
159 : IF (matrix%use_sp) THEN
160 : CALL spotri('U', my_n, a_sp(1, 1), SIZE(a_sp, 1), info)
161 : ELSE
162 : CALL dpotri('U', my_n, a(1, 1), SIZE(a, 1), info)
163 : END IF
164 :
165 : #endif
166 :
167 14206 : IF (PRESENT(info_out)) THEN
168 0 : info_out = info
169 : ELSE
170 14206 : IF (info /= 0) &
171 0 : CPABORT("Cholesky invert failed: the matrix is not positive definite or ill-conditioned.")
172 : END IF
173 :
174 14206 : CALL timestop(handle)
175 :
176 14206 : END SUBROUTINE cp_fm_cholesky_invert
177 :
178 : ! **************************************************************************************************
179 : !> \brief reduce a matrix pencil A,B to normal form
180 : !> B has to be cholesky decomposed with cp_fm_cholesky_decompose
181 : !> before calling this routine
182 : !> A,B -> inv(U^T)*A*inv(U),1
183 : !> (AX=BX -> inv(U^T)*A*inv(U)*U*X=U*X hence evecs U*X)
184 : !> \param matrix the symmetric matrix A
185 : !> \param matrixb the cholesky decomposition of matrix B
186 : !> \param itype ...
187 : !> \par History
188 : !> 05.2002 created [JVdV]
189 : !> \author Joost VandeVondele
190 : ! **************************************************************************************************
191 3650 : SUBROUTINE cp_fm_cholesky_reduce(matrix, matrixb, itype)
192 : TYPE(cp_fm_type), INTENT(IN) :: matrix, matrixb
193 : INTEGER, OPTIONAL :: itype
194 :
195 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_reduce'
196 3650 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b
197 : INTEGER :: info, handle
198 : INTEGER :: n, my_itype
199 : #if defined(__parallel)
200 : REAL(KIND=dp) :: scale
201 : INTEGER, DIMENSION(9) :: desca, descb
202 : #endif
203 :
204 3650 : CALL timeset(routineN, handle)
205 :
206 3650 : n = matrix%matrix_struct%nrow_global
207 :
208 3650 : my_itype = 1
209 3650 : IF (PRESENT(itype)) my_itype = itype
210 :
211 3650 : a => matrix%local_data
212 3650 : b => matrixb%local_data
213 :
214 : #if defined(__parallel)
215 :
216 36500 : desca(:) = matrix%matrix_struct%descriptor(:)
217 36500 : descb(:) = matrixb%matrix_struct%descriptor(:)
218 :
219 3650 : CALL pdsygst(my_itype, 'U', n, a(1, 1), 1, 1, desca, b(1, 1), 1, 1, descb, scale, info)
220 :
221 : ! this is supposed to be one in current version of lapack
222 : ! if not, eigenvalues have to be scaled by this number
223 3650 : IF (scale /= 1.0_dp) &
224 0 : CPABORT("scale not equal 1 (scale="//cp_to_string(scale)//")")
225 : #else
226 :
227 : CALL dsygst(my_itype, 'U', n, a(1, 1), n, b(1, 1), n, info)
228 :
229 : #endif
230 :
231 3650 : CPASSERT(info == 0)
232 :
233 3650 : CALL timestop(handle)
234 :
235 3650 : END SUBROUTINE cp_fm_cholesky_reduce
236 :
237 : !
238 : ! op can be "SOLVE" (out = U^-1 * in ) or "MULTIPLY" (out = U * in )
239 : ! pos can be "LEFT" or "RIGHT" (U at the left or at the right)
240 : !
241 : ! DEPRECATED, see cp_fm_basic_linalg:cp_fm_triangular_multiply
242 : !
243 : ! **************************************************************************************************
244 : !> \brief ...
245 : !> \param matrix ...
246 : !> \param neig ...
247 : !> \param matrixb ...
248 : !> \param matrixout ...
249 : !> \param op ...
250 : !> \param pos ...
251 : !> \param transa ...
252 : ! **************************************************************************************************
253 219261 : SUBROUTINE cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
254 : TYPE(cp_fm_type), INTENT(IN) :: matrix, matrixb, matrixout
255 : INTEGER, INTENT(IN) :: neig
256 : CHARACTER(LEN=*), INTENT(IN) :: op
257 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pos
258 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: transa
259 :
260 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_cholesky_restore'
261 219261 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, b, out
262 219261 : REAL(KIND=sp), DIMENSION(:, :), POINTER :: a_sp, b_sp, out_sp
263 : INTEGER :: itype, handle
264 : INTEGER :: n
265 : REAL(KIND=dp) :: alpha
266 : INTEGER :: myprow, mypcol
267 : TYPE(cp_blacs_env_type), POINTER :: context
268 : CHARACTER :: chol_pos, chol_transa
269 : #if defined(__parallel)
270 : INTEGER :: i
271 : INTEGER, DIMENSION(9) :: desca, descb, descout
272 : #endif
273 :
274 219261 : CALL timeset(routineN, handle)
275 :
276 219261 : context => matrix%matrix_struct%context
277 219261 : myprow = context%mepos(1)
278 219261 : mypcol = context%mepos(2)
279 219261 : n = matrix%matrix_struct%nrow_global
280 219261 : itype = 1
281 219261 : IF (op /= "SOLVE" .AND. op /= "MULTIPLY") &
282 0 : CPABORT("wrong argument op")
283 :
284 219261 : IF (PRESENT(pos)) THEN
285 71445 : SELECT CASE (pos)
286 : CASE ("LEFT")
287 71445 : chol_pos = 'L'
288 : CASE ("RIGHT")
289 70893 : chol_pos = 'R'
290 : CASE DEFAULT
291 142338 : CPABORT("wrong argument pos")
292 : END SELECT
293 : ELSE
294 76923 : chol_pos = 'L'
295 : END IF
296 :
297 219261 : chol_transa = 'N'
298 219261 : IF (PRESENT(transa)) chol_transa = transa
299 :
300 219261 : IF ((matrix%use_sp .NEQV. matrixb%use_sp) .OR. (matrix%use_sp .NEQV. matrixout%use_sp)) &
301 0 : CPABORT("not the same precision")
302 :
303 : ! notice b is the cholesky guy
304 219261 : a => matrix%local_data
305 219261 : b => matrixb%local_data
306 219261 : out => matrixout%local_data
307 219261 : a_sp => matrix%local_data_sp
308 219261 : b_sp => matrixb%local_data_sp
309 219261 : out_sp => matrixout%local_data_sp
310 :
311 : #if defined(__parallel)
312 :
313 2192610 : desca(:) = matrix%matrix_struct%descriptor(:)
314 2192610 : descb(:) = matrixb%matrix_struct%descriptor(:)
315 2192610 : descout(:) = matrixout%matrix_struct%descriptor(:)
316 219261 : alpha = 1.0_dp
317 4222601 : DO i = 1, neig
318 4222601 : IF (matrix%use_sp) THEN
319 0 : CALL pscopy(n, a_sp(1, 1), 1, i, desca, 1, out_sp(1, 1), 1, i, descout, 1)
320 : ELSE
321 4003340 : CALL pdcopy(n, a(1, 1), 1, i, desca, 1, out(1, 1), 1, i, descout, 1)
322 : END IF
323 : END DO
324 219261 : IF (op .EQ. "SOLVE") THEN
325 218247 : IF (matrix%use_sp) THEN
326 : CALL pstrsm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, &
327 0 : out_sp(1, 1), 1, 1, descout)
328 : ELSE
329 218247 : CALL pdtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
330 : END IF
331 : ELSE
332 1014 : IF (matrix%use_sp) THEN
333 : CALL pstrmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), 1, 1, descb, &
334 0 : out_sp(1, 1), 1, 1, descout)
335 : ELSE
336 1014 : CALL pdtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), 1, 1, descb, out(1, 1), 1, 1, descout)
337 : END IF
338 : END IF
339 : #else
340 :
341 : alpha = 1.0_dp
342 : IF (matrix%use_sp) THEN
343 : CALL scopy(neig*n, a_sp(1, 1), 1, out_sp(1, 1), 1)
344 : ELSE
345 : CALL dcopy(neig*n, a(1, 1), 1, out(1, 1), 1)
346 : END IF
347 : IF (op .EQ. "SOLVE") THEN
348 : IF (matrix%use_sp) THEN
349 : 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)
350 : ELSE
351 : CALL dtrsm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), SIZE(b, 1), out(1, 1), n)
352 : END IF
353 : ELSE
354 : IF (matrix%use_sp) THEN
355 : CALL strmm(chol_pos, 'U', chol_transa, 'N', n, neig, REAL(alpha, sp), b_sp(1, 1), n, out_sp(1, 1), n)
356 : ELSE
357 : CALL dtrmm(chol_pos, 'U', chol_transa, 'N', n, neig, alpha, b(1, 1), n, out(1, 1), n)
358 : END IF
359 : END IF
360 :
361 : #endif
362 :
363 219261 : CALL timestop(handle)
364 :
365 219261 : END SUBROUTINE cp_fm_cholesky_restore
366 :
367 : END MODULE cp_fm_cholesky
|