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
31 #include "../base/base_uses.f90"
63 nval_request, nrestarts, generalized_ev, iram)
65 TYPE(dbcsr_p_type),
DIMENSION(:) :: matrix
68 INTEGER :: selection_crit, nval_request, nrestarts
69 LOGICAL :: generalized_ev, iram
71 CALL setup_arnoldi_control(arnoldi_data, matrix, max_iter, threshold, selection_crit, &
72 nval_request, nrestarts, generalized_ev, iram)
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)
95 SUBROUTINE setup_arnoldi_control(arnoldi_data, matrix, max_iter, threshold, selection_crit, nval_request, &
96 nrestarts, generalized_ev, iram)
98 TYPE(dbcsr_p_type),
DIMENSION(:) :: matrix
101 INTEGER :: selection_crit, nval_request, nrestarts
102 LOGICAL :: generalized_ev, iram
104 INTEGER :: pcol_handle, group_handle
105 LOGICAL :: subgroups_defined
107 TYPE(dbcsr_distribution_type) :: distri
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)
119 CALL control%mp_group%set_handle(group_handle)
120 CALL control%pcol_group%set_handle(pcol_handle)
122 IF (.NOT. subgroups_defined) &
123 cpabort(
"arnoldi only with subgroups")
125 control%symmetric = .false.
127 IF (
SIZE(matrix) == 1) &
128 control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric
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.
139 control%nrestart = nrestarts
140 control%generalized_ev = generalized_ev
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')
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)')
160 ALLOCATE (control%selected_ind(max_iter))
163 END SUBROUTINE setup_arnoldi_control
175 TYPE(dbcsr_type) :: matrix, vector
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)
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)
195 DEALLOCATE (control%selected_ind)
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)
217 SUBROUTINE set_eval_selection(ar_data, itype)
224 control%selection_crit = itype
225 END SUBROUTINE set_eval_selection
239 nrestart = control%nrestart
248 FUNCTION get_nval_out(ar_data)
RESULT(nval_out)
255 nval_out = control%nval_out
257 END FUNCTION get_nval_out
265 FUNCTION get_subsp_size(ar_data)
RESULT(current_step)
267 INTEGER :: current_step
272 current_step = control%current_step
274 END FUNCTION get_subsp_size
288 converged = control%converged
301 COMPLEX(real_8) :: eval_out
303 COMPLEX(real_8),
DIMENSION(:),
POINTER :: evals_d
305 INTEGER,
DIMENSION(:),
POINTER :: selected_ind
307 IF (ind .GT. get_nval_out(ar_data)) &
308 cpabort(
'outside range of indexed evals')
311 ev_ind = selected_ind(ind)
313 evals_d =>
get_evals_d(ar_data); eval_out = evals_d(ev_ind)
315 evals_d =>
get_evals_z(ar_data); eval_out = evals_d(ev_ind)
326 SUBROUTINE get_all_selected_ritz_val(ar_data, eval_out)
328 COMPLEX(real_8),
DIMENSION(:) :: eval_out
330 COMPLEX(real_8),
DIMENSION(:),
POINTER :: evals_d
331 INTEGER :: ev_ind, ind
332 INTEGER,
DIMENSION(:),
POINTER :: selected_ind
335 IF (
SIZE(eval_out) .LT. get_nval_out(ar_data)) &
336 cpabort(
'array for eval output too small')
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)
347 END SUBROUTINE get_all_selected_ritz_val
356 TYPE(dbcsr_type) :: vector
361 control%has_initial_vector = .true.
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)
376 SUBROUTINE select_evals_d (arnoldi_data)
379 INTEGER :: my_crit, last_el, my_ind, i
380 REAL(
real_8) :: convergence
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)
394 CALL index_min_max_real_eval_d (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
397 CALL index_nmax_real_eval_d (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
400 CALL index_nmin_real_eval_d (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
402 cpabort(
"unknown selection index")
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)))
410 control%converged = convergence .LT. control%threshold
412 END SUBROUTINE select_evals_d
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
427 INTEGER,
DIMENSION(current_step) :: indexing
428 REAL(
real_8),
DIMENSION(current_step) :: tmp_array
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)
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)
450 END SUBROUTINE index_min_max_real_eval_d
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
466 INTEGER,
DIMENSION(current_step) :: indexing
467 REAL(
real_8),
DIMENSION(current_step) :: tmp_array
469 nlimit = neval; neval = 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)
477 IF (neval == nlimit)
EXIT
481 END SUBROUTINE index_nmax_real_eval_d
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
497 INTEGER,
DIMENSION(current_step) :: indexing
498 REAL(
real_8),
DIMENSION(current_step) :: tmp_array
500 nlimit = neval; neval = 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)
508 IF (neval == nlimit)
EXIT
512 END SUBROUTINE index_nmin_real_eval_d
520 SUBROUTINE setup_arnoldi_data_d (arnoldi_data, matrix, max_iter)
522 TYPE(dbcsr_p_type),
DIMENSION(:) :: matrix
525 INTEGER :: nrow_local
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))
535 ALLOCATE (ar_data%evals(max_iter))
536 ALLOCATE (ar_data%revec(max_iter, max_iter))
540 END SUBROUTINE setup_arnoldi_data_d
546 SUBROUTINE deallocate_arnoldi_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)
560 END SUBROUTINE deallocate_arnoldi_data_d
569 SUBROUTINE get_selected_ritz_vector_d (arnoldi_data, ind, matrix, vector)
572 TYPE(dbcsr_type) :: matrix
573 TYPE(dbcsr_type) :: vector
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
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)
591 CALL dbcsr_release(vector)
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)
597 data_vec => dbcsr_get_data_p(vector, select_data_type=0.0_real_8)
600 data_vec(1:vsize) = real(ritz_v(1:vsize), kind=
real_8)
605 END SUBROUTINE get_selected_ritz_vector_d
612 SUBROUTINE set_initial_vector_d (arnoldi_data, vector)
614 TYPE(dbcsr_type) :: vector
617 REAL(kind=
real_8),
DIMENSION(:),
POINTER :: data_vec
618 INTEGER :: nrow_local, ncol_local
623 CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
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)
628 END SUBROUTINE set_initial_vector_d
634 SUBROUTINE select_evals_z (arnoldi_data)
637 INTEGER :: my_crit, last_el, my_ind, i
638 REAL(
real_8) :: convergence
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)
652 CALL index_min_max_real_eval_z (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
655 CALL index_nmax_real_eval_z (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
658 CALL index_nmin_real_eval_z (ar_data%evals, control%current_step, control%selected_ind, control%nval_out)
660 cpabort(
"unknown selection index")
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)))
668 control%converged = convergence .LT. control%threshold
670 END SUBROUTINE select_evals_z
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
685 INTEGER,
DIMENSION(current_step) :: indexing
686 REAL(
real_8),
DIMENSION(current_step) :: tmp_array
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)
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)
708 END SUBROUTINE index_min_max_real_eval_z
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
724 INTEGER,
DIMENSION(current_step) :: indexing
725 REAL(
real_8),
DIMENSION(current_step) :: tmp_array
727 nlimit = neval; neval = 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)
735 IF (neval == nlimit)
EXIT
739 END SUBROUTINE index_nmax_real_eval_z
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
755 INTEGER,
DIMENSION(current_step) :: indexing
756 REAL(
real_8),
DIMENSION(current_step) :: tmp_array
758 nlimit = neval; neval = 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)
766 IF (neval == nlimit)
EXIT
770 END SUBROUTINE index_nmin_real_eval_z
778 SUBROUTINE setup_arnoldi_data_z (arnoldi_data, matrix, max_iter)
780 TYPE(dbcsr_p_type),
DIMENSION(:) :: matrix
783 INTEGER :: nrow_local
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))
793 ALLOCATE (ar_data%evals(max_iter))
794 ALLOCATE (ar_data%revec(max_iter, max_iter))
798 END SUBROUTINE setup_arnoldi_data_z
804 SUBROUTINE deallocate_arnoldi_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)
818 END SUBROUTINE deallocate_arnoldi_data_z
827 SUBROUTINE get_selected_ritz_vector_z (arnoldi_data, ind, matrix, vector)
830 TYPE(dbcsr_type) :: matrix
831 TYPE(dbcsr_type) :: vector
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
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)
849 CALL dbcsr_release(vector)
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)
855 data_vec => dbcsr_get_data_p(vector, select_data_type=cmplx(0.0, 0.0,
real_8))
858 data_vec(1:vsize) = cmplx(ritz_v(1:vsize), kind=
real_8)
863 END SUBROUTINE get_selected_ritz_vector_z
870 SUBROUTINE set_initial_vector_z (arnoldi_data, vector)
872 TYPE(dbcsr_type) :: vector
875 COMPLEX(kind=real_8),
DIMENSION(:),
POINTER :: data_vec
876 INTEGER :: nrow_local, ncol_local
881 CALL dbcsr_get_info(matrix=vector, nfullrows_local=nrow_local, nfullcols_local=ncol_local)
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)
886 END SUBROUTINE set_initial_vector_z
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
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
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...
Defines the basic variable types.
integer, parameter, public real_8
All kind of helpful little routines.