Skip to content
Permalink
Browse files

XAS_TDP| Added ENERGY_RANGE keyword in diagonalization to only compute

energies in a given window. Added SYEVR warning due to instability.
  • Loading branch information...
abussy authored and pseewald committed Sep 26, 2019
1 parent e01c96e commit cf644c8ede9cdba4409bc45f2a3954d567114ab5
@@ -87,6 +87,7 @@ MODULE cp_fm_diag
cp_fm_syevx, &
cp_fm_geeig, &
cp_fm_geeig_canon, &
cp_fm_geeig_syevr, &
diag_init, &
diag_finalize

@@ -720,6 +721,201 @@ 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 ...
@@ -1160,6 +1356,69 @@ 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)
@@ -8089,7 +8089,7 @@ SUBROUTINE create_xas_tdp_section(section)
CALL section_create(subsection, __LOCATION__, name="DIAGONALIZATION", &
description="Set the choice of method and parameters for solving the "// &
"XAS TDP eigenvalue problem.", &
n_keywords=2, &
n_keywords=3, &
n_subsections=0, &
repeats=.FALSE.)

@@ -8104,6 +8104,18 @@ SUBROUTINE create_xas_tdp_section(section)
CALL section_add_keyword(subsection, 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 "// &
"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", &
@@ -359,11 +359,6 @@ SUBROUTINE init_xas_atom_grid_harmo(xas_atom_env, grid_info, do_xc, qs_env)
grid_atom%ng_sphere = na
grid_atom%nr = nr

IF (llmax > la) THEN
CPWARN("A Lebedev grid for XAS TDP is built for a too low angular momentum l")
CPWARN("To fix this (l too low), increase the value of 'na' in the GRID keyword.")
END IF

! If this is an excited kind, create the harmonics with the RI_XAS basis, otherwise the ORB
IF (ANY(xas_atom_env%excited_kinds == ikind) .AND. do_xc) THEN
CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="RI_XAS")

0 comments on commit cf644c8

Please sign in to comment.
You can’t perform that action at this time.