1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief solves the preconditioner, contains to utility function for
8!>        fm<->dbcsr transfers, should be moved soon
9!> \par History
10!>      - [UB] 2009-05-13 Adding stable approximate inverse (full and sparse)
11!> \author Joost VandeVondele (09.2002)
12! **************************************************************************************************
13MODULE preconditioner_solvers
14   USE arnoldi_api,                     ONLY: arnoldi_data_type,&
15                                              arnoldi_ev,&
16                                              deallocate_arnoldi_data,&
17                                              get_selected_ritz_val,&
18                                              setup_arnoldi_data
19   USE bibliography,                    ONLY: Schiffmann2015,&
20                                              cite_reference
21   USE cp_blacs_env,                    ONLY: cp_blacs_env_type
22   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
23                                              copy_fm_to_dbcsr
24   USE cp_fm_basic_linalg,              ONLY: cp_fm_upper_to_full
25   USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose,&
26                                              cp_fm_cholesky_invert
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_set_all,&
33                                              cp_fm_type
34   USE cp_para_types,                   ONLY: cp_para_env_type
35   USE dbcsr_api,                       ONLY: &
36        dbcsr_create, dbcsr_filter, dbcsr_get_info, dbcsr_get_occupation, dbcsr_init_p, &
37        dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_real_8
38   USE input_constants,                 ONLY: ot_precond_solver_default,&
39                                              ot_precond_solver_direct,&
40                                              ot_precond_solver_inv_chol,&
41                                              ot_precond_solver_update
42   USE iterate_matrix,                  ONLY: invert_Hotelling
43   USE kinds,                           ONLY: dp
44   USE preconditioner_types,            ONLY: preconditioner_type
45#include "./base/base_uses.f90"
46
47   IMPLICIT NONE
48
49   PRIVATE
50
51   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_solvers'
52
53   PUBLIC :: solve_preconditioner, transfer_fm_to_dbcsr, transfer_dbcsr_to_fm
54
55CONTAINS
56
57! **************************************************************************************************
58!> \brief ...
59!> \param my_solver_type ...
60!> \param preconditioner_env ...
61!> \param matrix_s ...
62!> \param matrix_h ...
63! **************************************************************************************************
64   SUBROUTINE solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, &
65                                   matrix_h)
66      INTEGER                                            :: my_solver_type
67      TYPE(preconditioner_type)                          :: preconditioner_env
68      TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
69      TYPE(dbcsr_type), POINTER                          :: matrix_h
70
71      CHARACTER(len=*), PARAMETER :: routineN = 'solve_preconditioner', &
72         routineP = moduleN//':'//routineN
73
74      REAL(dp)                                           :: occ_matrix
75
76! here comes the solver
77
78      SELECT CASE (my_solver_type)
79      CASE (ot_precond_solver_inv_chol)
80         !
81         ! compute the full inverse
82         preconditioner_env%solver = ot_precond_solver_inv_chol
83         CALL make_full_inverse_cholesky(preconditioner_env, matrix_s)
84      CASE (ot_precond_solver_direct)
85         !
86         ! prepare for the direct solver
87         preconditioner_env%solver = ot_precond_solver_direct
88         CALL make_full_fact_cholesky(preconditioner_env, matrix_s)
89      CASE (ot_precond_solver_update)
90         !
91         ! uses an update of the full inverse (needs to be computed the first time)
92         ! make sure preconditioner_env is not destroyed in between
93         occ_matrix = 1.0_dp
94         IF (ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
95            IF (preconditioner_env%condition_num < 0.0_dp) &
96               CALL estimate_cond_num(preconditioner_env%sparse_matrix, preconditioner_env%condition_num)
97            CALL dbcsr_filter(preconditioner_env%sparse_matrix, &
98                              1.0_dp/preconditioner_env%condition_num*0.01_dp)
99            occ_matrix = dbcsr_get_occupation(preconditioner_env%sparse_matrix)
100         END IF
101         ! check whether we are in the first step and if it is a good idea to use cholesky (matrix sparsity)
102         IF (preconditioner_env%solver .NE. ot_precond_solver_update .AND. occ_matrix > 0.5_dp) THEN
103            preconditioner_env%solver = ot_precond_solver_update
104            CALL make_full_inverse_cholesky(preconditioner_env, matrix_s)
105         ELSE
106            preconditioner_env%solver = ot_precond_solver_update
107            CALL make_inverse_update(preconditioner_env, matrix_h)
108         END IF
109      CASE (ot_precond_solver_default)
110         preconditioner_env%solver = ot_precond_solver_default
111      CASE DEFAULT
112         !
113         CPABORT("Doesn't know this type of solver")
114      END SELECT
115
116   END SUBROUTINE solve_preconditioner
117
118! **************************************************************************************************
119!> \brief Compute the inverse using cholseky factorization
120!> \param preconditioner_env ...
121!> \param matrix_s ...
122! **************************************************************************************************
123   SUBROUTINE make_full_inverse_cholesky(preconditioner_env, matrix_s)
124
125      TYPE(preconditioner_type)                          :: preconditioner_env
126      TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
127
128      CHARACTER(len=*), PARAMETER :: routineN = 'make_full_inverse_cholesky', &
129         routineP = moduleN//':'//routineN
130
131      INTEGER                                            :: handle, info
132      TYPE(cp_fm_type), POINTER                          :: fm, fm_work
133
134      CALL timeset(routineN, handle)
135
136      ! Maybe we will get a sparse Cholesky at a given point then this can go,
137      ! if stuff was stored in fm anyway this simple returns
138      CALL transfer_dbcsr_to_fm(preconditioner_env%sparse_matrix, preconditioner_env%fm, &
139                                preconditioner_env%para_env, preconditioner_env%ctxt)
140      fm => preconditioner_env%fm
141
142      NULLIFY (fm_work)
143
144      CALL cp_fm_create(fm_work, fm%matrix_struct, name="fm_work")
145      !
146      ! compute the inverse of SPD matrix fm using the Cholesky factorization
147      CALL cp_fm_cholesky_decompose(fm, info_out=info)
148
149      !
150      ! if fm not SPD we go with the overlap matrix
151      IF (info /= 0) THEN
152         !
153         ! just the overlap matrix
154         IF (PRESENT(matrix_s)) THEN
155            CALL copy_dbcsr_to_fm(matrix_s, fm)
156            CALL cp_fm_cholesky_decompose(fm)
157         ELSE
158            CALL cp_fm_set_all(fm, alpha=0._dp, beta=1._dp)
159         ENDIF
160      ENDIF
161      CALL cp_fm_cholesky_invert(fm)
162
163      CALL cp_fm_upper_to_full(fm, fm_work)
164      CALL cp_fm_release(fm_work)
165
166      CALL timestop(handle)
167
168   END SUBROUTINE make_full_inverse_cholesky
169
170! **************************************************************************************************
171!> \brief Only perform the factorization, can be used later to solve the linear
172!>        system on the fly
173!> \param preconditioner_env ...
174!> \param matrix_s ...
175! **************************************************************************************************
176   SUBROUTINE make_full_fact_cholesky(preconditioner_env, matrix_s)
177
178      TYPE(preconditioner_type)                          :: preconditioner_env
179      TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
180
181      CHARACTER(len=*), PARAMETER :: routineN = 'make_full_fact_cholesky', &
182         routineP = moduleN//':'//routineN
183
184      INTEGER                                            :: handle, info_out
185      TYPE(cp_fm_type), POINTER                          :: fm
186
187      CALL timeset(routineN, handle)
188
189      ! Maybe we will get a sparse Cholesky at a given point then this can go,
190      ! if stuff was stored in fm anyway this simple returns
191      CALL transfer_dbcsr_to_fm(preconditioner_env%sparse_matrix, preconditioner_env%fm, &
192                                preconditioner_env%para_env, preconditioner_env%ctxt)
193
194      fm => preconditioner_env%fm
195      !
196      ! compute the inverse of SPD matrix fm using the Cholesky factorization
197      CALL cp_fm_cholesky_decompose(fm, info_out=info_out)
198      !
199      ! if fm not SPD we go with the overlap matrix
200      IF (info_out .NE. 0) THEN
201         !
202         ! just the overlap matrix
203         IF (PRESENT(matrix_s)) THEN
204            CALL copy_dbcsr_to_fm(matrix_s, fm)
205            CALL cp_fm_cholesky_decompose(fm)
206         ELSE
207            CALL cp_fm_set_all(fm, alpha=0._dp, beta=1._dp)
208         ENDIF
209      ENDIF
210
211      CALL timestop(handle)
212
213   END SUBROUTINE make_full_fact_cholesky
214
215! **************************************************************************************************
216!> \brief computes an approximate inverse using Hotelling iterations
217!> \param preconditioner_env ...
218!> \param matrix_h as S is not always present this is a safe template for the transfer
219! **************************************************************************************************
220   SUBROUTINE make_inverse_update(preconditioner_env, matrix_h)
221      TYPE(preconditioner_type)                          :: preconditioner_env
222      TYPE(dbcsr_type), POINTER                          :: matrix_h
223
224      CHARACTER(len=*), PARAMETER :: routineN = 'make_inverse_update', &
225         routineP = moduleN//':'//routineN
226
227      INTEGER                                            :: handle
228      LOGICAL                                            :: use_guess
229      REAL(KIND=dp)                                      :: filter_eps
230
231      CALL timeset(routineN, handle)
232      use_guess = .TRUE.
233      !
234      ! uses an update of the full inverse (needs to be computed the first time)
235      ! make sure preconditioner_env is not destroyed in between
236
237      CALL cite_reference(Schiffmann2015)
238
239      ! Maybe I gonna add a fm Hotelling, ... for now the same as above make sure we are dbcsr
240      CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%sparse_matrix, matrix_h)
241      IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN
242         use_guess = .FALSE.
243         CALL dbcsr_init_p(preconditioner_env%dbcsr_matrix)
244         CALL dbcsr_create(preconditioner_env%dbcsr_matrix, "prec_dbcsr", &
245                           template=matrix_h, matrix_type=dbcsr_type_no_symmetry)
246      END IF
247
248      ! Try to get a reasonbale guess for the filtering threshold
249      filter_eps = 1.0_dp/preconditioner_env%condition_num*0.1_dp
250
251      ! Agressive filtering on the initial guess is needed to avoid fill ins and retain sparsity
252      CALL dbcsr_filter(preconditioner_env%dbcsr_matrix, filter_eps*100.0_dp)
253      ! We don't need a high accuracy for the inverse so 0.4 is reasonable for convergence
254      CALL invert_Hotelling(preconditioner_env%dbcsr_matrix, preconditioner_env%sparse_matrix, filter_eps*10.0_dp, &
255                            use_inv_as_guess=use_guess, norm_convergence=0.4_dp, filter_eps=filter_eps)
256
257      CALL timestop(handle)
258
259   END SUBROUTINE make_inverse_update
260
261! **************************************************************************************************
262!> \brief computes an approximation to the condition number of a matrix using
263!>        arnoldi iterations
264!> \param matrix ...
265!> \param cond_num ...
266! **************************************************************************************************
267   SUBROUTINE estimate_cond_num(matrix, cond_num)
268      TYPE(dbcsr_type), POINTER                          :: matrix
269      REAL(KIND=dp)                                      :: cond_num
270
271      CHARACTER(len=*), PARAMETER :: routineN = 'estimate_cond_num', &
272         routineP = moduleN//':'//routineN
273
274      INTEGER                                            :: handle
275      REAL(KIND=dp)                                      :: max_ev, min_ev
276      TYPE(arnoldi_data_type)                            :: my_arnoldi
277      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrices
278
279      CALL timeset(routineN, handle)
280
281      ! its better to do 2 calculations as the maximum should quickly converge and the minimum won't need iram
282      ALLOCATE (matrices(1))
283      matrices(1)%matrix => matrix
284      ! compute the minimum ev
285      CALL setup_arnoldi_data(my_arnoldi, matrices, max_iter=20, threshold=5.0E-4_dp, selection_crit=2, &
286                              nval_request=1, nrestarts=15, generalized_ev=.FALSE., iram=.FALSE.)
287      CALL arnoldi_ev(matrices, my_arnoldi)
288      max_ev = REAL(get_selected_ritz_val(my_arnoldi, 1), dp)
289      CALL deallocate_arnoldi_data(my_arnoldi)
290
291      CALL setup_arnoldi_data(my_arnoldi, matrices, max_iter=20, threshold=5.0E-4_dp, selection_crit=3, &
292                              nval_request=1, nrestarts=15, generalized_ev=.FALSE., iram=.FALSE.)
293      CALL arnoldi_ev(matrices, my_arnoldi)
294      min_ev = REAL(get_selected_ritz_val(my_arnoldi, 1), dp)
295      CALL deallocate_arnoldi_data(my_arnoldi)
296
297      cond_num = max_ev/min_ev
298      DEALLOCATE (matrices)
299
300      CALL timestop(handle)
301   END SUBROUTINE estimate_cond_num
302
303! **************************************************************************************************
304!> \brief transfers a full matrix to a dbcsr
305!> \param fm_matrix a full matrix gets deallocated in the end
306!> \param dbcsr_matrix a dbcsr matrix, gets create from a template
307!> \param template_mat the template which is used for the structure
308! **************************************************************************************************
309   SUBROUTINE transfer_fm_to_dbcsr(fm_matrix, dbcsr_matrix, template_mat)
310
311      TYPE(cp_fm_type), POINTER                          :: fm_matrix
312      TYPE(dbcsr_type), POINTER                          :: dbcsr_matrix, template_mat
313
314      CHARACTER(len=*), PARAMETER :: routineN = 'transfer_fm_to_dbcsr', &
315         routineP = moduleN//':'//routineN
316
317      INTEGER                                            :: handle
318
319      CALL timeset(routineN, handle)
320      IF (ASSOCIATED(fm_matrix)) THEN
321         IF (.NOT. ASSOCIATED(dbcsr_matrix)) THEN
322            CALL dbcsr_init_p(dbcsr_matrix)
323            CALL dbcsr_create(dbcsr_matrix, template=template_mat, &
324                              name="preconditioner_env%dbcsr_matrix", &
325                              matrix_type=dbcsr_type_no_symmetry, &
326                              nze=0, data_type=dbcsr_type_real_8)
327         ENDIF
328!       CALL cp_fm_create(fm_tmp, matrix_struct=fm_matrix%matrix_struct)
329!       CALL cp_fm_upper_to_full(fm_matrix, fm_tmp)
330         CALL copy_fm_to_dbcsr(fm_matrix, dbcsr_matrix)
331!       CALL cp_fm_release(fm_tmp)
332         CALL cp_fm_release(fm_matrix)
333      END IF
334
335      CALL timestop(handle)
336
337   END SUBROUTINE transfer_fm_to_dbcsr
338
339! **************************************************************************************************
340!> \brief transfers a dbcsr to a full matrix
341!> \param dbcsr_matrix a dbcsr matrix, gets deallocated at the end
342!> \param fm_matrix a full matrix gets created if not yet done
343!> \param para_env the para_env
344!> \param context the blacs context
345! **************************************************************************************************
346   SUBROUTINE transfer_dbcsr_to_fm(dbcsr_matrix, fm_matrix, para_env, context)
347
348      TYPE(dbcsr_type), POINTER                          :: dbcsr_matrix
349      TYPE(cp_fm_type), POINTER                          :: fm_matrix
350      TYPE(cp_para_env_type), POINTER                    :: para_env
351      TYPE(cp_blacs_env_type), POINTER                   :: context
352
353      CHARACTER(len=*), PARAMETER :: routineN = 'transfer_dbcsr_to_fm', &
354         routineP = moduleN//':'//routineN
355
356      INTEGER                                            :: handle, n
357      TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
358
359      CALL timeset(routineN, handle)
360      IF (ASSOCIATED(dbcsr_matrix)) THEN
361         NULLIFY (fm_struct_tmp)
362
363         IF (ASSOCIATED(fm_matrix)) THEN
364            CALL cp_fm_release(fm_matrix)
365         ENDIF
366
367         CALL dbcsr_get_info(dbcsr_matrix, nfullrows_total=n)
368         CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
369                                  context=context, para_env=para_env)
370         CALL cp_fm_create(fm_matrix, fm_struct_tmp)
371         CALL cp_fm_struct_release(fm_struct_tmp)
372
373         CALL copy_dbcsr_to_fm(dbcsr_matrix, fm_matrix)
374         CALL dbcsr_release(dbcsr_matrix)
375         DEALLOCATE (dbcsr_matrix)
376      END IF
377
378      CALL timestop(handle)
379
380   END SUBROUTINE transfer_dbcsr_to_fm
381
382END MODULE preconditioner_solvers
383