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 137255 : 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 137255 : nval_request, nrestarts, generalized_ev, iram)
69 :
70 137255 : CALL setup_arnoldi_data(arnoldi_env, matrix, max_iter)
71 :
72 137255 : 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 137255 : 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 137255 : ALLOCATE (ar_data)
89 137255 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, nfullrows_local=nrow_local)
90 364127 : ALLOCATE (ar_data%f_vec(nrow_local))
91 226872 : ALLOCATE (ar_data%x_vec(nrow_local))
92 549020 : ALLOCATE (ar_data%Hessenberg(max_iter + 1, max_iter))
93 501382 : ALLOCATE (ar_data%local_history(nrow_local, max_iter))
94 :
95 411765 : ALLOCATE (ar_data%evals(max_iter))
96 549020 : ALLOCATE (ar_data%revec(max_iter, max_iter))
97 :
98 137255 : CALL set_data(arnoldi_env, ar_data)
99 :
100 137255 : 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 137255 : 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 137255 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, distribution=distri)
131 137255 : 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 137255 : pcol_group=pcol_handle)
137 :
138 137255 : CALL control%mp_group%set_handle(group_handle)
139 137255 : CALL control%pcol_group%set_handle(pcol_handle)
140 :
141 137255 : IF (.NOT. subgroups_defined) &
142 0 : CPABORT("arnoldi only with subgroups")
143 :
144 137255 : control%symmetric = .FALSE.
145 : ! Will need a fix for complex because there it has to be hermitian
146 137255 : IF (SIZE(matrix) == 1) &
147 133312 : control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
148 :
149 : ! Set the control parameters
150 137255 : control%max_iter = max_iter
151 137255 : control%current_step = 0
152 137255 : control%selection_crit = selection_crit
153 137255 : control%nval_req = nval_request
154 137255 : control%threshold = threshold
155 137255 : control%converged = .FALSE.
156 137255 : control%has_initial_vector = .FALSE.
157 137255 : control%iram = iram
158 137255 : control%nrestart = nrestarts
159 137255 : control%generalized_ev = generalized_ev
160 :
161 137255 : 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 137255 : IF (control%generalized_ev .AND. selection_crit == 1) &
167 : CALL cp_abort(__LOCATION__, &
168 0 : 'generalized ev can only highest OR lowest EV')
169 137255 : IF (control%generalized_ev .AND. nval_request /= 1) &
170 : CALL cp_abort(__LOCATION__, &
171 0 : 'generalized ev can only compute one EV at the time')
172 137255 : 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 137255 : IF (SIZE(matrix) /= 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 411765 : ALLOCATE (control%selected_ind(max_iter))
180 137255 : CALL set_control(arnoldi_env, control)
181 :
182 274510 : END SUBROUTINE setup_arnoldi_control
183 :
184 : ! **************************************************************************************************
185 : !> \brief ...
186 : !> \param arnoldi_env ...
187 : !> \param ind ...
188 : !> \param matrix ...
189 : !> \param vector ...
190 : ! **************************************************************************************************
191 141914 : 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 141914 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: ritz_v
197 : INTEGER :: i, myind, sspace_size, vsize
198 141914 : INTEGER, DIMENSION(:), POINTER :: selected_ind
199 141914 : REAL(kind=dp), DIMENSION(:), POINTER :: data_vec
200 : TYPE(arnoldi_control_type), POINTER :: control
201 : TYPE(arnoldi_data_type), POINTER :: ar_data
202 :
203 283828 : control => get_control(arnoldi_env)
204 141914 : selected_ind => get_sel_ind(arnoldi_env)
205 141914 : ar_data => get_data(arnoldi_env)
206 141914 : sspace_size = get_subsp_size(arnoldi_env)
207 141914 : vsize = SIZE(ar_data%f_vec)
208 141914 : myind = selected_ind(ind)
209 376397 : ALLOCATE (ritz_v(vsize))
210 1013840 : ritz_v = CMPLX(0.0, 0.0, dp)
211 :
212 141914 : CALL dbcsr_release(vector)
213 141914 : CALL create_col_vec_from_matrix(vector, matrix, 1)
214 141914 : IF (control%local_comp) THEN
215 660244 : DO i = 1, sspace_size
216 9721257 : ritz_v(:) = ritz_v(:) + ar_data%local_history(:, i)*ar_data%revec(i, myind)
217 : END DO
218 92569 : 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 964495 : data_vec(1:vsize) = REAL(ritz_v(1:vsize), KIND=dp)
222 : END IF
223 :
224 141914 : DEALLOCATE (ritz_v)
225 :
226 141914 : END SUBROUTINE get_selected_ritz_vector
227 :
228 : ! **************************************************************************************************
229 : !> \brief Deallocate the data in arnoldi_env
230 : !> \param arnoldi_env ...
231 : ! **************************************************************************************************
232 137255 : 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 137255 : ar_data => get_data(arnoldi_env)
239 137255 : IF (ASSOCIATED(ar_data%f_vec)) DEALLOCATE (ar_data%f_vec)
240 137255 : IF (ASSOCIATED(ar_data%x_vec)) DEALLOCATE (ar_data%x_vec)
241 137255 : IF (ASSOCIATED(ar_data%Hessenberg)) DEALLOCATE (ar_data%Hessenberg)
242 137255 : IF (ASSOCIATED(ar_data%local_history)) DEALLOCATE (ar_data%local_history)
243 137255 : IF (ASSOCIATED(ar_data%evals)) DEALLOCATE (ar_data%evals)
244 137255 : IF (ASSOCIATED(ar_data%revec)) DEALLOCATE (ar_data%revec)
245 137255 : DEALLOCATE (ar_data)
246 :
247 137255 : control => get_control(arnoldi_env)
248 137255 : DEALLOCATE (control%selected_ind)
249 137255 : DEALLOCATE (control)
250 :
251 137255 : 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 138267 : 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 138267 : control => get_control(arnoldi_env)
266 138267 : ar_data => get_data(arnoldi_env)
267 :
268 138267 : last_el = control%current_step
269 138267 : convergence = REAL(0.0, dp)
270 138267 : my_crit = control%selection_crit
271 138267 : control%nval_out = MIN(control%nval_req, control%current_step)
272 129045 : SELECT CASE (my_crit)
273 : ! minimum and maximum real eval
274 : CASE (1)
275 129045 : 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 4201 : 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 5021 : CALL index_nmin_real_eval(ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
282 : CASE DEFAULT
283 138267 : CPABORT("unknown selection index")
284 : END SELECT
285 : ! test whether we are converged
286 405579 : DO i = 1, control%nval_out
287 267312 : my_ind = control%selected_ind(i)
288 : convergence = MAX(convergence, &
289 405579 : ABS(ar_data%revec(last_el, my_ind)*ar_data%Hessenberg(last_el + 1, last_el)))
290 : END DO
291 138267 : control%converged = convergence < control%threshold
292 :
293 138267 : 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 137059 : 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 137059 : control => get_control(arnoldi_env)
322 137059 : nrestart = control%nrestart
323 :
324 137059 : 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 265844 : 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 265844 : control => get_control(arnoldi_env)
338 265844 : nval_out = control%nval_out
339 :
340 265844 : 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 141914 : 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 141914 : control => get_control(arnoldi_env)
354 141914 : current_step = control%current_step
355 :
356 141914 : 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 128981 : 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 128981 : control => get_control(arnoldi_env)
370 128981 : converged = control%converged
371 :
372 128981 : 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 265844 : 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 265844 : COMPLEX(dp), DIMENSION(:), POINTER :: evals
386 : INTEGER :: ev_ind
387 265844 : INTEGER, DIMENSION(:), POINTER :: selected_ind
388 :
389 265844 : IF (ind > get_nval_out(arnoldi_env)) &
390 0 : CPABORT('outside range of indexed evals')
391 :
392 265844 : selected_ind => get_sel_ind(arnoldi_env)
393 265844 : ev_ind = selected_ind(ind)
394 265844 : evals => get_evals(arnoldi_env)
395 265844 : eval_out = evals(ev_ind)
396 :
397 265844 : 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) < 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 9464 : 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 4732 : REAL(kind=dp), DIMENSION(:), POINTER :: data_vec
438 : TYPE(arnoldi_control_type), POINTER :: control
439 : TYPE(arnoldi_data_type), POINTER :: ar_data
440 :
441 9464 : control => get_control(arnoldi_env)
442 4732 : control%has_initial_vector = .TRUE.
443 4732 : ar_data => get_data(arnoldi_env)
444 :
445 4732 : CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
446 4732 : data_vec => dbcsr_get_data_p(vector)
447 160688 : IF (nrow_local*ncol_local > 0) ar_data%f_vec(1:nrow_local) = data_vec(1:nrow_local)
448 :
449 4732 : 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 129045 : 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 258090 : INTEGER, DIMENSION(current_step) :: indexing
470 258090 : REAL(dp), DIMENSION(current_step) :: tmp_array
471 :
472 129045 : neval = 0
473 4272693 : selected_ind = 0
474 744196 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
475 129045 : CALL sort(tmp_array, current_step, indexing)
476 129045 : DO i = 1, current_step
477 129045 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
478 129045 : selected_ind(1) = indexing(i)
479 129045 : neval = neval + 1
480 129045 : EXIT
481 : END IF
482 : END DO
483 129045 : DO i = current_step, 1, -1
484 129045 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
485 129045 : selected_ind(2) = indexing(i)
486 129045 : neval = neval + 1
487 129045 : EXIT
488 : END IF
489 : END DO
490 :
491 129045 : 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 4201 : 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 8402 : INTEGER, DIMENSION(current_step) :: indexing
508 8402 : REAL(dp), DIMENSION(current_step) :: tmp_array
509 :
510 4201 : nlimit = neval; neval = 0
511 88221 : selected_ind = 0
512 22412 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
513 4201 : CALL sort(tmp_array, current_step, indexing)
514 4201 : DO i = 1, current_step
515 8402 : IF (ABS(AIMAG(evals(indexing(current_step + 1 - i)))) < EPSILON(0.0_dp)) THEN
516 4201 : selected_ind(i) = indexing(current_step + 1 - i)
517 4201 : neval = neval + 1
518 4201 : IF (neval == nlimit) EXIT
519 : END IF
520 : END DO
521 :
522 4201 : 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 5021 : 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 10042 : INTEGER, DIMENSION(current_step) :: indexing
539 10042 : REAL(dp), DIMENSION(current_step) :: tmp_array
540 :
541 5021 : nlimit = neval; neval = 0
542 105441 : selected_ind = 0
543 85136 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
544 5021 : CALL sort(tmp_array, current_step, indexing)
545 5021 : DO i = 1, current_step
546 10042 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
547 5021 : selected_ind(i) = indexing(i)
548 5021 : neval = neval + 1
549 5021 : IF (neval == nlimit) EXIT
550 : END IF
551 : END DO
552 :
553 5021 : END SUBROUTINE index_nmin_real_eval
554 :
555 : END MODULE arnoldi_data_methods
|
