26 USE dbcsr_api,
ONLY: &
27 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_gershgorin_norm, &
28 dbcsr_get_block_diag, dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_iterator_blocks_left, &
29 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
30 dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
33 #include "./base/base_uses.f90"
41 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_condnum'
63 SUBROUTINE overlap_condnum(matrixkp_s, condnum, iunit, norml1, norml2, use_arnoldi, blacs_env)
65 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrixkp_s
66 REAL(kind=
dp),
DIMENSION(2),
INTENT(INOUT) :: condnum
67 INTEGER,
INTENT(IN) :: iunit
68 LOGICAL,
INTENT(IN) :: norml1, norml2, use_arnoldi
69 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
71 CHARACTER(len=*),
PARAMETER :: routinen =
'overlap_condnum'
73 INTEGER :: handle, ic, maxiter, nbas, ndep
75 REAL(kind=
dp) :: amnorm, anorm, eps_ev, max_ev, min_ev, &
77 REAL(kind=
dp),
DIMENSION(2) :: eigvals
78 TYPE(cp_fm_struct_type),
POINTER :: matrix_struct
79 TYPE(cp_fm_type) :: fmsmat, fmwork
80 TYPE(dbcsr_type) :: tempmat
81 TYPE(dbcsr_type),
POINTER :: smat
83 CALL timeset(routinen, handle)
87 IF (
SIZE(matrixkp_s, 2) == 1)
THEN
88 IF (iunit > 0)
WRITE (iunit,
'(/,T2,A)')
"OVERLAP MATRIX CONDITION NUMBER"
89 smat => matrixkp_s(1, 1)%matrix
91 IF (iunit > 0)
WRITE (iunit,
'(/,T2,A)')
"OVERLAP MATRIX CONDITION NUMBER AT GAMMA POINT"
93 CALL dbcsr_create(smat, template=matrixkp_s(1, 1)%matrix)
94 CALL dbcsr_copy(smat, matrixkp_s(1, 1)%matrix)
95 DO ic = 2,
SIZE(matrixkp_s, 2)
96 CALL dbcsr_add(smat, matrixkp_s(1, ic)%matrix, 1.0_dp, 1.0_dp)
100 IF (
ASSOCIATED(smat))
THEN
101 cpassert(dbcsr_get_matrix_type(smat) .EQ. dbcsr_type_symmetric)
104 anorm = dbcsr_gershgorin_norm(smat)
105 CALL estimate_norm_invmat(smat, amnorm)
107 WRITE (iunit,
'(T2,A)')
"1-Norm Condition Number (Estimate)"
108 WRITE (iunit,
'(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
109 "CN : |A|*|A^-1|: ", anorm,
" * ", amnorm,
"=", anorm*amnorm,
"Log(1-CN):", log10(anorm*amnorm)
111 condnum(1) = anorm*amnorm
116 CALL dbcsr_get_info(smat, nfullrows_total=nbas)
118 nrow_global=nbas, ncol_global=nbas)
122 CALL dbcsr_create(tempmat, template=smat, matrix_type=dbcsr_type_no_symmetry)
123 CALL dbcsr_desymmetrize(smat, tempmat)
128 CALL cp_fm_power(fmsmat, fmwork, -1.0_dp, eps_ev, ndep, eigvals=eigvals)
133 CALL dbcsr_release(tempmat)
134 CALL cp_fm_release(fmsmat)
135 CALL cp_fm_release(fmwork)
139 WRITE (iunit,
'(T2,A)')
"1-Norm and 2-Norm Condition Numbers using Diagonalization"
141 WRITE (iunit,
'(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
142 "CN : |A|*|A^-1|: ", anorm,
" * ", amnorm,
"=", anorm*amnorm,
"Log(1-CN):", log10(anorm*amnorm)
143 WRITE (iunit,
'(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
144 "CN : max/min ev: ", max_ev,
" / ", min_ev,
"=", max_ev/min_ev,
"Log(2-CN):", log10(max_ev/min_ev)
146 WRITE (iunit,
'(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
147 "CN : max/min EV: ", max_ev,
" / ", min_ev,
"Log(CN): infinity"
151 condnum(1) = anorm*amnorm
152 condnum(2) = max_ev/min_ev
154 condnum(1:2) = 0.0_dp
157 IF (use_arnoldi)
THEN
159 threshold = 1.0e-6_dp
163 threshold=threshold, max_iter=maxiter, converged=converged)
165 WRITE (iunit,
'(T2,A)')
"2-Norm Condition Number using Arnoldi iterations"
167 WRITE (iunit,
'(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
168 "CN : max/min ev: ", max_ev,
" / ", min_ev,
"=", max_ev/min_ev,
"Log(2-CN):", log10(max_ev/min_ev)
170 WRITE (iunit,
'(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
171 "CN : max/min ev: ", max_ev,
" / ", min_ev,
"Log(CN): infinity"
175 condnum(2) = max_ev/min_ev
180 IF (min_ev == 0)
THEN
181 cpwarn(
"Ill-conditioned S matrix: basis set is overcomplete.")
182 ELSE IF ((max_ev/min_ev) > eps_ev)
THEN
183 cpwarn(
"Ill-conditioned S matrix: basis set is overcomplete.")
186 cpwarn(
"Condition number estimate of overlap matrix is not reliable (not converged).")
190 IF (
SIZE(matrixkp_s, 2) == 1)
THEN
193 CALL dbcsr_release(smat)
197 CALL timestop(handle)
212 SUBROUTINE estimate_norm_invmat(amat, anorm)
213 TYPE(dbcsr_type),
POINTER :: amat
214 REAL(kind=
dp),
INTENT(OUT) :: anorm
216 INTEGER :: i, k, nbas
217 INTEGER,
DIMENSION(1) :: r
218 REAL(kind=
dp) :: g, gg
219 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: x, xsi
220 REAL(kind=
dp),
DIMENSION(2) :: work
221 REAL(kind=
dp),
EXTERNAL :: dlange
222 TYPE(dbcsr_type) :: pmat
225 CALL dbcsr_create(pmat, name=
"SMAT Preconditioner", template=amat)
227 CALL dbcsr_get_block_diag(amat, pmat)
229 CALL smat_precon_diag(pmat)
232 CALL dbcsr_get_info(amat, nfullrows_total=nbas)
233 ALLOCATE (x(nbas), xsi(nbas))
234 x(1:nbas) = 1._dp/real(nbas, kind=
dp)
235 CALL dbcsr_solve(amat, x, pmat)
236 g = dlange(
"1", nbas, 1, x, nbas, work)
237 xsi(1:nbas) = sign(1._dp, x(1:nbas))
238 x(1:nbas) = xsi(1:nbas)
239 CALL dbcsr_solve(amat, x, pmat)
245 CALL dbcsr_solve(amat, x, pmat)
247 g = dlange(
"1", nbas, 1, x, nbas, work)
249 x(1:nbas) = sign(1._dp, x(1:nbas))
250 IF (sum(abs(x - xsi)) == 0 .OR. sum(abs(x + xsi)) == 0)
EXIT
251 xsi(1:nbas) = x(1:nbas)
252 CALL dbcsr_solve(amat, x, pmat)
255 IF (sum(r) == sum(maxloc(abs(x))))
EXIT
260 x(i) = -1._dp**(i + 1)*(1._dp + real(i - 1,
dp)/real(nbas - 1,
dp))
265 CALL dbcsr_solve(amat, x, pmat)
266 gg = dlange(
"1", nbas, 1, x, nbas, work)
267 gg = 2._dp*gg/real(3*nbas,
dp)
270 CALL dbcsr_release(pmat)
272 END SUBROUTINE estimate_norm_invmat
280 SUBROUTINE dbcsr_solve(amat, x, pmat)
281 TYPE(dbcsr_type),
POINTER :: amat
282 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: x
283 TYPE(dbcsr_type) :: pmat
285 INTEGER :: max_iter, nbas
287 REAL(kind=
dp) :: threshold
289 CALL dbcsr_get_info(amat, nfullrows_total=nbas)
290 max_iter = min(1000, nbas)
294 END SUBROUTINE dbcsr_solve
300 SUBROUTINE smat_precon_diag(pmat)
301 TYPE(dbcsr_type) :: pmat
303 INTEGER :: iatom, info, jatom, n
304 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: sblock
305 TYPE(dbcsr_iterator_type) :: dbcsr_iter
307 CALL dbcsr_iterator_start(dbcsr_iter, pmat)
308 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
309 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, sblock)
310 cpassert(iatom == jatom)
312 CALL invmat(sblock(1:n, 1:n), info)
315 CALL dbcsr_iterator_stop(dbcsr_iter)
317 END SUBROUTINE smat_precon_diag
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...
subroutine, public arnoldi_conjugate_gradient(matrix_a, vec_x, matrix_p, converged, threshold, max_iter)
Wrapper for conjugated gradient algorithm for Ax=b.
methods related to the blacs parallel environment
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
basic linear algebra operations for full matrices
real(kind=dp) function, public cp_fm_norm(matrix, mode)
norm of matrix using (p)dlange
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Defines the basic variable types.
integer, parameter, public dp
Collection of simple mathematical functions and subroutines.
subroutine, public invmat(a, info)
returns inverse of matrix using the lapack routines DGETRF and DGETRI
Calculation of overlap matrix condition numbers.
subroutine, public overlap_condnum(matrixkp_s, condnum, iunit, norml1, norml2, use_arnoldi, blacs_env)
Calculation of the overlap matrix Condition Number.