(git:21ef868)
Loading...
Searching...
No Matches
pao_optimizer.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Optimizers used by pao_main.F
10!> \author Ole Schuett
11! **************************************************************************************************
12MODULE pao_optimizer
13 USE arnoldi_api, ONLY: arnoldi_extremal
14 USE cp_dbcsr_api, ONLY: &
15 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, dbcsr_release, &
16 dbcsr_scale, dbcsr_set, dbcsr_type
17 USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
18 dbcsr_dot,&
19 dbcsr_frobenius_norm,&
20 dbcsr_reserve_diag_blocks
21 USE kinds, ONLY: dp
22 USE pao_input, ONLY: pao_opt_bfgs,&
23 pao_opt_cg
24 USE pao_types, ONLY: pao_env_type
25#include "./base/base_uses.f90"
26
27 IMPLICIT NONE
28
29 PRIVATE
30
31 PUBLIC :: pao_opt_init, pao_opt_finalize, pao_opt_new_dir
32
33CONTAINS
34
35! **************************************************************************************************
36!> \brief Initialize the optimizer
37!> \param pao ...
38! **************************************************************************************************
39 SUBROUTINE pao_opt_init(pao)
40 TYPE(pao_env_type), POINTER :: pao
41
42 CALL dbcsr_copy(pao%matrix_D, pao%matrix_G)
43 CALL dbcsr_set(pao%matrix_D, 0.0_dp)
44
45 CALL dbcsr_copy(pao%matrix_G_prev, pao%matrix_D)
46
47 IF (pao%precondition) THEN
48 CALL dbcsr_copy(pao%matrix_D_preconed, pao%matrix_D)
49 END IF
50
51 IF (pao%optimizer == pao_opt_bfgs) THEN
52 CALL pao_opt_init_bfgs(pao)
53 END IF
54
55 END SUBROUTINE pao_opt_init
56
57! **************************************************************************************************
58!> \brief Initialize the BFGS optimizer
59!> \param pao ...
60! **************************************************************************************************
61 SUBROUTINE pao_opt_init_bfgs(pao)
62 TYPE(pao_env_type), POINTER :: pao
63
64 INTEGER, DIMENSION(:), POINTER :: nparams
65
66 CALL dbcsr_get_info(pao%matrix_X, row_blk_size=nparams)
67
68 CALL dbcsr_create(pao%matrix_BFGS, &
69 template=pao%matrix_X, &
70 row_blk_size=nparams, &
71 col_blk_size=nparams, &
72 name="PAO matrix_BFGS")
73
74 CALL dbcsr_reserve_diag_blocks(pao%matrix_BFGS)
75 CALL dbcsr_set(pao%matrix_BFGS, 0.0_dp)
76 CALL dbcsr_add_on_diag(pao%matrix_BFGS, 1.0_dp)
77
78 END SUBROUTINE pao_opt_init_bfgs
79
80! **************************************************************************************************
81!> \brief Finalize the optimizer
82!> \param pao ...
83! **************************************************************************************************
84 SUBROUTINE pao_opt_finalize(pao)
85 TYPE(pao_env_type), POINTER :: pao
86
87 CALL dbcsr_release(pao%matrix_D)
88 CALL dbcsr_release(pao%matrix_G_prev)
89 IF (pao%precondition) THEN
90 CALL dbcsr_release(pao%matrix_D_preconed)
91 END IF
92
93 IF (pao%optimizer == pao_opt_bfgs) THEN
94 CALL dbcsr_release(pao%matrix_BFGS)
95 END IF
96
97 END SUBROUTINE pao_opt_finalize
98
99! **************************************************************************************************
100!> \brief Calculates the new search direction.
101!> \param pao ...
102!> \param icycle ...
103! **************************************************************************************************
104 SUBROUTINE pao_opt_new_dir(pao, icycle)
105 TYPE(pao_env_type), POINTER :: pao
106 INTEGER, INTENT(IN) :: icycle
107
108 CHARACTER(len=*), PARAMETER :: routinen = 'pao_opt_new_dir'
109
110 INTEGER :: handle
111 TYPE(dbcsr_type) :: matrix_g_preconed
112
113 CALL timeset(routinen, handle)
114
115 IF (pao%precondition) THEN
116 ! We can't convert matrix_D for and back every time, the numeric noise would disturb the CG,
117 ! hence we keep matrix_D_preconed around.
118 CALL dbcsr_copy(matrix_g_preconed, pao%matrix_G)
119 CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_precon, pao%matrix_G, &
120 0.0_dp, matrix_g_preconed, retain_sparsity=.true.)
121 CALL pao_opt_new_dir_low(pao, icycle, matrix_g_preconed, pao%matrix_G_prev, pao%matrix_D_preconed)
122 CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_precon, pao%matrix_D_preconed, &
123 0.0_dp, pao%matrix_D, retain_sparsity=.true.)
124
125 ! store preconditioned gradient for next iteration
126 CALL dbcsr_copy(pao%matrix_G_prev, matrix_g_preconed)
127
128 pao%norm_G = dbcsr_frobenius_norm(matrix_g_preconed)
129 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| norm of preconditioned gradient:", pao%norm_G
130 CALL dbcsr_release(matrix_g_preconed)
131
132 ELSE
133 CALL pao_opt_new_dir_low(pao, icycle, pao%matrix_G, pao%matrix_G_prev, pao%matrix_D)
134 CALL dbcsr_copy(pao%matrix_G_prev, pao%matrix_G) ! store gradient for next iteration
135 pao%norm_G = dbcsr_frobenius_norm(pao%matrix_G)
136 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| norm of gradient:", pao%norm_G
137 END IF
138
139 CALL timestop(handle)
140
141 END SUBROUTINE pao_opt_new_dir
142
143! **************************************************************************************************
144!> \brief Calculates the new search direction.
145!> \param pao ...
146!> \param icycle ...
147!> \param matrix_G ...
148!> \param matrix_G_prev ...
149!> \param matrix_D ...
150! **************************************************************************************************
151 SUBROUTINE pao_opt_new_dir_low(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
152 TYPE(pao_env_type), POINTER :: pao
153 INTEGER, INTENT(IN) :: icycle
154 TYPE(dbcsr_type) :: matrix_g, matrix_g_prev, matrix_d
155
156 SELECT CASE (pao%optimizer)
157 CASE (pao_opt_cg)
158 CALL pao_opt_newdir_cg(pao, icycle, matrix_g, matrix_g_prev, matrix_d)
159 CASE (pao_opt_bfgs)
160 CALL pao_opt_newdir_bfgs(pao, icycle, matrix_g, matrix_g_prev, matrix_d)
161 CASE DEFAULT
162 cpabort("PAO: unknown optimizer")
163 END SELECT
164
165 END SUBROUTINE pao_opt_new_dir_low
166
167! **************************************************************************************************
168!> \brief Conjugate Gradient algorithm
169!> \param pao ...
170!> \param icycle ...
171!> \param matrix_G ...
172!> \param matrix_G_prev ...
173!> \param matrix_D ...
174! **************************************************************************************************
175 SUBROUTINE pao_opt_newdir_cg(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
176 TYPE(pao_env_type), POINTER :: pao
177 INTEGER, INTENT(IN) :: icycle
178 TYPE(dbcsr_type) :: matrix_g, matrix_g_prev, matrix_d
179
180 REAL(kind=dp) :: beta, change, trace_d, trace_d_gnew, &
181 trace_g_mix, trace_g_new, trace_g_prev
182
183 ! determine CG mixing factor
184 IF (icycle <= pao%cg_init_steps) THEN
185 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| warming up with steepest descent"
186 beta = 0.0_dp
187 ELSE
188 CALL dbcsr_dot(matrix_g, matrix_g, trace_g_new)
189 CALL dbcsr_dot(matrix_g_prev, matrix_g_prev, trace_g_prev)
190 CALL dbcsr_dot(matrix_g, matrix_g_prev, trace_g_mix)
191 CALL dbcsr_dot(matrix_d, matrix_g, trace_d_gnew)
192 CALL dbcsr_dot(matrix_d, matrix_d, trace_d)
193 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_new ", trace_g_new
194 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_prev ", trace_g_prev
195 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_mix ", trace_g_mix
196 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_D ", trace_d
197 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_D_Gnew", trace_d_gnew
198
199 IF (trace_g_prev /= 0.0_dp) THEN
200 beta = (trace_g_new - trace_g_mix)/trace_g_prev !Polak-Ribiere
201 END IF
202
203 IF (beta < 0.0_dp) THEN
204 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| resetting because beta < 0"
205 beta = 0.0_dp
206 END IF
207
208 change = trace_d_gnew**2/trace_d*trace_g_new
209 IF (change > pao%cg_reset_limit) THEN
210 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| resetting because change > CG_RESET_LIMIT"
211 beta = 0.0_dp
212 END IF
213
214 END IF
215
216 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| beta: ", beta
217
218 ! calculate new CG direction matrix_D
219 CALL dbcsr_add(matrix_d, matrix_g, beta, -1.0_dp)
220
221 END SUBROUTINE pao_opt_newdir_cg
222
223! **************************************************************************************************
224!> \brief Broyden-Fletcher-Goldfarb-Shanno algorithm
225!> \param pao ...
226!> \param icycle ...
227!> \param matrix_G ...
228!> \param matrix_G_prev ...
229!> \param matrix_D ...
230! **************************************************************************************************
231 SUBROUTINE pao_opt_newdir_bfgs(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
232 TYPE(pao_env_type), POINTER :: pao
233 INTEGER, INTENT(IN) :: icycle
234 TYPE(dbcsr_type) :: matrix_g, matrix_g_prev, matrix_d
235
236 CHARACTER(len=*), PARAMETER :: routinen = 'pao_opt_newdir_bfgs'
237
238 INTEGER :: handle
239 LOGICAL :: arnoldi_converged
240 REAL(dp) :: eval_max, eval_min, theta, trace_ry, &
241 trace_sy, trace_yhy, trace_yy
242 TYPE(dbcsr_type) :: matrix_hy, matrix_hyr, matrix_r, &
243 matrix_rr, matrix_ryh, matrix_ryhyr, &
244 matrix_s, matrix_y, matrix_yr
245
246 CALL timeset(routinen, handle)
247
248 !TODO add filtering?
249
250 ! Notation according to the book from Nocedal and Wright, see chapter 6.
251 IF (icycle > 1) THEN
252 ! y = G - G_prev
253 CALL dbcsr_copy(matrix_y, matrix_g)
254 CALL dbcsr_add(matrix_y, matrix_g_prev, 1.0_dp, -1.0_dp) ! dG
255
256 ! s = X - X_prev
257 CALL dbcsr_copy(matrix_s, matrix_d)
258 CALL dbcsr_scale(matrix_s, pao%linesearch%step_size) ! dX
259
260 ! sy = MATMUL(TRANPOSE(s), y)
261 CALL dbcsr_dot(matrix_s, matrix_y, trace_sy)
262
263 ! heuristic initialization
264 IF (icycle == 2) THEN
265 CALL dbcsr_dot(matrix_y, matrix_y, trace_yy)
266 CALL dbcsr_scale(pao%matrix_BFGS, trace_sy/trace_yy)
267 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| Initializing with:", trace_sy/trace_yy
268 END IF
269
270 ! Hy = MATMUL(H, y)
271 CALL dbcsr_create(matrix_hy, template=matrix_g, matrix_type="N")
272 CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_BFGS, matrix_y, 0.0_dp, matrix_hy)
273
274 ! yHy = MATMUL(TRANPOSE(y), Hy)
275 CALL dbcsr_dot(matrix_y, matrix_hy, trace_yhy)
276
277 ! Use damped BFGS algorithm to ensure H remains positive definite.
278 ! See chapter 18 in Nocedal and Wright's book for details.
279 ! The formulas were adopted to inverse Hessian algorithm.
280 IF (trace_sy < 0.2_dp*trace_yhy) THEN
281 theta = 0.8_dp*trace_yhy/(trace_yhy - trace_sy)
282 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| Dampening theta:", theta
283 ELSE
284 theta = 1.0
285 END IF
286
287 ! r = theta*s + (1-theta)*Hy
288 CALL dbcsr_copy(matrix_r, matrix_s)
289 CALL dbcsr_add(matrix_r, matrix_hy, theta, (1.0_dp - theta))
290
291 ! use t instead of y to update B matrix
292 CALL dbcsr_dot(matrix_r, matrix_y, trace_ry)
293 cpassert(trace_ry > 0.0_dp)
294
295 ! yr = MATMUL(y, TRANSPOSE(r))
296 CALL dbcsr_create(matrix_yr, template=pao%matrix_BFGS, matrix_type="N")
297 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_y, matrix_r, 0.0_dp, matrix_yr)
298
299 ! Hyr = MATMUL(H, yr)
300 CALL dbcsr_create(matrix_hyr, template=pao%matrix_BFGS, matrix_type="N")
301 CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_BFGS, matrix_yr, 0.0_dp, matrix_hyr)
302
303 ! ryH = MATMUL(TRANSPOSE(yr), H)
304 CALL dbcsr_create(matrix_ryh, template=pao%matrix_BFGS, matrix_type="N")
305 CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_yr, pao%matrix_BFGS, 0.0_dp, matrix_ryh)
306
307 ! ryHry = MATMUL(ryH,yr)
308 CALL dbcsr_create(matrix_ryhyr, template=pao%matrix_BFGS, matrix_type="N")
309 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ryh, matrix_yr, 0.0_dp, matrix_ryhyr)
310
311 ! rr = MATMUL(r,TRANSPOSE(r))
312 CALL dbcsr_create(matrix_rr, template=pao%matrix_BFGS, matrix_type="N")
313 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_r, matrix_r, 0.0_dp, matrix_rr)
314
315 ! H = H - Hyr/ry - ryH/ry + ryHyr/(ry**2) + rr/ry
316 CALL dbcsr_add(pao%matrix_BFGS, matrix_hyr, 1.0_dp, -1.0_dp/trace_ry)
317 CALL dbcsr_add(pao%matrix_BFGS, matrix_ryh, 1.0_dp, -1.0_dp/trace_ry)
318 CALL dbcsr_add(pao%matrix_BFGS, matrix_ryhyr, 1.0_dp, +1.0_dp/(trace_ry**2))
319 CALL dbcsr_add(pao%matrix_BFGS, matrix_rr, 1.0_dp, +1.0_dp/trace_ry)
320
321 ! clean up
322 CALL dbcsr_release(matrix_y)
323 CALL dbcsr_release(matrix_s)
324 CALL dbcsr_release(matrix_r)
325 CALL dbcsr_release(matrix_hy)
326 CALL dbcsr_release(matrix_yr)
327 CALL dbcsr_release(matrix_hyr)
328 CALL dbcsr_release(matrix_ryh)
329 CALL dbcsr_release(matrix_ryhyr)
330 CALL dbcsr_release(matrix_rr)
331 END IF
332
333 ! approximate condition of Hessian
334 !TODO: good setting for arnoldi?
335 CALL arnoldi_extremal(pao%matrix_BFGS, eval_max, eval_min, max_iter=100, &
336 threshold=1e-2_dp, converged=arnoldi_converged)
337 IF (arnoldi_converged) THEN
338 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| evals of inv. Hessian: min, max, max/min", &
339 eval_min, eval_max, eval_max/eval_min
340 ELSE
341 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| arnoldi of inv. Hessian did not converged."
342 END IF
343
344 ! calculate new direction
345 ! d = MATMUL(H, -g)
346 CALL dbcsr_multiply("N", "N", -1.0_dp, pao%matrix_BFGS, matrix_g, &
347 0.0_dp, matrix_d, retain_sparsity=.true.)
348
349 CALL timestop(handle)
350 END SUBROUTINE pao_opt_newdir_bfgs
351
352END MODULE pao_optimizer
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:16
arnoldi_api::arnoldi_extremal
subroutine, public arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface this hi...
Definition arnoldi_api.F:254
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_add_on_diag
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Definition cp_dbcsr_contrib.F:331
cp_dbcsr_contrib::dbcsr_frobenius_norm
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:127
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_contrib::dbcsr_reserve_diag_blocks
subroutine, public dbcsr_reserve_diag_blocks(matrix)
Reserves all diagonal blocks.
Definition cp_dbcsr_contrib.F:255
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
pao_input
Definition pao_input.F:8
pao_input::pao_opt_cg
integer, parameter, public pao_opt_cg
Definition pao_input.F:45
pao_input::pao_opt_bfgs
integer, parameter, public pao_opt_bfgs
Definition pao_input.F:45
pao_optimizer
Optimizers used by pao_main.F.
Definition pao_optimizer.F:12
pao_optimizer::pao_opt_init
subroutine, public pao_opt_init(pao)
Initialize the optimizer.
Definition pao_optimizer.F:40
pao_optimizer::pao_opt_finalize
subroutine, public pao_opt_finalize(pao)
Finalize the optimizer.
Definition pao_optimizer.F:85
pao_optimizer::pao_opt_new_dir
subroutine, public pao_opt_new_dir(pao, icycle)
Calculates the new search direction.
Definition pao_optimizer.F:105
pao_types
Types used by the PAO machinery.
Definition pao_types.F:12
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
pao_types::pao_env_type
Definition pao_types.F:148