LCOV - code coverage report
Current view: top level - src - qs_condnum.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 93.9 % 131 123
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 4 4

            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         1740 :          xsi(1:nbas) = x(1:nbas)
     254           36 :          CALL dbcsr_solve(amat, x, pmat)
     255           36 :          k = k + 1
     256           36 :          IF (k > 5) EXIT
     257         1780 :          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          196 :    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          196 :       CALL dbcsr_get_info(amat, nfullrows_total=nbas)
     292          196 :       max_iter = MIN(1000, nbas)
     293          196 :       threshold = 1.e-6_dp
     294          196 :       CALL arnoldi_conjugate_gradient(amat, x, pmat, converged=converged, threshold=threshold, max_iter=max_iter)
     295              : 
     296          196 :    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              : 
        

Generated by: LCOV version 2.0-1