Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 145343 : 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 145343 : nval_request, nrestarts, generalized_ev, iram)
69 :
70 145343 : CALL setup_arnoldi_data(arnoldi_env, matrix, max_iter)
71 :
72 145343 : 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 145343 : 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 145343 : ALLOCATE (ar_data)
89 145343 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, nfullrows_local=nrow_local)
90 384711 : ALLOCATE (ar_data%f_vec(nrow_local))
91 239368 : ALLOCATE (ar_data%x_vec(nrow_local))
92 581372 : ALLOCATE (ar_data%Hessenberg(max_iter + 1, max_iter))
93 530054 : ALLOCATE (ar_data%local_history(nrow_local, max_iter))
94 :
95 436029 : ALLOCATE (ar_data%evals(max_iter))
96 581372 : ALLOCATE (ar_data%revec(max_iter, max_iter))
97 :
98 145343 : CALL set_data(arnoldi_env, ar_data)
99 :
100 145343 : 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 145343 : 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 145343 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, distribution=distri)
131 145343 : 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 145343 : pcol_group=pcol_handle)
137 :
138 145343 : CALL control%mp_group%set_handle(group_handle)
139 145343 : CALL control%pcol_group%set_handle(pcol_handle)
140 :
141 145343 : IF (.NOT. subgroups_defined) THEN
142 0 : CPABORT("arnoldi only with subgroups")
143 : END IF
144 :
145 145343 : control%symmetric = .FALSE.
146 : ! Will need a fix for complex because there it has to be hermitian
147 145343 : IF (SIZE(matrix) == 1) THEN
148 141114 : control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
149 : END IF
150 :
151 : ! Set the control parameters
152 145343 : control%max_iter = max_iter
153 145343 : control%current_step = 0
154 145343 : control%selection_crit = selection_crit
155 145343 : control%nval_req = nval_request
156 145343 : control%threshold = threshold
157 145343 : control%converged = .FALSE.
158 145343 : control%has_initial_vector = .FALSE.
159 145343 : control%iram = iram
160 145343 : control%nrestart = nrestarts
161 145343 : control%generalized_ev = generalized_ev
162 :
163 145343 : IF (control%nval_req > 1 .AND. control%nrestart > 0 .AND. .NOT. control%iram) THEN
164 : CALL cp_abort(__LOCATION__, 'with more than one eigenvalue requested '// &
165 0 : 'internal restarting with a previous EVEC is a bad idea, set IRAM or nrestsart=0')
166 : END IF
167 :
168 : ! some checks for the generalized EV mode
169 145343 : IF (control%generalized_ev .AND. selection_crit == 1) THEN
170 : CALL cp_abort(__LOCATION__, &
171 0 : 'generalized ev can only highest OR lowest EV')
172 : END IF
173 145343 : IF (control%generalized_ev .AND. nval_request /= 1) THEN
174 : CALL cp_abort(__LOCATION__, &
175 0 : 'generalized ev can only compute one EV at the time')
176 : END IF
177 145343 : IF (control%generalized_ev .AND. control%nrestart == 0) THEN
178 : CALL cp_abort(__LOCATION__, &
179 0 : 'outer loops are mandatory for generalized EV, set nrestart appropriatly')
180 : END IF
181 145343 : IF (SIZE(matrix) /= 2 .AND. control%generalized_ev) THEN
182 : CALL cp_abort(__LOCATION__, &
183 0 : 'generalized ev needs exactly two matrices as input (2nd is the metric)')
184 : END IF
185 :
186 436029 : ALLOCATE (control%selected_ind(max_iter))
187 145343 : CALL set_control(arnoldi_env, control)
188 :
189 290686 : END SUBROUTINE setup_arnoldi_control
190 :
191 : ! **************************************************************************************************
192 : !> \brief ...
193 : !> \param arnoldi_env ...
194 : !> \param ind ...
195 : !> \param matrix ...
196 : !> \param vector ...
197 : ! **************************************************************************************************
198 150236 : SUBROUTINE get_selected_ritz_vector(arnoldi_env, ind, matrix, vector)
199 : TYPE(arnoldi_env_type) :: arnoldi_env
200 : INTEGER :: ind
201 : TYPE(dbcsr_type) :: matrix, vector
202 :
203 150236 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: ritz_v
204 : INTEGER :: i, myind, sspace_size, vsize
205 150236 : INTEGER, DIMENSION(:), POINTER :: selected_ind
206 150236 : REAL(kind=dp), DIMENSION(:), POINTER :: data_vec
207 : TYPE(arnoldi_control_type), POINTER :: control
208 : TYPE(arnoldi_data_type), POINTER :: ar_data
209 :
210 300472 : control => get_control(arnoldi_env)
211 150236 : selected_ind => get_sel_ind(arnoldi_env)
212 150236 : ar_data => get_data(arnoldi_env)
213 150236 : sspace_size = get_subsp_size(arnoldi_env)
214 150236 : vsize = SIZE(ar_data%f_vec)
215 150236 : myind = selected_ind(ind)
216 397531 : ALLOCATE (ritz_v(vsize))
217 150236 : ritz_v = CMPLX(0.0, 0.0, dp)
218 :
219 150236 : CALL dbcsr_release(vector)
220 150236 : CALL create_col_vec_from_matrix(vector, matrix, 1)
221 150236 : IF (control%local_comp) THEN
222 682938 : DO i = 1, sspace_size
223 9648258 : ritz_v(:) = ritz_v(:) + ar_data%local_history(:, i)*ar_data%revec(i, myind)
224 : END DO
225 97059 : data_vec => dbcsr_get_data_p(vector)
226 : ! is a bit odd but ritz_v is always complex and matrix type determines where it goes
227 : ! again I hope the user knows what is required
228 979864 : data_vec(1:vsize) = REAL(ritz_v(1:vsize), KIND=dp)
229 : END IF
230 :
231 150236 : DEALLOCATE (ritz_v)
232 :
233 150236 : END SUBROUTINE get_selected_ritz_vector
234 :
235 : ! **************************************************************************************************
236 : !> \brief Deallocate the data in arnoldi_env
237 : !> \param arnoldi_env ...
238 : ! **************************************************************************************************
239 145343 : SUBROUTINE deallocate_arnoldi_env(arnoldi_env)
240 : TYPE(arnoldi_env_type) :: arnoldi_env
241 :
242 : TYPE(arnoldi_control_type), POINTER :: control
243 : TYPE(arnoldi_data_type), POINTER :: ar_data
244 :
245 145343 : ar_data => get_data(arnoldi_env)
246 145343 : IF (ASSOCIATED(ar_data%f_vec)) DEALLOCATE (ar_data%f_vec)
247 145343 : IF (ASSOCIATED(ar_data%x_vec)) DEALLOCATE (ar_data%x_vec)
248 145343 : IF (ASSOCIATED(ar_data%Hessenberg)) DEALLOCATE (ar_data%Hessenberg)
249 145343 : IF (ASSOCIATED(ar_data%local_history)) DEALLOCATE (ar_data%local_history)
250 145343 : IF (ASSOCIATED(ar_data%evals)) DEALLOCATE (ar_data%evals)
251 145343 : IF (ASSOCIATED(ar_data%revec)) DEALLOCATE (ar_data%revec)
252 145343 : DEALLOCATE (ar_data)
253 :
254 145343 : control => get_control(arnoldi_env)
255 145343 : DEALLOCATE (control%selected_ind)
256 145343 : DEALLOCATE (control)
257 :
258 145343 : END SUBROUTINE deallocate_arnoldi_env
259 :
260 : ! **************************************************************************************************
261 : !> \brief perform the selection of eigenvalues, fills the selected_ind array
262 : !> \param arnoldi_env ...
263 : ! **************************************************************************************************
264 146371 : SUBROUTINE select_evals(arnoldi_env)
265 : TYPE(arnoldi_env_type) :: arnoldi_env
266 :
267 : INTEGER :: i, last_el, my_crit, my_ind
268 : REAL(dp) :: convergence
269 : TYPE(arnoldi_control_type), POINTER :: control
270 : TYPE(arnoldi_data_type), POINTER :: ar_data
271 :
272 146371 : control => get_control(arnoldi_env)
273 146371 : ar_data => get_data(arnoldi_env)
274 :
275 146371 : last_el = control%current_step
276 146371 : convergence = REAL(0.0, dp)
277 146371 : my_crit = control%selection_crit
278 146371 : control%nval_out = MIN(control%nval_req, control%current_step)
279 136583 : SELECT CASE (my_crit)
280 : ! minimum and maximum real eval
281 : CASE (1)
282 136583 : CALL index_min_max_real_eval(ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
283 : ! n maximum real eval
284 : CASE (2)
285 4457 : CALL index_nmax_real_eval(ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
286 : ! n minimum real eval
287 : CASE (3)
288 5331 : CALL index_nmin_real_eval(ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
289 : CASE DEFAULT
290 146371 : CPABORT("unknown selection index")
291 : END SELECT
292 : ! test whether we are converged
293 429325 : DO i = 1, control%nval_out
294 282954 : my_ind = control%selected_ind(i)
295 : convergence = MAX(convergence, &
296 429325 : ABS(ar_data%revec(last_el, my_ind)*ar_data%Hessenberg(last_el + 1, last_el)))
297 : END DO
298 146371 : control%converged = convergence < control%threshold
299 :
300 146371 : END SUBROUTINE select_evals
301 :
302 : ! **************************************************************************************************
303 : !> \brief set a new selection type, if you notice you didn't like the initial one
304 : !> \param arnoldi_env ...
305 : !> \param itype ...
306 : ! **************************************************************************************************
307 0 : SUBROUTINE set_eval_selection(arnoldi_env, itype)
308 : TYPE(arnoldi_env_type) :: arnoldi_env
309 : INTEGER :: itype
310 :
311 : TYPE(arnoldi_control_type), POINTER :: control
312 :
313 0 : control => get_control(arnoldi_env)
314 0 : control%selection_crit = itype
315 0 : END SUBROUTINE set_eval_selection
316 :
317 : ! **************************************************************************************************
318 : !> \brief returns the number of restarts allowed for arnoldi
319 : !> \param arnoldi_env ...
320 : !> \return ...
321 : ! **************************************************************************************************
322 145125 : FUNCTION get_nrestart(arnoldi_env) RESULT(nrestart)
323 : TYPE(arnoldi_env_type) :: arnoldi_env
324 : INTEGER :: nrestart
325 :
326 : TYPE(arnoldi_control_type), POINTER :: control
327 :
328 145125 : control => get_control(arnoldi_env)
329 145125 : nrestart = control%nrestart
330 :
331 145125 : END FUNCTION get_nrestart
332 :
333 : ! **************************************************************************************************
334 : !> \brief get the number of eigenvalues matching the search criterion
335 : !> \param arnoldi_env ...
336 : !> \return ...
337 : ! **************************************************************************************************
338 281448 : FUNCTION get_nval_out(arnoldi_env) RESULT(nval_out)
339 : TYPE(arnoldi_env_type) :: arnoldi_env
340 : INTEGER :: nval_out
341 :
342 : TYPE(arnoldi_control_type), POINTER :: control
343 :
344 281448 : control => get_control(arnoldi_env)
345 281448 : nval_out = control%nval_out
346 :
347 281448 : END FUNCTION get_nval_out
348 :
349 : ! **************************************************************************************************
350 : !> \brief get dimension of the krylov space. Can be less than max_iter if subspace converged early
351 : !> \param arnoldi_env ...
352 : !> \return ...
353 : ! **************************************************************************************************
354 150236 : FUNCTION get_subsp_size(arnoldi_env) RESULT(current_step)
355 : TYPE(arnoldi_env_type) :: arnoldi_env
356 : INTEGER :: current_step
357 :
358 : TYPE(arnoldi_control_type), POINTER :: control
359 :
360 150236 : control => get_control(arnoldi_env)
361 150236 : current_step = control%current_step
362 :
363 150236 : END FUNCTION get_subsp_size
364 :
365 : ! **************************************************************************************************
366 : !> \brief Find out whether the method with the current search criterion is converged
367 : !> \param arnoldi_env ...
368 : !> \return ...
369 : ! **************************************************************************************************
370 136541 : FUNCTION arnoldi_is_converged(arnoldi_env) RESULT(converged)
371 : TYPE(arnoldi_env_type) :: arnoldi_env
372 : LOGICAL :: converged
373 :
374 : TYPE(arnoldi_control_type), POINTER :: control
375 :
376 136541 : control => get_control(arnoldi_env)
377 136541 : converged = control%converged
378 :
379 136541 : END FUNCTION arnoldi_is_converged
380 :
381 : ! **************************************************************************************************
382 : !> \brief get a single specific Ritz value from the set of selected
383 : !> \param arnoldi_env ...
384 : !> \param ind ...
385 : !> \return ...
386 : ! **************************************************************************************************
387 281448 : FUNCTION get_selected_ritz_val(arnoldi_env, ind) RESULT(eval_out)
388 : TYPE(arnoldi_env_type) :: arnoldi_env
389 : INTEGER :: ind
390 : COMPLEX(dp) :: eval_out
391 :
392 281448 : COMPLEX(dp), DIMENSION(:), POINTER :: evals
393 : INTEGER :: ev_ind
394 281448 : INTEGER, DIMENSION(:), POINTER :: selected_ind
395 :
396 281448 : IF (ind > get_nval_out(arnoldi_env)) THEN
397 0 : CPABORT('outside range of indexed evals')
398 : END IF
399 :
400 281448 : selected_ind => get_sel_ind(arnoldi_env)
401 281448 : ev_ind = selected_ind(ind)
402 281448 : evals => get_evals(arnoldi_env)
403 281448 : eval_out = evals(ev_ind)
404 :
405 281448 : END FUNCTION get_selected_ritz_val
406 :
407 : ! **************************************************************************************************
408 : !> \brief Get all Ritz values of the selected set. eval_out has to be allocated
409 : !> at least the size of get_neval_out()
410 : !> \param arnoldi_env ...
411 : !> \param eval_out ...
412 : ! **************************************************************************************************
413 0 : SUBROUTINE get_all_selected_ritz_val(arnoldi_env, eval_out)
414 : TYPE(arnoldi_env_type) :: arnoldi_env
415 : COMPLEX(dp), DIMENSION(:) :: eval_out
416 :
417 0 : COMPLEX(dp), DIMENSION(:), POINTER :: evals
418 : INTEGER :: ev_ind, ind
419 0 : INTEGER, DIMENSION(:), POINTER :: selected_ind
420 :
421 0 : NULLIFY (evals)
422 0 : IF (SIZE(eval_out) < get_nval_out(arnoldi_env)) THEN
423 0 : CPABORT('array for eval output too small')
424 : END IF
425 0 : selected_ind => get_sel_ind(arnoldi_env)
426 :
427 0 : evals => get_evals(arnoldi_env)
428 :
429 0 : DO ind = 1, get_nval_out(arnoldi_env)
430 0 : ev_ind = selected_ind(ind)
431 0 : eval_out(ind) = evals(ev_ind)
432 : END DO
433 :
434 0 : END SUBROUTINE get_all_selected_ritz_val
435 :
436 : ! **************************************************************************************************
437 : !> \brief ...
438 : !> \param arnoldi_env ...
439 : !> \param vector ...
440 : ! **************************************************************************************************
441 9792 : SUBROUTINE set_arnoldi_initial_vector(arnoldi_env, vector)
442 : TYPE(arnoldi_env_type) :: arnoldi_env
443 : TYPE(dbcsr_type) :: vector
444 :
445 : INTEGER :: ncol_local, nrow_local
446 4896 : REAL(kind=dp), DIMENSION(:), POINTER :: data_vec
447 : TYPE(arnoldi_control_type), POINTER :: control
448 : TYPE(arnoldi_data_type), POINTER :: ar_data
449 :
450 9792 : control => get_control(arnoldi_env)
451 4896 : control%has_initial_vector = .TRUE.
452 4896 : ar_data => get_data(arnoldi_env)
453 :
454 4896 : CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
455 4896 : data_vec => dbcsr_get_data_p(vector)
456 151293 : IF (nrow_local*ncol_local > 0) ar_data%f_vec(1:nrow_local) = data_vec(1:nrow_local)
457 :
458 4896 : END SUBROUTINE set_arnoldi_initial_vector
459 :
460 : !!! Here come the methods handling the selection of eigenvalues and eigenvectors !!!
461 : !!! If you want a personal method, simply created a Subroutine returning the index
462 : !!! array selected ind which contains as the first nval_out entries the index of the evals
463 :
464 : ! **************************************************************************************************
465 : !> \brief ...
466 : !> \param evals ...
467 : !> \param current_step ...
468 : !> \param selected_ind ...
469 : !> \param neval ...
470 : ! **************************************************************************************************
471 136583 : SUBROUTINE index_min_max_real_eval(evals, current_step, selected_ind, neval)
472 : COMPLEX(dp), DIMENSION(:) :: evals
473 : INTEGER, INTENT(IN) :: current_step
474 : INTEGER, DIMENSION(:) :: selected_ind
475 : INTEGER :: neval
476 :
477 : INTEGER :: i
478 273166 : INTEGER, DIMENSION(current_step) :: indexing
479 273166 : REAL(dp), DIMENSION(current_step) :: tmp_array
480 :
481 136583 : neval = 0
482 4522855 : selected_ind = 0
483 776886 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
484 136583 : CALL sort(tmp_array, current_step, indexing)
485 136583 : DO i = 1, current_step
486 136583 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
487 136583 : selected_ind(1) = indexing(i)
488 136583 : neval = neval + 1
489 136583 : EXIT
490 : END IF
491 : END DO
492 136583 : DO i = current_step, 1, -1
493 136583 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
494 136583 : selected_ind(2) = indexing(i)
495 136583 : neval = neval + 1
496 136583 : EXIT
497 : END IF
498 : END DO
499 :
500 136583 : END SUBROUTINE index_min_max_real_eval
501 :
502 : ! **************************************************************************************************
503 : !> \brief ...
504 : !> \param evals ...
505 : !> \param current_step ...
506 : !> \param selected_ind ...
507 : !> \param neval ...
508 : ! **************************************************************************************************
509 4457 : SUBROUTINE index_nmax_real_eval(evals, current_step, selected_ind, neval)
510 : COMPLEX(dp), DIMENSION(:) :: evals
511 : INTEGER, INTENT(IN) :: current_step
512 : INTEGER, DIMENSION(:) :: selected_ind
513 : INTEGER :: neval
514 :
515 : INTEGER :: i, nlimit
516 8914 : INTEGER, DIMENSION(current_step) :: indexing
517 8914 : REAL(dp), DIMENSION(current_step) :: tmp_array
518 :
519 4457 : nlimit = neval; neval = 0
520 93597 : selected_ind = 0
521 22758 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
522 4457 : CALL sort(tmp_array, current_step, indexing)
523 4457 : DO i = 1, current_step
524 8914 : IF (ABS(AIMAG(evals(indexing(current_step + 1 - i)))) < EPSILON(0.0_dp)) THEN
525 4457 : selected_ind(i) = indexing(current_step + 1 - i)
526 4457 : neval = neval + 1
527 4457 : IF (neval == nlimit) EXIT
528 : END IF
529 : END DO
530 :
531 4457 : END SUBROUTINE index_nmax_real_eval
532 :
533 : ! **************************************************************************************************
534 : !> \brief ...
535 : !> \param evals ...
536 : !> \param current_step ...
537 : !> \param selected_ind ...
538 : !> \param neval ...
539 : ! **************************************************************************************************
540 5331 : SUBROUTINE index_nmin_real_eval(evals, current_step, selected_ind, neval)
541 : COMPLEX(dp), DIMENSION(:) :: evals
542 : INTEGER, INTENT(IN) :: current_step
543 : INTEGER, DIMENSION(:) :: selected_ind
544 : INTEGER :: neval
545 :
546 : INTEGER :: i, nlimit
547 10662 : INTEGER, DIMENSION(current_step) :: indexing
548 10662 : REAL(dp), DIMENSION(current_step) :: tmp_array
549 :
550 5331 : nlimit = neval; neval = 0
551 111951 : selected_ind = 0
552 90200 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
553 5331 : CALL sort(tmp_array, current_step, indexing)
554 5331 : DO i = 1, current_step
555 10662 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
556 5331 : selected_ind(i) = indexing(i)
557 5331 : neval = neval + 1
558 5331 : IF (neval == nlimit) EXIT
559 : END IF
560 : END DO
561 :
562 5331 : END SUBROUTINE index_nmin_real_eval
563 :
564 : END MODULE arnoldi_data_methods
|
