Skip to content

Commit

Permalink
XAS_TDP| Removed SYEVR diagonalization option (and all related code)
Browse files Browse the repository at this point in the history
  • Loading branch information
@abussy @pseewald
abussy authored and pseewald committed Dec 10, 2019
1 parent 73cb223 commit abc627f
Show file tree
Hide file tree
Showing 14 changed files with 57 additions and 451 deletions.
259 changes: 0 additions & 259 deletions src/fm/cp_fm_diag.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,6 @@ MODULE cp_fm_diag
cp_fm_syevx, &
cp_fm_geeig, &
cp_fm_geeig_canon, &
cp_fm_geeig_syevr, &
diag_init, &
diag_finalize

Expand Down Expand Up @@ -721,201 +720,6 @@ SUBROUTINE cp_fm_syevx(matrix, eigenvectors, eigenvalues, neig, work_syevx)

END SUBROUTINE cp_fm_syevx

! **************************************************************************************************
!> \brief computes selected eigenvalues and, optionally, eigenvectors of
!> a real symmetric matrix A distributed in 2D blockcyclic format by
!> calling the recommended sequence of ScaLAPACK routines.
!> WARNING: PDSEYVR is not always stable and might fail, especially for small matrices
!> The typical error message comes from PDORMTR with illegal value on entry
!> ScalaPACK version 2+ is required
!>
!> \param matrix ...
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \param ilow index of lowest desired eigenvalue
!> \param iup index of highest desired eigenvalue
!> \param vlow the lowest value from which to start calculating eigenvalues
!> \param vup the highest value up to which eigenvalues are computed
!> \note matrix is supposed to be in upper triangular form, and overwritten by this routine
!> subsets of eigenvalues/vectors can be selected by
!> specifying a range of values or a range of indices for the desired eigenvalues.
!> Reintroduced after removal because needed for XAS_TDP A. Bussy (07.2019)
! **************************************************************************************************
SUBROUTINE cp_fm_syevr(matrix, eigenvectors, eigenvalues, ilow, iup, vlow, vup)

TYPE(cp_fm_type), POINTER :: matrix
TYPE(cp_fm_type), POINTER, OPTIONAL :: eigenvectors
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
INTEGER, INTENT(IN), OPTIONAL :: ilow, iup
REAL(KIND=dp), INTENT(IN), OPTIONAL :: vlow, vup

CHARACTER(LEN=*), PARAMETER :: routineN = "cp_fm_syevr", &
routineP = moduleN//":"//routineN

CHARACTER(LEN=1) :: job_type, range_type
INTEGER :: handle, ilow_local, &
info, iup_local, &
lwork, liwork, mypcol, &
myprow, n, neig
REAL(KIND=dp) :: vlow_local, vup_local
INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
LOGICAL :: needs_evecs
REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: w, work
REAL(KIND=dp), DIMENSION(:, :), POINTER :: a, z

TYPE(cp_blacs_env_type), POINTER :: context

REAL(KIND=dp), EXTERNAL :: dlamch

#if defined(__SCALAPACK)
INTEGER, DIMENSION(9) :: desca, descz
INTEGER :: m, nz
#else
INTEGER :: m, nb
REAL(dp) :: abstol
INTEGER, DIMENSION(:), ALLOCATABLE :: ifail
INTEGER, EXTERNAL :: ilaenv
#endif

! by default all
n = matrix%matrix_struct%nrow_global
neig = n
iup_local = n
ilow_local = 1
range_type = "I"
IF (PRESENT(ilow) .AND. PRESENT(iup)) THEN
neig = iup - ilow + 1
iup_local = iup
ilow_local = ilow
ELSE IF (PRESENT(vlow) .AND. PRESENT(vup)) THEN
vlow_local = vlow
vup_local = vup
range_type = "V"
END IF
IF (neig <= 0) RETURN

CALL timeset(routineN, handle)

needs_evecs = PRESENT(eigenvectors)

n = matrix%matrix_struct%nrow_global

! set scalapack job type
IF (needs_evecs) THEN
job_type = "V"
ELSE
job_type = "N"
ENDIF

context => matrix%matrix_struct%context
myprow = context%mepos(1)
mypcol = context%mepos(2)

ALLOCATE (w(n))

eigenvalues(:) = 0.0_dp

#if defined(__SCALAPACK)

IF (matrix%matrix_struct%nrow_block /= matrix%matrix_struct%ncol_block) THEN
CPABORT("")
END IF

a => matrix%local_data
desca(:) = matrix%matrix_struct%descriptor(:)

IF (needs_evecs) THEN
z => eigenvectors%local_data
descz(:) = eigenvectors%matrix_struct%descriptor(:)
ELSE
! z will not be referenced
z => matrix%local_data
descz = desca
ENDIF

! First Call: Determine the needed work_space
lwork = -1
ALLOCATE (work(5*n))
ALLOCATE (iwork(6*n))

CALL pdsyevr(job_type, range_type, 'U', n, a, 1, 1, desca, vlow, vup, ilow_local, iup_local, &
m, nz, w(1), z, 1, 1, descz, work, lwork, iwork, liwork, info)

lwork = INT(work(1))
lwork = NINT(work(1) + 300000)
liwork = iwork(1)
IF (lwork > SIZE(work, 1)) THEN
DEALLOCATE (work)
ALLOCATE (work(lwork))
END IF
IF (liwork > SIZE(iwork, 1)) THEN
DEALLOCATE (iwork)
ALLOCATE (iwork(liwork))
END IF

!Second call: solve the eigenvalue problem
info = 0

CALL pdsyevr(job_type, range_type, 'U', n, a, 1, 1, desca, vlow, vup, ilow_local, iup_local, &
m, nz, w(1), z, 1, 1, descz, work, lwork, iwork, liwork, info)

IF (info > 0) THEN
WRITE (*, *) 'Processor ', myprow, mypcol, ': Error! INFO code = ', INFO
END IF
CPASSERT(info == 0)

! Release work storage
DEALLOCATE (iwork)
DEALLOCATE (work)

#else

a => matrix%local_data
IF (needs_evecs) THEN
z => eigenvectors%local_data
ELSE
! z will not be referenced
z => matrix%local_data
ENDIF

! Get the optimal work storage size

nb = MAX(ilaenv(1, "DSYTRD", "U", n, -1, -1, -1), &
ilaenv(1, "DORMTR", "U", n, -1, -1, -1))

lwork = MAX((nb + 3)*n, 8*n) + n ! sun bug fix
liwork = 5*n

ALLOCATE (ifail(n))
ifail = 0

ALLOCATE (iwork(liwork))
ALLOCATE (work(lwork))

! target the most accurate calculation of the eigenvalues
abstol = 2.0_dp*dlamch("S")

info = 0
CALL dsyevx(job_type, range_type, "U", n, a(1, 1), n, vlow, vup, ilow_local, iup_local, &
abstol, m, w, z(1, 1), n, work(1), lwork, iwork(1), ifail(1), info)

! Error handling
CPASSERT(info == 0)

! Release work storage
DEALLOCATE (iwork)
DEALLOCATE (work)

#endif

!put eigenvalues from index 1
eigenvalues(1:m) = w(1:m)
DEALLOCATE (w)

CALL timestop(handle)

END SUBROUTINE cp_fm_syevr

! **************************************************************************************************
!> \brief ...
!> \param matrix ...
Expand Down Expand Up @@ -1356,69 +1160,6 @@ SUBROUTINE cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
END SUBROUTINE cp_fm_geeig
! **************************************************************************************************
!> \brief A geeig solver using Cholesky decomposition and syevr so that only a couple of eigenvalues
!> can be computed. Can be considerably cheaper than the standard diagonalization for very
!> large matrices. WARNING: PDSYEVR is not always stable and might produce errors,
!> especially for smaller matrices (typically from inputs to PDORMTR)
!> amatrix ...
!> \param amatrix ...
!> \param bmatrix ...
!> \param eigenvectors ...
!> \param eigenvalues ...
!> \param work ...
!> \param ilow index of the first eigenvalue to compute
!> \param iup index of the last eigenvalue to compute
!> \param vlow lowest value from which eigenvalues are sought
!> \param vup upper value up to which eigenvalues are sought
!> \author A.Bussy (07.2019), originally for application in XAS_TDP
!> \note Blatantly copied from cp_fm_geeig. The eigenvalues are the first elements of the
!> `eigenvalues` array and the eigenvectors are in the first columns of `eigenvectors`
! **************************************************************************************************
SUBROUTINE cp_fm_geeig_syevr(amatrix, bmatrix, eigenvectors, eigenvalues, work, ilow, iup, &
vlow, vup)
TYPE(cp_fm_type), POINTER :: amatrix, bmatrix, eigenvectors
REAL(dp), DIMENSION(:) :: eigenvalues
TYPE(cp_fm_type), POINTER :: work
INTEGER, INTENT(IN), OPTIONAL :: ilow, iup
REAL(dp), INTENT(IN), OPTIONAL :: vlow, vup
CHARACTER(len=*), PARAMETER :: routineN = 'cp_fm_geeig_syevr', &
routineP = moduleN//':'//routineN
INTEGER :: nao
LOGICAL :: range_idx, range_val
range_idx = .FALSE.; range_val = .FALSE.
IF (PRESENT(ilow) .AND. PRESENT(iup)) THEN
range_idx = .TRUE.
ELSE IF (PRESENT(vlow) .AND. PRESENT(vup)) THEN
range_val = .TRUE.
END IF
CALL cp_fm_get_info(amatrix, nrow_global=nao)
! Cholesky decompose S=U(T)U
CALL cp_fm_cholesky_decompose(bmatrix)
! Invert to get U^(-1)
CALL cp_fm_triangular_invert(bmatrix)
! Reduce to get U^(-T) * H * U^(-1)
CALL cp_fm_triangular_multiply(bmatrix, amatrix, side="R")
CALL cp_fm_triangular_multiply(bmatrix, amatrix, transpose_tr=.TRUE.)
! Diagonalize
IF (range_idx) THEN
CALL cp_fm_syevr(amatrix, work, eigenvalues, ilow=ilow, iup=iup)
ELSE IF (range_val) THEN
CALL cp_fm_syevr(amatrix, work, eigenvalues, vlow=vlow, vup=vup)
ELSE
CALL cp_fm_syevr(amatrix, work, eigenvalues)
END IF
! Restore vectors C = U^(-1) * C*
CALL cp_fm_triangular_multiply(bmatrix, work)
CALL cp_fm_to_fm(work, eigenvectors, nao)
END SUBROUTINE cp_fm_geeig_syevr
! **************************************************************************************************
!> \brief General Eigenvalue Problem AX = BXE
!> Use canonical diagonalization : U*s**(-1/2)
Expand Down
4 changes: 1 addition & 3 deletions src/input_constants.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -632,9 +632,7 @@ MODULE input_constants
INTEGER, PARAMETER, PUBLIC :: xas_tdp_by_index = 1, &
xas_tdp_by_kind = 2, &
xas_tdp_roks = 1, &
xas_tdp_uks = 2, &
xas_tdp_diag_std = 1, &
xas_tdp_diag_syevr = 2
xas_tdp_uks = 2

! Form of dipole operator for TDDFPT oscillator strength calculation
INTEGER, PARAMETER, PUBLIC :: tddfpt_dipole_berry = 1, &
Expand Down
42 changes: 6 additions & 36 deletions src/input_cp2k_dft.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,8 +95,8 @@ MODULE input_cp2k_dft
wfi_use_prev_p_method_nr, wfi_use_prev_rho_r_method_nr, wfi_use_prev_wf_method_nr, &
xas_1s_type, xas_2p_type, xas_2s_type, xas_3d_type, xas_3p_type, xas_3s_type, xas_4d_type, &
xas_4f_type, xas_4p_type, xas_4s_type, xas_dip_len, xas_dip_vel, xas_dscf, xas_none, &
xas_not_excited, xas_tdp_by_index, xas_tdp_by_kind, xas_tdp_diag_std, xas_tdp_diag_syevr, &
xas_tp_fh, xas_tp_flex, xas_tp_hh, xas_tp_xfh, xas_tp_xhh, xes_tp_val
xas_not_excited, xas_tdp_by_index, xas_tdp_by_kind, xas_tp_fh, xas_tp_flex, xas_tp_hh, &
xas_tp_xfh, xas_tp_xhh, xes_tp_val
USE input_cp2k_almo, ONLY: create_almo_scf_section
USE input_cp2k_distribution, ONLY: create_distribution_section
USE input_cp2k_field, ONLY: create_efield_section,&
Expand Down Expand Up @@ -7891,7 +7891,7 @@ SUBROUTINE create_xas_tdp_section(section)
description="XAS simulations using linear-response TDDFT. Excitation from "// &
"specified core orbitals are considered one at a time. In case of high "// &
"symmetry structures, donor core orbitals should be localized.", &
n_keywords=13, n_subsections=4, repeats=.FALSE.)
n_keywords=15, n_subsections=3, repeats=.FALSE.)

NULLIFY (keyword, subsection, print_key)

Expand Down Expand Up @@ -8110,14 +8110,6 @@ SUBROUTINE create_xas_tdp_section(section)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

! The diagonalization subsection
CALL section_create(subsection, __LOCATION__, name="DIAGONALIZATION", &
description="Set the choice of method and parameters for solving the "// &
"XAS TDP eigenvalue problem.", &
n_keywords=3, &
n_subsections=0, &
repeats=.FALSE.)

CALL keyword_create(keyword, __LOCATION__, name="N_EXCITED", &
variants=(/"N_ROOTS"/), &
description="The number of excited states to compute per donor "// &
Expand All @@ -8126,43 +8118,21 @@ SUBROUTINE create_xas_tdp_section(section)
"If N_EXCITED is set to -1, all excitations are considered", &
usage="N_EXCITED {integer}", &
default_i_val=-1)
CALL section_add_keyword(subsection, keyword)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="ENERGY_RANGE", &
variants=s2a("E_RANGE"), &
description="The energy range in eV for which excitations are considered. "// &
"Only excitated states within the range of: first excitation "// &
"energy + ENERGY_RANGE are kept/computed. If ENERGY_RANGE "// &
"energy + ENERGY_RANGE are kept. If ENERGY_RANGE "// &
"and N_EXCITED are specified, the former has priority. "// &
"Negative values are ignored and N_EXCITED takes over.", &
usage="ENERGY_RANGE {real}", &
default_r_val=-1.0_dp)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)

CALL keyword_create(keyword, __LOCATION__, name="DIAGONALIZATION_TYPE", &
variants=s2a("DIAG_TYPE", "DIAG"), &
description="Which digonalization procedure to use", &
usage="DIAGONALIZATION_TYPE {string}", &
enum_c_vals=s2a("STD", "SYEVR"), &
enum_desc=s2a("The standard used everywhere else in CP2K (usually ELPA), "// &
"all eigenvectors are computed but only N_EXCITED are kept.", &
"Using the PDSYEVR routine of ScalaPACK which allows to "// &
"only solve for N_EXCITED eigenvectors. Can be considerably "// &
"faster than the standard for large problems, but might "// &
"crash for smaller problems due to a ScalaPACK bug (PDORMTR)."// &
" Requires version 2 of ScalaPACK. Final note: "// &
"ELPA is known to scale better with the number of MPI ranks."), &
enum_i_vals=(/xas_tdp_diag_std, xas_tdp_diag_syevr/), &
default_i_val=xas_tdp_diag_std)
CALL section_add_keyword(subsection, keyword)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
! End of the diagonalization subsection

! The KERNEL subsection
CALL section_create(subsection, __LOCATION__, name="KERNEL", &
description="Defines how the kernel is built in terms of functionals.", &
Expand Down
11 changes: 1 addition & 10 deletions src/xas_tdp_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ MODULE xas_tdp_methods
do_loc_none, do_potential_coulomb, do_potential_id, do_potential_short, &
do_potential_truncated, op_loc_berry, state_loc_list, tddfpt_singlet, tddfpt_spin_cons, &
tddfpt_spin_flip, tddfpt_triplet, xas_1s_type, xas_2p_type, xas_2s_type, xas_dip_len, &
xas_dip_vel, xas_not_excited, xas_tdp_by_index, xas_tdp_by_kind, xas_tdp_diag_syevr
xas_dip_vel, xas_not_excited, xas_tdp_by_index, xas_tdp_by_kind
USE input_cp2k_dft, ONLY: create_localize_section
USE input_section_types, ONLY: &
section_release, section_type, section_vals_create, section_vals_duplicate, &
Expand Down Expand Up @@ -408,15 +408,6 @@ SUBROUTINE xas_tdp_core(xas_tdp_section, qs_env)
! Do main calculations here
CALL setup_xas_tdp_prob(current_state, qs_env, xas_tdp_env, xas_tdp_control)
! Warn against possible failure of SYEVR if this is the solver of choice
IF (output_unit > 0 .AND. xas_tdp_control%diag_type == xas_tdp_diag_syevr) THEN
WRITE (UNIT=output_unit, FMT="(/T5,A,/,T5,A,/,T5,A)") &
"Diagonalizing using ScaLAPACK's PDSYEVR routine. This routine is not always", &
"stable and might lead to a crash. If it so happens, switch back to the", &
"default DIAG_TYPE STD keyword."
CALL m_flush(output_unit)
END IF
IF (xas_tdp_control%do_spin_cons) THEN
CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
ex_type=tddfpt_spin_cons)
Expand Down

0 comments on commit abc627f

Please sign in to comment.