(git:ed6f26b)
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-2025 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) &
52 CALL pao_opt_init_bfgs(pao)
53
54 END SUBROUTINE pao_opt_init
55
56! **************************************************************************************************
57!> \brief Initialize the BFGS optimizer
58!> \param pao ...
59! **************************************************************************************************
60 SUBROUTINE pao_opt_init_bfgs(pao)
61 TYPE(pao_env_type), POINTER :: pao
62
63 INTEGER, DIMENSION(:), POINTER :: nparams
64
65 CALL dbcsr_get_info(pao%matrix_X, row_blk_size=nparams)
66
67 CALL dbcsr_create(pao%matrix_BFGS, &
68 template=pao%matrix_X, &
69 row_blk_size=nparams, &
70 col_blk_size=nparams, &
71 name="PAO matrix_BFGS")
72
73 CALL dbcsr_reserve_diag_blocks(pao%matrix_BFGS)
74 CALL dbcsr_set(pao%matrix_BFGS, 0.0_dp)
75 CALL dbcsr_add_on_diag(pao%matrix_BFGS, 1.0_dp)
76
77 END SUBROUTINE pao_opt_init_bfgs
78
79! **************************************************************************************************
80!> \brief Finalize the optimizer
81!> \param pao ...
82! **************************************************************************************************
83 SUBROUTINE pao_opt_finalize(pao)
84 TYPE(pao_env_type), POINTER :: pao
85
86 CALL dbcsr_release(pao%matrix_D)
87 CALL dbcsr_release(pao%matrix_G_prev)
88 IF (pao%precondition) &
89 CALL dbcsr_release(pao%matrix_D_preconed)
90
91 IF (pao%optimizer == pao_opt_bfgs) &
92 CALL dbcsr_release(pao%matrix_BFGS)
93
94 END SUBROUTINE pao_opt_finalize
95
96! **************************************************************************************************
97!> \brief Calculates the new search direction.
98!> \param pao ...
99!> \param icycle ...
100! **************************************************************************************************
101 SUBROUTINE pao_opt_new_dir(pao, icycle)
102 TYPE(pao_env_type), POINTER :: pao
103 INTEGER, INTENT(IN) :: icycle
104
105 CHARACTER(len=*), PARAMETER :: routinen = 'pao_opt_new_dir'
106
107 INTEGER :: handle
108 TYPE(dbcsr_type) :: matrix_g_preconed
109
110 CALL timeset(routinen, handle)
111
112 IF (pao%precondition) THEN
113 ! We can't convert matrix_D for and back every time, the numeric noise would disturb the CG,
114 ! hence we keep matrix_D_preconed around.
115 CALL dbcsr_copy(matrix_g_preconed, pao%matrix_G)
116 CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_precon, pao%matrix_G, &
117 0.0_dp, matrix_g_preconed, retain_sparsity=.true.)
118 CALL pao_opt_new_dir_low(pao, icycle, matrix_g_preconed, pao%matrix_G_prev, pao%matrix_D_preconed)
119 CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_precon, pao%matrix_D_preconed, &
120 0.0_dp, pao%matrix_D, retain_sparsity=.true.)
121
122 ! store preconditioned gradient for next iteration
123 CALL dbcsr_copy(pao%matrix_G_prev, matrix_g_preconed)
124
125 pao%norm_G = dbcsr_frobenius_norm(matrix_g_preconed)
126 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| norm of preconditioned gradient:", pao%norm_G
127 CALL dbcsr_release(matrix_g_preconed)
128
129 ELSE
130 CALL pao_opt_new_dir_low(pao, icycle, pao%matrix_G, pao%matrix_G_prev, pao%matrix_D)
131 CALL dbcsr_copy(pao%matrix_G_prev, pao%matrix_G) ! store gradient for next iteration
132 pao%norm_G = dbcsr_frobenius_norm(pao%matrix_G)
133 IF (pao%iw > 0) WRITE (pao%iw, *) "PAO| norm of gradient:", pao%norm_G
134 END IF
135
136 CALL timestop(handle)
137
138 END SUBROUTINE pao_opt_new_dir
139
140! **************************************************************************************************
141!> \brief Calculates the new search direction.
142!> \param pao ...
143!> \param icycle ...
144!> \param matrix_G ...
145!> \param matrix_G_prev ...
146!> \param matrix_D ...
147! **************************************************************************************************
148 SUBROUTINE pao_opt_new_dir_low(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
149 TYPE(pao_env_type), POINTER :: pao
150 INTEGER, INTENT(IN) :: icycle
151 TYPE(dbcsr_type) :: matrix_g, matrix_g_prev, matrix_d
152
153 SELECT CASE (pao%optimizer)
154 CASE (pao_opt_cg)
155 CALL pao_opt_newdir_cg(pao, icycle, matrix_g, matrix_g_prev, matrix_d)
156 CASE (pao_opt_bfgs)
157 CALL pao_opt_newdir_bfgs(pao, icycle, matrix_g, matrix_g_prev, matrix_d)
158 CASE DEFAULT
159 cpabort("PAO: unknown optimizer")
160 END SELECT
161
162 END SUBROUTINE pao_opt_new_dir_low
163
164! **************************************************************************************************
165!> \brief Conjugate Gradient algorithm
166!> \param pao ...
167!> \param icycle ...
168!> \param matrix_G ...
169!> \param matrix_G_prev ...
170!> \param matrix_D ...
171! **************************************************************************************************
172 SUBROUTINE pao_opt_newdir_cg(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
173 TYPE(pao_env_type), POINTER :: pao
174 INTEGER, INTENT(IN) :: icycle
175 TYPE(dbcsr_type) :: matrix_g, matrix_g_prev, matrix_d
176
177 REAL(kind=dp) :: beta, change, trace_d, trace_d_gnew, &
178 trace_g_mix, trace_g_new, trace_g_prev
179
180 ! determine CG mixing factor
181 IF (icycle <= pao%cg_init_steps) THEN
182 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| warming up with steepest descent"
183 beta = 0.0_dp
184 ELSE
185 CALL dbcsr_dot(matrix_g, matrix_g, trace_g_new)
186 CALL dbcsr_dot(matrix_g_prev, matrix_g_prev, trace_g_prev)
187 CALL dbcsr_dot(matrix_g, matrix_g_prev, trace_g_mix)
188 CALL dbcsr_dot(matrix_d, matrix_g, trace_d_gnew)
189 CALL dbcsr_dot(matrix_d, matrix_d, trace_d)
190 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_new ", trace_g_new
191 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_prev ", trace_g_prev
192 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_G_mix ", trace_g_mix
193 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_D ", trace_d
194 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| trace_D_Gnew", trace_d_gnew
195
196 IF (trace_g_prev /= 0.0_dp) THEN
197 beta = (trace_g_new - trace_g_mix)/trace_g_prev !Polak-Ribiere
198 END IF
199
200 IF (beta < 0.0_dp) THEN
201 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| resetting because beta < 0"
202 beta = 0.0_dp
203 END IF
204
205 change = trace_d_gnew**2/trace_d*trace_g_new
206 IF (change > pao%cg_reset_limit) THEN
207 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| resetting because change > CG_RESET_LIMIT"
208 beta = 0.0_dp
209 END IF
210
211 END IF
212
213 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|CG| beta: ", beta
214
215 ! calculate new CG direction matrix_D
216 CALL dbcsr_add(matrix_d, matrix_g, beta, -1.0_dp)
217
218 END SUBROUTINE pao_opt_newdir_cg
219
220! **************************************************************************************************
221!> \brief Broyden-Fletcher-Goldfarb-Shanno algorithm
222!> \param pao ...
223!> \param icycle ...
224!> \param matrix_G ...
225!> \param matrix_G_prev ...
226!> \param matrix_D ...
227! **************************************************************************************************
228 SUBROUTINE pao_opt_newdir_bfgs(pao, icycle, matrix_G, matrix_G_prev, matrix_D)
229 TYPE(pao_env_type), POINTER :: pao
230 INTEGER, INTENT(IN) :: icycle
231 TYPE(dbcsr_type) :: matrix_g, matrix_g_prev, matrix_d
232
233 CHARACTER(len=*), PARAMETER :: routinen = 'pao_opt_newdir_bfgs'
234
235 INTEGER :: handle
236 LOGICAL :: arnoldi_converged
237 REAL(dp) :: eval_max, eval_min, theta, trace_ry, &
238 trace_sy, trace_yhy, trace_yy
239 TYPE(dbcsr_type) :: matrix_hy, matrix_hyr, matrix_r, &
240 matrix_rr, matrix_ryh, matrix_ryhyr, &
241 matrix_s, matrix_y, matrix_yr
242
243 CALL timeset(routinen, handle)
244
245 !TODO add filtering?
246
247 ! Notation according to the book from Nocedal and Wright, see chapter 6.
248 IF (icycle > 1) THEN
249 ! y = G - G_prev
250 CALL dbcsr_copy(matrix_y, matrix_g)
251 CALL dbcsr_add(matrix_y, matrix_g_prev, 1.0_dp, -1.0_dp) ! dG
252
253 ! s = X - X_prev
254 CALL dbcsr_copy(matrix_s, matrix_d)
255 CALL dbcsr_scale(matrix_s, pao%linesearch%step_size) ! dX
256
257 ! sy = MATMUL(TRANPOSE(s), y)
258 CALL dbcsr_dot(matrix_s, matrix_y, trace_sy)
259
260 ! heuristic initialization
261 IF (icycle == 2) THEN
262 CALL dbcsr_dot(matrix_y, matrix_y, trace_yy)
263 CALL dbcsr_scale(pao%matrix_BFGS, trace_sy/trace_yy)
264 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| Initializing with:", trace_sy/trace_yy
265 END IF
266
267 ! Hy = MATMUL(H, y)
268 CALL dbcsr_create(matrix_hy, template=matrix_g, matrix_type="N")
269 CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_BFGS, matrix_y, 0.0_dp, matrix_hy)
270
271 ! yHy = MATMUL(TRANPOSE(y), Hy)
272 CALL dbcsr_dot(matrix_y, matrix_hy, trace_yhy)
273
274 ! Use damped BFGS algorithm to ensure H remains positive definite.
275 ! See chapter 18 in Nocedal and Wright's book for details.
276 ! The formulas were adopted to inverse Hessian algorithm.
277 IF (trace_sy < 0.2_dp*trace_yhy) THEN
278 theta = 0.8_dp*trace_yhy/(trace_yhy - trace_sy)
279 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| Dampening theta:", theta
280 ELSE
281 theta = 1.0
282 END IF
283
284 ! r = theta*s + (1-theta)*Hy
285 CALL dbcsr_copy(matrix_r, matrix_s)
286 CALL dbcsr_add(matrix_r, matrix_hy, theta, (1.0_dp - theta))
287
288 ! use t instead of y to update B matrix
289 CALL dbcsr_dot(matrix_r, matrix_y, trace_ry)
290 cpassert(trace_ry > 0.0_dp)
291
292 ! yr = MATMUL(y, TRANSPOSE(r))
293 CALL dbcsr_create(matrix_yr, template=pao%matrix_BFGS, matrix_type="N")
294 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_y, matrix_r, 0.0_dp, matrix_yr)
295
296 ! Hyr = MATMUL(H, yr)
297 CALL dbcsr_create(matrix_hyr, template=pao%matrix_BFGS, matrix_type="N")
298 CALL dbcsr_multiply("N", "N", 1.0_dp, pao%matrix_BFGS, matrix_yr, 0.0_dp, matrix_hyr)
299
300 ! ryH = MATMUL(TRANSPOSE(yr), H)
301 CALL dbcsr_create(matrix_ryh, template=pao%matrix_BFGS, matrix_type="N")
302 CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_yr, pao%matrix_BFGS, 0.0_dp, matrix_ryh)
303
304 ! ryHry = MATMUL(ryH,yr)
305 CALL dbcsr_create(matrix_ryhyr, template=pao%matrix_BFGS, matrix_type="N")
306 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ryh, matrix_yr, 0.0_dp, matrix_ryhyr)
307
308 ! rr = MATMUL(r,TRANSPOSE(r))
309 CALL dbcsr_create(matrix_rr, template=pao%matrix_BFGS, matrix_type="N")
310 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_r, matrix_r, 0.0_dp, matrix_rr)
311
312 ! H = H - Hyr/ry - ryH/ry + ryHyr/(ry**2) + rr/ry
313 CALL dbcsr_add(pao%matrix_BFGS, matrix_hyr, 1.0_dp, -1.0_dp/trace_ry)
314 CALL dbcsr_add(pao%matrix_BFGS, matrix_ryh, 1.0_dp, -1.0_dp/trace_ry)
315 CALL dbcsr_add(pao%matrix_BFGS, matrix_ryhyr, 1.0_dp, +1.0_dp/(trace_ry**2))
316 CALL dbcsr_add(pao%matrix_BFGS, matrix_rr, 1.0_dp, +1.0_dp/trace_ry)
317
318 ! clean up
319 CALL dbcsr_release(matrix_y)
320 CALL dbcsr_release(matrix_s)
321 CALL dbcsr_release(matrix_r)
322 CALL dbcsr_release(matrix_hy)
323 CALL dbcsr_release(matrix_yr)
324 CALL dbcsr_release(matrix_hyr)
325 CALL dbcsr_release(matrix_ryh)
326 CALL dbcsr_release(matrix_ryhyr)
327 CALL dbcsr_release(matrix_rr)
328 END IF
329
330 ! approximate condition of Hessian
331 !TODO: good setting for arnoldi?
332 CALL arnoldi_extremal(pao%matrix_BFGS, eval_max, eval_min, max_iter=100, &
333 threshold=1e-2_dp, converged=arnoldi_converged)
334 IF (arnoldi_converged) THEN
335 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| evals of inv. Hessian: min, max, max/min", &
336 eval_min, eval_max, eval_max/eval_min
337 ELSE
338 IF (pao%iw_opt > 0) WRITE (pao%iw_opt, *) "PAO|BFGS| arnoldi of inv. Hessian did not converged."
339 END IF
340
341 ! calculate new direction
342 ! d = MATMUL(H, -g)
343 CALL dbcsr_multiply("N", "N", -1.0_dp, pao%matrix_BFGS, matrix_g, &
344 0.0_dp, matrix_d, retain_sparsity=.true.)
345
346 CALL timestop(handle)
347 END SUBROUTINE pao_opt_newdir_bfgs
348
349END MODULE pao_optimizer
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
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:1165
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
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:1073
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:794
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
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:84
pao_optimizer::pao_opt_new_dir
subroutine, public pao_opt_new_dir(pao, icycle)
Calculates the new search direction.
Definition pao_optimizer.F:102
pao_types
Types used by the PAO machinery.
Definition pao_types.F:12
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
pao_types::pao_env_type
Definition pao_types.F:147