Line data Source code
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 Small, matrix-free mathematical kernels used by ADIIS.
10 : !>
11 : !> The quadratic model uses the convention
12 : !>
13 : !> f(c) = 1/2 c^T hessian c + linear^T c,
14 : !>
15 : !> with c on the probability simplex. Only the symmetric part of the
16 : !> supplied Hessian contributes to the objective.
17 : ! **************************************************************************************************
18 : MODULE qs_scf_subspace_math
19 :
20 : USE ieee_arithmetic, ONLY: ieee_is_finite
21 : USE kinds, ONLY: dp
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 :
27 : INTEGER, PARAMETER, PUBLIC :: simplex_qp_success = 0
28 : INTEGER, PARAMETER, PUBLIC :: simplex_qp_invalid_shape = 1
29 : INTEGER, PARAMETER, PUBLIC :: simplex_qp_nonfinite_input = 2
30 : INTEGER, PARAMETER, PUBLIC :: simplex_qp_iteration_limit = 3
31 : INTEGER, PARAMETER, PUBLIC :: simplex_qp_pairwise_stationary = 4
32 :
33 : INTEGER, PARAMETER :: max_exact_dimension = 12
34 :
35 : PUBLIC :: qs_scf_subspace_build_adiis_model, &
36 : qs_scf_subspace_fifo_slot, &
37 : simplex_quadratic_minimize
38 :
39 : CONTAINS
40 :
41 : ! **************************************************************************************************
42 : !> \brief Minimize a quadratic model over the probability simplex.
43 : !> \param hessian Hessian in f(c) = 1/2 c^T hessian c + linear^T c.
44 : !> \param linear Linear part of the quadratic model.
45 : !> \param coeff Resulting non-negative coefficients, normalized to sum to one.
46 : !> \param objective Objective value at coeff.
47 : !> \param status Completion status. Invalid/nonfinite inputs use the preferred-point fallback;
48 : !> iteration limits return the best legal point found.
49 : !> \param preferred_index History entry preferred for fallbacks and numerical ties; defaults to n.
50 : !>
51 : !> SCF histories are normally very small. Up to max_exact_dimension, all
52 : !> active faces are enumerated. A global minimizer of a quadratic over a
53 : !> simplex is a stationary point in the relative interior of one of these
54 : !> faces (a singular stationary face always contains an equivalent boundary
55 : !> solution). Larger problems use multistart, simplex-preserving pairwise
56 : !> line searches and return simplex_qp_pairwise_stationary. This only proves
57 : !> that no improving two-coordinate transfer was found; for an indefinite
58 : !> quadratic it is a heuristic result, not a general local/global optimum.
59 : ! **************************************************************************************************
60 140 : SUBROUTINE simplex_quadratic_minimize(hessian, linear, coeff, objective, status, preferred_index)
61 :
62 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hessian
63 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: linear
64 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: coeff
65 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: objective
66 : INTEGER, INTENT(OUT), OPTIONAL :: status
67 : INTEGER, INTENT(IN), OPTIONAL :: preferred_index
68 :
69 : INTEGER :: local_status, n, preferred
70 : REAL(KIND=dp) :: best_value, model_scale
71 140 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: working_linear
72 140 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: symmetric_hessian, working_hessian
73 :
74 140 : n = SIZE(linear)
75 140 : preferred = n
76 140 : local_status = simplex_qp_success
77 : best_value = HUGE(0.0_dp)
78 :
79 866 : coeff = 0.0_dp
80 140 : IF (PRESENT(preferred_index)) THEN
81 140 : IF (preferred_index >= 1 .AND. preferred_index <= n) THEN
82 140 : preferred = preferred_index
83 : ELSE
84 : local_status = simplex_qp_invalid_shape
85 0 : IF (SIZE(coeff) > 0) coeff(SIZE(coeff)) = 1.0_dp
86 0 : IF (PRESENT(objective)) objective = best_value
87 0 : IF (PRESENT(status)) status = local_status
88 0 : RETURN
89 : END IF
90 : END IF
91 140 : IF (preferred >= 1 .AND. preferred <= SIZE(coeff)) coeff(preferred) = 1.0_dp
92 :
93 140 : IF (n < 1 .OR. SIZE(coeff) /= n .OR. SIZE(hessian, 1) /= n .OR. SIZE(hessian, 2) /= n) THEN
94 : local_status = simplex_qp_invalid_shape
95 0 : IF (PRESENT(objective)) objective = best_value
96 0 : IF (PRESENT(status)) status = local_status
97 : RETURN
98 : END IF
99 :
100 6726 : IF (.NOT. ALL(ieee_is_finite(hessian)) .OR. .NOT. ALL(ieee_is_finite(linear))) THEN
101 : local_status = simplex_qp_nonfinite_input
102 0 : IF (PRESENT(objective)) objective = best_value
103 0 : IF (PRESENT(status)) status = local_status
104 0 : RETURN
105 : END IF
106 :
107 1120 : ALLOCATE (symmetric_hessian(n, n), working_hessian(n, n), working_linear(n))
108 5860 : symmetric_hessian(:, :) = 0.5_dp*hessian + 0.5_dp*TRANSPOSE(hessian)
109 :
110 : ! Remove terms that are constant on the simplex, then normalize the remaining model.
111 : ! Besides improving the KKT conditioning, this makes coefficient selection invariant to
112 5860 : working_hessian(:, :) = symmetric_hessian
113 866 : working_linear(:) = linear
114 6726 : model_scale = MAX(MAXVAL(ABS(working_hessian)), MAXVAL(ABS(working_linear)))
115 : ! Scale first only at extreme magnitudes where subtracting opposite-signed finite
116 : ! values could overflow. At ordinary scales, removing the gauge first preserves
117 140 : IF (model_scale > 0.25_dp*HUGE(0.0_dp)) THEN
118 0 : working_hessian(:, :) = working_hessian/model_scale
119 0 : working_linear(:) = working_linear/model_scale
120 : END IF
121 : ! Subtracting one scalar from every Hessian element is a H_nn*1*1^T gauge
122 : ! shift and is therefore constant when sum(coeff)=1.
123 5860 : working_hessian(:, :) = working_hessian - working_hessian(n, n)
124 866 : working_linear(:) = working_linear - working_linear(n)
125 6726 : model_scale = MAX(MAXVAL(ABS(working_hessian)), MAXVAL(ABS(working_linear)))
126 140 : IF (model_scale > 0.0_dp) THEN
127 5812 : working_hessian(:, :) = working_hessian/model_scale
128 834 : working_linear(:) = working_linear/model_scale
129 : END IF
130 :
131 140 : CALL best_vertex(working_hessian, working_linear, preferred, coeff, best_value)
132 :
133 140 : IF (n <= max_exact_dimension) THEN
134 138 : CALL enumerate_active_faces(working_hessian, working_linear, preferred, coeff, best_value)
135 : ELSE
136 2 : CALL pairwise_minimize(working_hessian, working_linear, preferred, coeff, best_value, local_status)
137 : END IF
138 :
139 : ! Remove harmless round-off at active constraints before returning.
140 866 : coeff = MAX(coeff, 0.0_dp)
141 1592 : coeff = coeff/SUM(coeff)
142 140 : best_value = quadratic_value(symmetric_hessian, linear, coeff)
143 :
144 140 : IF (PRESENT(objective)) objective = best_value
145 140 : IF (PRESENT(status)) status = local_status
146 :
147 140 : END SUBROUTINE simplex_quadratic_minimize
148 :
149 : ! **************************************************************************************************
150 : !> \brief Construct the canonical ADIIS model relative to the newest history entry.
151 : !> \param pf_metric Matrix T_ij=Tr(P_i F_j), already summed over spin/k-points as needed.
152 : !> \param newest Index of the reference history entry.
153 : !> \param hessian Hessian of the canonical simplex quadratic model.
154 : !> \param linear Linear term of the canonical simplex quadratic model.
155 : !> \param valid Whether all dimensions and input values were valid.
156 : ! **************************************************************************************************
157 140 : SUBROUTINE qs_scf_subspace_build_adiis_model(pf_metric, newest, hessian, linear, valid)
158 :
159 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pf_metric
160 : INTEGER, INTENT(IN) :: newest
161 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: hessian
162 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: linear
163 : LOGICAL, INTENT(OUT), OPTIONAL :: valid
164 :
165 : INTEGER :: i, j, n
166 : LOGICAL :: local_valid
167 : REAL(KIND=dp) :: model_ij, model_ji, tnn
168 :
169 140 : n = SIZE(linear)
170 5860 : hessian = 0.0_dp
171 866 : linear = 0.0_dp
172 :
173 : local_valid = n > 0 .AND. newest >= 1 .AND. newest <= n .AND. &
174 : SIZE(pf_metric, 1) == n .AND. SIZE(pf_metric, 2) == n .AND. &
175 140 : SIZE(hessian, 1) == n .AND. SIZE(hessian, 2) == n
176 5860 : IF (local_valid) local_valid = ALL(ieee_is_finite(pf_metric))
177 140 : IF (.NOT. local_valid) THEN
178 0 : IF (PRESENT(valid)) valid = .FALSE.
179 0 : RETURN
180 : END IF
181 :
182 140 : tnn = pf_metric(newest, newest)
183 866 : DO i = 1, n
184 866 : linear(i) = 2.0_dp*(pf_metric(i, newest) - tnn)
185 : END DO
186 :
187 866 : DO j = 1, n
188 5860 : DO i = 1, n
189 : model_ij = (pf_metric(i, j) - pf_metric(i, newest)) + &
190 4994 : (tnn - pf_metric(newest, j))
191 : model_ji = (pf_metric(j, i) - pf_metric(j, newest)) + &
192 4994 : (tnn - pf_metric(newest, i))
193 5720 : hessian(i, j) = model_ij + model_ji
194 : END DO
195 : END DO
196 :
197 140 : IF (PRESENT(valid)) valid = .TRUE.
198 :
199 : END SUBROUTINE qs_scf_subspace_build_adiis_model
200 :
201 : ! **************************************************************************************************
202 : !> \brief Select the next physical slot for a strict FIFO history.
203 : !> \param generation Insertion generation for each physical slot; zero denotes an empty slot.
204 : !> \return First empty slot, otherwise the oldest occupied slot; zero for invalid input.
205 : ! **************************************************************************************************
206 140 : PURE FUNCTION qs_scf_subspace_fifo_slot(generation) RESULT(slot)
207 :
208 : INTEGER, DIMENSION(:), INTENT(IN) :: generation
209 : INTEGER :: slot
210 :
211 : INTEGER :: i
212 :
213 140 : slot = 0
214 140 : IF (SIZE(generation) < 1) RETURN
215 2380 : IF (ANY(generation < 0)) RETURN
216 :
217 726 : DO i = 1, SIZE(generation)
218 726 : IF (generation(i) == 0) THEN
219 140 : slot = i
220 : RETURN
221 : END IF
222 : END DO
223 0 : slot = MINLOC(generation, DIM=1)
224 :
225 0 : END FUNCTION qs_scf_subspace_fifo_slot
226 :
227 : ! **************************************************************************************************
228 : !> \brief Select the best simplex vertex, preferring the newest one on ties.
229 : !> \param hessian Symmetric quadratic-model Hessian.
230 : !> \param linear Quadratic-model linear term.
231 : !> \param preferred History index favored when vertex values are tied.
232 : !> \param coeff Coefficients of the selected vertex.
233 : !> \param best_value Objective value at the selected vertex.
234 : ! **************************************************************************************************
235 140 : SUBROUTINE best_vertex(hessian, linear, preferred, coeff, best_value)
236 :
237 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hessian
238 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: linear
239 : INTEGER, INTENT(IN) :: preferred
240 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: coeff
241 : REAL(KIND=dp), INTENT(OUT) :: best_value
242 :
243 : INTEGER :: i, n
244 : REAL(KIND=dp) :: value
245 :
246 140 : n = SIZE(linear)
247 866 : coeff = 0.0_dp
248 140 : coeff(preferred) = 1.0_dp
249 140 : best_value = 0.5_dp*hessian(preferred, preferred) + linear(preferred)
250 :
251 866 : DO i = n, 1, -1
252 726 : IF (i == preferred) CYCLE
253 586 : value = 0.5_dp*hessian(i, i) + linear(i)
254 726 : IF (strictly_better(value, best_value)) THEN
255 178 : coeff = 0.0_dp
256 30 : coeff(i) = 1.0_dp
257 30 : best_value = value
258 : END IF
259 : END DO
260 :
261 140 : END SUBROUTINE best_vertex
262 :
263 : ! **************************************************************************************************
264 : !> \brief Enumerate stationary points on every active face of a small simplex.
265 : !> \param hessian Symmetric quadratic-model Hessian.
266 : !> \param linear Quadratic-model linear term.
267 : !> \param preferred History index favored when objective values are tied.
268 : !> \param coeff Best feasible coefficients found so far.
269 : !> \param best_value Objective value corresponding to coeff.
270 : ! **************************************************************************************************
271 138 : SUBROUTINE enumerate_active_faces(hessian, linear, preferred, coeff, best_value)
272 :
273 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hessian
274 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: linear
275 : INTEGER, INTENT(IN) :: preferred
276 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coeff
277 : REAL(KIND=dp), INTENT(INOUT) :: best_value
278 :
279 : INTEGER :: i, j, k, last_mask, mask, n
280 138 : INTEGER, ALLOCATABLE, DIMENSION(:) :: active
281 : LOGICAL :: solved
282 : REAL(KIND=dp) :: feasibility_tolerance, scale, value
283 138 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: candidate, rhs, solution
284 138 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: kkt
285 :
286 138 : n = SIZE(linear)
287 138 : last_mask = ISHFT(1, n) - 1
288 690 : ALLOCATE (active(n), candidate(n))
289 :
290 28704 : DO mask = 1, last_mask
291 28566 : k = POPCNT(mask)
292 28566 : IF (k < 2) CYCLE
293 :
294 : j = 0
295 303310 : DO i = 1, n
296 303310 : IF (BTEST(mask, i - 1)) THEN
297 139700 : j = j + 1
298 139700 : active(j) = i
299 : END IF
300 : END DO
301 :
302 195062 : ALLOCATE (kkt(k + 1, k + 1), rhs(k + 1), solution(k + 1))
303 27866 : kkt = 0.0_dp
304 27866 : rhs = 0.0_dp
305 :
306 167566 : DO j = 1, k
307 139700 : rhs(j) = -linear(active(j))
308 139700 : kkt(j, k + 1) = 1.0_dp
309 139700 : kkt(k + 1, j) = 1.0_dp
310 953506 : DO i = 1, k
311 925640 : kkt(i, j) = hessian(active(i), active(j))
312 : END DO
313 : END DO
314 27866 : rhs(k + 1) = 1.0_dp
315 :
316 27866 : CALL solve_dense_linear(kkt, rhs, solution, solved)
317 27866 : IF (solved) THEN
318 121192 : scale = MAX(1.0_dp, MAXVAL(ABS(solution(1:k))))
319 22042 : feasibility_tolerance = 2048.0_dp*EPSILON(1.0_dp)*REAL(n, KIND=dp)*scale
320 121192 : IF (MINVAL(solution(1:k)) >= -feasibility_tolerance) THEN
321 532 : candidate = 0.0_dp
322 1596 : DO i = 1, k
323 1596 : candidate(active(i)) = MAX(solution(i), 0.0_dp)
324 : END DO
325 4070 : IF (SUM(candidate) > 0.0_dp) THEN
326 7608 : candidate = candidate/SUM(candidate)
327 532 : value = quadratic_value(hessian, linear, candidate)
328 532 : CALL consider_candidate(candidate, value, preferred, coeff, best_value)
329 : END IF
330 : END IF
331 : END IF
332 :
333 56570 : DEALLOCATE (kkt, rhs, solution)
334 : END DO
335 :
336 138 : END SUBROUTINE enumerate_active_faces
337 :
338 : ! **************************************************************************************************
339 : !> \brief Simplex-preserving pairwise coordinate descent for larger histories.
340 : !> \param hessian Symmetric quadratic-model Hessian.
341 : !> \param linear Quadratic-model linear term.
342 : !> \param preferred History index favored when objective values are tied.
343 : !> \param coeff Best multistart coefficients found.
344 : !> \param best_value Objective value corresponding to coeff.
345 : !> \param status Local-convergence or iteration-limit status.
346 : ! **************************************************************************************************
347 2 : SUBROUTINE pairwise_minimize(hessian, linear, preferred, coeff, best_value, status)
348 :
349 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hessian
350 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: linear
351 : INTEGER, INTENT(IN) :: preferred
352 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coeff
353 : REAL(KIND=dp), INTENT(INOUT) :: best_value
354 : INTEGER, INTENT(OUT) :: status
355 :
356 : INTEGER :: n, start, trial_status
357 : LOGICAL :: select_trial
358 : REAL(KIND=dp) :: trial_value
359 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: trial
360 :
361 2 : n = SIZE(linear)
362 6 : ALLOCATE (trial(n))
363 :
364 : ! Starting from every vertex avoids simple multi-coordinate traps while retaining
365 : ! deterministic newest-first tie breaking. A uniform start also covers broad interiors.
366 2 : best_value = HUGE(0.0_dp)
367 2 : status = simplex_qp_iteration_limit
368 28 : DO start = n, 1, -1
369 26 : trial = 0.0_dp
370 26 : trial(start) = 1.0_dp
371 26 : trial_value = quadratic_value(hessian, linear, trial)
372 26 : CALL pairwise_descent(hessian, linear, trial, trial_value, trial_status)
373 :
374 26 : select_trial = strictly_better(trial_value, best_value)
375 26 : IF (.NOT. select_trial .AND. numerically_equal(trial_value, best_value)) THEN
376 24 : select_trial = prefer_newer(trial, coeff, preferred)
377 : END IF
378 26 : IF (select_trial) THEN
379 28 : coeff = trial
380 2 : best_value = trial_value
381 2 : IF (trial_status == simplex_qp_success) THEN
382 2 : status = simplex_qp_pairwise_stationary
383 : ELSE
384 0 : status = trial_status
385 : END IF
386 : END IF
387 : END DO
388 :
389 28 : trial = 1.0_dp/REAL(n, KIND=dp)
390 2 : trial_value = quadratic_value(hessian, linear, trial)
391 2 : CALL pairwise_descent(hessian, linear, trial, trial_value, trial_status)
392 2 : IF (strictly_better(trial_value, best_value)) THEN
393 0 : coeff = trial
394 0 : best_value = trial_value
395 0 : IF (trial_status == simplex_qp_success) THEN
396 0 : status = simplex_qp_pairwise_stationary
397 : ELSE
398 0 : status = trial_status
399 : END IF
400 : END IF
401 :
402 2 : END SUBROUTINE pairwise_minimize
403 :
404 : ! **************************************************************************************************
405 : !> \brief Run pairwise coordinate descent from one legal simplex point.
406 : !> \param hessian Symmetric quadratic-model Hessian.
407 : !> \param linear Quadratic-model linear term.
408 : !> \param coeff Initial and final legal simplex coefficients.
409 : !> \param value Initial and final objective value.
410 : !> \param status Convergence or iteration-limit status.
411 : ! **************************************************************************************************
412 28 : SUBROUTINE pairwise_descent(hessian, linear, coeff, value, status)
413 :
414 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hessian
415 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: linear
416 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coeff
417 : REAL(KIND=dp), INTENT(INOUT) :: value
418 : INTEGER, INTENT(OUT) :: status
419 :
420 : INTEGER :: best_from, best_to, from, iter, &
421 : max_iter, n, to
422 : REAL(KIND=dp) :: alpha, alpha_max, best_alpha, &
423 : best_delta, curvature, delta, &
424 : derivative, objective_tolerance
425 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gradient
426 :
427 28 : n = SIZE(linear)
428 28 : max_iter = MAX(500, 50*n*n)
429 84 : ALLOCATE (gradient(n))
430 :
431 28 : status = simplex_qp_iteration_limit
432 76 : DO iter = 1, max_iter
433 14896 : gradient(:) = MATMUL(hessian, coeff) + linear
434 : best_delta = 0.0_dp
435 : best_alpha = 0.0_dp
436 : best_from = 0
437 : best_to = 0
438 :
439 1064 : DO from = 1, n
440 988 : alpha_max = coeff(from)
441 988 : IF (alpha_max <= 64.0_dp*EPSILON(1.0_dp)) CYCLE
442 3380 : DO to = 1, n
443 3068 : IF (to == from) CYCLE
444 :
445 2832 : derivative = gradient(to) - gradient(from)
446 2832 : curvature = hessian(to, to) + hessian(from, from) - 2.0_dp*hessian(to, from)
447 :
448 2832 : IF (curvature > 0.0_dp) THEN
449 : ! Compare before dividing so a tiny positive curvature cannot overflow
450 : ! -derivative/curvature under CP2K's Debug FPE traps.
451 2782 : IF (derivative >= 0.0_dp) THEN
452 : alpha = 0.0_dp
453 894 : ELSE IF (-derivative >= alpha_max*curvature) THEN
454 : alpha = alpha_max
455 : ELSE
456 48 : alpha = -derivative/curvature
457 : END IF
458 : ELSE
459 : alpha = alpha_max
460 : END IF
461 :
462 2832 : delta = alpha*derivative + 0.5_dp*alpha*alpha*curvature
463 3820 : IF (delta < best_delta) THEN
464 308 : best_delta = delta
465 308 : best_alpha = alpha
466 308 : best_from = from
467 308 : best_to = to
468 : END IF
469 : END DO
470 : END DO
471 :
472 76 : objective_tolerance = 16.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, ABS(value))
473 76 : IF (best_from == 0 .OR. best_delta >= -objective_tolerance) THEN
474 28 : status = simplex_qp_success
475 28 : EXIT
476 : END IF
477 :
478 48 : coeff(best_from) = coeff(best_from) - best_alpha
479 48 : coeff(best_to) = coeff(best_to) + best_alpha
480 672 : coeff = MAX(coeff, 0.0_dp)
481 1296 : coeff = coeff/SUM(coeff)
482 48 : value = quadratic_value(hessian, linear, coeff)
483 : END DO
484 :
485 28 : END SUBROUTINE pairwise_descent
486 :
487 : ! **************************************************************************************************
488 : !> \brief Dense Gaussian elimination with partial pivoting for tiny KKT systems.
489 : !> \param matrix Square coefficient matrix.
490 : !> \param rhs Right-hand side vector.
491 : !> \param solution Computed solution, or zeros when the system is rejected.
492 : !> \param solved Whether a finite nonsingular solution was obtained.
493 : ! **************************************************************************************************
494 27866 : SUBROUTINE solve_dense_linear(matrix, rhs, solution, solved)
495 :
496 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
497 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rhs
498 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: solution
499 : LOGICAL, INTENT(OUT) :: solved
500 :
501 : INTEGER :: i, k, n, pivot
502 : REAL(KIND=dp) :: factor, pivot_tolerance, scale, temp
503 27866 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work_rhs, work_row
504 27866 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work_matrix
505 :
506 27866 : n = SIZE(rhs)
507 27866 : solved = .FALSE.
508 195432 : solution = 0.0_dp
509 27866 : IF (SIZE(matrix, 1) /= n .OR. SIZE(matrix, 2) /= n .OR. SIZE(solution) /= n) RETURN
510 :
511 195062 : ALLOCATE (work_matrix(n, n), work_rhs(n), work_row(n))
512 1288638 : work_matrix(:, :) = matrix
513 195432 : work_rhs(:) = rhs
514 1288638 : scale = MAX(1.0_dp, MAXVAL(ABS(work_matrix)))
515 27866 : pivot_tolerance = 1024.0_dp*EPSILON(1.0_dp)*REAL(n, KIND=dp)*scale
516 :
517 181624 : DO k = 1, n
518 779472 : pivot = k - 1 + MAXLOC(ABS(work_matrix(k:n, k)), DIM=1)
519 159582 : IF (ABS(work_matrix(pivot, k)) <= pivot_tolerance) RETURN
520 :
521 153758 : IF (pivot /= k) THEN
522 860628 : work_row(:) = work_matrix(k, :)
523 860628 : work_matrix(k, :) = work_matrix(pivot, :)
524 860628 : work_matrix(pivot, :) = work_row
525 113590 : temp = work_rhs(k)
526 113590 : work_rhs(k) = work_rhs(pivot)
527 113590 : work_rhs(pivot) = temp
528 : END IF
529 :
530 628124 : DO i = k + 1, n
531 452324 : factor = work_matrix(i, k)/work_matrix(k, k)
532 452324 : work_matrix(i, k) = 0.0_dp
533 : work_matrix(i, k + 1:n) = work_matrix(i, k + 1:n) - &
534 2466784 : factor*work_matrix(k, k + 1:n)
535 606082 : work_rhs(i) = work_rhs(i) - factor*work_rhs(k)
536 : END DO
537 : END DO
538 :
539 143234 : DO i = n, 1, -1
540 : solution(i) = (work_rhs(i) - DOT_PRODUCT(work_matrix(i, i + 1:n), solution(i + 1:n)))/ &
541 439120 : work_matrix(i, i)
542 : END DO
543 :
544 143234 : solved = ALL(ieee_is_finite(solution))
545 :
546 27866 : END SUBROUTINE solve_dense_linear
547 :
548 : ! **************************************************************************************************
549 : !> \brief Compare and, when appropriate, retain a finite simplex candidate.
550 : !> \param candidate Feasible simplex candidate.
551 : !> \param value Objective value at candidate.
552 : !> \param preferred History index favored when objective values are tied.
553 : !> \param coeff Best coefficients found so far.
554 : !> \param best_value Objective value corresponding to coeff.
555 : ! **************************************************************************************************
556 532 : SUBROUTINE consider_candidate(candidate, value, preferred, coeff, best_value)
557 :
558 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: candidate
559 : REAL(KIND=dp), INTENT(IN) :: value
560 : INTEGER, INTENT(IN) :: preferred
561 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coeff
562 : REAL(KIND=dp), INTENT(INOUT) :: best_value
563 :
564 532 : IF (.NOT. ieee_is_finite(value)) RETURN
565 532 : IF (strictly_better(value, best_value)) THEN
566 734 : coeff = candidate
567 116 : best_value = value
568 416 : ELSE IF (numerically_equal(value, best_value) .AND. prefer_newer(candidate, coeff, preferred)) THEN
569 0 : coeff = candidate
570 0 : best_value = value
571 : END IF
572 :
573 : END SUBROUTINE consider_candidate
574 :
575 : ! **************************************************************************************************
576 : !> \brief Evaluate the canonical quadratic objective.
577 : !> \param hessian Symmetric quadratic-model Hessian.
578 : !> \param linear Quadratic-model linear term.
579 : !> \param coeff Simplex coefficient vector.
580 : !> \return Canonical quadratic objective value.
581 : ! **************************************************************************************************
582 748 : PURE FUNCTION quadratic_value(hessian, linear, coeff) RESULT(value)
583 :
584 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hessian
585 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: linear, coeff
586 : REAL(KIND=dp) :: value
587 :
588 62204 : value = 0.5_dp*DOT_PRODUCT(coeff, MATMUL(hessian, coeff)) + DOT_PRODUCT(linear, coeff)
589 :
590 748 : END FUNCTION quadratic_value
591 :
592 : ! **************************************************************************************************
593 : !> \brief Scale-aware strict comparison used for deterministic candidate selection.
594 : !> \param value Candidate objective value.
595 : !> \param reference Current best objective value.
596 : !> \return Whether value improves reference beyond round-off tolerance.
597 : ! **************************************************************************************************
598 1146 : PURE FUNCTION strictly_better(value, reference) RESULT(better)
599 :
600 : REAL(KIND=dp), INTENT(IN) :: value, reference
601 : LOGICAL :: better
602 :
603 : better = value < reference - 512.0_dp*EPSILON(1.0_dp)* &
604 1146 : MAX(1.0_dp, ABS(value), ABS(reference))
605 :
606 1146 : END FUNCTION strictly_better
607 :
608 : ! **************************************************************************************************
609 : !> \brief Scale-aware equality comparison.
610 : !> \param value Candidate objective value.
611 : !> \param reference Current best objective value.
612 : !> \return Whether both values agree within round-off tolerance.
613 : ! **************************************************************************************************
614 440 : PURE FUNCTION numerically_equal(value, reference) RESULT(equal)
615 :
616 : REAL(KIND=dp), INTENT(IN) :: value, reference
617 : LOGICAL :: equal
618 :
619 : equal = ABS(value - reference) <= 512.0_dp*EPSILON(1.0_dp)* &
620 440 : MAX(1.0_dp, ABS(value), ABS(reference))
621 :
622 440 : END FUNCTION numerically_equal
623 :
624 : ! **************************************************************************************************
625 : !> \brief Deterministic tie break that first favors the caller's preferred history entry.
626 : !> \param candidate Candidate coefficient vector.
627 : !> \param reference Current coefficient vector.
628 : !> \param preferred History index favored first.
629 : !> \return Whether candidate wins the deterministic tie break.
630 : ! **************************************************************************************************
631 26 : PURE FUNCTION prefer_newer(candidate, reference, preferred) RESULT(prefer)
632 :
633 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: candidate, reference
634 : INTEGER, INTENT(IN) :: preferred
635 : LOGICAL :: prefer
636 :
637 : INTEGER :: i
638 : REAL(KIND=dp) :: tolerance
639 :
640 26 : prefer = .FALSE.
641 26 : tolerance = 512.0_dp*EPSILON(1.0_dp)
642 :
643 26 : IF (candidate(preferred) > reference(preferred) + tolerance) THEN
644 26 : prefer = .TRUE.
645 : RETURN
646 26 : ELSE IF (candidate(preferred) < reference(preferred) - tolerance) THEN
647 : RETURN
648 : END IF
649 :
650 308 : DO i = SIZE(candidate), 1, -1
651 286 : IF (i == preferred) CYCLE
652 286 : IF (candidate(i) > reference(i) + tolerance) THEN
653 0 : prefer = .TRUE.
654 : RETURN
655 264 : ELSE IF (candidate(i) < reference(i) - tolerance) THEN
656 : RETURN
657 : END IF
658 : END DO
659 :
660 : END FUNCTION prefer_newer
661 :
662 748 : END MODULE qs_scf_subspace_math
|
