(git:0de0cc2)
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 
18  USE arnoldi_types, ONLY: &
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
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 
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  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 
889 END MODULE arnoldi_data_methods
890 
The methods which allow to analyze and manipulate the arnoldi procedure The main routine and this sho...
subroutine, public select_evals(arnoldi_data)
perform the selection of eigenvalues, fills the selected_ind array
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....
complex(real_8) function, public get_selected_ritz_val(ar_data, ind)
get a single specific Ritz value from the set of selected
subroutine, public get_selected_ritz_vector(arnoldi_data, ind, matrix, vector)
Deallocate the data in arnoldi_data.
subroutine, public deallocate_arnoldi_data(arnoldi_data)
Deallocate the data in arnoldi_data.
logical function, public arnoldi_is_converged(ar_data)
Find out whether the method with the current search criterion is converged.
subroutine, public set_arnoldi_initial_vector(ar_data, vector)
...
integer function, public get_nrestart(arnoldi_data)
returns the number of restarts allowed for arnoldi
collection of types used in arnoldi
Definition: arnoldi_types.F:15
subroutine, public set_data_s(ar_data, data_s)
...
subroutine, public set_control(ar_data, control)
...
type(arnoldi_data_c_type) function, pointer, public get_data_c(ar_data)
...
subroutine, public set_data_c(ar_data, data_c)
...
complex(real_4) function, dimension(:), pointer, public get_evals_s(ar_data)
...
type(arnoldi_data_s_type) function, pointer, public get_data_s(ar_data)
...
type(arnoldi_data_z_type) function, pointer, public get_data_z(ar_data)
...
subroutine, public set_data_d(ar_data, data_d)
...
logical function, public has_d_real(ar_data)
...
subroutine, public set_data_z(ar_data, data_z)
...
complex(real_4) function, dimension(:), pointer, public get_evals_c(ar_data)
...
type(arnoldi_data_d_type) function, pointer, public get_data_d(ar_data)
...
type(arnoldi_control_type) function, pointer, public get_control(ar_data)
...
elemental logical function, public has_d_cmplx(ar_data)
...
complex(real_8) function, dimension(:), pointer, public get_evals_d(ar_data)
...
complex(real_8) function, dimension(:), pointer, public get_evals_z(ar_data)
...
integer function, dimension(:), pointer, public get_sel_ind(ar_data)
...
operations for skinny matrices/vectors expressed in dbcsr form
Definition: dbcsr_vector.F:15
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public real_8
Definition: kinds.F:41
All kind of helpful little routines.
Definition: util.F:14