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 The methods which allow to analyze and manipulate the arnoldi procedure
10 : !> The main routine and this should eb the only public access point for the
11 : !> method
12 : !> \par History
13 : !> 2014.09 created [Florian Schiffmann]
14 : !> \author Florian Schiffmann
15 : ! **************************************************************************************************
16 :
17 : MODULE arnoldi_data_methods
18 : USE arnoldi_types, ONLY: &
19 : arnoldi_control_type, arnoldi_data_c_type, arnoldi_data_d_type, arnoldi_data_s_type, &
20 : arnoldi_data_type, arnoldi_data_z_type, get_control, get_data_c, get_data_d, get_data_s, &
21 : get_data_z, get_evals_c, get_evals_d, get_evals_s, get_evals_z, get_sel_ind, has_d_cmplx, &
22 : has_d_real, set_control, set_data_c, set_data_d, set_data_s, &
23 : set_data_z
24 : USE cp_dbcsr_api, ONLY: &
25 : dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_data_p, dbcsr_get_data_type, &
26 : dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_mp_grid_setup, dbcsr_p_type, dbcsr_release, &
27 : dbcsr_type, dbcsr_type_complex_8, dbcsr_type_real_8, dbcsr_type_symmetric
28 : USE dbcsr_vector, ONLY: create_col_vec_from_matrix
29 : USE kinds, ONLY: real_8
30 : USE util, ONLY: sort
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : PUBLIC :: select_evals, get_selected_ritz_val, arnoldi_is_converged, &
38 : arnoldi_data_type, get_nrestart, set_arnoldi_initial_vector, &
39 : setup_arnoldi_data, deallocate_arnoldi_data, get_selected_ritz_vector
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief This routine sets the environment for the arnoldi iteration and
45 : !> the krylov subspace creation. All simulation parameters have to be given
46 : !> at this stage so the rest can run fully automated
47 : !> In addition, this routine allocates the data necessary for
48 : !> \param arnoldi_data this type which gets filled with information and
49 : !> on output contains all information necessary to extract
50 : !> whatever the user desires
51 : !> \param matrix vector of matrices, only the first gets used to get some dimensions
52 : !> and parallel information needed later on
53 : !> \param max_iter maximum dimension of the krylov subspace
54 : !> \param threshold convergence threshold, this is used for both subspace and eigenval
55 : !> \param selection_crit integer defining according to which criterion the
56 : !> eigenvalues are selected for the subspace
57 : !> \param nval_request for some sel_crit useful, how many eV to select
58 : !> \param nrestarts ...
59 : !> \param generalized_ev ...
60 : !> \param iram ...
61 : ! **************************************************************************************************
62 128911 : SUBROUTINE setup_arnoldi_data(arnoldi_data, matrix, max_iter, threshold, selection_crit, &
63 : nval_request, nrestarts, generalized_ev, iram)
64 : TYPE(arnoldi_data_type) :: arnoldi_data
65 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
66 : INTEGER :: max_iter
67 : REAL(real_8) :: threshold
68 : INTEGER :: selection_crit, nval_request, nrestarts
69 : LOGICAL :: generalized_ev, iram
70 :
71 : CALL setup_arnoldi_control(arnoldi_data, matrix, max_iter, threshold, selection_crit, &
72 128911 : nval_request, nrestarts, generalized_ev, iram)
73 :
74 128911 : SELECT CASE (dbcsr_get_data_type(matrix(1)%matrix))
75 : CASE (dbcsr_type_real_8)
76 128911 : CALL setup_arnoldi_data_d(arnoldi_data, matrix, max_iter)
77 : CASE (dbcsr_type_complex_8)
78 128911 : CALL setup_arnoldi_data_z(arnoldi_data, matrix, max_iter)
79 : END SELECT
80 :
81 128911 : END SUBROUTINE setup_arnoldi_data
82 :
83 : ! **************************************************************************************************
84 : !> \brief Creates the control type for arnoldi, see above for details
85 : !> \param arnoldi_data ...
86 : !> \param matrix ...
87 : !> \param max_iter ...
88 : !> \param threshold ...
89 : !> \param selection_crit ...
90 : !> \param nval_request ...
91 : !> \param nrestarts ...
92 : !> \param generalized_ev ...
93 : !> \param iram ...
94 : ! **************************************************************************************************
95 128911 : SUBROUTINE setup_arnoldi_control(arnoldi_data, matrix, max_iter, threshold, selection_crit, nval_request, &
96 : nrestarts, generalized_ev, iram)
97 : TYPE(arnoldi_data_type) :: arnoldi_data
98 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
99 : INTEGER :: max_iter
100 : REAL(real_8) :: threshold
101 : INTEGER :: selection_crit, nval_request, nrestarts
102 : LOGICAL :: generalized_ev, iram
103 :
104 : INTEGER :: pcol_handle, group_handle
105 : LOGICAL :: subgroups_defined
106 : TYPE(arnoldi_control_type), POINTER :: control
107 : TYPE(dbcsr_distribution_type) :: distri
108 :
109 0 : ALLOCATE (control)
110 : ! Fill the information which will later on control the behavior of the arnoldi method and allow synchronization
111 128911 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, distribution=distri)
112 128911 : CALL dbcsr_mp_grid_setup(distri)
113 : CALL dbcsr_distribution_get(distri, &
114 : group=group_handle, &
115 : mynode=control%myproc, &
116 : subgroups_defined=subgroups_defined, &
117 128911 : pcol_group=pcol_handle)
118 :
119 128911 : CALL control%mp_group%set_handle(group_handle)
120 128911 : CALL control%pcol_group%set_handle(pcol_handle)
121 :
122 128911 : IF (.NOT. subgroups_defined) &
123 0 : CPABORT("arnoldi only with subgroups")
124 :
125 128911 : control%symmetric = .FALSE.
126 : ! Will need a fix for complex because there it has to be hermitian
127 128911 : IF (SIZE(matrix) == 1) &
128 124987 : control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
129 :
130 : ! Set the control parameters
131 128911 : control%max_iter = max_iter
132 128911 : control%current_step = 0
133 128911 : control%selection_crit = selection_crit
134 128911 : control%nval_req = nval_request
135 128911 : control%threshold = threshold
136 128911 : control%converged = .FALSE.
137 128911 : control%has_initial_vector = .FALSE.
138 128911 : control%iram = iram
139 128911 : control%nrestart = nrestarts
140 128911 : control%generalized_ev = generalized_ev
141 :
142 128911 : IF (control%nval_req > 1 .AND. control%nrestart > 0 .AND. .NOT. control%iram) &
143 : CALL cp_abort(__LOCATION__, 'with more than one eigenvalue requested '// &
144 0 : 'internal restarting with a previous EVEC is a bad idea, set IRAM or nrestsart=0')
145 :
146 : ! some checks for the generalized EV mode
147 128911 : IF (control%generalized_ev .AND. selection_crit == 1) &
148 : CALL cp_abort(__LOCATION__, &
149 0 : 'generalized ev can only highest OR lowest EV')
150 128911 : IF (control%generalized_ev .AND. nval_request .NE. 1) &
151 : CALL cp_abort(__LOCATION__, &
152 0 : 'generalized ev can only compute one EV at the time')
153 128911 : IF (control%generalized_ev .AND. control%nrestart == 0) &
154 : CALL cp_abort(__LOCATION__, &
155 0 : 'outer loops are mandatory for generalized EV, set nrestart appropriatly')
156 128911 : IF (SIZE(matrix) .NE. 2 .AND. control%generalized_ev) &
157 : CALL cp_abort(__LOCATION__, &
158 0 : 'generalized ev needs exactly two matrices as input (2nd is the metric)')
159 :
160 386733 : ALLOCATE (control%selected_ind(max_iter))
161 128911 : CALL set_control(arnoldi_data, control)
162 :
163 257822 : END SUBROUTINE setup_arnoldi_control
164 :
165 : ! **************************************************************************************************
166 : !> \brief Deallocate the data in arnoldi_data
167 : !> \param arnoldi_data ...
168 : !> \param ind ...
169 : !> \param matrix ...
170 : !> \param vector ...
171 : ! **************************************************************************************************
172 133455 : SUBROUTINE get_selected_ritz_vector(arnoldi_data, ind, matrix, vector)
173 : TYPE(arnoldi_data_type) :: arnoldi_data
174 : INTEGER :: ind
175 : TYPE(dbcsr_type) :: matrix, vector
176 :
177 133455 : IF (has_d_real(arnoldi_data)) CALL get_selected_ritz_vector_d(arnoldi_data, ind, matrix, vector)
178 133455 : IF (has_d_cmplx(arnoldi_data)) CALL get_selected_ritz_vector_z(arnoldi_data, ind, matrix, vector)
179 :
180 133455 : END SUBROUTINE get_selected_ritz_vector
181 :
182 : ! **************************************************************************************************
183 : !> \brief Deallocate the data in arnoldi_data
184 : !> \param arnoldi_data ...
185 : ! **************************************************************************************************
186 128911 : SUBROUTINE deallocate_arnoldi_data(arnoldi_data)
187 : TYPE(arnoldi_data_type) :: arnoldi_data
188 :
189 : TYPE(arnoldi_control_type), POINTER :: control
190 :
191 128911 : IF (has_d_real(arnoldi_data)) CALL deallocate_arnoldi_data_d(arnoldi_data)
192 128911 : IF (has_d_cmplx(arnoldi_data)) CALL deallocate_arnoldi_data_z(arnoldi_data)
193 :
194 128911 : control => get_control(arnoldi_data)
195 128911 : DEALLOCATE (control%selected_ind)
196 128911 : DEALLOCATE (control)
197 :
198 128911 : END SUBROUTINE deallocate_arnoldi_data
199 :
200 : ! **************************************************************************************************
201 : !> \brief perform the selection of eigenvalues, fills the selected_ind array
202 : !> \param arnoldi_data ...
203 : ! **************************************************************************************************
204 129914 : SUBROUTINE select_evals(arnoldi_data)
205 : TYPE(arnoldi_data_type) :: arnoldi_data
206 :
207 129914 : IF (has_d_real(arnoldi_data)) CALL select_evals_d(arnoldi_data)
208 129914 : IF (has_d_cmplx(arnoldi_data)) CALL select_evals_z(arnoldi_data)
209 :
210 129914 : END SUBROUTINE select_evals
211 :
212 : ! **************************************************************************************************
213 : !> \brief set a new selection type, if you notice you didn't like the initial one
214 : !> \param ar_data ...
215 : !> \param itype ...
216 : ! **************************************************************************************************
217 0 : SUBROUTINE set_eval_selection(ar_data, itype)
218 : TYPE(arnoldi_data_type) :: ar_data
219 : INTEGER :: itype
220 :
221 : TYPE(arnoldi_control_type), POINTER :: control
222 :
223 0 : control => get_control(ar_data)
224 0 : control%selection_crit = itype
225 0 : END SUBROUTINE set_eval_selection
226 :
227 : ! **************************************************************************************************
228 : !> \brief returns the number of restarts allowed for arnoldi
229 : !> \param arnoldi_data ...
230 : !> \return ...
231 : ! **************************************************************************************************
232 128715 : FUNCTION get_nrestart(arnoldi_data) RESULT(nrestart)
233 : TYPE(arnoldi_data_type) :: arnoldi_data
234 : INTEGER :: nrestart
235 :
236 : TYPE(arnoldi_control_type), POINTER :: control
237 :
238 128715 : control => get_control(arnoldi_data)
239 128715 : nrestart = control%nrestart
240 :
241 128715 : END FUNCTION get_nrestart
242 :
243 : ! **************************************************************************************************
244 : !> \brief get the number of eigenvalues matching the search criterion
245 : !> \param ar_data ...
246 : !> \return ...
247 : ! **************************************************************************************************
248 249280 : FUNCTION get_nval_out(ar_data) RESULT(nval_out)
249 : TYPE(arnoldi_data_type) :: ar_data
250 : INTEGER :: nval_out
251 :
252 : TYPE(arnoldi_control_type), POINTER :: control
253 :
254 249280 : control => get_control(ar_data)
255 249280 : nval_out = control%nval_out
256 :
257 249280 : END FUNCTION get_nval_out
258 :
259 : ! **************************************************************************************************
260 : !> \brief get the dimension of the krylov space. Can be less than max_iter
261 : !> if subspace converged early
262 : !> \param ar_data ...
263 : !> \return ...
264 : ! **************************************************************************************************
265 133455 : FUNCTION get_subsp_size(ar_data) RESULT(current_step)
266 : TYPE(arnoldi_data_type) :: ar_data
267 : INTEGER :: current_step
268 :
269 : TYPE(arnoldi_control_type), POINTER :: control
270 :
271 133455 : control => get_control(ar_data)
272 133455 : current_step = control%current_step
273 :
274 133455 : END FUNCTION get_subsp_size
275 :
276 : ! **************************************************************************************************
277 : !> \brief Find out whether the method with the current search criterion is converged
278 : !> \param ar_data ...
279 : !> \return ...
280 : ! **************************************************************************************************
281 120761 : FUNCTION arnoldi_is_converged(ar_data) RESULT(converged)
282 : TYPE(arnoldi_data_type) :: ar_data
283 : LOGICAL :: converged
284 :
285 : TYPE(arnoldi_control_type), POINTER :: control
286 :
287 120761 : control => get_control(ar_data)
288 120761 : converged = control%converged
289 :
290 120761 : END FUNCTION
291 :
292 : ! **************************************************************************************************
293 : !> \brief get a single specific Ritz value from the set of selected
294 : !> \param ar_data ...
295 : !> \param ind ...
296 : !> \return ...
297 : ! **************************************************************************************************
298 249280 : FUNCTION get_selected_ritz_val(ar_data, ind) RESULT(eval_out)
299 : TYPE(arnoldi_data_type) :: ar_data
300 : INTEGER :: ind
301 : COMPLEX(real_8) :: eval_out
302 :
303 249280 : COMPLEX(real_8), DIMENSION(:), POINTER :: evals_d
304 : INTEGER :: ev_ind
305 249280 : INTEGER, DIMENSION(:), POINTER :: selected_ind
306 :
307 249280 : IF (ind .GT. get_nval_out(ar_data)) &
308 0 : CPABORT('outside range of indexed evals')
309 :
310 249280 : selected_ind => get_sel_ind(ar_data)
311 249280 : ev_ind = selected_ind(ind)
312 249280 : IF (has_d_real(ar_data)) THEN
313 249280 : evals_d => get_evals_d(ar_data); eval_out = evals_d(ev_ind)
314 0 : ELSE IF (has_d_cmplx(ar_data)) THEN
315 0 : evals_d => get_evals_z(ar_data); eval_out = evals_d(ev_ind)
316 : END IF
317 :
318 249280 : END FUNCTION get_selected_ritz_val
319 :
320 : ! **************************************************************************************************
321 : !> \brief Get all Ritz values of the selected set. eval_out has to be allocated
322 : !> at least the size of get_neval_out()
323 : !> \param ar_data ...
324 : !> \param eval_out ...
325 : ! **************************************************************************************************
326 0 : SUBROUTINE get_all_selected_ritz_val(ar_data, eval_out)
327 : TYPE(arnoldi_data_type) :: ar_data
328 : COMPLEX(real_8), DIMENSION(:) :: eval_out
329 :
330 0 : COMPLEX(real_8), DIMENSION(:), POINTER :: evals_d
331 : INTEGER :: ev_ind, ind
332 0 : INTEGER, DIMENSION(:), POINTER :: selected_ind
333 :
334 0 : NULLIFY (evals_d)
335 0 : IF (SIZE(eval_out) .LT. get_nval_out(ar_data)) &
336 0 : CPABORT('array for eval output too small')
337 0 : selected_ind => get_sel_ind(ar_data)
338 :
339 0 : IF (has_d_real(ar_data)) evals_d => get_evals_d(ar_data)
340 0 : IF (has_d_cmplx(ar_data)) evals_d => get_evals_d(ar_data)
341 :
342 0 : DO ind = 1, get_nval_out(ar_data)
343 0 : ev_ind = selected_ind(ind)
344 0 : IF (ASSOCIATED(evals_d)) eval_out(ind) = evals_d(ev_ind)
345 : END DO
346 :
347 0 : END SUBROUTINE get_all_selected_ritz_val
348 :
349 : ! **************************************************************************************************
350 : !> \brief ...
351 : !> \param ar_data ...
352 : !> \param vector ...
353 : ! **************************************************************************************************
354 4786 : SUBROUTINE set_arnoldi_initial_vector(ar_data, vector)
355 : TYPE(arnoldi_data_type) :: ar_data
356 : TYPE(dbcsr_type) :: vector
357 :
358 : TYPE(arnoldi_control_type), POINTER :: control
359 :
360 4786 : control => get_control(ar_data)
361 4786 : control%has_initial_vector = .TRUE.
362 :
363 4786 : IF (has_d_real(ar_data)) CALL set_initial_vector_d(ar_data, vector)
364 4786 : IF (has_d_cmplx(ar_data)) CALL set_initial_vector_z(ar_data, vector)
365 :
366 4786 : END SUBROUTINE set_arnoldi_initial_vector
367 :
368 : !!! Here come the methods handling the selection of eigenvalues and eigenvectors !!!
369 : !!! If you want a personal method, simply created a Subroutine returning the index
370 : !!! array selected ind which contains as the first nval_out entries the index of the evals
371 :
372 : #:include 'arnoldi.fypp'
373 : #:for nametype1, type_prec, real_zero, nametype_zero, type_nametype1, vartype in inst_params_1
374 : ! **************************************************************************************************
375 : !> \brief ...
376 : !> \param arnoldi_data ...
377 : ! **************************************************************************************************
378 129914 : SUBROUTINE select_evals_${nametype1}$ (arnoldi_data)
379 : TYPE(arnoldi_data_type) :: arnoldi_data
380 :
381 : INTEGER :: my_crit, last_el, my_ind, i
382 : REAL(${type_prec}$) :: convergence
383 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
384 : TYPE(arnoldi_control_type), POINTER :: control
385 :
386 129914 : control => get_control(arnoldi_data)
387 129914 : ar_data => get_data_${nametype1}$ (arnoldi_data)
388 :
389 129914 : last_el = control%current_step
390 129914 : convergence = REAL(0.0, ${type_prec}$)
391 129914 : my_crit = control%selection_crit
392 129914 : control%nval_out = MIN(control%nval_req, control%current_step)
393 120815 : SELECT CASE (my_crit)
394 : ! minimum and maximum real eval
395 : CASE (1)
396 120815 : CALL index_min_max_real_eval_${nametype1}$ (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
397 : ! n maximum real eval
398 : CASE (2)
399 4236 : CALL index_nmax_real_eval_${nametype1}$ (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
400 : ! n minimum real eval
401 : CASE (3)
402 4863 : CALL index_nmin_real_eval_${nametype1}$ (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
403 : CASE DEFAULT
404 129914 : CPABORT("unknown selection index")
405 : END SELECT
406 : ! test whether we are converged
407 380643 : DO i = 1, control%nval_out
408 250729 : my_ind = control%selected_ind(i)
409 : convergence = MAX(convergence, &
410 380643 : ABS(ar_data%revec(last_el, my_ind)*ar_data%Hessenberg(last_el + 1, last_el)))
411 : END DO
412 129914 : control%converged = convergence .LT. control%threshold
413 :
414 129914 : END SUBROUTINE select_evals_${nametype1}$
415 :
416 : ! **************************************************************************************************
417 : !> \brief ...
418 : !> \param evals ...
419 : !> \param current_step ...
420 : !> \param selected_ind ...
421 : !> \param neval ...
422 : ! **************************************************************************************************
423 120815 : SUBROUTINE index_min_max_real_eval_${nametype1}$ (evals, current_step, selected_ind, neval)
424 : COMPLEX(${type_prec}$), DIMENSION(:) :: evals
425 : INTEGER, INTENT(IN) :: current_step
426 : INTEGER, DIMENSION(:) :: selected_ind
427 : INTEGER :: neval
428 :
429 241630 : INTEGER, DIMENSION(current_step) :: indexing
430 241630 : REAL(${type_prec}$), DIMENSION(current_step) :: tmp_array
431 : INTEGER :: i
432 :
433 120815 : neval = 0
434 4001103 : selected_ind = 0
435 693228 : tmp_array(1:current_step) = REAL(evals(1:current_step), ${type_prec}$)
436 120815 : CALL sort(tmp_array, current_step, indexing)
437 120815 : DO i = 1, current_step
438 120815 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(${real_zero}$)) THEN
439 120815 : selected_ind(1) = indexing(i)
440 120815 : neval = neval + 1
441 120815 : EXIT
442 : END IF
443 : END DO
444 120815 : DO i = current_step, 1, -1
445 120815 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(${real_zero}$)) THEN
446 120815 : selected_ind(2) = indexing(i)
447 120815 : neval = neval + 1
448 120815 : EXIT
449 : END IF
450 : END DO
451 :
452 120815 : END SUBROUTINE index_min_max_real_eval_${nametype1}$
453 :
454 : ! **************************************************************************************************
455 : !> \brief ...
456 : !> \param evals ...
457 : !> \param current_step ...
458 : !> \param selected_ind ...
459 : !> \param neval ...
460 : ! **************************************************************************************************
461 4236 : SUBROUTINE index_nmax_real_eval_${nametype1}$ (evals, current_step, selected_ind, neval)
462 : COMPLEX(${type_prec}$), DIMENSION(:) :: evals
463 : INTEGER, INTENT(IN) :: current_step
464 : INTEGER, DIMENSION(:) :: selected_ind
465 : INTEGER :: neval
466 :
467 : INTEGER :: i, nlimit
468 8472 : INTEGER, DIMENSION(current_step) :: indexing
469 8472 : REAL(${type_prec}$), DIMENSION(current_step) :: tmp_array
470 :
471 4236 : nlimit = neval; neval = 0
472 88956 : selected_ind = 0
473 24683 : tmp_array(1:current_step) = REAL(evals(1:current_step), ${type_prec}$)
474 4236 : CALL sort(tmp_array, current_step, indexing)
475 4236 : DO i = 1, current_step
476 8472 : IF (ABS(AIMAG(evals(indexing(current_step + 1 - i)))) < EPSILON(${real_zero}$)) THEN
477 4236 : selected_ind(i) = indexing(current_step + 1 - i)
478 4236 : neval = neval + 1
479 4236 : IF (neval == nlimit) EXIT
480 : END IF
481 : END DO
482 :
483 4236 : END SUBROUTINE index_nmax_real_eval_${nametype1}$
484 :
485 : ! **************************************************************************************************
486 : !> \brief ...
487 : !> \param evals ...
488 : !> \param current_step ...
489 : !> \param selected_ind ...
490 : !> \param neval ...
491 : ! **************************************************************************************************
492 4863 : SUBROUTINE index_nmin_real_eval_${nametype1}$ (evals, current_step, selected_ind, neval)
493 : COMPLEX(${type_prec}$), DIMENSION(:) :: evals
494 : INTEGER, INTENT(IN) :: current_step
495 : INTEGER, DIMENSION(:) :: selected_ind
496 : INTEGER :: neval, nlimit
497 :
498 : INTEGER :: i
499 9726 : INTEGER, DIMENSION(current_step) :: indexing
500 9726 : REAL(${type_prec}$), DIMENSION(current_step) :: tmp_array
501 :
502 4863 : nlimit = neval; neval = 0
503 102123 : selected_ind = 0
504 81884 : tmp_array(1:current_step) = REAL(evals(1:current_step), ${type_prec}$)
505 4863 : CALL sort(tmp_array, current_step, indexing)
506 4863 : DO i = 1, current_step
507 9726 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(${real_zero}$)) THEN
508 4863 : selected_ind(i) = indexing(i)
509 4863 : neval = neval + 1
510 4863 : IF (neval == nlimit) EXIT
511 : END IF
512 : END DO
513 :
514 4863 : END SUBROUTINE index_nmin_real_eval_${nametype1}$
515 :
516 : ! **************************************************************************************************
517 : !> \brief ...
518 : !> \param arnoldi_data ...
519 : !> \param matrix ...
520 : !> \param max_iter ...
521 : ! **************************************************************************************************
522 128911 : SUBROUTINE setup_arnoldi_data_${nametype1}$ (arnoldi_data, matrix, max_iter)
523 : TYPE(arnoldi_data_type) :: arnoldi_data
524 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
525 : INTEGER :: max_iter
526 :
527 : INTEGER :: nrow_local
528 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
529 :
530 128911 : ALLOCATE (ar_data)
531 128911 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, nfullrows_local=nrow_local)
532 343037 : ALLOCATE (ar_data%f_vec(nrow_local))
533 214126 : ALLOCATE (ar_data%x_vec(nrow_local))
534 515644 : ALLOCATE (ar_data%Hessenberg(max_iter + 1, max_iter))
535 471948 : ALLOCATE (ar_data%local_history(nrow_local, max_iter))
536 :
537 386733 : ALLOCATE (ar_data%evals(max_iter))
538 515644 : ALLOCATE (ar_data%revec(max_iter, max_iter))
539 :
540 128911 : CALL set_data_${nametype1}$ (arnoldi_data, ar_data)
541 :
542 128911 : END SUBROUTINE setup_arnoldi_data_${nametype1}$
543 :
544 : ! **************************************************************************************************
545 : !> \brief ...
546 : !> \param arnoldi_data ...
547 : ! **************************************************************************************************
548 128911 : SUBROUTINE deallocate_arnoldi_data_${nametype1}$ (arnoldi_data)
549 : TYPE(arnoldi_data_type) :: arnoldi_data
550 :
551 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
552 :
553 128911 : ar_data => get_data_${nametype1}$ (arnoldi_data)
554 128911 : IF (ASSOCIATED(ar_data%f_vec)) DEALLOCATE (ar_data%f_vec)
555 128911 : IF (ASSOCIATED(ar_data%x_vec)) DEALLOCATE (ar_data%x_vec)
556 128911 : IF (ASSOCIATED(ar_data%Hessenberg)) DEALLOCATE (ar_data%Hessenberg)
557 128911 : IF (ASSOCIATED(ar_data%local_history)) DEALLOCATE (ar_data%local_history)
558 128911 : IF (ASSOCIATED(ar_data%evals)) DEALLOCATE (ar_data%evals)
559 128911 : IF (ASSOCIATED(ar_data%revec)) DEALLOCATE (ar_data%revec)
560 128911 : DEALLOCATE (ar_data)
561 :
562 128911 : END SUBROUTINE deallocate_arnoldi_data_${nametype1}$
563 :
564 : ! **************************************************************************************************
565 : !> \brief ...
566 : !> \param arnoldi_data ...
567 : !> \param ind ...
568 : !> \param matrix ...
569 : !> \param vector ...
570 : ! **************************************************************************************************
571 133455 : SUBROUTINE get_selected_ritz_vector_${nametype1}$ (arnoldi_data, ind, matrix, vector)
572 : TYPE(arnoldi_data_type) :: arnoldi_data
573 : INTEGER :: ind
574 : TYPE(dbcsr_type) :: matrix
575 : TYPE(dbcsr_type) :: vector
576 :
577 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
578 : INTEGER :: vsize, myind, sspace_size, i
579 133455 : INTEGER, DIMENSION(:), POINTER :: selected_ind
580 133455 : COMPLEX(${type_prec}$), DIMENSION(:), ALLOCATABLE :: ritz_v
581 133455 : ${type_nametype1}$, DIMENSION(:), POINTER :: data_vec
582 : TYPE(arnoldi_control_type), POINTER :: control
583 :
584 266910 : control => get_control(arnoldi_data)
585 133455 : selected_ind => get_sel_ind(arnoldi_data)
586 133455 : ar_data => get_data_${nametype1}$ (arnoldi_data)
587 133455 : sspace_size = get_subsp_size(arnoldi_data)
588 133455 : vsize = SIZE(ar_data%f_vec)
589 133455 : myind = selected_ind(ind)
590 355018 : ALLOCATE (ritz_v(vsize))
591 961587 : ritz_v = CMPLX(0.0, 0.0, ${type_prec}$)
592 :
593 133455 : CALL dbcsr_release(vector)
594 133455 : CALL create_col_vec_from_matrix(vector, matrix, 1)
595 133455 : IF (control%local_comp) THEN
596 627421 : DO i = 1, sspace_size
597 9066946 : ritz_v(:) = ritz_v(:) + ar_data%local_history(:, i)*ar_data%revec(i, myind)
598 : END DO
599 88108 : data_vec => dbcsr_get_data_p(vector, select_data_type=${nametype_zero}$)
600 : ! is a bit odd but ritz_v is always complex and matrix type determines where it goes
601 : ! again I hope the user knows what is required
602 916240 : data_vec(1:vsize) = ${vartype}$ (ritz_v(1:vsize), KIND=${type_prec}$)
603 : END IF
604 :
605 133455 : DEALLOCATE (ritz_v)
606 :
607 133455 : END SUBROUTINE get_selected_ritz_vector_${nametype1}$
608 :
609 : ! **************************************************************************************************
610 : !> \brief ...
611 : !> \param arnoldi_data ...
612 : !> \param vector ...
613 : ! **************************************************************************************************
614 9572 : SUBROUTINE set_initial_vector_${nametype1}$ (arnoldi_data, vector)
615 : TYPE(arnoldi_data_type) :: arnoldi_data
616 : TYPE(dbcsr_type) :: vector
617 :
618 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
619 4786 : ${type_nametype1}$, DIMENSION(:), POINTER :: data_vec
620 : INTEGER :: nrow_local, ncol_local
621 : TYPE(arnoldi_control_type), POINTER :: control
622 :
623 4786 : control => get_control(arnoldi_data)
624 :
625 4786 : CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
626 4786 : ar_data => get_data_${nametype1}$ (arnoldi_data)
627 4786 : data_vec => dbcsr_get_data_p(vector, select_data_type=${nametype_zero}$)
628 166932 : IF (nrow_local*ncol_local > 0) ar_data%f_vec(1:nrow_local) = data_vec(1:nrow_local)
629 :
630 4786 : END SUBROUTINE set_initial_vector_${nametype1}$
631 :
632 : #:endfor
633 :
634 : END MODULE arnoldi_data_methods
635 :
|
