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