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