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, &
22 #include "./base/base_uses.f90"
37 TYPE(pao_env_type),
POINTER :: pao
39 CALL dbcsr_copy(pao%matrix_D, pao%matrix_G)
40 CALL dbcsr_set(pao%matrix_D, 0.0_dp)
42 CALL dbcsr_copy(pao%matrix_G_prev, pao%matrix_D)
44 IF (pao%precondition)
THEN
45 CALL dbcsr_copy(pao%matrix_D_preconed, pao%matrix_D)
49 CALL pao_opt_init_bfgs(pao)
57 SUBROUTINE pao_opt_init_bfgs(pao)
58 TYPE(pao_env_type),
POINTER :: pao
60 INTEGER,
DIMENSION(:),
POINTER :: nparams
62 CALL dbcsr_get_info(pao%matrix_X, row_blk_size=nparams)
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")
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)
74 END SUBROUTINE pao_opt_init_bfgs
81 TYPE(pao_env_type),
POINTER :: pao
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)
89 CALL dbcsr_release(pao%matrix_BFGS)
99 TYPE(pao_env_type),
POINTER :: pao
100 INTEGER,
INTENT(IN) :: icycle
102 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_opt_new_dir'
105 TYPE(dbcsr_type) :: matrix_g_preconed
107 CALL timeset(routinen, handle)
109 IF (pao%precondition)
THEN
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.)
120 CALL dbcsr_copy(pao%matrix_G_prev, matrix_g_preconed)
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)
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)
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
133 CALL timestop(handle)
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
150 SELECT CASE (pao%optimizer)
152 CALL pao_opt_newdir_cg(pao, icycle, matrix_g, matrix_g_prev, matrix_d)
154 CALL pao_opt_newdir_bfgs(pao, icycle, matrix_g, matrix_g_prev, matrix_d)
156 cpabort(
"PAO: unknown optimizer")
159 END SUBROUTINE pao_opt_new_dir_low
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
174 REAL(kind=
dp) :: beta, change, trace_d, trace_d_gnew, &
175 trace_g_mix, trace_g_new, trace_g_prev
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"
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
193 IF (trace_g_prev /= 0.0_dp)
THEN
194 beta = (trace_g_new - trace_g_mix)/trace_g_prev
197 IF (beta < 0.0_dp)
THEN
198 IF (pao%iw_opt > 0)
WRITE (pao%iw_opt, *)
"PAO|CG| resetting because beta < 0"
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"
210 IF (pao%iw_opt > 0)
WRITE (pao%iw_opt, *)
"PAO|CG| beta: ", beta
213 CALL dbcsr_add(matrix_d, matrix_g, beta, -1.0_dp)
215 END SUBROUTINE pao_opt_newdir_cg
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
230 CHARACTER(len=*),
PARAMETER :: routinen =
'pao_opt_newdir_bfgs'
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
240 CALL timeset(routinen, handle)
247 CALL dbcsr_copy(matrix_y, matrix_g)
248 CALL dbcsr_add(matrix_y, matrix_g_prev, 1.0_dp, -1.0_dp)
251 CALL dbcsr_copy(matrix_s, matrix_d)
252 CALL dbcsr_scale(matrix_s, pao%linesearch%step_size)
255 CALL dbcsr_dot(matrix_s, matrix_y, trace_sy)
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
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)
269 CALL dbcsr_dot(matrix_y, matrix_hy, trace_yhy)
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
282 CALL dbcsr_copy(matrix_r, matrix_s)
283 CALL dbcsr_add(matrix_r, matrix_hy, theta, (1.0_dp - theta))
286 CALL dbcsr_dot(matrix_r, matrix_y, trace_ry)
287 cpassert(trace_ry > 0.0_dp)
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)
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)
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)
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)
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)
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)
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)
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
335 IF (pao%iw_opt > 0)
WRITE (pao%iw_opt, *)
"PAO|BFGS| arnoldi of inv. Hessian did not converged."
340 CALL dbcsr_multiply(
"N",
"N", -1.0_dp, pao%matrix_BFGS, matrix_g, &
341 0.0_dp, matrix_d, retain_sparsity=.true.)
343 CALL timestop(handle)
344 END SUBROUTINE pao_opt_newdir_bfgs
arnoldi iteration using dbcsr
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...
Defines the basic variable types.
integer, parameter, public dp
Optimizers used by pao_main.F.
subroutine, public pao_opt_init(pao)
Initialize the optimizer.
subroutine, public pao_opt_finalize(pao)
Finalize the optimizer.
subroutine, public pao_opt_new_dir(pao, icycle)
Calculates the new search direction.
Types used by the PAO machinery.