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 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 method
11 : !> \par History
12 : !> 2014.09 created [Florian Schiffmann]
13 : !> 2023.12 Removed support for single-precision [Ole Schuett]
14 : !> 2024.12 Removed support for complex input matrices [Ole Schuett]
15 : !> \author Florian Schiffmann
16 : ! **************************************************************************************************
17 : MODULE arnoldi_data_methods
18 : USE arnoldi_types, ONLY: &
19 : arnoldi_control_type, arnoldi_data_type, arnoldi_env_type, get_control, get_data, &
20 : get_evals, get_sel_ind, set_control, set_data
21 : USE arnoldi_vector, ONLY: create_col_vec_from_matrix
22 : USE cp_dbcsr_api, ONLY: &
23 : dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_data_p, dbcsr_get_info, &
24 : dbcsr_get_matrix_type, dbcsr_mp_grid_setup, dbcsr_p_type, dbcsr_release, dbcsr_type, &
25 : dbcsr_type_symmetric
26 : USE kinds, ONLY: dp
27 : USE util, ONLY: sort
28 : #include "../base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 :
34 : PUBLIC :: select_evals, get_selected_ritz_val, arnoldi_is_converged, &
35 : arnoldi_env_type, get_nrestart, set_arnoldi_initial_vector, &
36 : setup_arnoldi_env, deallocate_arnoldi_env, get_selected_ritz_vector
37 :
38 : CONTAINS
39 :
40 : ! **************************************************************************************************
41 : !> \brief This routine sets the environment for the arnoldi iteration and
42 : !> the krylov subspace creation. All simulation parameters have to be given
43 : !> at this stage so the rest can run fully automated
44 : !> In addition, this routine allocates the data necessary for
45 : !> \param arnoldi_env this type which gets filled with information and on output contains all
46 : !> information necessary to extract whatever the user desires
47 : !> \param matrix vector of matrices, only the first gets used to get some dimensions
48 : !> and parallel information needed later on
49 : !> \param max_iter maximum dimension of the krylov subspace
50 : !> \param threshold convergence threshold, this is used for both subspace and eigenval
51 : !> \param selection_crit integer defining according to which criterion the
52 : !> eigenvalues are selected for the subspace
53 : !> \param nval_request for some sel_crit useful, how many eV to select
54 : !> \param nrestarts ...
55 : !> \param generalized_ev ...
56 : !> \param iram ...
57 : ! **************************************************************************************************
58 134339 : SUBROUTINE setup_arnoldi_env(arnoldi_env, matrix, max_iter, threshold, selection_crit, &
59 : nval_request, nrestarts, generalized_ev, iram)
60 : TYPE(arnoldi_env_type) :: arnoldi_env
61 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
62 : INTEGER :: max_iter
63 : REAL(dp) :: threshold
64 : INTEGER :: selection_crit, nval_request, nrestarts
65 : LOGICAL :: generalized_ev, iram
66 :
67 : CALL setup_arnoldi_control(arnoldi_env, matrix, max_iter, threshold, selection_crit, &
68 134339 : nval_request, nrestarts, generalized_ev, iram)
69 :
70 134339 : CALL setup_arnoldi_data(arnoldi_env, matrix, max_iter)
71 :
72 134339 : END SUBROUTINE setup_arnoldi_env
73 :
74 : ! **************************************************************************************************
75 : !> \brief Creates the data type for arnoldi, see above for details
76 : !> \param arnoldi_env ...
77 : !> \param matrix ...
78 : !> \param max_iter ...
79 : ! **************************************************************************************************
80 134339 : SUBROUTINE setup_arnoldi_data(arnoldi_env, matrix, max_iter)
81 : TYPE(arnoldi_env_type) :: arnoldi_env
82 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
83 : INTEGER :: max_iter
84 :
85 : INTEGER :: nrow_local
86 : TYPE(arnoldi_data_type), POINTER :: ar_data
87 :
88 134339 : ALLOCATE (ar_data)
89 134339 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, nfullrows_local=nrow_local)
90 356799 : ALLOCATE (ar_data%f_vec(nrow_local))
91 222460 : ALLOCATE (ar_data%x_vec(nrow_local))
92 537356 : ALLOCATE (ar_data%Hessenberg(max_iter + 1, max_iter))
93 491138 : ALLOCATE (ar_data%local_history(nrow_local, max_iter))
94 :
95 403017 : ALLOCATE (ar_data%evals(max_iter))
96 537356 : ALLOCATE (ar_data%revec(max_iter, max_iter))
97 :
98 134339 : CALL set_data(arnoldi_env, ar_data)
99 :
100 134339 : END SUBROUTINE setup_arnoldi_data
101 :
102 : ! **************************************************************************************************
103 : !> \brief Creates the control type for arnoldi, see above for details
104 : !> \param arnoldi_env ...
105 : !> \param matrix ...
106 : !> \param max_iter ...
107 : !> \param threshold ...
108 : !> \param selection_crit ...
109 : !> \param nval_request ...
110 : !> \param nrestarts ...
111 : !> \param generalized_ev ...
112 : !> \param iram ...
113 : ! **************************************************************************************************
114 134339 : SUBROUTINE setup_arnoldi_control(arnoldi_env, matrix, max_iter, threshold, selection_crit, &
115 : nval_request, nrestarts, generalized_ev, iram)
116 : TYPE(arnoldi_env_type) :: arnoldi_env
117 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
118 : INTEGER :: max_iter
119 : REAL(dp) :: threshold
120 : INTEGER :: selection_crit, nval_request, nrestarts
121 : LOGICAL :: generalized_ev, iram
122 :
123 : INTEGER :: group_handle, pcol_handle
124 : LOGICAL :: subgroups_defined
125 : TYPE(arnoldi_control_type), POINTER :: control
126 : TYPE(dbcsr_distribution_type) :: distri
127 :
128 0 : ALLOCATE (control)
129 : ! Fill the information which will later control the arnoldi method and allow synchronization.
130 134339 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, distribution=distri)
131 134339 : CALL dbcsr_mp_grid_setup(distri)
132 : CALL dbcsr_distribution_get(distri, &
133 : group=group_handle, &
134 : mynode=control%myproc, &
135 : subgroups_defined=subgroups_defined, &
136 134339 : pcol_group=pcol_handle)
137 :
138 134339 : CALL control%mp_group%set_handle(group_handle)
139 134339 : CALL control%pcol_group%set_handle(pcol_handle)
140 :
141 134339 : IF (.NOT. subgroups_defined) &
142 0 : CPABORT("arnoldi only with subgroups")
143 :
144 134339 : control%symmetric = .FALSE.
145 : ! Will need a fix for complex because there it has to be hermitian
146 134339 : IF (SIZE(matrix) == 1) &
147 130496 : control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
148 :
149 : ! Set the control parameters
150 134339 : control%max_iter = max_iter
151 134339 : control%current_step = 0
152 134339 : control%selection_crit = selection_crit
153 134339 : control%nval_req = nval_request
154 134339 : control%threshold = threshold
155 134339 : control%converged = .FALSE.
156 134339 : control%has_initial_vector = .FALSE.
157 134339 : control%iram = iram
158 134339 : control%nrestart = nrestarts
159 134339 : control%generalized_ev = generalized_ev
160 :
161 134339 : IF (control%nval_req > 1 .AND. control%nrestart > 0 .AND. .NOT. control%iram) &
162 : CALL cp_abort(__LOCATION__, 'with more than one eigenvalue requested '// &
163 0 : 'internal restarting with a previous EVEC is a bad idea, set IRAM or nrestsart=0')
164 :
165 : ! some checks for the generalized EV mode
166 134339 : IF (control%generalized_ev .AND. selection_crit == 1) &
167 : CALL cp_abort(__LOCATION__, &
168 0 : 'generalized ev can only highest OR lowest EV')
169 134339 : IF (control%generalized_ev .AND. nval_request .NE. 1) &
170 : CALL cp_abort(__LOCATION__, &
171 0 : 'generalized ev can only compute one EV at the time')
172 134339 : IF (control%generalized_ev .AND. control%nrestart == 0) &
173 : CALL cp_abort(__LOCATION__, &
174 0 : 'outer loops are mandatory for generalized EV, set nrestart appropriatly')
175 134339 : IF (SIZE(matrix) .NE. 2 .AND. control%generalized_ev) &
176 : CALL cp_abort(__LOCATION__, &
177 0 : 'generalized ev needs exactly two matrices as input (2nd is the metric)')
178 :
179 403017 : ALLOCATE (control%selected_ind(max_iter))
180 134339 : CALL set_control(arnoldi_env, control)
181 :
182 268678 : END SUBROUTINE setup_arnoldi_control
183 :
184 : ! **************************************************************************************************
185 : !> \brief ...
186 : !> \param arnoldi_env ...
187 : !> \param ind ...
188 : !> \param matrix ...
189 : !> \param vector ...
190 : ! **************************************************************************************************
191 138884 : SUBROUTINE get_selected_ritz_vector(arnoldi_env, ind, matrix, vector)
192 : TYPE(arnoldi_env_type) :: arnoldi_env
193 : INTEGER :: ind
194 : TYPE(dbcsr_type) :: matrix, vector
195 :
196 138884 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: ritz_v
197 : INTEGER :: i, myind, sspace_size, vsize
198 138884 : INTEGER, DIMENSION(:), POINTER :: selected_ind
199 138884 : REAL(kind=dp), DIMENSION(:), POINTER :: data_vec
200 : TYPE(arnoldi_control_type), POINTER :: control
201 : TYPE(arnoldi_data_type), POINTER :: ar_data
202 :
203 277768 : control => get_control(arnoldi_env)
204 138884 : selected_ind => get_sel_ind(arnoldi_env)
205 138884 : ar_data => get_data(arnoldi_env)
206 138884 : sspace_size = get_subsp_size(arnoldi_env)
207 138884 : vsize = SIZE(ar_data%f_vec)
208 138884 : myind = selected_ind(ind)
209 368779 : ALLOCATE (ritz_v(vsize))
210 991312 : ritz_v = CMPLX(0.0, 0.0, dp)
211 :
212 138884 : CALL dbcsr_release(vector)
213 138884 : CALL create_col_vec_from_matrix(vector, matrix, 1)
214 138884 : IF (control%local_comp) THEN
215 651211 : DO i = 1, sspace_size
216 9560777 : ritz_v(:) = ritz_v(:) + ar_data%local_history(:, i)*ar_data%revec(i, myind)
217 : END DO
218 91011 : data_vec => dbcsr_get_data_p(vector)
219 : ! is a bit odd but ritz_v is always complex and matrix type determines where it goes
220 : ! again I hope the user knows what is required
221 943439 : data_vec(1:vsize) = REAL(ritz_v(1:vsize), KIND=dp)
222 : END IF
223 :
224 138884 : DEALLOCATE (ritz_v)
225 :
226 138884 : END SUBROUTINE get_selected_ritz_vector
227 :
228 : ! **************************************************************************************************
229 : !> \brief Deallocate the data in arnoldi_env
230 : !> \param arnoldi_env ...
231 : ! **************************************************************************************************
232 134339 : SUBROUTINE deallocate_arnoldi_env(arnoldi_env)
233 : TYPE(arnoldi_env_type) :: arnoldi_env
234 :
235 : TYPE(arnoldi_control_type), POINTER :: control
236 : TYPE(arnoldi_data_type), POINTER :: ar_data
237 :
238 134339 : ar_data => get_data(arnoldi_env)
239 134339 : IF (ASSOCIATED(ar_data%f_vec)) DEALLOCATE (ar_data%f_vec)
240 134339 : IF (ASSOCIATED(ar_data%x_vec)) DEALLOCATE (ar_data%x_vec)
241 134339 : IF (ASSOCIATED(ar_data%Hessenberg)) DEALLOCATE (ar_data%Hessenberg)
242 134339 : IF (ASSOCIATED(ar_data%local_history)) DEALLOCATE (ar_data%local_history)
243 134339 : IF (ASSOCIATED(ar_data%evals)) DEALLOCATE (ar_data%evals)
244 134339 : IF (ASSOCIATED(ar_data%revec)) DEALLOCATE (ar_data%revec)
245 134339 : DEALLOCATE (ar_data)
246 :
247 134339 : control => get_control(arnoldi_env)
248 134339 : DEALLOCATE (control%selected_ind)
249 134339 : DEALLOCATE (control)
250 :
251 134339 : END SUBROUTINE deallocate_arnoldi_env
252 :
253 : ! **************************************************************************************************
254 : !> \brief perform the selection of eigenvalues, fills the selected_ind array
255 : !> \param arnoldi_env ...
256 : ! **************************************************************************************************
257 135309 : SUBROUTINE select_evals(arnoldi_env)
258 : TYPE(arnoldi_env_type) :: arnoldi_env
259 :
260 : INTEGER :: i, last_el, my_crit, my_ind
261 : REAL(dp) :: convergence
262 : TYPE(arnoldi_control_type), POINTER :: control
263 : TYPE(arnoldi_data_type), POINTER :: ar_data
264 :
265 135309 : control => get_control(arnoldi_env)
266 135309 : ar_data => get_data(arnoldi_env)
267 :
268 135309 : last_el = control%current_step
269 135309 : convergence = REAL(0.0, dp)
270 135309 : my_crit = control%selection_crit
271 135309 : control%nval_out = MIN(control%nval_req, control%current_step)
272 126321 : SELECT CASE (my_crit)
273 : ! minimum and maximum real eval
274 : CASE (1)
275 126321 : CALL index_min_max_real_eval(ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
276 : ! n maximum real eval
277 : CASE (2)
278 4099 : CALL index_nmax_real_eval(ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
279 : ! n minimum real eval
280 : CASE (3)
281 4889 : CALL index_nmin_real_eval(ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
282 : CASE DEFAULT
283 135309 : CPABORT("unknown selection index")
284 : END SELECT
285 : ! test whether we are converged
286 396939 : DO i = 1, control%nval_out
287 261630 : my_ind = control%selected_ind(i)
288 : convergence = MAX(convergence, &
289 396939 : ABS(ar_data%revec(last_el, my_ind)*ar_data%Hessenberg(last_el + 1, last_el)))
290 : END DO
291 135309 : control%converged = convergence .LT. control%threshold
292 :
293 135309 : END SUBROUTINE select_evals
294 :
295 : ! **************************************************************************************************
296 : !> \brief set a new selection type, if you notice you didn't like the initial one
297 : !> \param arnoldi_env ...
298 : !> \param itype ...
299 : ! **************************************************************************************************
300 0 : SUBROUTINE set_eval_selection(arnoldi_env, itype)
301 : TYPE(arnoldi_env_type) :: arnoldi_env
302 : INTEGER :: itype
303 :
304 : TYPE(arnoldi_control_type), POINTER :: control
305 :
306 0 : control => get_control(arnoldi_env)
307 0 : control%selection_crit = itype
308 0 : END SUBROUTINE set_eval_selection
309 :
310 : ! **************************************************************************************************
311 : !> \brief returns the number of restarts allowed for arnoldi
312 : !> \param arnoldi_env ...
313 : !> \return ...
314 : ! **************************************************************************************************
315 134141 : FUNCTION get_nrestart(arnoldi_env) RESULT(nrestart)
316 : TYPE(arnoldi_env_type) :: arnoldi_env
317 : INTEGER :: nrestart
318 :
319 : TYPE(arnoldi_control_type), POINTER :: control
320 :
321 134141 : control => get_control(arnoldi_env)
322 134141 : nrestart = control%nrestart
323 :
324 134141 : END FUNCTION get_nrestart
325 :
326 : ! **************************************************************************************************
327 : !> \brief get the number of eigenvalues matching the search criterion
328 : !> \param arnoldi_env ...
329 : !> \return ...
330 : ! **************************************************************************************************
331 260212 : FUNCTION get_nval_out(arnoldi_env) RESULT(nval_out)
332 : TYPE(arnoldi_env_type) :: arnoldi_env
333 : INTEGER :: nval_out
334 :
335 : TYPE(arnoldi_control_type), POINTER :: control
336 :
337 260212 : control => get_control(arnoldi_env)
338 260212 : nval_out = control%nval_out
339 :
340 260212 : END FUNCTION get_nval_out
341 :
342 : ! **************************************************************************************************
343 : !> \brief get dimension of the krylov space. Can be less than max_iter if subspace converged early
344 : !> \param arnoldi_env ...
345 : !> \return ...
346 : ! **************************************************************************************************
347 138884 : FUNCTION get_subsp_size(arnoldi_env) RESULT(current_step)
348 : TYPE(arnoldi_env_type) :: arnoldi_env
349 : INTEGER :: current_step
350 :
351 : TYPE(arnoldi_control_type), POINTER :: control
352 :
353 138884 : control => get_control(arnoldi_env)
354 138884 : current_step = control%current_step
355 :
356 138884 : END FUNCTION get_subsp_size
357 :
358 : ! **************************************************************************************************
359 : !> \brief Find out whether the method with the current search criterion is converged
360 : !> \param arnoldi_env ...
361 : !> \return ...
362 : ! **************************************************************************************************
363 126269 : FUNCTION arnoldi_is_converged(arnoldi_env) RESULT(converged)
364 : TYPE(arnoldi_env_type) :: arnoldi_env
365 : LOGICAL :: converged
366 :
367 : TYPE(arnoldi_control_type), POINTER :: control
368 :
369 126269 : control => get_control(arnoldi_env)
370 126269 : converged = control%converged
371 :
372 126269 : END FUNCTION arnoldi_is_converged
373 :
374 : ! **************************************************************************************************
375 : !> \brief get a single specific Ritz value from the set of selected
376 : !> \param arnoldi_env ...
377 : !> \param ind ...
378 : !> \return ...
379 : ! **************************************************************************************************
380 260212 : FUNCTION get_selected_ritz_val(arnoldi_env, ind) RESULT(eval_out)
381 : TYPE(arnoldi_env_type) :: arnoldi_env
382 : INTEGER :: ind
383 : COMPLEX(dp) :: eval_out
384 :
385 260212 : COMPLEX(dp), DIMENSION(:), POINTER :: evals
386 : INTEGER :: ev_ind
387 260212 : INTEGER, DIMENSION(:), POINTER :: selected_ind
388 :
389 260212 : IF (ind .GT. get_nval_out(arnoldi_env)) &
390 0 : CPABORT('outside range of indexed evals')
391 :
392 260212 : selected_ind => get_sel_ind(arnoldi_env)
393 260212 : ev_ind = selected_ind(ind)
394 260212 : evals => get_evals(arnoldi_env)
395 260212 : eval_out = evals(ev_ind)
396 :
397 260212 : END FUNCTION get_selected_ritz_val
398 :
399 : ! **************************************************************************************************
400 : !> \brief Get all Ritz values of the selected set. eval_out has to be allocated
401 : !> at least the size of get_neval_out()
402 : !> \param arnoldi_env ...
403 : !> \param eval_out ...
404 : ! **************************************************************************************************
405 0 : SUBROUTINE get_all_selected_ritz_val(arnoldi_env, eval_out)
406 : TYPE(arnoldi_env_type) :: arnoldi_env
407 : COMPLEX(dp), DIMENSION(:) :: eval_out
408 :
409 0 : COMPLEX(dp), DIMENSION(:), POINTER :: evals
410 : INTEGER :: ev_ind, ind
411 0 : INTEGER, DIMENSION(:), POINTER :: selected_ind
412 :
413 0 : NULLIFY (evals)
414 0 : IF (SIZE(eval_out) .LT. get_nval_out(arnoldi_env)) &
415 0 : CPABORT('array for eval output too small')
416 0 : selected_ind => get_sel_ind(arnoldi_env)
417 :
418 0 : evals => get_evals(arnoldi_env)
419 :
420 0 : DO ind = 1, get_nval_out(arnoldi_env)
421 0 : ev_ind = selected_ind(ind)
422 0 : eval_out(ind) = evals(ev_ind)
423 : END DO
424 :
425 0 : END SUBROUTINE get_all_selected_ritz_val
426 :
427 : ! **************************************************************************************************
428 : !> \brief ...
429 : !> \param arnoldi_env ...
430 : !> \param vector ...
431 : ! **************************************************************************************************
432 9344 : SUBROUTINE set_arnoldi_initial_vector(arnoldi_env, vector)
433 : TYPE(arnoldi_env_type) :: arnoldi_env
434 : TYPE(dbcsr_type) :: vector
435 :
436 : INTEGER :: ncol_local, nrow_local
437 4672 : REAL(kind=dp), DIMENSION(:), POINTER :: data_vec
438 : TYPE(arnoldi_control_type), POINTER :: control
439 : TYPE(arnoldi_data_type), POINTER :: ar_data
440 :
441 9344 : control => get_control(arnoldi_env)
442 4672 : control%has_initial_vector = .TRUE.
443 4672 : ar_data => get_data(arnoldi_env)
444 :
445 4672 : CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
446 4672 : data_vec => dbcsr_get_data_p(vector)
447 159977 : IF (nrow_local*ncol_local > 0) ar_data%f_vec(1:nrow_local) = data_vec(1:nrow_local)
448 :
449 4672 : END SUBROUTINE set_arnoldi_initial_vector
450 :
451 : !!! Here come the methods handling the selection of eigenvalues and eigenvectors !!!
452 : !!! If you want a personal method, simply created a Subroutine returning the index
453 : !!! array selected ind which contains as the first nval_out entries the index of the evals
454 :
455 : ! **************************************************************************************************
456 : !> \brief ...
457 : !> \param evals ...
458 : !> \param current_step ...
459 : !> \param selected_ind ...
460 : !> \param neval ...
461 : ! **************************************************************************************************
462 126321 : SUBROUTINE index_min_max_real_eval(evals, current_step, selected_ind, neval)
463 : COMPLEX(dp), DIMENSION(:) :: evals
464 : INTEGER, INTENT(IN) :: current_step
465 : INTEGER, DIMENSION(:) :: selected_ind
466 : INTEGER :: neval
467 :
468 : INTEGER :: i
469 252642 : INTEGER, DIMENSION(current_step) :: indexing
470 252642 : REAL(dp), DIMENSION(current_step) :: tmp_array
471 :
472 126321 : neval = 0
473 4182801 : selected_ind = 0
474 731068 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
475 126321 : CALL sort(tmp_array, current_step, indexing)
476 126321 : DO i = 1, current_step
477 126321 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
478 126321 : selected_ind(1) = indexing(i)
479 126321 : neval = neval + 1
480 126321 : EXIT
481 : END IF
482 : END DO
483 126321 : DO i = current_step, 1, -1
484 126321 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
485 126321 : selected_ind(2) = indexing(i)
486 126321 : neval = neval + 1
487 126321 : EXIT
488 : END IF
489 : END DO
490 :
491 126321 : END SUBROUTINE index_min_max_real_eval
492 :
493 : ! **************************************************************************************************
494 : !> \brief ...
495 : !> \param evals ...
496 : !> \param current_step ...
497 : !> \param selected_ind ...
498 : !> \param neval ...
499 : ! **************************************************************************************************
500 4099 : SUBROUTINE index_nmax_real_eval(evals, current_step, selected_ind, neval)
501 : COMPLEX(dp), DIMENSION(:) :: evals
502 : INTEGER, INTENT(IN) :: current_step
503 : INTEGER, DIMENSION(:) :: selected_ind
504 : INTEGER :: neval
505 :
506 : INTEGER :: i, nlimit
507 8198 : INTEGER, DIMENSION(current_step) :: indexing
508 8198 : REAL(dp), DIMENSION(current_step) :: tmp_array
509 :
510 4099 : nlimit = neval; neval = 0
511 86079 : selected_ind = 0
512 21906 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
513 4099 : CALL sort(tmp_array, current_step, indexing)
514 4099 : DO i = 1, current_step
515 8198 : IF (ABS(AIMAG(evals(indexing(current_step + 1 - i)))) < EPSILON(0.0_dp)) THEN
516 4099 : selected_ind(i) = indexing(current_step + 1 - i)
517 4099 : neval = neval + 1
518 4099 : IF (neval == nlimit) EXIT
519 : END IF
520 : END DO
521 :
522 4099 : END SUBROUTINE index_nmax_real_eval
523 :
524 : ! **************************************************************************************************
525 : !> \brief ...
526 : !> \param evals ...
527 : !> \param current_step ...
528 : !> \param selected_ind ...
529 : !> \param neval ...
530 : ! **************************************************************************************************
531 4889 : SUBROUTINE index_nmin_real_eval(evals, current_step, selected_ind, neval)
532 : COMPLEX(dp), DIMENSION(:) :: evals
533 : INTEGER, INTENT(IN) :: current_step
534 : INTEGER, DIMENSION(:) :: selected_ind
535 : INTEGER :: neval
536 :
537 : INTEGER :: i, nlimit
538 9778 : INTEGER, DIMENSION(current_step) :: indexing
539 9778 : REAL(dp), DIMENSION(current_step) :: tmp_array
540 :
541 4889 : nlimit = neval; neval = 0
542 102669 : selected_ind = 0
543 82714 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
544 4889 : CALL sort(tmp_array, current_step, indexing)
545 4889 : DO i = 1, current_step
546 9778 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
547 4889 : selected_ind(i) = indexing(i)
548 4889 : neval = neval + 1
549 4889 : IF (neval == nlimit) EXIT
550 : END IF
551 : END DO
552 :
553 4889 : END SUBROUTINE index_nmin_real_eval
554 :
555 : END MODULE arnoldi_data_methods
|
