eigensolver.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module eigensolver_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use eigen_cg_oct_m
  use eigen_chebyshev_oct_m
  use eigen_evolution_oct_m
  use eigen_rmmdiis_oct_m
  use electron_space_oct_m
  use exponential_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_elec_base_oct_m
  use interaction_partner_oct_m
  use, intrinsic :: iso_fortran_env
  use ks_potential_oct_m
  use lalg_adv_oct_m
  use lalg_basic_oct_m
  use loct_oct_m
  use mesh_oct_m
  use mesh_batch_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use mpi_lib_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use potential_interpolation_oct_m
  use preconditioners_oct_m
  use profiling_oct_m
  use smear_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_calc_oct_m
  use states_elec_dim_oct_m
  use states_elec_parallel_oct_m
  use subspace_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use wfs_elec_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::            &
    eigensolver_t,     &
    eigensolver_init,  &
    eigensolver_end
 
  type eigensolver_t
    private
    integer, public :: es_type
 
    real(real64),   public :: tolerance
    integer, public :: es_maxiter
 
    real(real64) :: tau
    real(real64) :: tau0
    real(real64) :: time
    integer :: it_propagator
    type(potential_interpolation_t) :: vks_old
    real(real64), allocatable :: normalization_energies(:, :)
    !                                                               !! evolution eigensolver
    real(real64), allocatable :: normalization_energies_prev(:, :)
    logical :: variable_timestep
 
    real(real64), allocatable,   public :: diff(:, :)
    integer,              public :: matvec
    integer, allocatable, public :: converged(:)
 
    type(preconditioner_t), public :: pre
 
    type(subspace_t) :: sdiag
 
    integer :: rmmdiis_minimization_iter
 
    logical, public :: folded_spectrum
 
    ! cg options
    logical, public :: orthogonalize_to_all
    integer, public :: conjugate_direction
    logical, public :: additional_terms
    real(real64),   public :: energy_change_threshold
 
    ! Chebyshev filtering options
    type(eigen_chebyshev_t), public :: cheby_params
 
    type(exponential_t) :: exponential_operator
  contains
    procedure :: run =>  eigensolver_run
  end type eigensolver_t
 
 
  integer, public, parameter :: &
    RS_PLAN      = 11,          &
    rs_cg        =  5,          &
    rs_evo       =  9,          &
    rs_rmmdiis   = 10,          &
    rs_chebyshev = 12
 
contains
 
  ! ---------------------------------------------------------
  subroutine eigensolver_init(eigens, namespace, gr, st, hm, mc, space, deactivate_oracle)
    type(eigensolver_t),      intent(out)   :: eigens
    type(namespace_t),        intent(in)    :: namespace
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(in)    :: st
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(multicomm_t),        intent(in)    :: mc
    class(space_t),           intent(in)    :: space
    logical, optional,        intent(in)    :: deactivate_oracle
 
    integer :: default_iter, default_es
    real(real64)   :: default_tol
    real(real64)   :: mem
 
    push_sub(eigensolver_init)
 
    !%Variable Eigensolver
    !%Type integer
    !%Section SCF::Eigensolver
    !%Description
    !% Which eigensolver to use to obtain the lowest eigenvalues and
    !% eigenfunctions of the Kohn-Sham Hamiltonian. The default is
    !% conjugate gradients (<tt>cg</tt>), except that when parallelization in states is
    !% enabled, the default is <tt>chebyshev_filter</tt>.
    !%Option cg 5
    !% Conjugate-gradients algorithm.
    !%Option plan 11
    !% Preconditioned Lanczos scheme. Ref: Y. Saad, A. Stathopoulos, J. Chelikowsky, K. Wu and S. Ogut,
    !% "Solution of Large Eigenvalue Problems in Electronic Structure Calculations", <i>BIT</i> <b>36</b>, 1 (1996).
    !%Option evolution 9
    !% (Experimental) Propagation in imaginary time.
    !% This eigensolver uses a time propagation in imaginary time (see, e.g., Aichinger, M. & Krotscheck, E.,
    !% Computational Materials Science 34, 188–212 (2005), 10.1016/j.commatsci.2004.11.002).
    !% The timestep is set using <tt>EigensolverImaginaryTime</tt>. The convergence typically depends
    !% on the total imaginary time in a propagation which means that it should be faster with larger
    !% timesteps. However, simulations typically diverge above a certain timestep that depends on the system.
    !% In general, it is advisable to use an exponential method that always converges, especially for large
    !% timesteps, i.e., set <tt>TDExponentialMethod</tt> to <tt>lanczos</tt> or <tt>chebyshev</tt>.
    !% The propagation method can be set with <tt>ImaginaryTimePropagator</tt> and variable timestepping
    !% can be enabled with <tt>ImaginaryTimeVariableTimestep</tt>.
    !% Moreover, this method usually converges quite slowly, so you might need to increase <tt>MaximumIter</tt>.
    !% It is incompatible with smearing for the occupations.
    !%Option rmmdiis 10
    !% Residual minimization scheme, direct inversion in the
    !% iterative subspace eigensolver, based on the implementation of
    !% Kresse and Furthm&uuml;ller [<i>Phys. Rev. B</i> <b>54</b>, 11169
    !% (1996)]. This eigensolver requires almost no orthogonalization
    !% so it can be considerably faster than the other options for
    !% large systems. To improve its performance a large number of <tt>ExtraStates</tt>
    !% are required (around 10-20% of the number of occupied states).
    !% Note: with <tt>unocc</tt>, you will need to stop the calculation
    !% by hand, since the highest states will probably never converge.
    !% Usage with more than one block of states per node is experimental, unfortunately.
    !%Option chebyshev_filter 12
    !% Chebyshev polynomials are used to construct a spectral filter that is iteratively applied to a trial subspace,
    !% amplifying a subset of the lowest eigenvalues of the Hamiltonian, which effectively isolates the invariant
    !% subspace associated with the these states. The filtered subspace is projected to define a reduced eigenvalue problem,
    !% which is then solved using dense diagonalisation. Finally, the eigenstates of the original Hamiltonian are recovered
    !% via subspace rotation. For further details, see [Zhou et. al.](http://dx.doi.org/10.1016/j.jcp.2014.06.056)
    !%End
 
    if(st%parallel_in_states) then
      default_es = rs_chebyshev
    else
      default_es = rs_cg
    end if
 
    call parse_variable(namespace, 'Eigensolver', default_es, eigens%es_type)
 
    if(st%parallel_in_states .and. .not. eigensolver_parallel_in_states(eigens)) then
      message(1) = "The selected eigensolver is not parallel in states."
      message(2) = "Please use the rmmdiis or Chebyshev filtering eigensolvers."
      call messages_fatal(2, namespace=namespace)
    end if
 
    call messages_obsolete_variable(namespace, 'EigensolverVerbose')
    call messages_obsolete_variable(namespace, 'EigensolverSubspaceDiag', 'SubspaceDiagonalization')
 
    default_iter = 25
 
    select case(eigens%es_type)
    case(rs_cg)
      !%Variable CGOrthogonalizeAll
      !%Type logical
      !%Default no
      !%Section SCF::Eigensolver
      !%Description
      !% Used by the cg solver only.
      !% During the cg iterations, the current band can be orthogonalized
      !% against all other bands or only against the lower bands. Orthogonalizing
      !% against all other bands can improve convergence properties, whereas
      !% orthogonalizing against lower bands needs less operations.
      !% Moreover, orthogonalizing against all bands can make converging
      !% the highest band or unoccupied bands more difficult.
      !%End
      call parse_variable(namespace, 'CGOrthogonalizeAll', .false., eigens%orthogonalize_to_all)
 
      !%Variable CGDirection
      !%Type integer
      !%Section SCF::Eigensolver
      !%Description
      !% Used by the cg solver only.
      !% The conjugate direction is updated using a certain coefficient to the previous
      !% direction. This coeffiction can be computed in different ways. The default is
      !% to use Fletcher-Reeves (FR), an alternative is Polak-Ribiere (PR).
      !%Option fletcher 1
      !% The coefficient for Fletcher-Reeves consists of the current norm of the
      !% steepest descent vector divided by that of the previous iteration.
      !%Option polak 2
      !% For the Polak-Ribiere scheme, a product of the current with the previous
      !% steepest descent vector is subtracted in the nominator.
      !%End
      call parse_variable(namespace, 'CGDirection', option__cgdirection__fletcher, eigens%conjugate_direction)
 
      !%Variable CGAdditionalTerms
      !%Type logical
      !%Section SCF::Eigensolver
      !%Default no
      !%Description
      !% Used by the cg solver only.
      !% Add additional terms during the line minimization, see PTA92, eq. 5.31ff.
      !% These terms can improve convergence for some systems, but they are quite costly.
      !% If you experience convergence problems, you might try out this option.
      !% This feature is still experimental.
      !%End
      call parse_variable(namespace, 'CGAdditionalTerms', .false., eigens%additional_terms)
      if(eigens%additional_terms) then
        call messages_experimental("The additional terms for the CG eigensolver are not tested for all cases.")
      end if
 
      !%Variable CGEnergyChangeThreshold
      !%Type float
      !%Section SCF::Eigensolver
      !%Default 0.1
      !%Description
      !% Used by the cg solver only.
      !% For each band, the CG iterations are stopped when the change in energy is smaller than the
      !% change in the first iteration multiplied by this factor. This limits the number of CG
      !% iterations for each band, while still showing good convergence for the SCF cycle. The criterion
      !% is discussed in Sec. V.B.6 of Payne et al. (1992), Rev. Mod. Phys. 64, 4.
      !% The default value is 0.1, which is usually a good choice for LDA and GGA potentials. If you
      !% are solving the OEP equation, you might want to set this value to 1e-3 or smaller. In general,
      !% smaller values might help if you experience convergence problems.
      !% For very small convergence tolerances, choose 0 to disable this criterion.
      !%End
      call parse_variable(namespace, 'CGEnergyChangeThreshold', 0.1_real64, eigens%energy_change_threshold)
 
    case(rs_plan)
    case(rs_evo)
      call messages_experimental("imaginary-time evolution eigensolver")
 
      !%Variable EigensolverImaginaryTime
      !%Type float
      !%Default 0.1
      !%Section SCF::Eigensolver
      !%Description
      !% The imaginary-time step that is used in the imaginary-time evolution
      !% method (<tt>Eigensolver = evolution</tt>) to obtain the lowest eigenvalues/eigenvectors.
      !% It must satisfy <tt>EigensolverImaginaryTime > 0</tt>.
      !% Increasing this value can make the propagation faster, but could lead to unstable propagations.
      !%End
      call parse_variable(namespace, 'EigensolverImaginaryTime', 0.1_real64, eigens%tau)
      if(eigens%tau <= m_zero) call messages_input_error(namespace, 'EigensolverImaginaryTime')
      ! save initial timestep
      eigens%tau0 = eigens%tau
 
      call exponential_init(eigens%exponential_operator, namespace)
 
      if(st%smear%method /= smear_semiconductor .and. st%smear%method /= smear_fixed_occ) then
        message(1) = "Smearing of occupations is incompatible with imaginary time evolution."
        call messages_fatal(1, namespace=namespace)
      end if
 
      !%Variable ImaginaryTimePropagator
      !%Type integer
      !%Section SCF::Eigensolver
      !%Default it_forward_euler
      !%Description
      !% The propagation method for imaginary-time evolution.
      !%Option it_forward_euler 1
      !% Approximate the evolution operator by one forward Euler step with the exponential:
      !% <math>
      !%   U_{\rm FE}(t+\delta t, t) = \exp \left( -\tau H_{t}\right)\,.
      !% </math>
      !% This scheme is first order.
      !%Option it_expmid 2
      !% Similar to forward Euler, but evaluate the Hamiltonian at the half-timestep by
      !% extrapolation:
      !% <math>
      !%   U_{\rm EM}(t+\delta t, t) = \exp \left( -\tau H_{t+\delta t/2}\right)\,.
      !% </math>
      !% This scheme is second order.
      !%End
      call parse_variable(namespace, 'ImaginaryTimePropagator', it_forward_euler, eigens%it_propagator)
      ! initialize potential interpolation
      if (eigens%it_propagator == it_expmid) then
        call hm%ks_pot%init_interpolation(eigens%vks_old)
        call hm%ks_pot%run_zero_iter(eigens%vks_old)
      end if
 
      !%Variable ImaginaryTimeVariableTimestep
      !%Type logical
      !%Default no
      !%Section SCF::Eigensolver
      !%Description
      !% Enable variable timesteps for the imaginary time propagation.
      !%End
      call parse_variable(namespace, 'ImaginaryTimeVariableTimestep', .false., eigens%variable_timestep)
 
      eigens%time = m_zero
 
      safe_allocate(eigens%normalization_energies(1:st%nst, 1:st%nik))
      safe_allocate(eigens%normalization_energies_prev(1:st%nst, 1:st%nik))
 
    case(rs_rmmdiis)
      default_iter = 5
 
      !%Variable EigensolverMinimizationIter
      !%Type integer
      !%Default 0
      !%Section SCF::Eigensolver
      !%Description
      !% During the first iterations, the RMMDIIS eigensolver requires
      !% some steepest-descent minimizations to improve
      !% convergence. This variable determines the number of those
      !% minimizations.
      !%End
 
      call parse_variable(namespace, 'EigensolverMinimizationIter', 0, eigens%rmmdiis_minimization_iter)
 
      if(gr%use_curvilinear) call messages_experimental("RMMDIIS eigensolver for curvilinear coordinates")
 
    case(rs_chebyshev)
      !%Variable ChebyshevFilterLanczosOrder
      !%Type integer
      !%Default 5
      !%Section SCF::Eigensolver
      !%Description
      !% Used by the Chebyshev filter only.
      !% The number of Lanczos iterations used to construct the tridiagonal matrix,
      !% from which the upper bound of H is estimated.
      !% A value in the range 4 <= ChebyshevFilterLanczosOrder <= 10 is reasonable.
      !% Values greater than 10 will raise an assertion.
      !%End
      call parse_variable(namespace, 'ChebyshevFilterLanczosOrder', default_chebyshev_params%n_lanczos, &
        eigens%cheby_params%n_lanczos)
 
      ! TODO Rewrite this, having changed the behaviour
      !%Variable ChebyshevFilterDegree
      !%Type integer
      !%Default 25
      !%Section SCF::Eigensolver
      !%Description
      !% Used by the Chebyshev filter only.
      !% The degree of the Chebyshev polynomial used to dampen the interval of eigenstates one is not interested in.
      !% If used in conjunction with "OptimizeChebyshevFilterDegree", which is the default, "ChebyshevFilterDegree" defines the
      !% the maximum Chebyshev polynomial degree Octopus will consider when determining an approximate, optimal degree.
      !%
      !% If not used in conjunction with "OptimizeChebyshevFilterDegree", this defines the polynomial degree used for all SCF steps.
      !% The larger the degree, the larger the separation between the subspaces, which will reduce the number of SCF iterations
      !% required to reach convergence. However, ChebyshevFilterDegree also directly correlates with the number of matrix-vector
      !% products performed on the Hamiltonian per SCF step, and so larger values result in slower SCF cycles.
      !% A value in the range 8 <= ChebyshevFilterDegree <= 20 is reasonable.
      !%End
      call parse_variable(namespace, 'ChebyshevFilterDegree', default_chebyshev_params%degree, eigens%cheby_params%degree)
 
      !%Variable OptimizeChebyshevFilterDegree
      !%Type logical
      !%Default yes
      !%Section SCF::Eigensolver
      !%Description
      !% Used by the Chebyshev filter only.
      !% Octopus generates a best-estimate for the degree of the Chebyshev polynomial used to filter the subspace.
      !% A separate estimate is generated for each state block, per SCF iteration. Note that if running parallelism
      !% over states, the block/batch containing the largest eigenstates will converge more slowly and will
      !% therefore use a larger degree relative to all other batches. One should therefore avoid setting "ChebyshevFilterDegree"
      !% to an excessive value (for example > 50). For more details regarding how the degree is estimated, one can
      !% refer to Section 5.4 in "Computer Physics Communications 187 (2015) 98–105"
      !% (http://dx.doi.org/10.1016/j.cpc.2014.10.015).
      !%
      !% This is deactivated by default for spinors as this is not found to be advantageous in test systems.
      !%End
      if (st%d%ispin == spinors .or. optional_default(deactivate_oracle, .false.)) then
        default_chebyshev_params%optimize_degree = .false.
      end if
      call parse_variable(namespace, 'OptimizeChebyshevFilterDegree', default_chebyshev_params%optimize_degree, &
        eigens%cheby_params%optimize_degree)
 
      !%Variable ChebyshevFilterBoundMixing
      !%Type float
      !%Default 0.5
      !%Section SCF::Eigensolver
      !%Description
      !% In the first application of the filter, for the first SCF step, the initial lower bound estimate
      !% is defined as a linear combination of the smallest and largest eigenvalues of the Hamiltonian.
      !% The bound mixing determines the proportion of linear mixing, beta:
      !% $b_{lower} = beta min{\lambda} + (\beta - 1) max{\lambda}$
      !% of the smallest and largest eigenvalues.
      !%End
      call parse_variable(namespace, 'ChebyshevFilterBoundMixing', default_chebyshev_params%bound_mixing, &
        eigens%cheby_params%bound_mixing)
 
      !%Variable ChebyshevFilterNIter
      !%Type integer
      !%Default 5
      !%Section SCF::Eigensolver
      !%Description
      !% The max number of iterations in the first SCF step of the Chebyshev solver. In practice this
      !% does not need to exceed 10, as the eigenstates will vary alot between the first and second
      !% SCF steps.
      !%End
      call parse_variable(namespace, 'ChebyshevFilterNIter', default_chebyshev_params%n_iter, &
        eigens%cheby_params%n_iter)
 
    case default
      call messages_input_error(namespace, 'Eigensolver')
    end select
 
    call messages_print_with_emphasis(msg='Eigensolver', namespace=namespace)
 
    call messages_print_var_option("Eigensolver", eigens%es_type, namespace=namespace)
 
    call messages_obsolete_variable(namespace, 'EigensolverInitTolerance', 'EigensolverTolerance')
    call messages_obsolete_variable(namespace, 'EigensolverFinalTolerance', 'EigensolverTolerance')
    call messages_obsolete_variable(namespace, 'EigensolverFinalToleranceIteration')
 
    ! this is an internal option that makes the solver use the
    ! folded operator (H-shift)^2 to converge first eigenvalues around
    ! the values of shift
    ! c.f. L. W. Wang and A. Zunger
    ! JCP 100, 2394 (1994); doi: http://dx.doi.org/10.1063/1.466486
    eigens%folded_spectrum = .false.
 
    !%Variable EigensolverTolerance
    !%Type float
    !%Section SCF::Eigensolver
    !%Description
    !% This is the tolerance for the eigenvectors. The default is 1e-7.
    !% However, if <tt>ConvRelDens</tt> is defined, the default is set to a tenth of its value.
    !%End
    default_tol = 1e-7_real64
    if (parse_is_defined(namespace, 'ConvRelDens')) then
      call parse_variable(namespace, 'ConvRelDens', 1e-6_real64, default_tol)
      default_tol = default_tol / 10.0_real64
    end if
    call parse_variable(namespace, 'EigensolverTolerance', default_tol, eigens%tolerance)
 
    !%Variable EigensolverMaxIter
    !%Type integer
    !%Section SCF::Eigensolver
    !%Description
    !% Determines the maximum number of iterations that the
    !% eigensolver will perform if the desired tolerance is not
    !% achieved. The default is 25 iterations for all eigensolvers
    !% except for <tt>rmdiis</tt>, which performs only 5 iterations.
    !% Increasing this value for <tt>rmdiis</tt> increases the convergence speed,
    !% at the cost of an increased memory footprint.
    !%End
    call parse_variable(namespace, 'EigensolverMaxIter', default_iter, eigens%es_maxiter)
    if(eigens%es_maxiter < 1) call messages_input_error(namespace, 'EigensolverMaxIter')
 
    if(eigens%es_maxiter > default_iter) then
      call messages_write('You have specified a large number of eigensolver iterations (')
      call messages_write(eigens%es_maxiter)
      call messages_write(').', new_line = .true.)
      call messages_write('This is not a good idea as it might slow down convergence, even for', new_line = .true.)
      call messages_write('independent particles, as subspace diagonalization will not be used', new_line = .true.)
      call messages_write('often enough.')
      call messages_warning(namespace=namespace)
    end if
 
    if (any(eigens%es_type == (/rs_plan, rs_cg, rs_rmmdiis/))) then
      call preconditioner_init(eigens%pre, namespace, gr, mc, space)
    end if
 
    safe_allocate(eigens%diff(1:st%nst, 1:st%nik))
    eigens%diff(1:st%nst, 1:st%nik) = 0
 
    safe_allocate(eigens%converged(1:st%nik))
    eigens%converged(1:st%nik) = 0
    eigens%matvec = 0
 
    call eigens%sdiag%init(namespace, st)
 
    ! print memory requirements
    select case(eigens%es_type)
    case(rs_rmmdiis)
      call messages_write('Info: The rmmdiis eigensolver requires ')
      mem = (m_two*eigens%es_maxiter - m_one)*st%block_size*real(gr%np_part, real64)
      if(states_are_real(st)) then
        mem = mem*8.0_real64
      else
        mem = mem*16.0_real64
      end if
      call messages_write(mem, units = unit_megabytes, fmt = '(f9.1)')
      call messages_write(' of additional')
      call messages_new_line()
      call messages_write('      memory.  This amount can be reduced by decreasing the value')
      call messages_new_line()
      call messages_write('      of the variable StatesBlockSize (currently set to ')
      call messages_write(st%block_size)
      call messages_write(').')
      call messages_info()
    end select
 
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(eigensolver_init)
  end subroutine eigensolver_init
 
 
  ! ---------------------------------------------------------
  subroutine eigensolver_end(eigens)
    type(eigensolver_t), intent(inout) :: eigens
 
    push_sub(eigensolver_end)
 
    select case(eigens%es_type)
    case(rs_plan, rs_cg, rs_rmmdiis)
      call preconditioner_end(eigens%pre)
    end select
 
    safe_deallocate_a(eigens%converged)
    safe_deallocate_a(eigens%diff)
 
    safe_deallocate_a(eigens%normalization_energies)
    safe_deallocate_a(eigens%normalization_energies_prev)
 
    pop_sub(eigensolver_end)
  end subroutine eigensolver_end
 
 
  ! ---------------------------------------------------------
  subroutine eigensolver_run(eigens, namespace, gr, st, hm, space, ext_partners, iter, conv, nstconv)
    class(eigensolver_t),     intent(inout) :: eigens
    type(namespace_t),        intent(in)    :: namespace
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    class(space_t),           intent(in)    :: space
    type(partner_list_t),     intent(in)    :: ext_partners
    integer,                  intent(in)    :: iter
    logical,        optional, intent(out)   :: conv
    integer,        optional, intent(in)    :: nstconv
    !                                                  !< the convergence criteria
 
    integer :: nstconv_
#ifdef HAVE_MPI
    logical :: conv_reduced
    integer :: ist, outcount, lmatvec
    real(real64), allocatable :: ldiff(:), leigenval(:)
    real(real64), allocatable :: ldiff_out(:), leigenval_out(:)
    integer, allocatable :: lconv(:)
#endif
 
    call profiling_in("EIGEN_SOLVER")
    push_sub(eigensolver_run)
 
    if(present(conv)) conv = .false.
    if(present(nstconv)) then
      nstconv_ = nstconv
    else
      nstconv_ = st%nst
    end if
 
    eigens%matvec = 0
 
    if(mpi_grp_is_root(mpi_world) .and. eigensolver_has_progress_bar(eigens) .and. .not. debug%info) then
      call loct_progress_bar(-1, st%lnst*st%d%kpt%nlocal)
    end if
 
    if (states_are_real(st)) then
      call deigensolver_run(eigens, namespace, gr, st, hm, space, ext_partners, iter)
    else
      call zeigensolver_run(eigens, namespace, gr, st, hm, space, ext_partners, iter)
    end if
 
    if(mpi_grp_is_root(mpi_world) .and. eigensolver_has_progress_bar(eigens) .and. .not. debug%info) then
      write(stdout, '(1x)')
    end if
 
    if(present(conv)) conv = all(eigens%converged(st%d%kpt%start:st%d%kpt%end) >= nstconv_)
 
#ifdef HAVE_MPI
    if (st%d%kpt%parallel) then
      if (present(conv)) then
        call st%d%kpt%mpi_grp%allreduce(conv, conv_reduced, 1, mpi_logical, mpi_land)
        conv = conv_reduced
      end if
 
      lmatvec = eigens%matvec
      call st%d%kpt%mpi_grp%allreduce(lmatvec, eigens%matvec, 1, mpi_integer, mpi_sum)
 
      safe_allocate(lconv(1:st%d%kpt%nlocal))
      lconv(1:st%d%kpt%nlocal) = eigens%converged(st%d%kpt%start:st%d%kpt%end)
      call lmpi_gen_allgatherv(st%d%kpt%nlocal, lconv, outcount, eigens%converged, st%d%kpt%mpi_grp)
      assert(outcount == st%nik)
      safe_deallocate_a(lconv)
 
      ! every node needs to know all eigenvalues (and diff)
      safe_allocate(ldiff(1:st%d%kpt%nlocal))
      safe_allocate(leigenval(1:st%d%kpt%nlocal))
      safe_allocate(ldiff_out(1:st%nik))
      safe_allocate(leigenval_out(1:st%nik))
      do ist = st%st_start, st%st_end
        ldiff(1:st%d%kpt%nlocal) = eigens%diff(ist, st%d%kpt%start:st%d%kpt%end)
        leigenval(1:st%d%kpt%nlocal) = st%eigenval(ist, st%d%kpt%start:st%d%kpt%end)
        call lmpi_gen_allgatherv(st%d%kpt%nlocal, ldiff, outcount, ldiff_out, st%d%kpt%mpi_grp)
        eigens%diff(ist, :) = ldiff_out
        assert(outcount == st%nik)
        call lmpi_gen_allgatherv(st%d%kpt%nlocal, leigenval, outcount, leigenval_out, st%d%kpt%mpi_grp)
        st%eigenval(ist, :) = leigenval_out
        assert(outcount == st%nik)
      end do
      safe_deallocate_a(ldiff)
      safe_deallocate_a(ldiff_out)
      safe_deallocate_a(leigenval)
      safe_deallocate_a(leigenval_out)
    end if
#endif
 
    pop_sub(eigensolver_run)
    call profiling_out("EIGEN_SOLVER")
  end subroutine eigensolver_run
 
 
  ! ---------------------------------------------------------
  logical function eigensolver_parallel_in_states(this) result(par_stat)
    type(eigensolver_t), intent(in) :: this
 
    push_sub(eigensolver_parallel_in_states)
 
    par_stat = .false.
 
    select case(this%es_type)
    case(rs_rmmdiis, rs_evo, rs_chebyshev)
      par_stat = .true.
    end select
 
    pop_sub(eigensolver_parallel_in_states)
  end function eigensolver_parallel_in_states
 
 
  ! ---------------------------------------------------------
  logical function eigensolver_has_progress_bar(this) result(has)
    type(eigensolver_t), intent(in) :: this
 
    push_sub(eigensolver_has_progress_bar)
 
    has = .false.
 
    select case(this%es_type)
    case(rs_rmmdiis, rs_cg)
      has = .true.
    end select
 
    pop_sub(eigensolver_has_progress_bar)
  end function eigensolver_has_progress_bar
 
 
#include "undef.F90"
#include "real.F90"
#include "eigensolver_inc.F90"
#include "eigen_plan_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "eigensolver_inc.F90"
#include "eigen_plan_inc.F90"
 
end module eigensolver_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: