Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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
11 : !> method
12 : !> \par History
13 : !> 2014.09 created [Florian Schiffmann]
14 : !> \author Florian Schiffmann
15 : ! **************************************************************************************************
16 :
17 : MODULE arnoldi_data_methods
18 : USE arnoldi_types, ONLY: &
19 : arnoldi_control_type, arnoldi_data_c_type, arnoldi_data_d_type, arnoldi_data_s_type, &
20 : arnoldi_data_type, arnoldi_data_z_type, get_control, get_data_c, get_data_d, get_data_s, &
21 : get_data_z, get_evals_c, get_evals_d, get_evals_s, get_evals_z, get_sel_ind, has_d_cmplx, &
22 : has_d_real, set_control, set_data_c, set_data_d, set_data_s, &
23 : set_data_z
24 : USE dbcsr_api, ONLY: &
25 : dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_data_p, dbcsr_get_data_type, &
26 : dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_mp_grid_setup, dbcsr_p_type, dbcsr_release, &
27 : dbcsr_type, dbcsr_type_complex_8, dbcsr_type_real_8, dbcsr_type_symmetric
28 : USE dbcsr_vector, ONLY: create_col_vec_from_matrix
29 : USE kinds, ONLY: real_8
30 : USE util, ONLY: sort
31 : #include "../base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : PUBLIC :: select_evals, get_selected_ritz_val, arnoldi_is_converged, &
38 : arnoldi_data_type, get_nrestart, set_arnoldi_initial_vector, &
39 : setup_arnoldi_data, deallocate_arnoldi_data, get_selected_ritz_vector
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief This routine sets the environment for the arnoldi iteration and
45 : !> the krylov subspace creation. All simulation parameters have to be given
46 : !> at this stage so the rest can run fully automated
47 : !> In addition, this routine allocates the data necessary for
48 : !> \param arnoldi_data this type which gets filled with information and
49 : !> on output contains all information necessary to extract
50 : !> whatever the user desires
51 : !> \param matrix vector of matrices, only the first gets used to get some dimensions
52 : !> and parallel information needed later on
53 : !> \param max_iter maximum dimension of the krylov subspace
54 : !> \param threshold convergence threshold, this is used for both subspace and eigenval
55 : !> \param selection_crit integer defining according to which criterion the
56 : !> eigenvalues are selected for the subspace
57 : !> \param nval_request for some sel_crit useful, how many eV to select
58 : !> \param nrestarts ...
59 : !> \param generalized_ev ...
60 : !> \param iram ...
61 : ! **************************************************************************************************
62 130727 : SUBROUTINE setup_arnoldi_data(arnoldi_data, matrix, max_iter, threshold, selection_crit, &
63 : nval_request, nrestarts, generalized_ev, iram)
64 : TYPE(arnoldi_data_type) :: arnoldi_data
65 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
66 : INTEGER :: max_iter
67 : REAL(real_8) :: threshold
68 : INTEGER :: selection_crit, nval_request, nrestarts
69 : LOGICAL :: generalized_ev, iram
70 :
71 : CALL setup_arnoldi_control(arnoldi_data, matrix, max_iter, threshold, selection_crit, &
72 130727 : nval_request, nrestarts, generalized_ev, iram)
73 :
74 130727 : SELECT CASE (dbcsr_get_data_type(matrix(1)%matrix))
75 : CASE (dbcsr_type_real_8)
76 130727 : CALL setup_arnoldi_data_d(arnoldi_data, matrix, max_iter)
77 : CASE (dbcsr_type_complex_8)
78 130727 : CALL setup_arnoldi_data_z(arnoldi_data, matrix, max_iter)
79 : END SELECT
80 :
81 130727 : END SUBROUTINE setup_arnoldi_data
82 :
83 : ! **************************************************************************************************
84 : !> \brief Creates the control type for arnoldi, see above for details
85 : !> \param arnoldi_data ...
86 : !> \param matrix ...
87 : !> \param max_iter ...
88 : !> \param threshold ...
89 : !> \param selection_crit ...
90 : !> \param nval_request ...
91 : !> \param nrestarts ...
92 : !> \param generalized_ev ...
93 : !> \param iram ...
94 : ! **************************************************************************************************
95 130727 : SUBROUTINE setup_arnoldi_control(arnoldi_data, matrix, max_iter, threshold, selection_crit, nval_request, &
96 : nrestarts, generalized_ev, iram)
97 : TYPE(arnoldi_data_type) :: arnoldi_data
98 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
99 : INTEGER :: max_iter
100 : REAL(real_8) :: threshold
101 : INTEGER :: selection_crit, nval_request, nrestarts
102 : LOGICAL :: generalized_ev, iram
103 :
104 : INTEGER :: pcol_handle, group_handle
105 : LOGICAL :: subgroups_defined
106 : TYPE(arnoldi_control_type), POINTER :: control
107 : TYPE(dbcsr_distribution_type) :: distri
108 :
109 0 : ALLOCATE (control)
110 : ! Fill the information which will later on control the behavior of the arnoldi method and allow synchronization
111 130727 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, distribution=distri)
112 130727 : CALL dbcsr_mp_grid_setup(distri)
113 : CALL dbcsr_distribution_get(distri, &
114 : group=group_handle, &
115 : mynode=control%myproc, &
116 : subgroups_defined=subgroups_defined, &
117 130727 : pcol_group=pcol_handle)
118 :
119 130727 : CALL control%mp_group%set_handle(group_handle)
120 130727 : CALL control%pcol_group%set_handle(pcol_handle)
121 :
122 130727 : IF (.NOT. subgroups_defined) &
123 0 : CPABORT("arnoldi only with subgroups")
124 :
125 130727 : control%symmetric = .FALSE.
126 : ! Will need a fix for complex because there it has to be hermitian
127 130727 : IF (SIZE(matrix) == 1) &
128 126837 : control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
129 :
130 : ! Set the control parameters
131 130727 : control%max_iter = max_iter
132 130727 : control%current_step = 0
133 130727 : control%selection_crit = selection_crit
134 130727 : control%nval_req = nval_request
135 130727 : control%threshold = threshold
136 130727 : control%converged = .FALSE.
137 130727 : control%has_initial_vector = .FALSE.
138 130727 : control%iram = iram
139 130727 : control%nrestart = nrestarts
140 130727 : control%generalized_ev = generalized_ev
141 :
142 130727 : IF (control%nval_req > 1 .AND. control%nrestart > 0 .AND. .NOT. control%iram) &
143 : CALL cp_abort(__LOCATION__, 'with more than one eigenvalue requested '// &
144 0 : 'internal restarting with a previous EVEC is a bad idea, set IRAM or nrestsart=0')
145 :
146 : ! some checks for the generalized EV mode
147 130727 : IF (control%generalized_ev .AND. selection_crit == 1) &
148 : CALL cp_abort(__LOCATION__, &
149 0 : 'generalized ev can only highest OR lowest EV')
150 130727 : IF (control%generalized_ev .AND. nval_request .NE. 1) &
151 : CALL cp_abort(__LOCATION__, &
152 0 : 'generalized ev can only compute one EV at the time')
153 130727 : IF (control%generalized_ev .AND. control%nrestart == 0) &
154 : CALL cp_abort(__LOCATION__, &
155 0 : 'outer loops are mandatory for generalized EV, set nrestart appropriatly')
156 130727 : IF (SIZE(matrix) .NE. 2 .AND. control%generalized_ev) &
157 : CALL cp_abort(__LOCATION__, &
158 0 : 'generalized ev needs exactly two matrices as input (2nd is the metric)')
159 :
160 392181 : ALLOCATE (control%selected_ind(max_iter))
161 130727 : CALL set_control(arnoldi_data, control)
162 :
163 261454 : END SUBROUTINE setup_arnoldi_control
164 :
165 : ! **************************************************************************************************
166 : !> \brief Deallocate the data in arnoldi_data
167 : !> \param arnoldi_data ...
168 : !> \param ind ...
169 : !> \param matrix ...
170 : !> \param vector ...
171 : ! **************************************************************************************************
172 135227 : SUBROUTINE get_selected_ritz_vector(arnoldi_data, ind, matrix, vector)
173 : TYPE(arnoldi_data_type) :: arnoldi_data
174 : INTEGER :: ind
175 : TYPE(dbcsr_type) :: matrix, vector
176 :
177 135227 : IF (has_d_real(arnoldi_data)) CALL get_selected_ritz_vector_d(arnoldi_data, ind, matrix, vector)
178 135227 : IF (has_d_cmplx(arnoldi_data)) CALL get_selected_ritz_vector_z(arnoldi_data, ind, matrix, vector)
179 :
180 135227 : END SUBROUTINE get_selected_ritz_vector
181 :
182 : ! **************************************************************************************************
183 : !> \brief Deallocate the data in arnoldi_data
184 : !> \param arnoldi_data ...
185 : ! **************************************************************************************************
186 130727 : SUBROUTINE deallocate_arnoldi_data(arnoldi_data)
187 : TYPE(arnoldi_data_type) :: arnoldi_data
188 :
189 : TYPE(arnoldi_control_type), POINTER :: control
190 :
191 130727 : IF (has_d_real(arnoldi_data)) CALL deallocate_arnoldi_data_d(arnoldi_data)
192 130727 : IF (has_d_cmplx(arnoldi_data)) CALL deallocate_arnoldi_data_z(arnoldi_data)
193 :
194 130727 : control => get_control(arnoldi_data)
195 130727 : DEALLOCATE (control%selected_ind)
196 130727 : DEALLOCATE (control)
197 :
198 130727 : END SUBROUTINE deallocate_arnoldi_data
199 :
200 : ! **************************************************************************************************
201 : !> \brief perform the selection of eigenvalues, fills the selected_ind array
202 : !> \param arnoldi_data ...
203 : ! **************************************************************************************************
204 131716 : SUBROUTINE select_evals(arnoldi_data)
205 : TYPE(arnoldi_data_type) :: arnoldi_data
206 :
207 131716 : IF (has_d_real(arnoldi_data)) CALL select_evals_d(arnoldi_data)
208 131716 : IF (has_d_cmplx(arnoldi_data)) CALL select_evals_z(arnoldi_data)
209 :
210 131716 : END SUBROUTINE select_evals
211 :
212 : ! **************************************************************************************************
213 : !> \brief set a new selection type, if you notice you didn't like the initial one
214 : !> \param ar_data ...
215 : !> \param itype ...
216 : ! **************************************************************************************************
217 0 : SUBROUTINE set_eval_selection(ar_data, itype)
218 : TYPE(arnoldi_data_type) :: ar_data
219 : INTEGER :: itype
220 :
221 : TYPE(arnoldi_control_type), POINTER :: control
222 :
223 0 : control => get_control(ar_data)
224 0 : control%selection_crit = itype
225 0 : END SUBROUTINE set_eval_selection
226 :
227 : ! **************************************************************************************************
228 : !> \brief returns the number of restarts allowed for arnoldi
229 : !> \param arnoldi_data ...
230 : !> \return ...
231 : ! **************************************************************************************************
232 130529 : FUNCTION get_nrestart(arnoldi_data) RESULT(nrestart)
233 : TYPE(arnoldi_data_type) :: arnoldi_data
234 : INTEGER :: nrestart
235 :
236 : TYPE(arnoldi_control_type), POINTER :: control
237 :
238 130529 : control => get_control(arnoldi_data)
239 130529 : nrestart = control%nrestart
240 :
241 130529 : END FUNCTION get_nrestart
242 :
243 : ! **************************************************************************************************
244 : !> \brief get the number of eigenvalues matching the search criterion
245 : !> \param ar_data ...
246 : !> \return ...
247 : ! **************************************************************************************************
248 252976 : FUNCTION get_nval_out(ar_data) RESULT(nval_out)
249 : TYPE(arnoldi_data_type) :: ar_data
250 : INTEGER :: nval_out
251 :
252 : TYPE(arnoldi_control_type), POINTER :: control
253 :
254 252976 : control => get_control(ar_data)
255 252976 : nval_out = control%nval_out
256 :
257 252976 : END FUNCTION get_nval_out
258 :
259 : ! **************************************************************************************************
260 : !> \brief get the dimension of the krylov space. Can be less than max_iter
261 : !> if subspace converged early
262 : !> \param ar_data ...
263 : !> \return ...
264 : ! **************************************************************************************************
265 135227 : FUNCTION get_subsp_size(ar_data) RESULT(current_step)
266 : TYPE(arnoldi_data_type) :: ar_data
267 : INTEGER :: current_step
268 :
269 : TYPE(arnoldi_control_type), POINTER :: control
270 :
271 135227 : control => get_control(ar_data)
272 135227 : current_step = control%current_step
273 :
274 135227 : END FUNCTION get_subsp_size
275 :
276 : ! **************************************************************************************************
277 : !> \brief Find out whether the method with the current search criterion is converged
278 : !> \param ar_data ...
279 : !> \return ...
280 : ! **************************************************************************************************
281 122645 : FUNCTION arnoldi_is_converged(ar_data) RESULT(converged)
282 : TYPE(arnoldi_data_type) :: ar_data
283 : LOGICAL :: converged
284 :
285 : TYPE(arnoldi_control_type), POINTER :: control
286 :
287 122645 : control => get_control(ar_data)
288 122645 : converged = control%converged
289 :
290 122645 : END FUNCTION
291 :
292 : ! **************************************************************************************************
293 : !> \brief get a single specific Ritz value from the set of selected
294 : !> \param ar_data ...
295 : !> \param ind ...
296 : !> \return ...
297 : ! **************************************************************************************************
298 252976 : FUNCTION get_selected_ritz_val(ar_data, ind) RESULT(eval_out)
299 : TYPE(arnoldi_data_type) :: ar_data
300 : INTEGER :: ind
301 : COMPLEX(real_8) :: eval_out
302 :
303 252976 : COMPLEX(real_8), DIMENSION(:), POINTER :: evals_d
304 : INTEGER :: ev_ind
305 252976 : INTEGER, DIMENSION(:), POINTER :: selected_ind
306 :
307 252976 : IF (ind .GT. get_nval_out(ar_data)) &
308 0 : CPABORT('outside range of indexed evals')
309 :
310 252976 : selected_ind => get_sel_ind(ar_data)
311 252976 : ev_ind = selected_ind(ind)
312 252976 : IF (has_d_real(ar_data)) THEN
313 252976 : evals_d => get_evals_d(ar_data); eval_out = evals_d(ev_ind)
314 0 : ELSE IF (has_d_cmplx(ar_data)) THEN
315 0 : evals_d => get_evals_z(ar_data); eval_out = evals_d(ev_ind)
316 : END IF
317 :
318 252976 : END FUNCTION get_selected_ritz_val
319 :
320 : ! **************************************************************************************************
321 : !> \brief Get all Ritz values of the selected set. eval_out has to be allocated
322 : !> at least the size of get_neval_out()
323 : !> \param ar_data ...
324 : !> \param eval_out ...
325 : ! **************************************************************************************************
326 0 : SUBROUTINE get_all_selected_ritz_val(ar_data, eval_out)
327 : TYPE(arnoldi_data_type) :: ar_data
328 : COMPLEX(real_8), DIMENSION(:) :: eval_out
329 :
330 0 : COMPLEX(real_8), DIMENSION(:), POINTER :: evals_d
331 : INTEGER :: ev_ind, ind
332 0 : INTEGER, DIMENSION(:), POINTER :: selected_ind
333 :
334 0 : NULLIFY (evals_d)
335 0 : IF (SIZE(eval_out) .LT. get_nval_out(ar_data)) &
336 0 : CPABORT('array for eval output too small')
337 0 : selected_ind => get_sel_ind(ar_data)
338 :
339 0 : IF (has_d_real(ar_data)) evals_d => get_evals_d(ar_data)
340 0 : IF (has_d_cmplx(ar_data)) evals_d => get_evals_d(ar_data)
341 :
342 0 : DO ind = 1, get_nval_out(ar_data)
343 0 : ev_ind = selected_ind(ind)
344 0 : IF (ASSOCIATED(evals_d)) eval_out(ind) = evals_d(ev_ind)
345 : END DO
346 :
347 0 : END SUBROUTINE get_all_selected_ritz_val
348 :
349 : ! **************************************************************************************************
350 : !> \brief ...
351 : !> \param ar_data ...
352 : !> \param vector ...
353 : ! **************************************************************************************************
354 4770 : SUBROUTINE set_arnoldi_initial_vector(ar_data, vector)
355 : TYPE(arnoldi_data_type) :: ar_data
356 : TYPE(dbcsr_type) :: vector
357 :
358 : TYPE(arnoldi_control_type), POINTER :: control
359 :
360 4770 : control => get_control(ar_data)
361 4770 : control%has_initial_vector = .TRUE.
362 :
363 4770 : IF (has_d_real(ar_data)) CALL set_initial_vector_d(ar_data, vector)
364 4770 : IF (has_d_cmplx(ar_data)) CALL set_initial_vector_z(ar_data, vector)
365 :
366 4770 : END SUBROUTINE set_arnoldi_initial_vector
367 :
368 : !!! Here come the methods handling the selection of eigenvalues and eigenvectors !!!
369 : !!! If you want a personal method, simply created a Subroutine returning the index
370 : !!! array selected ind which contains as the first nval_out entries the index of the evals
371 :
372 : #:include 'arnoldi.fypp'
373 : #:for nametype1, type_prec, real_zero, nametype_zero, type_nametype1, vartype in inst_params_1
374 : ! **************************************************************************************************
375 : !> \brief ...
376 : !> \param arnoldi_data ...
377 : ! **************************************************************************************************
378 131716 : SUBROUTINE select_evals_${nametype1}$ (arnoldi_data)
379 : TYPE(arnoldi_data_type) :: arnoldi_data
380 :
381 : INTEGER :: my_crit, last_el, my_ind, i
382 : REAL(${type_prec}$) :: convergence
383 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
384 : TYPE(arnoldi_control_type), POINTER :: control
385 :
386 131716 : control => get_control(arnoldi_data)
387 131716 : ar_data => get_data_${nametype1}$ (arnoldi_data)
388 :
389 131716 : last_el = control%current_step
390 131716 : convergence = REAL(0.0, ${type_prec}$)
391 131716 : my_crit = control%selection_crit
392 131716 : control%nval_out = MIN(control%nval_req, control%current_step)
393 122687 : SELECT CASE (my_crit)
394 : ! minimum and maximum real eval
395 : CASE (1)
396 122687 : CALL index_min_max_real_eval_${nametype1}$ (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
397 : ! n maximum real eval
398 : CASE (2)
399 4202 : CALL index_nmax_real_eval_${nametype1}$ (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
400 : ! n minimum real eval
401 : CASE (3)
402 4827 : CALL index_nmin_real_eval_${nametype1}$ (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
403 : CASE DEFAULT
404 131716 : CPABORT("unknown selection index")
405 : END SELECT
406 : ! test whether we are converged
407 386119 : DO i = 1, control%nval_out
408 254403 : my_ind = control%selected_ind(i)
409 : convergence = MAX(convergence, &
410 386119 : ABS(ar_data%revec(last_el, my_ind)*ar_data%Hessenberg(last_el + 1, last_el)))
411 : END DO
412 131716 : control%converged = convergence .LT. control%threshold
413 :
414 131716 : END SUBROUTINE select_evals_${nametype1}$
415 :
416 : ! **************************************************************************************************
417 : !> \brief ...
418 : !> \param evals ...
419 : !> \param current_step ...
420 : !> \param selected_ind ...
421 : !> \param neval ...
422 : ! **************************************************************************************************
423 122687 : SUBROUTINE index_min_max_real_eval_${nametype1}$ (evals, current_step, selected_ind, neval)
424 : COMPLEX(${type_prec}$), DIMENSION(:) :: evals
425 : INTEGER, INTENT(IN) :: current_step
426 : INTEGER, DIMENSION(:) :: selected_ind
427 : INTEGER :: neval
428 :
429 245374 : INTEGER, DIMENSION(current_step) :: indexing
430 245374 : REAL(${type_prec}$), DIMENSION(current_step) :: tmp_array
431 : INTEGER :: i
432 :
433 122687 : neval = 0
434 4062623 : selected_ind = 0
435 694050 : tmp_array(1:current_step) = REAL(evals(1:current_step), ${type_prec}$)
436 122687 : CALL sort(tmp_array, current_step, indexing)
437 122687 : DO i = 1, current_step
438 122687 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(${real_zero}$)) THEN
439 122687 : selected_ind(1) = indexing(i)
440 122687 : neval = neval + 1
441 122687 : EXIT
442 : END IF
443 : END DO
444 122687 : DO i = current_step, 1, -1
445 122687 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(${real_zero}$)) THEN
446 122687 : selected_ind(2) = indexing(i)
447 122687 : neval = neval + 1
448 122687 : EXIT
449 : END IF
450 : END DO
451 :
452 122687 : END SUBROUTINE index_min_max_real_eval_${nametype1}$
453 :
454 : ! **************************************************************************************************
455 : !> \brief ...
456 : !> \param evals ...
457 : !> \param current_step ...
458 : !> \param selected_ind ...
459 : !> \param neval ...
460 : ! **************************************************************************************************
461 4202 : SUBROUTINE index_nmax_real_eval_${nametype1}$ (evals, current_step, selected_ind, neval)
462 : COMPLEX(${type_prec}$), DIMENSION(:) :: evals
463 : INTEGER, INTENT(IN) :: current_step
464 : INTEGER, DIMENSION(:) :: selected_ind
465 : INTEGER :: neval
466 :
467 : INTEGER :: i, nlimit
468 8404 : INTEGER, DIMENSION(current_step) :: indexing
469 8404 : REAL(${type_prec}$), DIMENSION(current_step) :: tmp_array
470 :
471 4202 : nlimit = neval; neval = 0
472 88242 : selected_ind = 0
473 24513 : tmp_array(1:current_step) = REAL(evals(1:current_step), ${type_prec}$)
474 4202 : CALL sort(tmp_array, current_step, indexing)
475 4202 : DO i = 1, current_step
476 8404 : IF (ABS(AIMAG(evals(indexing(current_step + 1 - i)))) < EPSILON(${real_zero}$)) THEN
477 4202 : selected_ind(i) = indexing(current_step + 1 - i)
478 4202 : neval = neval + 1
479 4202 : IF (neval == nlimit) EXIT
480 : END IF
481 : END DO
482 :
483 4202 : END SUBROUTINE index_nmax_real_eval_${nametype1}$
484 :
485 : ! **************************************************************************************************
486 : !> \brief ...
487 : !> \param evals ...
488 : !> \param current_step ...
489 : !> \param selected_ind ...
490 : !> \param neval ...
491 : ! **************************************************************************************************
492 4827 : SUBROUTINE index_nmin_real_eval_${nametype1}$ (evals, current_step, selected_ind, neval)
493 : COMPLEX(${type_prec}$), DIMENSION(:) :: evals
494 : INTEGER, INTENT(IN) :: current_step
495 : INTEGER, DIMENSION(:) :: selected_ind
496 : INTEGER :: neval, nlimit
497 :
498 : INTEGER :: i
499 9654 : INTEGER, DIMENSION(current_step) :: indexing
500 9654 : REAL(${type_prec}$), DIMENSION(current_step) :: tmp_array
501 :
502 4827 : nlimit = neval; neval = 0
503 101367 : selected_ind = 0
504 81578 : tmp_array(1:current_step) = REAL(evals(1:current_step), ${type_prec}$)
505 4827 : CALL sort(tmp_array, current_step, indexing)
506 4827 : DO i = 1, current_step
507 9654 : IF (ABS(AIMAG(evals(indexing(i)))) < EPSILON(${real_zero}$)) THEN
508 4827 : selected_ind(i) = indexing(i)
509 4827 : neval = neval + 1
510 4827 : IF (neval == nlimit) EXIT
511 : END IF
512 : END DO
513 :
514 4827 : END SUBROUTINE index_nmin_real_eval_${nametype1}$
515 :
516 : ! **************************************************************************************************
517 : !> \brief ...
518 : !> \param arnoldi_data ...
519 : !> \param matrix ...
520 : !> \param max_iter ...
521 : ! **************************************************************************************************
522 130727 : SUBROUTINE setup_arnoldi_data_${nametype1}$ (arnoldi_data, matrix, max_iter)
523 : TYPE(arnoldi_data_type) :: arnoldi_data
524 : TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
525 : INTEGER :: max_iter
526 :
527 : INTEGER :: nrow_local
528 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
529 :
530 130727 : ALLOCATE (ar_data)
531 130727 : CALL dbcsr_get_info(matrix=matrix(1)%matrix, nfullrows_local=nrow_local)
532 347455 : ALLOCATE (ar_data%f_vec(nrow_local))
533 216728 : ALLOCATE (ar_data%x_vec(nrow_local))
534 522908 : ALLOCATE (ar_data%Hessenberg(max_iter + 1, max_iter))
535 478182 : ALLOCATE (ar_data%local_history(nrow_local, max_iter))
536 :
537 392181 : ALLOCATE (ar_data%evals(max_iter))
538 522908 : ALLOCATE (ar_data%revec(max_iter, max_iter))
539 :
540 130727 : CALL set_data_${nametype1}$ (arnoldi_data, ar_data)
541 :
542 130727 : END SUBROUTINE setup_arnoldi_data_${nametype1}$
543 :
544 : ! **************************************************************************************************
545 : !> \brief ...
546 : !> \param arnoldi_data ...
547 : ! **************************************************************************************************
548 130727 : SUBROUTINE deallocate_arnoldi_data_${nametype1}$ (arnoldi_data)
549 : TYPE(arnoldi_data_type) :: arnoldi_data
550 :
551 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
552 :
553 130727 : ar_data => get_data_${nametype1}$ (arnoldi_data)
554 130727 : IF (ASSOCIATED(ar_data%f_vec)) DEALLOCATE (ar_data%f_vec)
555 130727 : IF (ASSOCIATED(ar_data%x_vec)) DEALLOCATE (ar_data%x_vec)
556 130727 : IF (ASSOCIATED(ar_data%Hessenberg)) DEALLOCATE (ar_data%Hessenberg)
557 130727 : IF (ASSOCIATED(ar_data%local_history)) DEALLOCATE (ar_data%local_history)
558 130727 : IF (ASSOCIATED(ar_data%evals)) DEALLOCATE (ar_data%evals)
559 130727 : IF (ASSOCIATED(ar_data%revec)) DEALLOCATE (ar_data%revec)
560 130727 : DEALLOCATE (ar_data)
561 :
562 130727 : END SUBROUTINE deallocate_arnoldi_data_${nametype1}$
563 :
564 : ! **************************************************************************************************
565 : !> \brief ...
566 : !> \param arnoldi_data ...
567 : !> \param ind ...
568 : !> \param matrix ...
569 : !> \param vector ...
570 : ! **************************************************************************************************
571 135227 : SUBROUTINE get_selected_ritz_vector_${nametype1}$ (arnoldi_data, ind, matrix, vector)
572 : TYPE(arnoldi_data_type) :: arnoldi_data
573 : INTEGER :: ind
574 : TYPE(dbcsr_type) :: matrix
575 : TYPE(dbcsr_type) :: vector
576 :
577 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
578 : INTEGER :: vsize, myind, sspace_size, i
579 135227 : INTEGER, DIMENSION(:), POINTER :: selected_ind
580 135227 : COMPLEX(${type_prec}$), DIMENSION(:), ALLOCATABLE :: ritz_v
581 135227 : ${type_nametype1}$, DIMENSION(:), POINTER :: data_vec
582 : TYPE(arnoldi_control_type), POINTER :: control
583 :
584 270454 : control => get_control(arnoldi_data)
585 135227 : selected_ind => get_sel_ind(arnoldi_data)
586 135227 : ar_data => get_data_${nametype1}$ (arnoldi_data)
587 135227 : sspace_size = get_subsp_size(arnoldi_data)
588 135227 : vsize = SIZE(ar_data%f_vec)
589 135227 : myind = selected_ind(ind)
590 359322 : ALLOCATE (ritz_v(vsize))
591 970530 : ritz_v = CMPLX(0.0, 0.0, ${type_prec}$)
592 :
593 135227 : CALL dbcsr_release(vector)
594 135227 : CALL create_col_vec_from_matrix(vector, matrix, 1)
595 135227 : IF (control%local_comp) THEN
596 625519 : DO i = 1, sspace_size
597 9003042 : ritz_v(:) = ritz_v(:) + ar_data%local_history(:, i)*ar_data%revec(i, myind)
598 : END DO
599 88868 : data_vec => dbcsr_get_data_p(vector, select_data_type=${nametype_zero}$)
600 : ! is a bit odd but ritz_v is always complex and matrix type determines where it goes
601 : ! again I hope the user knows what is required
602 924171 : data_vec(1:vsize) = ${vartype}$ (ritz_v(1:vsize), KIND=${type_prec}$)
603 : END IF
604 :
605 135227 : DEALLOCATE (ritz_v)
606 :
607 135227 : END SUBROUTINE get_selected_ritz_vector_${nametype1}$
608 :
609 : ! **************************************************************************************************
610 : !> \brief ...
611 : !> \param arnoldi_data ...
612 : !> \param vector ...
613 : ! **************************************************************************************************
614 9540 : SUBROUTINE set_initial_vector_${nametype1}$ (arnoldi_data, vector)
615 : TYPE(arnoldi_data_type) :: arnoldi_data
616 : TYPE(dbcsr_type) :: vector
617 :
618 : TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data
619 4770 : ${type_nametype1}$, DIMENSION(:), POINTER :: data_vec
620 : INTEGER :: nrow_local, ncol_local
621 : TYPE(arnoldi_control_type), POINTER :: control
622 :
623 4770 : control => get_control(arnoldi_data)
624 :
625 4770 : CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
626 4770 : ar_data => get_data_${nametype1}$ (arnoldi_data)
627 4770 : data_vec => dbcsr_get_data_p(vector, select_data_type=${nametype_zero}$)
628 166913 : IF (nrow_local*ncol_local > 0) ar_data%f_vec(1:nrow_local) = data_vec(1:nrow_local)
629 :
630 4770 : END SUBROUTINE set_initial_vector_${nametype1}$
631 :
632 : #:endfor
633 :
634 : END MODULE arnoldi_data_methods
635 :
|
