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