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 148009 : 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 148009 : nval_request, nrestarts, generalized_ev, iram)
69 :
70 148009 : CALL setup_arnoldi_data(arnoldi_env, matrix, max_iter)
71 :
72 148009 : 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 148009 : 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 148009 : ALLOCATE (ar_data)
89 148009 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, nfullrows_local=nrow_local)
90 391062 : ALLOCATE (ar_data%f_vec(nrow_local))
91 243053 : ALLOCATE (ar_data%x_vec(nrow_local))
92 592036 : ALLOCATE (ar_data%Hessenberg(max_iter + 1, max_iter))
93 539071 : ALLOCATE (ar_data%local_history(nrow_local, max_iter))
94 :
95 444027 : ALLOCATE (ar_data%evals(max_iter))
96 592036 : ALLOCATE (ar_data%revec(max_iter, max_iter))
97 :
98 148009 : CALL set_data(arnoldi_env, ar_data)
99 :
100 148009 : 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 148009 : 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 148009 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, distribution=distri)
131 148009 : 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 148009 : pcol_group=pcol_handle)
137 :
138 148009 : CALL control%mp_group%set_handle(group_handle)
139 148009 : CALL control%pcol_group%set_handle(pcol_handle)
140 :
141 148009 : IF (.NOT. subgroups_defined) THEN
142 0 : CPABORT("arnoldi only with subgroups")
143 : END IF
144 :
145 148009 : control%symmetric = .FALSE.
146 : ! Will need a fix for complex because there it has to be hermitian
147 148009 : IF (SIZE(matrix) == 1) THEN
148 143652 : control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
149 : END IF
150 :
151 : ! Set the control parameters
152 148009 : control%max_iter = max_iter
153 148009 : control%current_step = 0
154 148009 : control%selection_crit = selection_crit
155 148009 : control%nval_req = nval_request
156 148009 : control%threshold = threshold
157 148009 : control%converged = .FALSE.
158 148009 : control%has_initial_vector = .FALSE.
159 148009 : control%iram = iram
160 148009 : control%nrestart = nrestarts
161 148009 : control%generalized_ev = generalized_ev
162 :
163 148009 : 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 148009 : 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 148009 : 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 148009 : 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 148009 : 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 444027 : ALLOCATE (control%selected_ind(max_iter))
187 148009 : CALL set_control(arnoldi_env, control)
188 :
189 296018 : END SUBROUTINE setup_arnoldi_control
190 :
191 : ! **************************************************************************************************
192 : !> \brief ...
193 : !> \param arnoldi_env ...
194 : !> \param ind ...
195 : !> \param matrix ...
196 : !> \param vector ...
197 : ! **************************************************************************************************
198 153030 : 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 153030 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:) :: ritz_v
204 : INTEGER :: i, myind, sspace_size, vsize
205 153030 : INTEGER, DIMENSION(:), POINTER :: selected_ind
206 153030 : REAL(kind=dp), DIMENSION(:), POINTER :: data_vec
207 : TYPE(arnoldi_control_type), POINTER :: control
208 : TYPE(arnoldi_data_type), POINTER :: ar_data
209 :
210 306060 : control => get_control(arnoldi_env)
211 153030 : selected_ind => get_sel_ind(arnoldi_env)
212 153030 : ar_data => get_data(arnoldi_env)
213 153030 : sspace_size = get_subsp_size(arnoldi_env)
214 153030 : vsize = SIZE(ar_data%f_vec)
215 153030 : myind = selected_ind(ind)
216 404170 : ALLOCATE (ritz_v(vsize))
217 153030 : ritz_v = CMPLX(0.0, 0.0, dp)
218 :
219 153030 : CALL dbcsr_release(vector)
220 153030 : CALL create_col_vec_from_matrix(vector, matrix, 1)
221 153030 : IF (control%local_comp) THEN
222 686100 : DO i = 1, sspace_size
223 9655570 : ritz_v(:) = ritz_v(:) + ar_data%local_history(:, i)*ar_data%revec(i, myind)
224 : END DO
225 98110 : 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 983119 : data_vec(1:vsize) = REAL(ritz_v(1:vsize), KIND=dp)
229 : END IF
230 :
231 153030 : DEALLOCATE (ritz_v)
232 :
233 153030 : END SUBROUTINE get_selected_ritz_vector
234 :
235 : ! **************************************************************************************************
236 : !> \brief Deallocate the data in arnoldi_env
237 : !> \param arnoldi_env ...
238 : ! **************************************************************************************************
239 148009 : 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 148009 : ar_data => get_data(arnoldi_env)
246 148009 : IF (ASSOCIATED(ar_data%f_vec)) DEALLOCATE (ar_data%f_vec)
247 148009 : IF (ASSOCIATED(ar_data%x_vec)) DEALLOCATE (ar_data%x_vec)
248 148009 : IF (ASSOCIATED(ar_data%Hessenberg)) DEALLOCATE (ar_data%Hessenberg)
249 148009 : IF (ASSOCIATED(ar_data%local_history)) DEALLOCATE (ar_data%local_history)
250 148009 : IF (ASSOCIATED(ar_data%evals)) DEALLOCATE (ar_data%evals)
251 148009 : IF (ASSOCIATED(ar_data%revec)) DEALLOCATE (ar_data%revec)
252 148009 : DEALLOCATE (ar_data)
253 :
254 148009 : control => get_control(arnoldi_env)
255 148009 : DEALLOCATE (control%selected_ind)
256 148009 : DEALLOCATE (control)
257 :
258 148009 : 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 149039 : 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 149039 : control => get_control(arnoldi_env)
273 149039 : ar_data => get_data(arnoldi_env)
274 :
275 149039 : last_el = control%current_step
276 149039 : convergence = REAL(0.0, dp)
277 149039 : my_crit = control%selection_crit
278 149039 : control%nval_out = MIN(control%nval_req, control%current_step)
279 138993 : SELECT CASE (my_crit)
280 : ! minimum and maximum real eval
281 : CASE (1)
282 138993 : 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 4585 : 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 5461 : CALL index_nmin_real_eval(ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
289 : CASE DEFAULT
290 149039 : CPABORT("unknown selection index")
291 : END SELECT
292 : ! test whether we are converged
293 437071 : DO i = 1, control%nval_out
294 288032 : my_ind = control%selected_ind(i)
295 : convergence = MAX(convergence, &
296 437071 : ABS(ar_data%revec(last_el, my_ind)*ar_data%Hessenberg(last_el + 1, last_el)))
297 : END DO
298 149039 : control%converged = convergence < control%threshold
299 :
300 149039 : 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 147791 : 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 147791 : control => get_control(arnoldi_env)
329 147791 : nrestart = control%nrestart
330 :
331 147791 : 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 286524 : 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 286524 : control => get_control(arnoldi_env)
345 286524 : nval_out = control%nval_out
346 :
347 286524 : 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 153030 : 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 153030 : control => get_control(arnoldi_env)
361 153030 : current_step = control%current_step
362 :
363 153030 : 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 138951 : 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 138951 : control => get_control(arnoldi_env)
377 138951 : converged = control%converged
378 :
379 138951 : 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 286524 : 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 286524 : COMPLEX(dp), DIMENSION(:), POINTER :: evals
393 : INTEGER :: ev_ind
394 286524 : INTEGER, DIMENSION(:), POINTER :: selected_ind
395 :
396 286524 : IF (ind > get_nval_out(arnoldi_env)) THEN
397 0 : CPABORT('outside range of indexed evals')
398 : END IF
399 :
400 286524 : selected_ind => get_sel_ind(arnoldi_env)
401 286524 : ev_ind = selected_ind(ind)
402 286524 : evals => get_evals(arnoldi_env)
403 286524 : eval_out = evals(ev_ind)
404 :
405 286524 : 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 10272 : 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 5136 : REAL(kind=dp), DIMENSION(:), POINTER :: data_vec
447 : TYPE(arnoldi_control_type), POINTER :: control
448 : TYPE(arnoldi_data_type), POINTER :: ar_data
449 :
450 10272 : control => get_control(arnoldi_env)
451 5136 : control%has_initial_vector = .TRUE.
452 5136 : ar_data => get_data(arnoldi_env)
453 :
454 5136 : CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
455 5136 : data_vec => dbcsr_get_data_p(vector)
456 152223 : IF (nrow_local*ncol_local > 0) ar_data%f_vec(1:nrow_local) = data_vec(1:nrow_local)
457 :
458 5136 : 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 138993 : 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 277986 : INTEGER, DIMENSION(current_step) :: indexing
479 277986 : REAL(dp), DIMENSION(current_step) :: tmp_array
480 :
481 138993 : neval = 0
482 4602385 : selected_ind = 0
483 783058 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
484 138993 : CALL sort(tmp_array, current_step, indexing)
485 138993 : DO i = 1, current_step
486 138993 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
487 138993 : selected_ind(1) = indexing(i)
488 138993 : neval = neval + 1
489 138993 : EXIT
490 : END IF
491 : END DO
492 138993 : DO i = current_step, 1, -1
493 138993 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
494 138993 : selected_ind(2) = indexing(i)
495 138993 : neval = neval + 1
496 138993 : EXIT
497 : END IF
498 : END DO
499 :
500 138993 : 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 4585 : 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 9170 : INTEGER, DIMENSION(current_step) :: indexing
517 9170 : REAL(dp), DIMENSION(current_step) :: tmp_array
518 :
519 4585 : nlimit = neval; neval = 0
520 96285 : selected_ind = 0
521 23014 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
522 4585 : CALL sort(tmp_array, current_step, indexing)
523 4585 : DO i = 1, current_step
524 9170 : IF (ABS(AIMAG(evals(indexing(current_step + 1 - i)))) < EPSILON(0.0_dp)) THEN
525 4585 : selected_ind(i) = indexing(current_step + 1 - i)
526 4585 : neval = neval + 1
527 4585 : IF (neval == nlimit) EXIT
528 : END IF
529 : END DO
530 :
531 4585 : 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 5461 : 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 10922 : INTEGER, DIMENSION(current_step) :: indexing
548 10922 : REAL(dp), DIMENSION(current_step) :: tmp_array
549 :
550 5461 : nlimit = neval; neval = 0
551 114681 : selected_ind = 0
552 90880 : tmp_array(1:current_step) = REAL(evals(1:current_step), dp)
553 5461 : CALL sort(tmp_array, current_step, indexing)
554 5461 : DO i = 1, current_step
555 10922 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(0.0_dp)) THEN
556 5461 : selected_ind(i) = indexing(i)
557 5461 : neval = neval + 1
558 5461 : IF (neval == nlimit) EXIT
559 : END IF
560 : END DO
561 :
562 5461 : END SUBROUTINE index_nmin_real_eval
563 :
564 : END MODULE arnoldi_data_methods
|
