(git:374b731)
Loading...
Searching...
No Matches
arnoldi_data_methods.F
Go to the documentation of this file.
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
17MODULE 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
41CONTAINS
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 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 nval_request, nrestarts, generalized_ev, iram)
73
74 SELECT CASE (dbcsr_get_data_type(matrix(1)%matrix))
75 CASE (dbcsr_type_real_8)
76 CALL setup_arnoldi_data_d(arnoldi_data, matrix, max_iter)
77 CASE (dbcsr_type_complex_8)
78 CALL setup_arnoldi_data_z(arnoldi_data, matrix, max_iter)
79 END SELECT
80
81 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 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 ALLOCATE (control)
110! Fill the information which will later on control the behavior of the arnoldi method and allow synchronization
111 CALL dbcsr_get_info(matrix=matrix(1)%matrix, distribution=distri)
112 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 pcol_group=pcol_handle)
118
119 CALL control%mp_group%set_handle(group_handle)
120 CALL control%pcol_group%set_handle(pcol_handle)
121
122 IF (.NOT. subgroups_defined) &
123 cpabort("arnoldi only with subgroups")
124
125 control%symmetric = .false.
126! Will need a fix for complex because there it has to be hermitian
127 IF (SIZE(matrix) == 1) &
128 control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
129
130! Set the control parameters
131 control%max_iter = max_iter
132 control%current_step = 0
133 control%selection_crit = selection_crit
134 control%nval_req = nval_request
135 control%threshold = threshold
136 control%converged = .false.
137 control%has_initial_vector = .false.
138 control%iram = iram
139 control%nrestart = nrestarts
140 control%generalized_ev = generalized_ev
141
142 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 '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 IF (control%generalized_ev .AND. selection_crit == 1) &
148 CALL cp_abort(__location__, &
149 'generalized ev can only highest OR lowest EV')
150 IF (control%generalized_ev .AND. nval_request .NE. 1) &
151 CALL cp_abort(__location__, &
152 'generalized ev can only compute one EV at the time')
153 IF (control%generalized_ev .AND. control%nrestart == 0) &
154 CALL cp_abort(__location__, &
155 'outer loops are mandatory for generalized EV, set nrestart appropriatly')
156 IF (SIZE(matrix) .NE. 2 .AND. control%generalized_ev) &
157 CALL cp_abort(__location__, &
158 'generalized ev needs exactly two matrices as input (2nd is the metric)')
159
160 ALLOCATE (control%selected_ind(max_iter))
161 CALL set_control(arnoldi_data, control)
162
163 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 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 IF (has_d_real(arnoldi_data)) CALL get_selected_ritz_vector_d(arnoldi_data, ind, matrix, vector)
178 IF (has_d_cmplx(arnoldi_data)) CALL get_selected_ritz_vector_z(arnoldi_data, ind, matrix, vector)
179
180 END SUBROUTINE get_selected_ritz_vector
181
182! **************************************************************************************************
183!> \brief Deallocate the data in arnoldi_data
184!> \param arnoldi_data ...
185! **************************************************************************************************
186 SUBROUTINE deallocate_arnoldi_data(arnoldi_data)
187 TYPE(arnoldi_data_type) :: arnoldi_data
188
189 TYPE(arnoldi_control_type), POINTER :: control
190
191 IF (has_d_real(arnoldi_data)) CALL deallocate_arnoldi_data_d(arnoldi_data)
192 IF (has_d_cmplx(arnoldi_data)) CALL deallocate_arnoldi_data_z(arnoldi_data)
193
194 control => get_control(arnoldi_data)
195 DEALLOCATE (control%selected_ind)
196 DEALLOCATE (control)
197
198 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 SUBROUTINE select_evals(arnoldi_data)
205 TYPE(arnoldi_data_type) :: arnoldi_data
206
207 IF (has_d_real(arnoldi_data)) CALL select_evals_d(arnoldi_data)
208 IF (has_d_cmplx(arnoldi_data)) CALL select_evals_z(arnoldi_data)
209
210 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 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 control => get_control(ar_data)
224 control%selection_crit = itype
225 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 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 control => get_control(arnoldi_data)
239 nrestart = control%nrestart
240
241 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 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 control => get_control(ar_data)
255 nval_out = control%nval_out
256
257 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 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 control => get_control(ar_data)
272 current_step = control%current_step
273
274 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 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 control => get_control(ar_data)
288 converged = control%converged
289
290 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 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 COMPLEX(real_8), DIMENSION(:), POINTER :: evals_d
304 INTEGER :: ev_ind
305 INTEGER, DIMENSION(:), POINTER :: selected_ind
306
307 IF (ind .GT. get_nval_out(ar_data)) &
308 cpabort('outside range of indexed evals')
309
310 selected_ind => get_sel_ind(ar_data)
311 ev_ind = selected_ind(ind)
312 IF (has_d_real(ar_data)) THEN
313 evals_d => get_evals_d(ar_data); eval_out = evals_d(ev_ind)
314 ELSE IF (has_d_cmplx(ar_data)) THEN
315 evals_d => get_evals_z(ar_data); eval_out = evals_d(ev_ind)
316 END IF
317
318 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 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 COMPLEX(real_8), DIMENSION(:), POINTER :: evals_d
331 INTEGER :: ev_ind, ind
332 INTEGER, DIMENSION(:), POINTER :: selected_ind
333
334 NULLIFY (evals_d)
335 IF (SIZE(eval_out) .LT. get_nval_out(ar_data)) &
336 cpabort('array for eval output too small')
337 selected_ind => get_sel_ind(ar_data)
338
339 IF (has_d_real(ar_data)) evals_d => get_evals_d(ar_data)
340 IF (has_d_cmplx(ar_data)) evals_d => get_evals_d(ar_data)
341
342 DO ind = 1, get_nval_out(ar_data)
343 ev_ind = selected_ind(ind)
344 IF (ASSOCIATED(evals_d)) eval_out(ind) = evals_d(ev_ind)
345 END DO
346
347 END SUBROUTINE get_all_selected_ritz_val
348
349! **************************************************************************************************
350!> \brief ...
351!> \param ar_data ...
352!> \param vector ...
353! **************************************************************************************************
354 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 control => get_control(ar_data)
361 control%has_initial_vector = .true.
362
363 IF (has_d_real(ar_data)) CALL set_initial_vector_d(ar_data, vector)
364 IF (has_d_cmplx(ar_data)) CALL set_initial_vector_z(ar_data, vector)
365
366 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! **************************************************************************************************
373!> \brief ...
374!> \param arnoldi_data ...
375! **************************************************************************************************
376 SUBROUTINE select_evals_d (arnoldi_data)
377 TYPE(arnoldi_data_type) :: arnoldi_data
378
379 INTEGER :: my_crit, last_el, my_ind, i
380 REAL(real_8) :: convergence
381 TYPE(arnoldi_data_d_type), POINTER :: ar_data
382 TYPE(arnoldi_control_type), POINTER :: control
383
384 control => get_control(arnoldi_data)
385 ar_data => get_data_d(arnoldi_data)
386
387 last_el = control%current_step
388 convergence = real(0.0, real_8)
389 my_crit = control%selection_crit
390 control%nval_out = min(control%nval_req, control%current_step)
391 SELECT CASE (my_crit)
392 ! minimum and maximum real eval
393 CASE (1)
394 CALL index_min_max_real_eval_d (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
395 ! n maximum real eval
396 CASE (2)
397 CALL index_nmax_real_eval_d (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
398 ! n minimum real eval
399 CASE (3)
400 CALL index_nmin_real_eval_d (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
401 CASE DEFAULT
402 cpabort("unknown selection index")
403 END SELECT
404 ! test whether we are converged
405 DO i = 1, control%nval_out
406 my_ind = control%selected_ind(i)
407 convergence = max(convergence, &
408 abs(ar_data%revec(last_el, my_ind)*ar_data%Hessenberg(last_el + 1, last_el)))
409 END DO
410 control%converged = convergence .LT. control%threshold
411
412 END SUBROUTINE select_evals_d
413
414! **************************************************************************************************
415!> \brief ...
416!> \param evals ...
417!> \param current_step ...
418!> \param selected_ind ...
419!> \param neval ...
420! **************************************************************************************************
421 SUBROUTINE index_min_max_real_eval_d (evals, current_step, selected_ind, neval)
422 COMPLEX(real_8), DIMENSION(:) :: evals
423 INTEGER, INTENT(IN) :: current_step
424 INTEGER, DIMENSION(:) :: selected_ind
425 INTEGER :: neval
426
427 INTEGER, DIMENSION(current_step) :: indexing
428 REAL(real_8), DIMENSION(current_step) :: tmp_array
429 INTEGER :: i
430
431 neval = 0
432 selected_ind = 0
433 tmp_array(1:current_step) = real(evals(1:current_step), real_8)
434 CALL sort(tmp_array, current_step, indexing)
435 DO i = 1, current_step
436 IF (abs(aimag(evals(indexing(i)))) < epsilon(0.0_real_8)) THEN
437 selected_ind(1) = indexing(i)
438 neval = neval + 1
439 EXIT
440 END IF
441 END DO
442 DO i = current_step, 1, -1
443 IF (abs(aimag(evals(indexing(i)))) < epsilon(0.0_real_8)) THEN
444 selected_ind(2) = indexing(i)
445 neval = neval + 1
446 EXIT
447 END IF
448 END DO
449
450 END SUBROUTINE index_min_max_real_eval_d
451
452! **************************************************************************************************
453!> \brief ...
454!> \param evals ...
455!> \param current_step ...
456!> \param selected_ind ...
457!> \param neval ...
458! **************************************************************************************************
459 SUBROUTINE index_nmax_real_eval_d (evals, current_step, selected_ind, neval)
460 COMPLEX(real_8), DIMENSION(:) :: evals
461 INTEGER, INTENT(IN) :: current_step
462 INTEGER, DIMENSION(:) :: selected_ind
463 INTEGER :: neval
464
465 INTEGER :: i, nlimit
466 INTEGER, DIMENSION(current_step) :: indexing
467 REAL(real_8), DIMENSION(current_step) :: tmp_array
468
469 nlimit = neval; neval = 0
470 selected_ind = 0
471 tmp_array(1:current_step) = real(evals(1:current_step), real_8)
472 CALL sort(tmp_array, current_step, indexing)
473 DO i = 1, current_step
474 IF (abs(aimag(evals(indexing(current_step + 1 - i)))) < epsilon(0.0_real_8)) THEN
475 selected_ind(i) = indexing(current_step + 1 - i)
476 neval = neval + 1
477 IF (neval == nlimit) EXIT
478 END IF
479 END DO
480
481 END SUBROUTINE index_nmax_real_eval_d
482
483! **************************************************************************************************
484!> \brief ...
485!> \param evals ...
486!> \param current_step ...
487!> \param selected_ind ...
488!> \param neval ...
489! **************************************************************************************************
490 SUBROUTINE index_nmin_real_eval_d (evals, current_step, selected_ind, neval)
491 COMPLEX(real_8), DIMENSION(:) :: evals
492 INTEGER, INTENT(IN) :: current_step
493 INTEGER, DIMENSION(:) :: selected_ind
494 INTEGER :: neval, nlimit
495
496 INTEGER :: i
497 INTEGER, DIMENSION(current_step) :: indexing
498 REAL(real_8), DIMENSION(current_step) :: tmp_array
499
500 nlimit = neval; neval = 0
501 selected_ind = 0
502 tmp_array(1:current_step) = real(evals(1:current_step), real_8)
503 CALL sort(tmp_array, current_step, indexing)
504 DO i = 1, current_step
505 IF (abs(aimag(evals(indexing(i)))) < epsilon(0.0_real_8)) THEN
506 selected_ind(i) = indexing(i)
507 neval = neval + 1
508 IF (neval == nlimit) EXIT
509 END IF
510 END DO
511
512 END SUBROUTINE index_nmin_real_eval_d
513
514! **************************************************************************************************
515!> \brief ...
516!> \param arnoldi_data ...
517!> \param matrix ...
518!> \param max_iter ...
519! **************************************************************************************************
520 SUBROUTINE setup_arnoldi_data_d (arnoldi_data, matrix, max_iter)
521 TYPE(arnoldi_data_type) :: arnoldi_data
522 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
523 INTEGER :: max_iter
524
525 INTEGER :: nrow_local
526 TYPE(arnoldi_data_d_type), POINTER :: ar_data
527
528 ALLOCATE (ar_data)
529 CALL dbcsr_get_info(matrix=matrix(1)%matrix, nfullrows_local=nrow_local)
530 ALLOCATE (ar_data%f_vec(nrow_local))
531 ALLOCATE (ar_data%x_vec(nrow_local))
532 ALLOCATE (ar_data%Hessenberg(max_iter + 1, max_iter))
533 ALLOCATE (ar_data%local_history(nrow_local, max_iter))
534
535 ALLOCATE (ar_data%evals(max_iter))
536 ALLOCATE (ar_data%revec(max_iter, max_iter))
537
538 CALL set_data_d (arnoldi_data, ar_data)
539
540 END SUBROUTINE setup_arnoldi_data_d
541
542! **************************************************************************************************
543!> \brief ...
544!> \param arnoldi_data ...
545! **************************************************************************************************
546 SUBROUTINE deallocate_arnoldi_data_d (arnoldi_data)
547 TYPE(arnoldi_data_type) :: arnoldi_data
548
549 TYPE(arnoldi_data_d_type), POINTER :: ar_data
550
551 ar_data => get_data_d(arnoldi_data)
552 IF (ASSOCIATED(ar_data%f_vec)) DEALLOCATE (ar_data%f_vec)
553 IF (ASSOCIATED(ar_data%x_vec)) DEALLOCATE (ar_data%x_vec)
554 IF (ASSOCIATED(ar_data%Hessenberg)) DEALLOCATE (ar_data%Hessenberg)
555 IF (ASSOCIATED(ar_data%local_history)) DEALLOCATE (ar_data%local_history)
556 IF (ASSOCIATED(ar_data%evals)) DEALLOCATE (ar_data%evals)
557 IF (ASSOCIATED(ar_data%revec)) DEALLOCATE (ar_data%revec)
558 DEALLOCATE (ar_data)
559
560 END SUBROUTINE deallocate_arnoldi_data_d
561
562! **************************************************************************************************
563!> \brief ...
564!> \param arnoldi_data ...
565!> \param ind ...
566!> \param matrix ...
567!> \param vector ...
568! **************************************************************************************************
569 SUBROUTINE get_selected_ritz_vector_d (arnoldi_data, ind, matrix, vector)
570 TYPE(arnoldi_data_type) :: arnoldi_data
571 INTEGER :: ind
572 TYPE(dbcsr_type) :: matrix
573 TYPE(dbcsr_type) :: vector
574
575 TYPE(arnoldi_data_d_type), POINTER :: ar_data
576 INTEGER :: vsize, myind, sspace_size, i
577 INTEGER, DIMENSION(:), POINTER :: selected_ind
578 COMPLEX(real_8), DIMENSION(:), ALLOCATABLE :: ritz_v
579 REAL(kind=real_8), DIMENSION(:), POINTER :: data_vec
580 TYPE(arnoldi_control_type), POINTER :: control
581
582 control => get_control(arnoldi_data)
583 selected_ind => get_sel_ind(arnoldi_data)
584 ar_data => get_data_d(arnoldi_data)
585 sspace_size = get_subsp_size(arnoldi_data)
586 vsize = SIZE(ar_data%f_vec)
587 myind = selected_ind(ind)
588 ALLOCATE (ritz_v(vsize))
589 ritz_v = cmplx(0.0, 0.0, real_8)
590
591 CALL dbcsr_release(vector)
592 CALL create_col_vec_from_matrix(vector, matrix, 1)
593 IF (control%local_comp) THEN
594 DO i = 1, sspace_size
595 ritz_v(:) = ritz_v(:) + ar_data%local_history(:, i)*ar_data%revec(i, myind)
596 END DO
597 data_vec => dbcsr_get_data_p(vector, select_data_type=0.0_real_8)
598 ! is a bit odd but ritz_v is always complex and matrix type determines where it goes
599 ! again I hope the user knows what is required
600 data_vec(1:vsize) = real(ritz_v(1:vsize), kind=real_8)
601 END IF
602
603 DEALLOCATE (ritz_v)
604
605 END SUBROUTINE get_selected_ritz_vector_d
606
607! **************************************************************************************************
608!> \brief ...
609!> \param arnoldi_data ...
610!> \param vector ...
611! **************************************************************************************************
612 SUBROUTINE set_initial_vector_d (arnoldi_data, vector)
613 TYPE(arnoldi_data_type) :: arnoldi_data
614 TYPE(dbcsr_type) :: vector
615
616 TYPE(arnoldi_data_d_type), POINTER :: ar_data
617 REAL(kind=real_8), DIMENSION(:), POINTER :: data_vec
618 INTEGER :: nrow_local, ncol_local
619 TYPE(arnoldi_control_type), POINTER :: control
620
621 control => get_control(arnoldi_data)
622
623 CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
624 ar_data => get_data_d(arnoldi_data)
625 data_vec => dbcsr_get_data_p(vector, select_data_type=0.0_real_8)
626 IF (nrow_local*ncol_local > 0) ar_data%f_vec(1:nrow_local) = data_vec(1:nrow_local)
627
628 END SUBROUTINE set_initial_vector_d
629
630! **************************************************************************************************
631!> \brief ...
632!> \param arnoldi_data ...
633! **************************************************************************************************
634 SUBROUTINE select_evals_z (arnoldi_data)
635 TYPE(arnoldi_data_type) :: arnoldi_data
636
637 INTEGER :: my_crit, last_el, my_ind, i
638 REAL(real_8) :: convergence
639 TYPE(arnoldi_data_z_type), POINTER :: ar_data
640 TYPE(arnoldi_control_type), POINTER :: control
641
642 control => get_control(arnoldi_data)
643 ar_data => get_data_z(arnoldi_data)
644
645 last_el = control%current_step
646 convergence = real(0.0, real_8)
647 my_crit = control%selection_crit
648 control%nval_out = min(control%nval_req, control%current_step)
649 SELECT CASE (my_crit)
650 ! minimum and maximum real eval
651 CASE (1)
652 CALL index_min_max_real_eval_z (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
653 ! n maximum real eval
654 CASE (2)
655 CALL index_nmax_real_eval_z (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
656 ! n minimum real eval
657 CASE (3)
658 CALL index_nmin_real_eval_z (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
659 CASE DEFAULT
660 cpabort("unknown selection index")
661 END SELECT
662 ! test whether we are converged
663 DO i = 1, control%nval_out
664 my_ind = control%selected_ind(i)
665 convergence = max(convergence, &
666 abs(ar_data%revec(last_el, my_ind)*ar_data%Hessenberg(last_el + 1, last_el)))
667 END DO
668 control%converged = convergence .LT. control%threshold
669
670 END SUBROUTINE select_evals_z
671
672! **************************************************************************************************
673!> \brief ...
674!> \param evals ...
675!> \param current_step ...
676!> \param selected_ind ...
677!> \param neval ...
678! **************************************************************************************************
679 SUBROUTINE index_min_max_real_eval_z (evals, current_step, selected_ind, neval)
680 COMPLEX(real_8), DIMENSION(:) :: evals
681 INTEGER, INTENT(IN) :: current_step
682 INTEGER, DIMENSION(:) :: selected_ind
683 INTEGER :: neval
684
685 INTEGER, DIMENSION(current_step) :: indexing
686 REAL(real_8), DIMENSION(current_step) :: tmp_array
687 INTEGER :: i
688
689 neval = 0
690 selected_ind = 0
691 tmp_array(1:current_step) = real(evals(1:current_step), real_8)
692 CALL sort(tmp_array, current_step, indexing)
693 DO i = 1, current_step
694 IF (abs(aimag(evals(indexing(i)))) < epsilon(0.0_real_8)) THEN
695 selected_ind(1) = indexing(i)
696 neval = neval + 1
697 EXIT
698 END IF
699 END DO
700 DO i = current_step, 1, -1
701 IF (abs(aimag(evals(indexing(i)))) < epsilon(0.0_real_8)) THEN
702 selected_ind(2) = indexing(i)
703 neval = neval + 1
704 EXIT
705 END IF
706 END DO
707
708 END SUBROUTINE index_min_max_real_eval_z
709
710! **************************************************************************************************
711!> \brief ...
712!> \param evals ...
713!> \param current_step ...
714!> \param selected_ind ...
715!> \param neval ...
716! **************************************************************************************************
717 SUBROUTINE index_nmax_real_eval_z (evals, current_step, selected_ind, neval)
718 COMPLEX(real_8), DIMENSION(:) :: evals
719 INTEGER, INTENT(IN) :: current_step
720 INTEGER, DIMENSION(:) :: selected_ind
721 INTEGER :: neval
722
723 INTEGER :: i, nlimit
724 INTEGER, DIMENSION(current_step) :: indexing
725 REAL(real_8), DIMENSION(current_step) :: tmp_array
726
727 nlimit = neval; neval = 0
728 selected_ind = 0
729 tmp_array(1:current_step) = real(evals(1:current_step), real_8)
730 CALL sort(tmp_array, current_step, indexing)
731 DO i = 1, current_step
732 IF (abs(aimag(evals(indexing(current_step + 1 - i)))) < epsilon(0.0_real_8)) THEN
733 selected_ind(i) = indexing(current_step + 1 - i)
734 neval = neval + 1
735 IF (neval == nlimit) EXIT
736 END IF
737 END DO
738
739 END SUBROUTINE index_nmax_real_eval_z
740
741! **************************************************************************************************
742!> \brief ...
743!> \param evals ...
744!> \param current_step ...
745!> \param selected_ind ...
746!> \param neval ...
747! **************************************************************************************************
748 SUBROUTINE index_nmin_real_eval_z (evals, current_step, selected_ind, neval)
749 COMPLEX(real_8), DIMENSION(:) :: evals
750 INTEGER, INTENT(IN) :: current_step
751 INTEGER, DIMENSION(:) :: selected_ind
752 INTEGER :: neval, nlimit
753
754 INTEGER :: i
755 INTEGER, DIMENSION(current_step) :: indexing
756 REAL(real_8), DIMENSION(current_step) :: tmp_array
757
758 nlimit = neval; neval = 0
759 selected_ind = 0
760 tmp_array(1:current_step) = real(evals(1:current_step), real_8)
761 CALL sort(tmp_array, current_step, indexing)
762 DO i = 1, current_step
763 IF (abs(aimag(evals(indexing(i)))) < epsilon(0.0_real_8)) THEN
764 selected_ind(i) = indexing(i)
765 neval = neval + 1
766 IF (neval == nlimit) EXIT
767 END IF
768 END DO
769
770 END SUBROUTINE index_nmin_real_eval_z
771
772! **************************************************************************************************
773!> \brief ...
774!> \param arnoldi_data ...
775!> \param matrix ...
776!> \param max_iter ...
777! **************************************************************************************************
778 SUBROUTINE setup_arnoldi_data_z (arnoldi_data, matrix, max_iter)
779 TYPE(arnoldi_data_type) :: arnoldi_data
780 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix
781 INTEGER :: max_iter
782
783 INTEGER :: nrow_local
784 TYPE(arnoldi_data_z_type), POINTER :: ar_data
785
786 ALLOCATE (ar_data)
787 CALL dbcsr_get_info(matrix=matrix(1)%matrix, nfullrows_local=nrow_local)
788 ALLOCATE (ar_data%f_vec(nrow_local))
789 ALLOCATE (ar_data%x_vec(nrow_local))
790 ALLOCATE (ar_data%Hessenberg(max_iter + 1, max_iter))
791 ALLOCATE (ar_data%local_history(nrow_local, max_iter))
792
793 ALLOCATE (ar_data%evals(max_iter))
794 ALLOCATE (ar_data%revec(max_iter, max_iter))
795
796 CALL set_data_z (arnoldi_data, ar_data)
797
798 END SUBROUTINE setup_arnoldi_data_z
799
800! **************************************************************************************************
801!> \brief ...
802!> \param arnoldi_data ...
803! **************************************************************************************************
804 SUBROUTINE deallocate_arnoldi_data_z (arnoldi_data)
805 TYPE(arnoldi_data_type) :: arnoldi_data
806
807 TYPE(arnoldi_data_z_type), POINTER :: ar_data
808
809 ar_data => get_data_z(arnoldi_data)
810 IF (ASSOCIATED(ar_data%f_vec)) DEALLOCATE (ar_data%f_vec)
811 IF (ASSOCIATED(ar_data%x_vec)) DEALLOCATE (ar_data%x_vec)
812 IF (ASSOCIATED(ar_data%Hessenberg)) DEALLOCATE (ar_data%Hessenberg)
813 IF (ASSOCIATED(ar_data%local_history)) DEALLOCATE (ar_data%local_history)
814 IF (ASSOCIATED(ar_data%evals)) DEALLOCATE (ar_data%evals)
815 IF (ASSOCIATED(ar_data%revec)) DEALLOCATE (ar_data%revec)
816 DEALLOCATE (ar_data)
817
818 END SUBROUTINE deallocate_arnoldi_data_z
819
820! **************************************************************************************************
821!> \brief ...
822!> \param arnoldi_data ...
823!> \param ind ...
824!> \param matrix ...
825!> \param vector ...
826! **************************************************************************************************
827 SUBROUTINE get_selected_ritz_vector_z (arnoldi_data, ind, matrix, vector)
828 TYPE(arnoldi_data_type) :: arnoldi_data
829 INTEGER :: ind
830 TYPE(dbcsr_type) :: matrix
831 TYPE(dbcsr_type) :: vector
832
833 TYPE(arnoldi_data_z_type), POINTER :: ar_data
834 INTEGER :: vsize, myind, sspace_size, i
835 INTEGER, DIMENSION(:), POINTER :: selected_ind
836 COMPLEX(real_8), DIMENSION(:), ALLOCATABLE :: ritz_v
837 COMPLEX(kind=real_8), DIMENSION(:), POINTER :: data_vec
838 TYPE(arnoldi_control_type), POINTER :: control
839
840 control => get_control(arnoldi_data)
841 selected_ind => get_sel_ind(arnoldi_data)
842 ar_data => get_data_z(arnoldi_data)
843 sspace_size = get_subsp_size(arnoldi_data)
844 vsize = SIZE(ar_data%f_vec)
845 myind = selected_ind(ind)
846 ALLOCATE (ritz_v(vsize))
847 ritz_v = cmplx(0.0, 0.0, real_8)
848
849 CALL dbcsr_release(vector)
850 CALL create_col_vec_from_matrix(vector, matrix, 1)
851 IF (control%local_comp) THEN
852 DO i = 1, sspace_size
853 ritz_v(:) = ritz_v(:) + ar_data%local_history(:, i)*ar_data%revec(i, myind)
854 END DO
855 data_vec => dbcsr_get_data_p(vector, select_data_type=cmplx(0.0, 0.0, real_8))
856 ! is a bit odd but ritz_v is always complex and matrix type determines where it goes
857 ! again I hope the user knows what is required
858 data_vec(1:vsize) = cmplx(ritz_v(1:vsize), kind=real_8)
859 END IF
860
861 DEALLOCATE (ritz_v)
862
863 END SUBROUTINE get_selected_ritz_vector_z
864
865! **************************************************************************************************
866!> \brief ...
867!> \param arnoldi_data ...
868!> \param vector ...
869! **************************************************************************************************
870 SUBROUTINE set_initial_vector_z (arnoldi_data, vector)
871 TYPE(arnoldi_data_type) :: arnoldi_data
872 TYPE(dbcsr_type) :: vector
873
874 TYPE(arnoldi_data_z_type), POINTER :: ar_data
875 COMPLEX(kind=real_8), DIMENSION(:), POINTER :: data_vec
876 INTEGER :: nrow_local, ncol_local
877 TYPE(arnoldi_control_type), POINTER :: control
878
879 control => get_control(arnoldi_data)
880
881 CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
882 ar_data => get_data_z(arnoldi_data)
883 data_vec => dbcsr_get_data_p(vector, select_data_type=cmplx(0.0, 0.0, real_8))
884 IF (nrow_local*ncol_local > 0) ar_data%f_vec(1:nrow_local) = data_vec(1:nrow_local)
885
886 END SUBROUTINE set_initial_vector_z
887
888
889END MODULE arnoldi_data_methods
890
arnoldi_data_methods
The methods which allow to analyze and manipulate the arnoldi procedure The main routine and this sho...
Definition arnoldi_data_methods.F:17
arnoldi_data_methods::select_evals
subroutine, public select_evals(arnoldi_data)
perform the selection of eigenvalues, fills the selected_ind array
Definition arnoldi_data_methods.F:205
arnoldi_data_methods::setup_arnoldi_data
subroutine, public setup_arnoldi_data(arnoldi_data, matrix, max_iter, threshold, selection_crit, nval_request, nrestarts, generalized_ev, iram)
This routine sets the environment for the arnoldi iteration and the krylov subspace creation....
Definition arnoldi_data_methods.F:64
arnoldi_data_methods::get_selected_ritz_val
complex(real_8) function, public get_selected_ritz_val(ar_data, ind)
get a single specific Ritz value from the set of selected
Definition arnoldi_data_methods.F:299
arnoldi_data_methods::get_selected_ritz_vector
subroutine, public get_selected_ritz_vector(arnoldi_data, ind, matrix, vector)
Deallocate the data in arnoldi_data.
Definition arnoldi_data_methods.F:173
arnoldi_data_methods::deallocate_arnoldi_data
subroutine, public deallocate_arnoldi_data(arnoldi_data)
Deallocate the data in arnoldi_data.
Definition arnoldi_data_methods.F:187
arnoldi_data_methods::arnoldi_is_converged
logical function, public arnoldi_is_converged(ar_data)
Find out whether the method with the current search criterion is converged.
Definition arnoldi_data_methods.F:282
arnoldi_data_methods::set_arnoldi_initial_vector
subroutine, public set_arnoldi_initial_vector(ar_data, vector)
...
Definition arnoldi_data_methods.F:355
arnoldi_data_methods::get_nrestart
integer function, public get_nrestart(arnoldi_data)
returns the number of restarts allowed for arnoldi
Definition arnoldi_data_methods.F:233
arnoldi_types
collection of types used in arnoldi
Definition arnoldi_types.F:15
arnoldi_types::set_data_s
subroutine, public set_data_s(ar_data, data_s)
...
Definition arnoldi_types.F:206
arnoldi_types::set_control
subroutine, public set_control(ar_data, control)
...
Definition arnoldi_types.F:116
arnoldi_types::get_data_c
type(arnoldi_data_c_type) function, pointer, public get_data_c(ar_data)
...
Definition arnoldi_types.F:180
arnoldi_types::set_data_c
subroutine, public set_data_c(ar_data, data_c)
...
Definition arnoldi_types.F:219
arnoldi_types::get_evals_s
complex(real_4) function, dimension(:), pointer, public get_evals_s(ar_data)
...
Definition arnoldi_types.F:323
arnoldi_types::get_data_s
type(arnoldi_data_s_type) function, pointer, public get_data_s(ar_data)
...
Definition arnoldi_types.F:154
arnoldi_types::get_data_z
type(arnoldi_data_z_type) function, pointer, public get_data_z(ar_data)
...
Definition arnoldi_types.F:167
arnoldi_types::set_data_d
subroutine, public set_data_d(ar_data, data_d)
...
Definition arnoldi_types.F:193
arnoldi_types::has_d_real
logical function, public has_d_real(ar_data)
...
Definition arnoldi_types.F:258
arnoldi_types::set_data_z
subroutine, public set_data_z(ar_data, data_z)
...
Definition arnoldi_types.F:232
arnoldi_types::get_evals_c
complex(real_4) function, dimension(:), pointer, public get_evals_c(ar_data)
...
Definition arnoldi_types.F:349
arnoldi_types::get_data_d
type(arnoldi_data_d_type) function, pointer, public get_data_d(ar_data)
...
Definition arnoldi_types.F:141
arnoldi_types::get_control
type(arnoldi_control_type) function, pointer, public get_control(ar_data)
...
Definition arnoldi_types.F:245
arnoldi_types::has_d_cmplx
elemental logical function, public has_d_cmplx(ar_data)
...
Definition arnoldi_types.F:284
arnoldi_types::get_evals_d
complex(real_8) function, dimension(:), pointer, public get_evals_d(ar_data)
...
Definition arnoldi_types.F:310
arnoldi_types::get_evals_z
complex(real_8) function, dimension(:), pointer, public get_evals_z(ar_data)
...
Definition arnoldi_types.F:336
arnoldi_types::get_sel_ind
integer function, dimension(:), pointer, public get_sel_ind(ar_data)
...
Definition arnoldi_types.F:128
dbcsr_vector
operations for skinny matrices/vectors expressed in dbcsr form
Definition dbcsr_vector.F:15
dbcsr_vector::create_col_vec_from_matrix
subroutine, public create_col_vec_from_matrix(dbcsr_vec, matrix, ncol)
creates a dbcsr col vector like object which lives on proc_col 0 and has the same row dist as the tem...
Definition dbcsr_vector.F:109
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::real_8
integer, parameter, public real_8
Definition kinds.F:41
util
All kind of helpful little routines.
Definition util.F:14
arnoldi_types::arnoldi_control_type
Definition arnoldi_types.F:26
arnoldi_types::arnoldi_data_c_type
Definition arnoldi_types.F:74
arnoldi_types::arnoldi_data_d_type
Definition arnoldi_types.F:44
arnoldi_types::arnoldi_data_s_type
Definition arnoldi_types.F:54
arnoldi_types::arnoldi_data_type
Definition arnoldi_types.F:84
arnoldi_types::arnoldi_data_z_type
Definition arnoldi_types.F:64