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 Interface to (sca)lapack for the Cholesky based procedures
10 : !> \author VW
11 : !> \date 2009-09-08
12 : !> \version 0.8
13 : !>
14 : !> <b>Modification history:</b>
15 : !> - Created 2009-09-08
16 : ! **************************************************************************************************
17 : MODULE cp_dbcsr_cholesky
18 : USE cp_blacs_env, ONLY: cp_blacs_env_type
19 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
20 : copy_fm_to_dbcsr
21 : USE cp_fm_basic_linalg, ONLY: cp_fm_cholesky_restore,&
22 : cp_fm_potrf,&
23 : cp_fm_potri,&
24 : cp_fm_upper_to_full
25 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
26 : cp_fm_struct_release,&
27 : cp_fm_struct_type
28 : USE cp_fm_types, ONLY: cp_fm_create,&
29 : cp_fm_release,&
30 : cp_fm_type
31 : USE dbcsr_api, ONLY: dbcsr_get_info,&
32 : dbcsr_type
33 : USE message_passing, ONLY: mp_para_env_type
34 : #include "base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_dbcsr_cholesky'
39 :
40 : PUBLIC :: cp_dbcsr_cholesky_decompose, cp_dbcsr_cholesky_invert, &
41 : cp_dbcsr_cholesky_restore
42 :
43 : PRIVATE
44 :
45 : CONTAINS
46 :
47 : ! **************************************************************************************************
48 : !> \brief used to replace a symmetric positive def. matrix M with its cholesky
49 : !> decomposition U: M = U^T * U, with U upper triangular
50 : !> \param matrix the matrix to replace with its cholesky decomposition
51 : !> \param n the number of row (and columns) of the matrix &
52 : !> (defaults to the min(size(matrix)))
53 : !> \param para_env ...
54 : !> \param blacs_env ...
55 : !> \par History
56 : !> 05.2002 created [JVdV]
57 : !> 12.2002 updated, added n optional parm [fawzi]
58 : !> \author Joost
59 : ! **************************************************************************************************
60 20178 : SUBROUTINE cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
61 : TYPE(dbcsr_type) :: matrix
62 : INTEGER, INTENT(in), OPTIONAL :: n
63 : TYPE(mp_para_env_type), POINTER :: para_env
64 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
65 :
66 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_dbcsr_cholesky_decompose'
67 :
68 : INTEGER :: handle, my_n, nfullcols_total, &
69 : nfullrows_total
70 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
71 : TYPE(cp_fm_type) :: fm_matrix
72 :
73 10089 : CALL timeset(routineN, handle)
74 :
75 10089 : NULLIFY (fm_struct)
76 10089 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
77 :
78 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
79 10089 : ncol_global=nfullcols_total, para_env=para_env)
80 10089 : CALL cp_fm_create(fm_matrix, fm_struct, name="fm_matrix")
81 10089 : CALL cp_fm_struct_release(fm_struct)
82 :
83 10089 : CALL copy_dbcsr_to_fm(matrix, fm_matrix)
84 :
85 : my_n = MIN(fm_matrix%matrix_struct%nrow_global, &
86 10089 : fm_matrix%matrix_struct%ncol_global)
87 10089 : IF (PRESENT(n)) THEN
88 728 : CPASSERT(n <= my_n)
89 728 : my_n = n
90 : END IF
91 :
92 10089 : CALL cp_fm_potrf(fm_matrix, my_n)
93 :
94 10089 : CALL copy_fm_to_dbcsr(fm_matrix, matrix)
95 :
96 10089 : CALL cp_fm_release(fm_matrix)
97 :
98 10089 : CALL timestop(handle)
99 :
100 10089 : END SUBROUTINE cp_dbcsr_cholesky_decompose
101 :
102 : ! **************************************************************************************************
103 : !> \brief used to replace the cholesky decomposition by the inverse
104 : !> \param matrix the matrix to invert (must be an upper triangular matrix)
105 : !> \param n size of the matrix to invert (defaults to the min(size(matrix)))
106 : !> \param para_env ...
107 : !> \param blacs_env ...
108 : !> \param upper_to_full ...
109 : !> \par History
110 : !> 05.2002 created [JVdV]
111 : !> \author Joost VandeVondele
112 : ! **************************************************************************************************
113 18626 : SUBROUTINE cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, upper_to_full)
114 : TYPE(dbcsr_type) :: matrix
115 : INTEGER, INTENT(in), OPTIONAL :: n
116 : TYPE(mp_para_env_type), POINTER :: para_env
117 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
118 : LOGICAL, INTENT(IN) :: upper_to_full
119 :
120 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_dbcsr_cholesky_invert'
121 :
122 : INTEGER :: handle, my_n, nfullcols_total, &
123 : nfullrows_total
124 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
125 : TYPE(cp_fm_type) :: fm_matrix, fm_matrix_tmp
126 :
127 9313 : CALL timeset(routineN, handle)
128 :
129 9313 : NULLIFY (fm_struct)
130 9313 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
131 :
132 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
133 9313 : ncol_global=nfullrows_total, para_env=para_env)
134 9313 : CALL cp_fm_create(fm_matrix, fm_struct, name="fm_matrix")
135 9313 : CALL cp_fm_struct_release(fm_struct)
136 :
137 9313 : CALL copy_dbcsr_to_fm(matrix, fm_matrix)
138 :
139 : my_n = MIN(fm_matrix%matrix_struct%nrow_global, &
140 9313 : fm_matrix%matrix_struct%ncol_global)
141 9313 : IF (PRESENT(n)) THEN
142 0 : CPASSERT(n <= my_n)
143 0 : my_n = n
144 : END IF
145 :
146 9313 : CALL cp_fm_potri(fm_matrix, my_n)
147 :
148 9313 : IF (upper_to_full) THEN
149 9313 : CALL cp_fm_create(fm_matrix_tmp, fm_matrix%matrix_struct, name="fm_matrix_tmp")
150 9313 : CALL cp_fm_upper_to_full(fm_matrix, fm_matrix_tmp)
151 9313 : CALL cp_fm_release(fm_matrix_tmp)
152 : END IF
153 :
154 9313 : CALL copy_fm_to_dbcsr(fm_matrix, matrix)
155 :
156 9313 : CALL cp_fm_release(fm_matrix)
157 :
158 9313 : CALL timestop(handle)
159 :
160 9313 : END SUBROUTINE cp_dbcsr_cholesky_invert
161 :
162 : ! **************************************************************************************************
163 : !> \brief ...
164 : !> \param matrix ...
165 : !> \param neig ...
166 : !> \param matrixb ...
167 : !> \param matrixout ...
168 : !> \param op ...
169 : !> \param pos ...
170 : !> \param transa ...
171 : !> \param para_env ...
172 : !> \param blacs_env ...
173 : ! **************************************************************************************************
174 4736 : SUBROUTINE cp_dbcsr_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa, &
175 : para_env, blacs_env)
176 : TYPE(dbcsr_type) :: matrix
177 : INTEGER, INTENT(IN) :: neig
178 : TYPE(dbcsr_type) :: matrixb, matrixout
179 : CHARACTER(LEN=*), INTENT(IN) :: op
180 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pos, transa
181 : TYPE(mp_para_env_type), POINTER :: para_env
182 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
183 :
184 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_dbcsr_cholesky_restore'
185 :
186 : CHARACTER :: chol_pos, chol_transa
187 : INTEGER :: handle, nfullcols_total, nfullrows_total
188 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
189 : TYPE(cp_fm_type) :: fm_matrix, fm_matrixb, fm_matrixout
190 :
191 1184 : CALL timeset(routineN, handle)
192 :
193 1184 : NULLIFY (fm_struct)
194 :
195 1184 : CALL dbcsr_get_info(matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
196 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
197 1184 : ncol_global=nfullcols_total, para_env=para_env)
198 1184 : CALL cp_fm_create(fm_matrix, fm_struct, name="fm_matrix")
199 1184 : CALL cp_fm_struct_release(fm_struct)
200 :
201 1184 : CALL dbcsr_get_info(matrixb, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
202 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
203 1184 : ncol_global=nfullcols_total, para_env=para_env)
204 1184 : CALL cp_fm_create(fm_matrixb, fm_struct, name="fm_matrixb")
205 1184 : CALL cp_fm_struct_release(fm_struct)
206 :
207 1184 : CALL dbcsr_get_info(matrixout, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
208 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
209 1184 : ncol_global=nfullcols_total, para_env=para_env)
210 1184 : CALL cp_fm_create(fm_matrixout, fm_struct, name="fm_matrixout")
211 1184 : CALL cp_fm_struct_release(fm_struct)
212 :
213 1184 : CALL copy_dbcsr_to_fm(matrix, fm_matrix)
214 1184 : CALL copy_dbcsr_to_fm(matrixb, fm_matrixb)
215 : !CALL copy_dbcsr_to_fm(matrixout, fm_matrixout)
216 :
217 1184 : IF (op /= "SOLVE" .AND. op /= "MULTIPLY") &
218 0 : CPABORT("wrong argument op")
219 :
220 1184 : IF (PRESENT(pos)) THEN
221 0 : SELECT CASE (pos)
222 : CASE ("LEFT")
223 0 : chol_pos = 'L'
224 : CASE ("RIGHT")
225 1184 : chol_pos = 'R'
226 : CASE DEFAULT
227 1184 : CPABORT("wrong argument pos")
228 : END SELECT
229 : ELSE
230 0 : chol_pos = 'L'
231 : END IF
232 :
233 1184 : chol_transa = 'N'
234 1184 : IF (PRESENT(transa)) chol_transa = transa
235 :
236 1184 : IF ((fm_matrix%use_sp .NEQV. fm_matrixb%use_sp) .OR. (fm_matrix%use_sp .NEQV. fm_matrixout%use_sp)) &
237 0 : CPABORT("not the same precision")
238 :
239 1184 : CALL cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, chol_pos, chol_transa)
240 :
241 1184 : CALL copy_fm_to_dbcsr(fm_matrixout, matrixout)
242 :
243 1184 : CALL cp_fm_release(fm_matrix)
244 1184 : CALL cp_fm_release(fm_matrixb)
245 1184 : CALL cp_fm_release(fm_matrixout)
246 :
247 1184 : CALL timestop(handle)
248 :
249 1184 : END SUBROUTINE cp_dbcsr_cholesky_restore
250 :
251 : END MODULE cp_dbcsr_cholesky
252 :
|
