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 2 : PROGRAM ot_covariant_preconditioner_unittest
8 2 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
9 : cp_blacs_env_release,&
10 : cp_blacs_env_type
11 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
12 : cp_fm_struct_release,&
13 : cp_fm_struct_type
14 : USE cp_fm_types, ONLY: cp_fm_create,&
15 : cp_fm_get_submatrix,&
16 : cp_fm_release,&
17 : cp_fm_set_submatrix,&
18 : cp_fm_type
19 : USE ieee_arithmetic, ONLY: ieee_is_finite
20 : USE input_constants, ONLY: ot_low_rank_base_overlap,&
21 : ot_precond_fermi_low_rank,&
22 : ot_precond_full_all_covariant
23 : USE iso_fortran_env, ONLY: error_unit
24 : USE kinds, ONLY: dp
25 : USE message_passing, ONLY: mp_para_env_type,&
26 : mp_world_finalize,&
27 : mp_world_init
28 : USE preconditioner_apply, ONLY: apply_preconditioner_fm
29 : USE preconditioner_types, ONLY: destroy_preconditioner,&
30 : init_preconditioner,&
31 : preconditioner_type
32 :
33 : IMPLICIT NONE
34 :
35 : INTEGER, PARAMETER :: k = 3, n = 5
36 : REAL(KIND=dp), PARAMETER :: eps = 2.0E-11_dp
37 :
38 : INTEGER :: i, j, rank_used
39 : REAL(KIND=dp), DIMENSION(n, k) :: gradient, gradient_rotated, hessian_x, &
40 : output, output_reference, &
41 : output_rotated, x_reference
42 : REAL(KIND=dp), DIMENSION(n) :: full_evals
43 : REAL(KIND=dp), DIMENSION(k) :: occ_evals
44 : REAL(KIND=dp), DIMENSION(k, k) :: occupied_h, q, rotation, rotation_2, &
45 : rotation_gauge
46 : REAL(KIND=dp), DIMENSION(n, n) :: inverse, overlap_inverse, spectral_vectors
47 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
48 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
49 : TYPE(cp_fm_type) :: matrix_in, matrix_out
50 : TYPE(mp_para_env_type), POINTER :: para_env
51 8 : TYPE(preconditioner_type) :: preconditioner_env
52 :
53 2 : NULLIFY (blacs_env, fm_struct, para_env)
54 2 : ALLOCATE (para_env)
55 2 : CALL mp_world_init(para_env)
56 2 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
57 2 : CALL init_preconditioner(preconditioner_env, para_env, blacs_env)
58 :
59 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=k, &
60 2 : context=blacs_env, para_env=para_env)
61 2 : CALL cp_fm_create(matrix_in, fm_struct, name="covariant test input")
62 2 : CALL cp_fm_create(matrix_out, fm_struct, name="covariant test output")
63 2 : CALL cp_fm_struct_release(fm_struct)
64 :
65 : CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=n, &
66 2 : context=blacs_env, para_env=para_env)
67 2 : ALLOCATE (preconditioner_env%fm)
68 2 : CALL cp_fm_create(preconditioner_env%fm, fm_struct, name="spectral eigenvectors")
69 2 : CALL cp_fm_struct_release(fm_struct)
70 :
71 : CALL cp_fm_struct_create(fm_struct, nrow_global=k, ncol_global=k, &
72 2 : context=blacs_env, para_env=para_env)
73 2 : ALLOCATE (preconditioner_env%occ_rotation)
74 2 : CALL cp_fm_create(preconditioner_env%occ_rotation, fm_struct, name="occupied rotation")
75 2 : CALL cp_fm_struct_release(fm_struct)
76 :
77 2 : preconditioner_env%in_use = ot_precond_full_all_covariant
78 2 : preconditioner_env%energy_gap = 0.10_dp
79 2 : ALLOCATE (preconditioner_env%full_evals(n), preconditioner_env%occ_evals(k))
80 2 : full_evals = [0.8_dp, 1.3_dp, 2.1_dp, 3.0_dp, 4.4_dp]
81 2 : occ_evals = [-0.9_dp, -0.2_dp, 0.35_dp]
82 12 : preconditioner_env%full_evals = full_evals
83 8 : preconditioner_env%occ_evals = occ_evals
84 :
85 2 : spectral_vectors = 0.0_dp
86 12 : DO i = 1, n
87 12 : spectral_vectors(i, i) = 1.0_dp
88 : END DO
89 2 : CALL cp_fm_set_submatrix(preconditioner_env%fm, spectral_vectors)
90 2 : CALL plane_rotation(0.43_dp, 1, 2, rotation)
91 2 : CALL plane_rotation(-0.31_dp, 2, 3, q)
92 80 : rotation = MATMUL(rotation, q)
93 2 : CALL cp_fm_set_submatrix(preconditioner_env%occ_rotation, rotation)
94 :
95 8 : DO j = 1, k
96 38 : DO i = 1, n
97 : gradient(i, j) = SIN(0.37_dp*REAL(2*i + j, dp)) + &
98 30 : 0.2_dp*COS(0.19_dp*REAL(i - 3*j, dp))
99 : x_reference(i, j) = COS(0.23_dp*REAL(i + 2*j, dp)) - &
100 36 : 0.1_dp*SIN(0.41_dp*REAL(3*i - j, dp))
101 : END DO
102 : END DO
103 :
104 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient, output_reference)
105 :
106 : ! Invert the frozen-H Sylvester operator A X - X B exactly.
107 2 : occupied_h = 0.0_dp
108 8 : DO j = 1, k
109 80 : occupied_h = occupied_h + occ_evals(j)*outer_product(rotation(:, j), rotation(:, j))
110 : END DO
111 152 : hessian_x = -MATMUL(x_reference, occupied_h)
112 12 : DO i = 1, n
113 42 : hessian_x(i, :) = hessian_x(i, :) + full_evals(i)*x_reference(i, :)
114 : END DO
115 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, hessian_x, output)
116 2 : CALL assert_close(output, x_reference, eps, "Frozen-H Sylvester inverse")
117 :
118 : ! Rotate the occupied gauge and verify P_(C R)(G R) = P_C(G) R.
119 2 : CALL plane_rotation(-0.37_dp, 1, 2, rotation_gauge)
120 2 : CALL plane_rotation(0.28_dp, 2, 3, q)
121 80 : rotation_gauge = MATMUL(rotation_gauge, q)
122 80 : rotation_2 = MATMUL(TRANSPOSE(rotation_gauge), rotation)
123 116 : gradient_rotated = MATMUL(gradient, rotation_gauge)
124 2 : CALL cp_fm_set_submatrix(preconditioner_env%occ_rotation, rotation_2)
125 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient_rotated, output_rotated)
126 : CALL assert_close(output_rotated, MATMUL(output_reference, rotation_gauge), eps, &
127 116 : "Occupied-gauge covariance")
128 :
129 : ! The gap floor must leave the inverse positive definite.
130 2 : preconditioner_env%energy_gap = 0.75_dp
131 8 : preconditioner_env%occ_evals = [0.7_dp, 0.9_dp, 1.1_dp]
132 2 : CALL cp_fm_set_submatrix(preconditioner_env%occ_rotation, rotation)
133 2 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient, output)
134 38 : IF (SUM(gradient*output) <= 0.0_dp) THEN
135 0 : ERROR STOP "Gap-floor inverse is not positive definite"
136 : END IF
137 :
138 : ! Exercise the production low-rank FM path with a nontrivial overlap base,
139 : ! including zero rank, truncation, and a gap-clipped complementary eigenvalue.
140 2 : preconditioner_env%in_use = ot_precond_fermi_low_rank
141 2 : preconditioner_env%low_rank_base = ot_low_rank_base_overlap
142 2 : preconditioner_env%energy_gap = 0.1_dp
143 2 : preconditioner_env%spectral_reference = full_evals(k + 1) - 0.04_dp
144 2 : preconditioner_env%spectral_base_scale = 1.25_dp
145 2 : CALL plane_rotation(0.29_dp, 2, 4, spectral_vectors)
146 12 : spectral_vectors(:, 2) = 1.2_dp*spectral_vectors(:, 2)
147 12 : spectral_vectors(:, 4) = 0.8_dp*spectral_vectors(:, 4)
148 312 : overlap_inverse = MATMUL(spectral_vectors, TRANSPOSE(spectral_vectors))
149 2 : CALL cp_fm_set_submatrix(preconditioner_env%fm, spectral_vectors)
150 2 : ALLOCATE (preconditioner_env%base_fm)
151 2 : CALL cp_fm_create(preconditioner_env%base_fm, preconditioner_env%fm%matrix_struct)
152 2 : CALL cp_fm_set_submatrix(preconditioner_env%base_fm, overlap_inverse)
153 8 : DO rank_used = 0, n - k
154 6 : preconditioner_env%spectral_rank = rank_used
155 186 : inverse = preconditioner_env%spectral_base_scale*overlap_inverse
156 12 : DO i = k + 1, k + rank_used
157 : inverse = inverse + (1.0_dp/MAX(preconditioner_env%energy_gap, &
158 : full_evals(i) - preconditioner_env%spectral_reference) - &
159 : preconditioner_env%spectral_base_scale)* &
160 192 : outer_product(spectral_vectors(:, i), spectral_vectors(:, i))
161 : END DO
162 564 : output_reference = MATMUL(inverse, gradient)
163 6 : CALL apply_fm(preconditioner_env, matrix_in, matrix_out, gradient, output)
164 8 : CALL assert_close(output, output_reference, eps, "Low-rank spectral inverse")
165 : END DO
166 :
167 2 : CALL cp_fm_release(matrix_in)
168 2 : CALL cp_fm_release(matrix_out)
169 2 : CALL destroy_preconditioner(preconditioner_env)
170 2 : CALL cp_blacs_env_release(blacs_env)
171 2 : CALL mp_world_finalize()
172 4 : DEALLOCATE (para_env)
173 :
174 : CONTAINS
175 :
176 : ! **************************************************************************************************
177 : !> \brief Apply the production FM path to a replicated test matrix.
178 : !> \param preconditioner_env ...
179 : !> \param matrix_in ...
180 : !> \param matrix_out ...
181 : !> \param input ...
182 : !> \param RESULT ...
183 : ! **************************************************************************************************
184 14 : SUBROUTINE apply_fm(preconditioner_env, matrix_in, matrix_out, input, RESULT)
185 :
186 : TYPE(preconditioner_type) :: preconditioner_env
187 : TYPE(cp_fm_type), INTENT(IN) :: matrix_in, matrix_out
188 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: input
189 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: result
190 :
191 14 : CALL cp_fm_set_submatrix(matrix_in, input)
192 14 : CALL apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
193 14 : CALL cp_fm_get_submatrix(matrix_out, RESULT)
194 :
195 2 : END SUBROUTINE apply_fm
196 :
197 : ! **************************************************************************************************
198 : !> \brief Create an orthogonal plane rotation.
199 : !> \param angle ...
200 : !> \param axis_1 ...
201 : !> \param axis_2 ...
202 : !> \param rotation ...
203 : ! **************************************************************************************************
204 10 : SUBROUTINE plane_rotation(angle, axis_1, axis_2, rotation)
205 :
206 : REAL(KIND=dp), INTENT(IN) :: angle
207 : INTEGER, INTENT(IN) :: axis_1, axis_2
208 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: rotation
209 :
210 : INTEGER :: i
211 :
212 166 : rotation = 0.0_dp
213 44 : DO i = 1, SIZE(rotation, 1)
214 44 : rotation(i, i) = 1.0_dp
215 : END DO
216 10 : rotation(axis_1, axis_1) = COS(angle)
217 10 : rotation(axis_2, axis_2) = COS(angle)
218 10 : rotation(axis_1, axis_2) = -SIN(angle)
219 10 : rotation(axis_2, axis_1) = SIN(angle)
220 :
221 10 : END SUBROUTINE plane_rotation
222 :
223 : ! **************************************************************************************************
224 : !> \brief Form a real outer product.
225 : !> \param left ...
226 : !> \param right ...
227 : !> \return ...
228 : ! **************************************************************************************************
229 12 : PURE FUNCTION outer_product(left, right) RESULT(product)
230 :
231 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: left, right
232 : REAL(KIND=dp), DIMENSION(SIZE(left), SIZE(right)) :: product
233 :
234 : INTEGER :: i, j
235 :
236 60 : DO j = 1, SIZE(right)
237 264 : DO i = 1, SIZE(left)
238 252 : PRODUCT(i, j) = left(i)*right(j)
239 : END DO
240 : END DO
241 :
242 12 : END FUNCTION outer_product
243 :
244 : ! **************************************************************************************************
245 : !> \brief Abort when two matrices differ beyond the requested tolerance.
246 : !> \param actual ...
247 : !> \param reference ...
248 : !> \param tolerance ...
249 : !> \param label ...
250 : ! **************************************************************************************************
251 10 : SUBROUTINE assert_close(actual, reference, tolerance, label)
252 :
253 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: actual, reference
254 : REAL(KIND=dp), INTENT(IN) :: tolerance
255 : CHARACTER(LEN=*), INTENT(IN) :: label
256 :
257 : INTEGER :: io_unit
258 :
259 : IF (.NOT. ALL(ieee_is_finite(actual)) .OR. &
260 570 : .NOT. ALL(ieee_is_finite(reference)) .OR. MAXVAL(ABS(actual - reference)) > tolerance) THEN
261 0 : io_unit = error_unit
262 0 : WRITE (io_unit, '(A,2(1X,ES13.6))') TRIM(label)//": error/tolerance", &
263 0 : MAXVAL(ABS(actual - reference)), tolerance
264 0 : ERROR STOP "Matrix comparison failed"
265 : END IF
266 :
267 10 : END SUBROUTINE assert_close
268 :
269 : END PROGRAM ot_covariant_preconditioner_unittest
|
