LCOV - code coverage report
Current view: top level - src - qs_condnum.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 123 131 93.9 %
Date: 2025-02-18 08:24:35 Functions: 4 4 100.0 %

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

Generated by: LCOV version 1.15