Line data Source code
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 Calculation of overlap matrix condition numbers
10 : !> \par History
11 : !> \author JGH (14.11.2016)
12 : ! **************************************************************************************************
13 : MODULE qs_condnum
14 : USE arnoldi_api, ONLY: arnoldi_conjugate_gradient,&
15 : arnoldi_extremal
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_dbcsr_api, ONLY: &
18 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_get_info, &
19 : dbcsr_get_matrix_type, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
20 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
21 : dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
22 : USE cp_dbcsr_contrib, ONLY: dbcsr_gershgorin_norm,&
23 : dbcsr_get_block_diag
24 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
25 : USE cp_fm_basic_linalg, ONLY: cp_fm_norm
26 : USE cp_fm_diag, ONLY: cp_fm_power
27 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
28 : cp_fm_struct_release,&
29 : cp_fm_struct_type
30 : USE cp_fm_types, ONLY: cp_fm_create,&
31 : cp_fm_release,&
32 : cp_fm_type
33 : USE kinds, ONLY: dp
34 : USE mathlib, ONLY: invmat
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : ! *** Global parameters ***
42 :
43 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_condnum'
44 :
45 : ! *** Public subroutines ***
46 :
47 : PUBLIC :: overlap_condnum
48 :
49 : CONTAINS
50 :
51 : ! **************************************************************************************************
52 : !> \brief Calculation of the overlap matrix Condition Number
53 : !> \param matrixkp_s The overlap matrices to be calculated (kpoints, optional)
54 : !> \param condnum Condition numbers for 1 and 2 norm
55 : !> \param iunit output unit
56 : !> \param norml1 logical: calculate estimate to 1-norm
57 : !> \param norml2 logical: calculate estimate to 1-norm and 2-norm condition number
58 : !> \param use_arnoldi logical: use Arnoldi iteration to estimate 2-norm condition number
59 : !> \param blacs_env ...
60 : !> \date 07.11.2016
61 : !> \par History
62 : !> \author JHU
63 : !> \version 1.0
64 : ! **************************************************************************************************
65 274 : SUBROUTINE overlap_condnum(matrixkp_s, condnum, iunit, norml1, norml2, use_arnoldi, blacs_env)
66 :
67 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_s
68 : REAL(KIND=dp), DIMENSION(2), INTENT(INOUT) :: condnum
69 : INTEGER, INTENT(IN) :: iunit
70 : LOGICAL, INTENT(IN) :: norml1, norml2, use_arnoldi
71 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
72 :
73 : CHARACTER(len=*), PARAMETER :: routineN = 'overlap_condnum'
74 :
75 : INTEGER :: handle, ic, maxiter, nbas, ndep
76 : LOGICAL :: converged
77 : REAL(KIND=dp) :: amnorm, anorm, eps_ev, max_ev, min_ev, &
78 : threshold
79 : REAL(KIND=dp), DIMENSION(2) :: eigvals
80 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
81 : TYPE(cp_fm_type) :: fmsmat, fmwork
82 : TYPE(dbcsr_type) :: tempmat
83 : TYPE(dbcsr_type), POINTER :: smat
84 :
85 274 : CALL timeset(routineN, handle)
86 :
87 274 : condnum(1:2) = 0.0_dp
88 274 : NULLIFY (smat)
89 274 : IF (SIZE(matrixkp_s, 2) == 1) THEN
90 268 : IF (iunit > 0) WRITE (iunit, '(/,T2,A)') "OVERLAP MATRIX CONDITION NUMBER"
91 268 : smat => matrixkp_s(1, 1)%matrix
92 : ELSE
93 6 : IF (iunit > 0) WRITE (iunit, '(/,T2,A)') "OVERLAP MATRIX CONDITION NUMBER AT GAMMA POINT"
94 6 : ALLOCATE (smat)
95 6 : CALL dbcsr_create(smat, template=matrixkp_s(1, 1)%matrix)
96 6 : CALL dbcsr_copy(smat, matrixkp_s(1, 1)%matrix)
97 2050 : DO ic = 2, SIZE(matrixkp_s, 2)
98 2050 : CALL dbcsr_add(smat, matrixkp_s(1, ic)%matrix, 1.0_dp, 1.0_dp)
99 : END DO
100 : END IF
101 : !
102 274 : IF (ASSOCIATED(smat)) THEN
103 274 : CPASSERT(dbcsr_get_matrix_type(smat) .EQ. dbcsr_type_symmetric)
104 274 : IF (norml1) THEN
105 : ! norm of S
106 40 : anorm = dbcsr_gershgorin_norm(smat)
107 40 : CALL estimate_norm_invmat(smat, amnorm)
108 40 : IF (iunit > 0) THEN
109 20 : WRITE (iunit, '(T2,A)') "1-Norm Condition Number (Estimate)"
110 : WRITE (iunit, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
111 20 : "CN : |A|*|A^-1|: ", anorm, " * ", amnorm, "=", anorm*amnorm, "Log(1-CN):", LOG10(anorm*amnorm)
112 : END IF
113 40 : condnum(1) = anorm*amnorm
114 : END IF
115 274 : IF (norml2) THEN
116 246 : eps_ev = 1.0E-14_dp
117 : ! diagonalization
118 246 : CALL dbcsr_get_info(smat, nfullrows_total=nbas)
119 : CALL cp_fm_struct_create(fmstruct=matrix_struct, context=blacs_env, &
120 246 : nrow_global=nbas, ncol_global=nbas)
121 246 : CALL cp_fm_create(fmsmat, matrix_struct)
122 246 : CALL cp_fm_create(fmwork, matrix_struct)
123 : ! transfer to FM
124 246 : CALL dbcsr_create(tempmat, template=smat, matrix_type=dbcsr_type_no_symmetry)
125 246 : CALL dbcsr_desymmetrize(smat, tempmat)
126 246 : CALL copy_dbcsr_to_fm(tempmat, fmsmat)
127 :
128 : ! diagonalize
129 246 : anorm = cp_fm_norm(fmsmat, "1")
130 246 : CALL cp_fm_power(fmsmat, fmwork, -1.0_dp, eps_ev, ndep, eigvals=eigvals)
131 246 : min_ev = eigvals(1)
132 246 : max_ev = eigvals(2)
133 246 : amnorm = cp_fm_norm(fmsmat, "1")
134 :
135 246 : CALL dbcsr_release(tempmat)
136 246 : CALL cp_fm_release(fmsmat)
137 246 : CALL cp_fm_release(fmwork)
138 246 : CALL cp_fm_struct_release(matrix_struct)
139 :
140 246 : IF (iunit > 0) THEN
141 6 : WRITE (iunit, '(T2,A)') "1-Norm and 2-Norm Condition Numbers using Diagonalization"
142 6 : IF (min_ev > 0) THEN
143 : WRITE (iunit, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
144 6 : "CN : |A|*|A^-1|: ", anorm, " * ", amnorm, "=", anorm*amnorm, "Log(1-CN):", LOG10(anorm*amnorm)
145 : WRITE (iunit, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
146 6 : "CN : max/min ev: ", max_ev, " / ", min_ev, "=", max_ev/min_ev, "Log(2-CN):", LOG10(max_ev/min_ev)
147 : ELSE
148 : WRITE (iunit, '(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
149 0 : "CN : max/min EV: ", max_ev, " / ", min_ev, "Log(CN): infinity"
150 : END IF
151 : END IF
152 246 : IF (min_ev > 0) THEN
153 246 : condnum(1) = anorm*amnorm
154 246 : condnum(2) = max_ev/min_ev
155 : ELSE
156 0 : condnum(1:2) = 0.0_dp
157 : END IF
158 : END IF
159 274 : IF (use_arnoldi) THEN
160 : ! parameters for matrix condition test
161 12 : threshold = 1.0E-6_dp
162 12 : maxiter = 1000
163 12 : eps_ev = 1.0E8_dp
164 : CALL arnoldi_extremal(smat, max_ev, min_ev, &
165 12 : threshold=threshold, max_iter=maxiter, converged=converged)
166 12 : IF (iunit > 0) THEN
167 6 : WRITE (iunit, '(T2,A)') "2-Norm Condition Number using Arnoldi iterations"
168 6 : IF (min_ev > 0) THEN
169 : WRITE (iunit, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
170 6 : "CN : max/min ev: ", max_ev, " / ", min_ev, "=", max_ev/min_ev, "Log(2-CN):", LOG10(max_ev/min_ev)
171 : ELSE
172 : WRITE (iunit, '(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
173 0 : "CN : max/min ev: ", max_ev, " / ", min_ev, "Log(CN): infinity"
174 : END IF
175 : END IF
176 12 : IF (min_ev > 0) THEN
177 12 : condnum(2) = max_ev/min_ev
178 : ELSE
179 0 : condnum(2) = 0.0_dp
180 : END IF
181 12 : IF (converged) THEN
182 12 : IF (min_ev == 0) THEN
183 0 : CPWARN("Ill-conditioned S matrix: basis set is overcomplete.")
184 12 : ELSE IF ((max_ev/min_ev) > eps_ev) THEN
185 0 : CPWARN("Ill-conditioned S matrix: basis set is overcomplete.")
186 : END IF
187 : ELSE
188 0 : CPWARN("Condition number estimate of overlap matrix is not reliable (not converged).")
189 : END IF
190 : END IF
191 : END IF
192 274 : IF (SIZE(matrixkp_s, 2) == 1) THEN
193 : NULLIFY (smat)
194 : ELSE
195 6 : CALL dbcsr_release(smat)
196 6 : DEALLOCATE (smat)
197 : END IF
198 :
199 274 : CALL timestop(handle)
200 :
201 274 : END SUBROUTINE overlap_condnum
202 :
203 : ! **************************************************************************************************
204 : !> \brief Calculates an estimate of the 1-norm of the inverse of a matrix
205 : !> Uses LAPACK norm estimator algorithm
206 : !> NJ Higham, Function of Matrices, Algorithm 3.21, page 66
207 : !> \param amat Sparse, symmetric matrix
208 : !> \param anorm Estimate of ||INV(A)||
209 : !> \date 15.11.2016
210 : !> \par History
211 : !> \author JHU
212 : !> \version 1.0
213 : ! **************************************************************************************************
214 40 : SUBROUTINE estimate_norm_invmat(amat, anorm)
215 : TYPE(dbcsr_type), POINTER :: amat
216 : REAL(KIND=dp), INTENT(OUT) :: anorm
217 :
218 : INTEGER :: i, k, nbas
219 : INTEGER, DIMENSION(1) :: r
220 : REAL(KIND=dp) :: g, gg
221 40 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x, xsi
222 : REAL(KIND=dp), DIMENSION(2) :: work
223 : REAL(KIND=dp), EXTERNAL :: dlange
224 : TYPE(dbcsr_type) :: pmat
225 :
226 : ! generate a block diagonal preconditioner
227 40 : CALL dbcsr_create(pmat, name="SMAT Preconditioner", template=amat)
228 : ! replicate the diagonal blocks of the overlap matrix
229 40 : CALL dbcsr_get_block_diag(amat, pmat)
230 : ! invert preconditioner
231 40 : CALL smat_precon_diag(pmat)
232 :
233 40 : anorm = 1.0_dp
234 40 : CALL dbcsr_get_info(amat, nfullrows_total=nbas)
235 160 : ALLOCATE (x(nbas), xsi(nbas))
236 1796 : x(1:nbas) = 1._dp/REAL(nbas, KIND=dp)
237 40 : CALL dbcsr_solve(amat, x, pmat)
238 40 : g = dlange("1", nbas, 1, x, nbas, work)
239 1796 : xsi(1:nbas) = SIGN(1._dp, x(1:nbas))
240 1796 : x(1:nbas) = xsi(1:nbas)
241 40 : CALL dbcsr_solve(amat, x, pmat)
242 40 : k = 2
243 : DO
244 1836 : r = MAXLOC(ABS(x))
245 1796 : x(1:nbas) = 0._dp
246 80 : x(r) = 1._dp
247 40 : CALL dbcsr_solve(amat, x, pmat)
248 40 : gg = g
249 40 : g = dlange("1", nbas, 1, x, nbas, work)
250 40 : IF (g <= gg) EXIT
251 1796 : x(1:nbas) = SIGN(1._dp, x(1:nbas))
252 3552 : IF (SUM(ABS(x - xsi)) == 0 .OR. SUM(ABS(x + xsi)) == 0) EXIT
253 1768 : xsi(1:nbas) = x(1:nbas)
254 38 : CALL dbcsr_solve(amat, x, pmat)
255 38 : k = k + 1
256 38 : IF (k > 5) EXIT
257 1808 : IF (SUM(r) == SUM(MAXLOC(ABS(x)))) EXIT
258 : END DO
259 : !
260 40 : IF (nbas > 1) THEN
261 1796 : DO i = 1, nbas
262 1796 : x(i) = -1._dp**(i + 1)*(1._dp + REAL(i - 1, dp)/REAL(nbas - 1, dp))
263 : END DO
264 : ELSE
265 0 : x(1) = 1.0_dp
266 : END IF
267 40 : CALL dbcsr_solve(amat, x, pmat)
268 40 : gg = dlange("1", nbas, 1, x, nbas, work)
269 40 : gg = 2._dp*gg/REAL(3*nbas, dp)
270 40 : anorm = MAX(g, gg)
271 40 : DEALLOCATE (x, xsi)
272 40 : CALL dbcsr_release(pmat)
273 :
274 80 : END SUBROUTINE estimate_norm_invmat
275 :
276 : ! **************************************************************************************************
277 : !> \brief ...
278 : !> \param amat ...
279 : !> \param x ...
280 : !> \param pmat ...
281 : ! **************************************************************************************************
282 198 : SUBROUTINE dbcsr_solve(amat, x, pmat)
283 : TYPE(dbcsr_type), POINTER :: amat
284 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: x
285 : TYPE(dbcsr_type) :: pmat
286 :
287 : INTEGER :: max_iter, nbas
288 : LOGICAL :: converged
289 : REAL(KIND=dp) :: threshold
290 :
291 198 : CALL dbcsr_get_info(amat, nfullrows_total=nbas)
292 198 : max_iter = MIN(1000, nbas)
293 198 : threshold = 1.e-6_dp
294 198 : CALL arnoldi_conjugate_gradient(amat, x, pmat, converged=converged, threshold=threshold, max_iter=max_iter)
295 :
296 198 : END SUBROUTINE dbcsr_solve
297 :
298 : ! **************************************************************************************************
299 : !> \brief ...
300 : !> \param pmat ...
301 : ! **************************************************************************************************
302 40 : SUBROUTINE smat_precon_diag(pmat)
303 : TYPE(dbcsr_type) :: pmat
304 :
305 : INTEGER :: iatom, info, jatom, n
306 40 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sblock
307 : TYPE(dbcsr_iterator_type) :: dbcsr_iter
308 :
309 40 : CALL dbcsr_iterator_start(dbcsr_iter, pmat)
310 122 : DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
311 82 : CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, sblock)
312 82 : CPASSERT(iatom == jatom)
313 82 : n = SIZE(sblock, 1)
314 82 : CALL invmat(sblock(1:n, 1:n), info)
315 122 : CPASSERT(info == 0)
316 : END DO
317 40 : CALL dbcsr_iterator_stop(dbcsr_iter)
318 :
319 40 : END SUBROUTINE smat_precon_diag
320 :
321 : END MODULE qs_condnum
322 :
|
