Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief computes preconditioners, and implements methods to apply them
      10             : !>      currently used in qs_ot
      11             : !> \par History
      12             : !>      - [UB] 2009-05-13 Adding stable approximate inverse (full and sparse)
      13             : !> \author Joost VandeVondele (09.2002)
      14             : ! **************************************************************************************************
      15             : MODULE preconditioner_apply
      16             :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_restore
      17             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      18             :                                               cp_fm_get_info,&
      19             :                                               cp_fm_release,&
      20             :                                               cp_fm_type
      21             :    USE dbcsr_api,                       ONLY: &
      22             :         dbcsr_copy, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
      23             :         dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_release, dbcsr_type
      24             :    USE input_constants,                 ONLY: ot_precond_full_all,&
      25             :                                               ot_precond_full_kinetic,&
      26             :                                               ot_precond_full_single,&
      27             :                                               ot_precond_full_single_inverse,&
      28             :                                               ot_precond_s_inverse,&
      29             :                                               ot_precond_solver_direct,&
      30             :                                               ot_precond_solver_inv_chol,&
      31             :                                               ot_precond_solver_update
      32             :    USE kinds,                           ONLY: dp
      33             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      34             :    USE preconditioner_types,            ONLY: preconditioner_type
      35             : #include "./base/base_uses.f90"
      36             : 
      37             :    IMPLICIT NONE
      38             :    PRIVATE
      39             : 
      40             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_apply'
      41             : 
      42             :    PUBLIC :: apply_preconditioner_fm, apply_preconditioner_dbcsr
      43             : 
      44             : CONTAINS
      45             : 
      46             : ! **************************************************************************************************
      47             : !> \brief applies a previously created preconditioner to a full matrix
      48             : !> \param preconditioner_env ...
      49             : !> \param matrix_in ...
      50             : !> \param matrix_out ...
      51             : ! **************************************************************************************************
      52       47906 :    SUBROUTINE apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out)
      53             : 
      54             :       TYPE(preconditioner_type)                          :: preconditioner_env
      55             :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_in, matrix_out
      56             : 
      57             :       CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_fm'
      58             : 
      59             :       INTEGER                                            :: handle
      60             : 
      61       47906 :       CALL timeset(routineN, handle)
      62             : 
      63       47906 :       SELECT CASE (preconditioner_env%in_use)
      64             :       CASE (0)
      65           0 :          CPABORT("No preconditioner in use")
      66             :       CASE (ot_precond_full_single)
      67        1352 :          CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
      68             :       CASE (ot_precond_full_all)
      69       27476 :          CALL apply_full_all(preconditioner_env, matrix_in, matrix_out)
      70             :       CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
      71       38156 :          SELECT CASE (preconditioner_env%solver)
      72             :          CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
      73       19078 :             CALL apply_full_single(preconditioner_env, matrix_in, matrix_out)
      74             :          CASE (ot_precond_solver_direct)
      75           0 :             CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
      76             :          CASE DEFAULT
      77       19078 :             CPABORT("Solver not implemented")
      78             :          END SELECT
      79             :       CASE DEFAULT
      80       47906 :          CPABORT("Unknown preconditioner")
      81             :       END SELECT
      82             : 
      83       47906 :       CALL timestop(handle)
      84             : 
      85       47906 :    END SUBROUTINE apply_preconditioner_fm
      86             : 
      87             : ! **************************************************************************************************
      88             : !> \brief ...
      89             : !> \param preconditioner_env ...
      90             : !> \param matrix_in ...
      91             : !> \param matrix_out ...
      92             : ! **************************************************************************************************
      93       61017 :    SUBROUTINE apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out)
      94             : 
      95             :       TYPE(preconditioner_type)                          :: preconditioner_env
      96             :       TYPE(dbcsr_type)                                   :: matrix_in, matrix_out
      97             : 
      98             :       CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr'
      99             : 
     100             :       INTEGER                                            :: handle
     101             : 
     102       61017 :       CALL timeset(routineN, handle)
     103             : 
     104       61017 :       SELECT CASE (preconditioner_env%in_use)
     105             :       CASE (0)
     106           0 :          CPABORT("No preconditioner in use")
     107             :       CASE (ot_precond_full_single)
     108         226 :          CALL apply_single(preconditioner_env, matrix_in, matrix_out)
     109             :       CASE (ot_precond_full_all)
     110       16106 :          CALL apply_all(preconditioner_env, matrix_in, matrix_out)
     111             :       CASE (ot_precond_full_kinetic, ot_precond_full_single_inverse, ot_precond_s_inverse)
     112       89370 :          SELECT CASE (preconditioner_env%solver)
     113             :          CASE (ot_precond_solver_inv_chol, ot_precond_solver_update)
     114       44685 :             CALL apply_single(preconditioner_env, matrix_in, matrix_out)
     115             :          CASE (ot_precond_solver_direct)
     116           0 :             CPABORT("Apply_full_direct not supported with ot")
     117             :             !CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out)
     118             :          CASE DEFAULT
     119       44685 :             CPABORT("Wrong solver")
     120             :          END SELECT
     121             :       CASE DEFAULT
     122       61017 :          CPABORT("Wrong preconditioner")
     123             :       END SELECT
     124             : 
     125       61017 :       CALL timestop(handle)
     126             : 
     127       61017 :    END SUBROUTINE apply_preconditioner_dbcsr
     128             : 
     129             : ! **************************************************************************************************
     130             : !> \brief apply to full matrix, complete inversion has already been done
     131             : !> \param preconditioner_env ...
     132             : !> \param matrix_in ...
     133             : !> \param matrix_out ...
     134             : ! **************************************************************************************************
     135       40860 :    SUBROUTINE apply_full_single(preconditioner_env, matrix_in, matrix_out)
     136             : 
     137             :       TYPE(preconditioner_type)                          :: preconditioner_env
     138             :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_in, matrix_out
     139             : 
     140             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'apply_full_single'
     141             : 
     142             :       INTEGER                                            :: handle, k, n
     143             : 
     144       20430 :       CALL timeset(routineN, handle)
     145             : 
     146       20430 :       CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
     147             :       CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
     148       20430 :                          matrix_in, 0.0_dp, matrix_out)
     149       20430 :       CALL timestop(handle)
     150             : 
     151       20430 :    END SUBROUTINE apply_full_single
     152             : 
     153             : ! **************************************************************************************************
     154             : !> \brief apply to dbcsr matrix, complete inversion has already been done
     155             : !> \param preconditioner_env ...
     156             : !> \param matrix_in ...
     157             : !> \param matrix_out ...
     158             : ! **************************************************************************************************
     159       44911 :    SUBROUTINE apply_single(preconditioner_env, matrix_in, matrix_out)
     160             : 
     161             :       TYPE(preconditioner_type)                          :: preconditioner_env
     162             :       TYPE(dbcsr_type)                                   :: matrix_in, matrix_out
     163             : 
     164             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'apply_single'
     165             : 
     166             :       INTEGER                                            :: handle
     167             : 
     168       44911 :       CALL timeset(routineN, handle)
     169             : 
     170       44911 :       IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) &
     171           0 :          CPABORT("NOT ASSOCIATED preconditioner_env%dbcsr_matrix")
     172             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, matrix_in, &
     173       44911 :                           0.0_dp, matrix_out)
     174             : 
     175       44911 :       CALL timestop(handle)
     176             : 
     177       44911 :    END SUBROUTINE apply_single
     178             : 
     179             : ! **************************************************************************************************
     180             : !> \brief preconditioner contains the factorization, application done by
     181             : !>        solving the linear system
     182             : !> \param preconditioner_env ...
     183             : !> \param matrix_in ...
     184             : !> \param matrix_out ...
     185             : ! **************************************************************************************************
     186           0 :    SUBROUTINE apply_full_direct(preconditioner_env, matrix_in, matrix_out)
     187             : 
     188             :       TYPE(preconditioner_type)                          :: preconditioner_env
     189             :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_in, matrix_out
     190             : 
     191             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'apply_full_direct'
     192             : 
     193             :       INTEGER                                            :: handle, k, n
     194             :       TYPE(cp_fm_type)                                   :: work
     195             : 
     196           0 :       CALL timeset(routineN, handle)
     197             : 
     198           0 :       CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
     199             :       CALL cp_fm_create(work, matrix_in%matrix_struct, name="apply_full_single", &
     200           0 :                         use_sp=matrix_in%use_sp)
     201             :       CALL cp_fm_cholesky_restore(matrix_in, k, preconditioner_env%fm, work,&
     202           0 :            &                      "SOLVE", transa="T")
     203             :       CALL cp_fm_cholesky_restore(work, k, preconditioner_env%fm, matrix_out,&
     204           0 :            &                      "SOLVE", transa="N")
     205           0 :       CALL cp_fm_release(work)
     206             : 
     207           0 :       CALL timestop(handle)
     208             : 
     209           0 :    END SUBROUTINE apply_full_direct
     210             : 
     211             : ! **************************************************************************************************
     212             : !> \brief full all to a full matrix
     213             : !> \param preconditioner_env ...
     214             : !> \param matrix_in ...
     215             : !> \param matrix_out ...
     216             : ! **************************************************************************************************
     217       82428 :    SUBROUTINE apply_full_all(preconditioner_env, matrix_in, matrix_out)
     218             : 
     219             :       TYPE(preconditioner_type)                          :: preconditioner_env
     220             :       TYPE(cp_fm_type), INTENT(IN)                       :: matrix_in, matrix_out
     221             : 
     222             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'apply_full_all'
     223             : 
     224             :       INTEGER                                            :: handle, i, j, k, n, ncol_local, &
     225             :                                                             nrow_local
     226       27476 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     227             :       REAL(KIND=dp)                                      :: dum
     228             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
     229       27476 :          POINTER                                         :: local_data
     230             :       TYPE(cp_fm_type)                                   :: matrix_tmp
     231             : 
     232       27476 :       CALL timeset(routineN, handle)
     233             : 
     234       27476 :       CALL cp_fm_get_info(matrix_in, nrow_global=n, ncol_global=k)
     235             : 
     236       27476 :       CALL cp_fm_create(matrix_tmp, matrix_in%matrix_struct, name="apply_full_all")
     237             :       CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, &
     238       27476 :                           row_indices=row_indices, col_indices=col_indices, local_data=local_data)
     239             : 
     240             :       !
     241             :       CALL parallel_gemm('T', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
     242       27476 :                          matrix_in, 0.0_dp, matrix_tmp)
     243             : 
     244             :       ! do the right scaling
     245      218734 :       DO j = 1, ncol_local
     246     2726918 :       DO i = 1, nrow_local
     247             :          dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
     248     2508184 :                           preconditioner_env%full_evals(row_indices(i)) - preconditioner_env%occ_evals(col_indices(j)))
     249     2699442 :          local_data(i, j) = local_data(i, j)*dum
     250             :       END DO
     251             :       END DO
     252             : 
     253             :       ! mult back
     254             :       CALL parallel_gemm('N', 'N', n, k, n, 1.0_dp, preconditioner_env%fm, &
     255       27476 :                          matrix_tmp, 0.0_dp, matrix_out)
     256             : 
     257       27476 :       CALL cp_fm_release(matrix_tmp)
     258             : 
     259       27476 :       CALL timestop(handle)
     260             : 
     261       27476 :    END SUBROUTINE apply_full_all
     262             : 
     263             : ! **************************************************************************************************
     264             : !> \brief full all to a dbcsr matrix
     265             : !> \param preconditioner_env ...
     266             : !> \param matrix_in ...
     267             : !> \param matrix_out ...
     268             : ! **************************************************************************************************
     269       16106 :    SUBROUTINE apply_all(preconditioner_env, matrix_in, matrix_out)
     270             : 
     271             :       TYPE(preconditioner_type)                          :: preconditioner_env
     272             :       TYPE(dbcsr_type)                                   :: matrix_in, matrix_out
     273             : 
     274             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'apply_all'
     275             : 
     276             :       INTEGER                                            :: col, col_offset, col_size, handle, i, j, &
     277             :                                                             row, row_offset, row_size
     278             :       REAL(KIND=dp)                                      :: dum
     279       16106 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: DATA
     280             :       TYPE(dbcsr_iterator_type)                          :: iter
     281             :       TYPE(dbcsr_type)                                   :: matrix_tmp
     282             : 
     283       16106 :       CALL timeset(routineN, handle)
     284             : 
     285       16106 :       CALL dbcsr_copy(matrix_tmp, matrix_in, name="apply_full_all")
     286             :       CALL dbcsr_multiply('T', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
     287       16106 :                           matrix_in, 0.0_dp, matrix_tmp)
     288             :       ! do the right scaling
     289       16106 :       CALL dbcsr_iterator_start(iter, matrix_tmp)
     290       43881 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     291             :          CALL dbcsr_iterator_next_block(iter, row, col, DATA, &
     292             :                                         row_size=row_size, col_size=col_size, &
     293       27775 :                                         row_offset=row_offset, col_offset=col_offset)
     294      250990 :          DO j = 1, col_size
     295     2058285 :          DO i = 1, row_size
     296             :             dum = 1.0_dp/MAX(preconditioner_env%energy_gap, &
     297             :                              preconditioner_env%full_evals(row_offset + i - 1) &
     298     1823401 :                              - preconditioner_env%occ_evals(col_offset + j - 1))
     299     2030510 :             DATA(i, j) = DATA(i, j)*dum
     300             :          END DO
     301             :          END DO
     302             :       END DO
     303       16106 :       CALL dbcsr_iterator_stop(iter)
     304             : 
     305             :       ! mult back
     306             :       CALL dbcsr_multiply('N', 'N', 1.0_dp, preconditioner_env%dbcsr_matrix, &
     307       16106 :                           matrix_tmp, 0.0_dp, matrix_out)
     308       16106 :       CALL dbcsr_release(matrix_tmp)
     309       16106 :       CALL timestop(handle)
     310             : 
     311       16106 :    END SUBROUTINE apply_all
     312             : 
     313             : END MODULE preconditioner_apply