quadratic_pseudospectrum.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Metric-aware quadratic pseudospectra for finite AO operators.
!> \author Thomas D. Kuehne, tkuehne@cp2k.org
! **************************************************************************************************
MODULE quadratic_pseudospectrum
   USE ieee_arithmetic, ONLY: ieee_is_finite
#if defined(__HAS_IEEE_EXCEPTIONS)
   USE ieee_exceptions, ONLY: ieee_all, ieee_get_halting_mode, ieee_set_halting_mode
#endif
   USE kinds, ONLY: dp
   IMPLICIT NONE
   PRIVATE
   PUBLIC :: quadratic_dense, quadratic_solve, quadratic_result_type
 
   TYPE quadratic_result_type
      REAL(dp), ALLOCATABLE :: eigenvalues(:), residuals(:), errors(:, :), means(:, :)
      REAL(dp), ALLOCATABLE :: projected_errors(:, :), leakage(:, :)
      COMPLEX(dp), ALLOCATABLE :: states(:, :), complex_means(:, :)
      INTEGER :: iterations = 0
   END TYPE quadratic_result_type
 
   INTERFACE quadratic_dense
      MODULE PROCEDURE quadratic_dense_real, quadratic_dense_complex
   END INTERFACE
   INTERFACE quadratic_solve
      MODULE PROCEDURE quadratic_solve_real, quadratic_solve_complex
   END INTERFACE
 
   abstract INTERFACE
! **************************************************************************************************
!> \brief Apply a covariant AO observable or the overlap to a block of coefficients.
!> \param OPERATOR Zero for overlap, otherwise the requested observable index
!> \param v Input coefficients
!> \param y Covariant output vectors
!> \param status Zero on success
! **************************************************************************************************
      SUBROUTINE operator_apply(OPERATOR, v, y, status)
         IMPORT dp
         INTEGER, INTENT(IN) :: operator
         COMPLEX(dp), INTENT(IN) :: v(:, :)
         COMPLEX(dp), INTENT(OUT) :: y(:, :)
         INTEGER, INTENT(OUT) :: status
      END SUBROUTINE operator_apply
! **************************************************************************************************
!> \brief Solve the positive AO metric for a block of covariant right hand sides.
!> \param v Right hand sides
!> \param y Contravariant solutions
!> \param status Zero on success
! **************************************************************************************************
      SUBROUTINE metric_solve(v, y, status)
         IMPORT dp
         COMPLEX(dp), INTENT(IN) :: v(:, :)
         COMPLEX(dp), INTENT(OUT) :: y(:, :)
         INTEGER, INTENT(OUT) :: status
      END SUBROUTINE metric_solve
   END INTERFACE
 
   TYPE, ABSTRACT :: quadratic_operator_type
   CONTAINS
      PROCEDURE(operator_method), DEFERRED :: apply
      PROCEDURE(operator_method), DEFERRED :: adjoint
      PROCEDURE(metric_method), DEFERRED :: solve
   END TYPE quadratic_operator_type
 
   abstract INTERFACE
! **************************************************************************************************
!> \brief ...
!> \param this ...
!> \param OPERATOR ...
!> \param v ...
!> \param y ...
!> \param status ...
! **************************************************************************************************
      SUBROUTINE operator_method(this, OPERATOR, v, y, status)
         IMPORT dp, quadratic_operator_type
         CLASS(quadratic_operator_type), INTENT(INOUT) :: this
         INTEGER, INTENT(IN) :: operator
         COMPLEX(dp), INTENT(IN) :: v(:, :)
         COMPLEX(dp), INTENT(OUT) :: y(:, :)
         INTEGER, INTENT(OUT) :: status
      END SUBROUTINE operator_method
 
! **************************************************************************************************
!> \brief ...
!> \param this ...
!> \param v ...
!> \param y ...
!> \param status ...
! **************************************************************************************************
      SUBROUTINE metric_method(this, v, y, status)
         IMPORT dp, quadratic_operator_type
         CLASS(quadratic_operator_type), INTENT(INOUT) :: this
         COMPLEX(dp), INTENT(IN) :: v(:, :)
         COMPLEX(dp), INTENT(OUT) :: y(:, :)
         INTEGER, INTENT(OUT) :: status
      END SUBROUTINE metric_method
   END INTERFACE
 
   TYPE, EXTENDS(quadratic_operator_type) :: dense_operator_type
      COMPLEX(dp), POINTER :: operators(:, :, :) => null(), metric(:, :) => null()
      COMPLEX(dp), ALLOCATABLE :: factor(:, :)
   CONTAINS
      PROCEDURE :: apply => dense_apply
      PROCEDURE :: adjoint => dense_adjoint
      PROCEDURE :: solve => dense_solve
   END TYPE dense_operator_type
 
   TYPE, EXTENDS(quadratic_operator_type) :: callback_operator_type
      PROCEDURE(operator_apply), POINTER, NOPASS :: apply_callback => null()
      PROCEDURE(operator_apply), POINTER, NOPASS :: adjoint_callback => null()
      PROCEDURE(metric_solve), POINTER, NOPASS :: solve_callback => null()
   CONTAINS
      PROCEDURE :: apply => callback_apply
      PROCEDURE :: adjoint => callback_adjoint
      PROCEDURE :: solve => callback_solve
   END TYPE callback_operator_type
 
CONTAINS
 
! **************************************************************************************************
!> \brief Dense AO oracle, also providing a dense-metric reference for the iterative solver.
!> \param operators Hermitian observables in a covariant AO basis
!> \param metric Positive AO overlap
!> \param shift Query value for each observable
!> \param weights Nonnegative observable scales
!> \param nstates Number of lowest eigenpairs
!> \param iterative Use the matrix-free Ritz iteration rather than full diagonalization
!> \param max_subspace Maximum trial-space dimension
!> \param max_iter Maximum Ritz iterations
!> \param tolerance Relative eigenpair residual threshold
!> \param metric_tolerance Minimum overlap eigenvalue ratio
!> \param RESULT Eigenpairs, separate observable residuals and expectations
!> \param status 0 success, 1 unconverged, -1 input, -2 metric, -4 numerical failure
! **************************************************************************************************
   SUBROUTINE quadratic_dense_real(operators, metric, shift, weights, nstates, iterative, max_subspace, &
                                   max_iter, tolerance, metric_tolerance, RESULT, status)
      COMPLEX(dp), INTENT(IN)                            :: operators(:, :, :), metric(:, :)
      REAL(dp), INTENT(IN)                               :: shift(:), weights(:)
      INTEGER, INTENT(IN)                                :: nstates
      LOGICAL, INTENT(IN)                                :: iterative
      INTEGER, INTENT(IN)                                :: max_subspace, max_iter
      REAL(dp), INTENT(IN)                               :: tolerance, metric_tolerance
      TYPE(quadratic_result_type), INTENT(OUT)           :: result
      INTEGER, INTENT(OUT)                               :: status
 
      INTEGER                                            :: j
 
      status = -1
      IF (SIZE(operators, 1) < 1 .OR. SIZE(operators, 1) /= SIZE(operators, 2)) RETURN
      DO j = 1, SIZE(operators, 3)
         IF (.NOT. hermitian(operators(:, :, j))) RETURN
      END DO
      CALL quadratic_dense_complex(operators, metric, cmplx(shift, 0.0_dp, dp), weights, nstates, iterative, &
                                   max_subspace, max_iter, tolerance, metric_tolerance, result, status)
   END SUBROUTINE quadratic_dense_real
 
! **************************************************************************************************
!> \brief Dense oracle for general complex queries and possibly non-Hermitian AO operators.
!> \param operators Covariant AO operators, not necessarily normal or Hermitian
!> \param metric Positive AO overlap
!> \param shift Complex query value for each operator
!> \param weights Nonnegative residual scales
!> \param nstates Number of lowest eigenpairs
!> \param iterative Whether to use the bounded Ritz iteration
!> \param max_subspace Maximum trial-space dimension
!> \param max_iter Maximum Ritz iterations
!> \param tolerance Relative eigenpair residual threshold
!> \param metric_tolerance Minimum overlap eigenvalue ratio
!> \param RESULT Eigenpairs, residuals and complex expectations; means holds their real parts
!> \param status 0 success, 1 unconverged, -1 input, -2 metric, -4 numerical failure
!> \param unitary Optional mask for compressions of continuum unitary operators, whose full normal Gram is S
! **************************************************************************************************
   SUBROUTINE quadratic_dense_complex(operators, metric, shift, weights, nstates, iterative, max_subspace, &
                                      max_iter, tolerance, metric_tolerance, RESULT, status, unitary)
      COMPLEX(dp), INTENT(IN), TARGET :: operators(:, :, :), metric(:, :)
      COMPLEX(dp), INTENT(IN) :: shift(:)
      REAL(dp), INTENT(IN) :: weights(:), tolerance, metric_tolerance
      INTEGER, INTENT(IN) :: nstates, max_subspace, max_iter
      LOGICAL, INTENT(IN) :: iterative
      TYPE(quadratic_result_type), INTENT(OUT) :: result
      INTEGER, INTENT(OUT) :: status
      LOGICAL, INTENT(IN), OPTIONAL :: unitary(:)
      TYPE(dense_operator_type) :: operator_backend
      REAL(dp), ALLOCATABLE :: eigenvalues(:)
      INTEGER :: n
#if defined(__HAS_IEEE_EXCEPTIONS)
      LOGICAL :: halt(SIZE(IEEE_ALL))
#endif
 
      status = -1
      n = SIZE(metric, 1)
      IF (n < 1 .OR. SIZE(metric, 2) /= n) RETURN
      IF (SIZE(shift) < 1 .OR. SIZE(weights) /= SIZE(shift)) RETURN
      IF (any(shape(operators) /= [n, n, SIZE(shift)])) RETURN
      IF (.NOT. ieee_is_finite(metric_tolerance)) RETURN
      IF (metric_tolerance <= 0.0_dp .OR. metric_tolerance >= 1.0_dp) RETURN
      IF (.NOT. hermitian(metric)) RETURN
      IF (.NOT. all(ieee_is_finite(real(operators, dp))) .OR. &
          .NOT. all(ieee_is_finite(aimag(operators)))) RETURN
      operator_backend%operators => operators
      operator_backend%metric => metric
      ALLOCATE (operator_backend%factor(n, n), eigenvalues(n))
      operator_backend%factor(:, :) = metric
      CALL diagonalize(operator_backend%factor, eigenvalues, status)
      IF (status /= 0) RETURN
      status = -2
      IF (eigenvalues(n) <= 0.0_dp) RETURN
      IF (eigenvalues(1) <= metric_tolerance*eigenvalues(n)) RETURN
      operator_backend%factor(:, :) = metric
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_get_halting_mode(ieee_all, halt)
      CALL ieee_set_halting_mode(ieee_all, .false.)
#endif
      CALL zpotrf('L', n, operator_backend%factor, n, status)
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_set_halting_mode(ieee_all, halt)
#endif
      IF (status /= 0) THEN
         status = -2
         RETURN
      END IF
      CALL quadratic_solve_core(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
                                operator_backend, result, status, unitary)
   END SUBROUTINE quadratic_dense_complex
 
! **************************************************************************************************
!> \brief Apply a dense covariant operator or the metric.
!> \param this Dense operator backend
!> \param OPERATOR Zero for metric, otherwise an observable index
!> \param v Input block
!> \param y Output block
!> \param status Return status
! **************************************************************************************************
   SUBROUTINE dense_apply(this, OPERATOR, v, y, status)
      CLASS(dense_operator_type), INTENT(INOUT) :: this
      INTEGER, INTENT(IN) :: operator
      COMPLEX(dp), INTENT(IN) :: v(:, :)
      COMPLEX(dp), INTENT(OUT) :: y(:, :)
      INTEGER, INTENT(OUT) :: status
 
      status = 0
      IF (operator == 0) THEN
         y = matmul(this%metric, v)
      ELSE
         y = matmul(this%operators(:, :, operator), v)
      END IF
   END SUBROUTINE dense_apply
 
! **************************************************************************************************
!> \brief Apply the actual dense covariant adjoint, without assuming a normal operator.
!> \param this Dense operator backend
!> \param OPERATOR Observable index
!> \param v Input block
!> \param y Adjoint output block
!> \param status Return status
! **************************************************************************************************
   SUBROUTINE dense_adjoint(this, OPERATOR, v, y, status)
      CLASS(dense_operator_type), INTENT(INOUT) :: this
      INTEGER, INTENT(IN) :: operator
      COMPLEX(dp), INTENT(IN) :: v(:, :)
      COMPLEX(dp), INTENT(OUT) :: y(:, :)
      INTEGER, INTENT(OUT) :: status
 
      y = matmul(conjg(transpose(this%operators(:, :, operator))), v)
      status = 0
   END SUBROUTINE dense_adjoint
 
! **************************************************************************************************
!> \brief Reuse one dense Cholesky factorization for AO-metric solutions.
!> \param this Dense operator backend
!> \param v Right hand sides
!> \param y Solutions
!> \param status Return status
! **************************************************************************************************
   SUBROUTINE dense_solve(this, v, y, status)
      CLASS(dense_operator_type), INTENT(INOUT) :: this
      COMPLEX(dp), INTENT(IN) :: v(:, :)
      COMPLEX(dp), INTENT(OUT) :: y(:, :)
      INTEGER, INTENT(OUT) :: status
 
      INTEGER :: n
#if defined(__HAS_IEEE_EXCEPTIONS)
      LOGICAL :: halt(SIZE(IEEE_ALL))
#endif
 
      n = SIZE(this%factor, 1)
      y = v
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_get_halting_mode(ieee_all, halt)
      CALL ieee_set_halting_mode(ieee_all, .false.)
#endif
      CALL zpotrs('L', n, SIZE(v, 2), this%factor, n, y, n, status)
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_set_halting_mode(ieee_all, halt)
#endif
      IF (status /= 0) status = -4
   END SUBROUTINE dense_solve
 
! **************************************************************************************************
!> \brief Solve sum_j w_j^2 A_j S^-1 A_j c = lambda S c without squaring sparse operators.
!> \param n AO/spinor dimension
!> \param shift Query values for the observables
!> \param weights Nonnegative observable scales
!> \param nstates Number of lowest eigenpairs
!> \param iterative Whether to use a bounded restarted block Ritz space
!> \param max_subspace Maximum trial-space size
!> \param max_iter Maximum Ritz iterations
!> \param tolerance Relative S^-1 eigenpair residual threshold
!> \param apply Apply covariant observables or metric (operator zero)
!> \param solve Apply inverse metric using already validated factors
!> \param RESULT Normalized eigenstates and separated physical residuals
!> \param status 0 success, 1 unconverged, -1 invalid input, -4 numerical failure
! **************************************************************************************************
   SUBROUTINE quadratic_solve_real(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
                                   apply, solve, RESULT, status)
      INTEGER, INTENT(IN)                                :: n
      REAL(dp), INTENT(IN)                               :: shift(:), weights(:)
      INTEGER, INTENT(IN)                                :: nstates
      LOGICAL, INTENT(IN)                                :: iterative
      INTEGER, INTENT(IN)                                :: max_subspace, max_iter
      REAL(dp), INTENT(IN)                               :: tolerance
 
      PROCEDURE(operator_apply) :: apply
      PROCEDURE(metric_solve) :: solve
      TYPE(quadratic_result_type), INTENT(OUT) :: result
      INTEGER, INTENT(OUT) :: status
 
      CALL quadratic_solve_complex(n, cmplx(shift, 0.0_dp, dp), weights, nstates, iterative, max_subspace, &
                                   max_iter, tolerance, apply, solve, result, status, apply)
   END SUBROUTINE quadratic_solve_real
 
! **************************************************************************************************
!> \brief Solve sum_j w_j^2 (A_j-z_j S)^dagger S^-1 (A_j-z_j S) c = lambda S c.
!> \param n AO/spinor dimension
!> \param shift Complex queries for arbitrary, possibly nonnormal operators
!> \param weights Nonnegative residual scales
!> \param nstates Number of lowest eigenpairs
!> \param iterative Whether to use a bounded restarted block Ritz space
!> \param max_subspace Maximum trial-space size
!> \param max_iter Maximum Ritz iterations
!> \param tolerance Relative S^-1 eigenpair residual threshold
!> \param apply Apply covariant operators or the positive metric (operator zero)
!> \param solve Solve with the already validated positive metric
!> \param RESULT Eigenpairs, residuals and complex expectations; means holds their real parts
!> \param status 0 success, 1 unconverged, -1 input, -4 numerical failure
!> \param adjoint Apply A_j^dagger, not the coefficient-space S-adjoint; positive operator indices
!> \param unitary True for compressions of continuum unitary operators; retain full norms and projection leakage
! **************************************************************************************************
   SUBROUTINE quadratic_solve_complex(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
                                      apply, solve, RESULT, status, adjoint, unitary)
      INTEGER, INTENT(IN)                                :: n
      COMPLEX(dp), INTENT(IN)                            :: shift(:)
      REAL(dp), INTENT(IN)                               :: weights(:)
      INTEGER, INTENT(IN)                                :: nstates
      LOGICAL, INTENT(IN)                                :: iterative
      INTEGER, INTENT(IN)                                :: max_subspace, max_iter
      REAL(dp), INTENT(IN)                               :: tolerance
 
      PROCEDURE(operator_apply) :: apply, adjoint
      PROCEDURE(metric_solve) :: solve
      TYPE(quadratic_result_type), INTENT(OUT) :: result
      INTEGER, INTENT(OUT) :: status
      LOGICAL, INTENT(IN), OPTIONAL :: unitary(:)
 
      TYPE(callback_operator_type) :: operator_backend
 
      operator_backend%apply_callback => apply
      operator_backend%adjoint_callback => adjoint
      operator_backend%solve_callback => solve
      CALL quadratic_solve_core(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
                                operator_backend, result, status, unitary)
   END SUBROUTINE quadratic_solve_complex
 
! **************************************************************************************************
!> \brief Shared quadratic solver using an explicit operator backend.
!> \param n AO/spinor dimension
!> \param shift Complex query values
!> \param weights Nonnegative observable scales
!> \param nstates Number of lowest eigenpairs
!> \param iterative Whether to use a bounded restarted block Ritz space
!> \param max_subspace Maximum trial-space size
!> \param max_iter Maximum Ritz iterations
!> \param tolerance Relative S^-1 eigenpair residual threshold
!> \param operator_backend Operator application and metric-solve backend
!> \param RESULT Eigenpairs, residuals and expectations
!> \param status 0 success, 1 unconverged, -1 input, -4 numerical failure
!> \param unitary Optional continuum-unitary mask
! **************************************************************************************************
   SUBROUTINE quadratic_solve_core(n, shift, weights, nstates, iterative, max_subspace, max_iter, tolerance, &
                                   operator_backend, RESULT, status, unitary)
      INTEGER, INTENT(IN) :: n, nstates, max_subspace, max_iter
      COMPLEX(dp), INTENT(IN) :: shift(:)
      REAL(dp), INTENT(IN) :: weights(:), tolerance
      LOGICAL, INTENT(IN) :: iterative
      CLASS(quadratic_operator_type), INTENT(INOUT) :: operator_backend
      TYPE(quadratic_result_type), INTENT(OUT) :: result
      INTEGER, INTENT(OUT) :: status
      LOGICAL, INTENT(IN), OPTIONAL :: unitary(:)
      COMPLEX(dp), ALLOCATABLE :: v(:, :), sv(:, :), qv(:, :), projected(:, :), sc(:, :), qc(:, :), &
                                  r(:, :), z(:, :), trial(:, :), strial(:, :), work(:, :), invwork(:, :), &
                                  restarted(:, :)
      REAL(dp), ALLOCATABLE :: eigenvalues(:)
      REAL(dp) :: square, scale, norm_before, norm_after, full_square, lost_square
      COMPLEX(dp) :: coefficient
      INTEGER :: capacity, m, old, i, j, p, iteration, seed, initial, pass, keep
      LOGICAL :: accepted
      LOGICAL, ALLOCATABLE :: continuum(:)
 
      status = -1
      IF (n < 1 .OR. nstates < 1 .OR. nstates > n .OR. max_iter < 1) RETURN
      IF (SIZE(shift) < 1 .OR. SIZE(weights) /= SIZE(shift)) RETURN
      ALLOCATE (continuum(SIZE(shift)))
      continuum = .false.
      IF (PRESENT(unitary)) THEN
         IF (SIZE(unitary) /= SIZE(shift)) RETURN
         continuum(:) = unitary
      END IF
      IF (.NOT. all(ieee_is_finite([real(shift, dp), aimag(shift), weights, tolerance]))) RETURN
      IF (any(weights < 0.0_dp) .OR. maxval(weights) <= 0.0_dp .OR. tolerance <= 0.0_dp) RETURN
      capacity = n
      IF (iterative) capacity = min(n, max_subspace)
      IF (capacity < nstates .OR. (capacity < n .AND. capacity < 2*nstates)) RETURN
      ALLOCATE (v(n, capacity), sv(n, capacity), qv(n, capacity), sc(n, nstates), qc(n, nstates), &
                r(n, nstates), z(n, nstates), trial(n, 1), strial(n, 1), work(n, nstates), invwork(n, nstates))
      ALLOCATE (result%states(n, nstates), result%eigenvalues(nstates), result%residuals(nstates), &
                result%errors(SIZE(shift), nstates), result%means(SIZE(shift), nstates), &
                result%complex_means(SIZE(shift), nstates), &
                result%projected_errors(SIZE(shift), nstates), result%leakage(SIZE(shift), nstates))
      result%leakage = 0.0_dp
      m = 0
      seed = 0
      initial = n
      IF (iterative) initial = min(capacity, max(nstates + 2, 2*nstates))
      DO j = 1, initial
         trial = cmplx(0.0_dp, 0.0_dp, dp)
         IF (.NOT. iterative) THEN
            trial(j, 1) = 1.0_dp
         ELSE
            seed = seed + 1
            CALL seed_vector()
         END IF
         CALL append_vector(accepted)
         IF (status /= 0) RETURN
         IF (.NOT. accepted) THEN
            status = -4
            RETURN
         END IF
      END DO
      CALL apply_q(v(:, :m), qv(:, :m), status)
      IF (status /= 0) RETURN
      DO iteration = 1, max_iter
         result%iterations = iteration
         ALLOCATE (projected(m, m), eigenvalues(m))
         projected(:, :) = matmul(conjg(transpose(v(:, :m))), qv(:, :m))
         scale = max(1.0_dp, maxval(abs(projected)))
         IF (maxval(abs(projected - conjg(transpose(projected)))) > 1.e-8_dp*scale) THEN
            status = -4
            RETURN
         END IF
         projected(:, :) = 0.5_dp*(projected + conjg(transpose(projected)))
         CALL diagonalize(projected, eigenvalues, status)
         IF (status /= 0) RETURN
         result%eigenvalues(:) = eigenvalues(:nstates)
         result%states(:, :) = matmul(v(:, :m), projected(:, :nstates))
         sc(:, :) = matmul(sv(:, :m), projected(:, :nstates))
         qc(:, :) = matmul(qv(:, :m), projected(:, :nstates))
         DO pass = 1, 2
            IF (pass == 2) THEN
               ! Cached metric/operator images can drift through repeated Ritz restarts.
               CALL operator_backend%apply(0, result%states, sc, status)
               IF (status /= 0) RETURN
               DO j = 1, nstates
                  square = real(dot_product(result%states(:, j), sc(:, j)), dp)
                  IF (.NOT. ieee_is_finite(square) .OR. square <= 0.0_dp) THEN
                     status = -4
                     RETURN
                  END IF
                  result%states(:, j) = result%states(:, j)/sqrt(square)
               END DO
               CALL operator_backend%apply(0, result%states, sc, status)
               IF (status /= 0) RETURN
               CALL apply_q(result%states, qc, status)
               IF (status /= 0) RETURN
            END IF
            DO j = 1, nstates
               r(:, j) = qc(:, j) - result%eigenvalues(j)*sc(:, j)
            END DO
            CALL operator_backend%solve(r, z, status)
            IF (status /= 0) RETURN
            DO j = 1, nstates
               square = real(dot_product(r(:, j), z(:, j)), dp)
               IF (.NOT. ieee_is_finite(square) .OR. square < -epsilon(square)*scale) THEN
                  status = -4
                  RETURN
               END IF
               result%residuals(j) = sqrt(max(0.0_dp, square))
            END DO
            IF (any(result%residuals > tolerance*max(1.0_dp, abs(result%eigenvalues)))) EXIT
         END DO
         IF (all(result%residuals <= tolerance*max(1.0_dp, abs(result%eigenvalues)))) THEN
            DEALLOCATE (projected, eigenvalues)
            EXIT
         END IF
         IF (.NOT. iterative .OR. m == n .OR. iteration == max_iter) THEN
            status = 1
            RETURN
         END IF
         ! Fill remaining slots with a partial residual block before restarting.
         IF (m == capacity) THEN
            ! Retain nearby Ritz vectors across restarts, important for clustered spectra.
            keep = min(m, max(nstates, capacity/2))
            ALLOCATE (restarted(n, keep))
            restarted(:, :) = matmul(v(:, :m), projected(:, :keep))
            v(:, :keep) = restarted
            restarted(:, :) = matmul(sv(:, :m), projected(:, :keep))
            sv(:, :keep) = restarted
            restarted(:, :) = matmul(qv(:, :m), projected(:, :keep))
            qv(:, :keep) = restarted
            DEALLOCATE (restarted)
            m = keep
         END IF
         DEALLOCATE (projected, eigenvalues)
         old = m
         DO j = 1, nstates
            IF (m == capacity) EXIT
            trial(:, 1) = z(:, j)
            CALL append_vector(accepted)
            IF (status /= 0) RETURN
         END DO
         IF (m == old) THEN
            seed = seed + 1
            CALL seed_vector()
            CALL append_vector(accepted)
            IF (status /= 0) RETURN
            IF (.NOT. accepted) THEN
               status = 1
               RETURN
            END IF
         END IF
         CALL apply_q(v(:, old + 1:m), qv(:, old + 1:m), status)
         IF (status /= 0) RETURN
      END DO
      IF (any(result%eigenvalues < -tolerance*max(1.0_dp, maxval(abs(result%eigenvalues))))) THEN
         status = -4
         RETURN
      END IF
      result%eigenvalues = max(0.0_dp, result%eigenvalues)
      DO p = 1, SIZE(shift)
         CALL operator_backend%apply(p, result%states, work, status)
         IF (status /= 0) RETURN
         DO j = 1, nstates
            result%complex_means(p, j) = dot_product(result%states(:, j), work(:, j))
            result%means(p, j) = real(result%complex_means(p, j), dp)
         END DO
         work(:, :) = work - shift(p)*sc
         CALL operator_backend%solve(work, invwork, status)
         IF (status /= 0) RETURN
         DO j = 1, nstates
            square = real(dot_product(work(:, j), invwork(:, j)), dp)
            IF (.NOT. ieee_is_finite(square) .OR. square < -tolerance) THEN
               status = -4
               RETURN
            END IF
            result%projected_errors(p, j) = sqrt(max(0.0_dp, square))
            IF (continuum(p)) THEN
               full_square = (1.0_dp + abs(shift(p))**2)*real(dot_product(result%states(:, j), sc(:, j)), dp) &
                             - 2.0_dp*real(conjg(shift(p))*result%complex_means(p, j), dp)
               lost_square = full_square - square
               scale = max(1.0_dp, abs(full_square), abs(square))
               IF (.NOT. ieee_is_finite(full_square) .OR. full_square < -tolerance*scale .OR. &
                   lost_square < -tolerance*scale) THEN
                  status = -4
                  RETURN
               END IF
               result%leakage(p, j) = sqrt(max(0.0_dp, lost_square))
               square = full_square
            END IF
            result%errors(p, j) = sqrt(max(0.0_dp, square))
         END DO
      END DO
      DO j = 1, nstates
         square = sum((weights*result%errors(:, j))**2)
         IF (abs(square - result%eigenvalues(j)) > 10.0_dp*tolerance*max(1.0_dp, square)) THEN
            status = -4
            RETURN
         END IF
      END DO
      status = 0
   CONTAINS
! **************************************************************************************************
!> \brief Deterministic complex trial block, independent of compiler random-number generators.
! **************************************************************************************************
      SUBROUTINE seed_vector()
      INTEGER                                            :: k
 
         DO k = 1, n
            trial(k, 1) = cmplx(sin(real(k, dp)*real(seed + 1, dp)*1.4142135623730951_dp), &
                                cos(real(k, dp)*real(seed + 2, dp)*1.7320508075688772_dp), dp)
         END DO
      END SUBROUTINE seed_vector
! **************************************************************************************************
!> \brief Twice-reorthogonalized S-metric expansion, rejecting linearly dependent vectors.
!> \param accepted Whether a new basis vector was appended
! **************************************************************************************************
      SUBROUTINE append_vector(accepted)
      LOGICAL, INTENT(OUT)                               :: accepted
 
         accepted = .false.
         CALL operator_backend%apply(0, trial, strial, status)
         IF (status /= 0) RETURN
         norm_before = real(dot_product(trial(:, 1), strial(:, 1)), dp)
         IF (.NOT. ieee_is_finite(norm_before) .OR. norm_before < 0.0_dp) THEN
            status = -4
            RETURN
         END IF
         DO pass = 1, 2
            DO i = 1, m
               coefficient = dot_product(v(:, i), strial(:, 1))
               trial(:, 1) = trial(:, 1) - coefficient*v(:, i)
               strial(:, 1) = strial(:, 1) - coefficient*sv(:, i)
            END DO
         END DO
         norm_after = real(dot_product(trial(:, 1), strial(:, 1)), dp)
         IF (.NOT. ieee_is_finite(norm_after)) THEN
            status = -4
            RETURN
         END IF
         IF (norm_after <= 1.e-24_dp*max(tiny(1.0_dp), norm_before)) RETURN
         m = m + 1
         v(:, m) = trial(:, 1)/sqrt(norm_after)
         sv(:, m) = strial(:, 1)/sqrt(norm_after)
         accepted = .true.
      END SUBROUTINE append_vector
! **************************************************************************************************
!> \brief Matrix-free covariant normal products, reusing the positive metric factorization.
!> \param block Input AO coefficients
!> \param image Covariant Q times the input
!> \param status Return status
! **************************************************************************************************
      SUBROUTINE apply_q(block, image, status)
      COMPLEX(dp), INTENT(IN)                            :: block(:, :)
      COMPLEX(dp), INTENT(OUT)                           :: image(:, :)
      INTEGER, INTENT(OUT)                               :: status
 
      COMPLEX(dp), ALLOCATABLE                           :: ab(:, :), sb(:, :), second(:, :), &
                                                            solved(:, :), ss(:, :)
      INTEGER                                            :: nrhs, op
 
         nrhs = SIZE(block, 2)
         ALLOCATE (sb(n, nrhs), ab(n, nrhs), solved(n, nrhs), second(n, nrhs), ss(n, nrhs))
         CALL operator_backend%apply(0, block, sb, status)
         IF (status /= 0) RETURN
         image = cmplx(0.0_dp, 0.0_dp, dp)
         DO op = 1, SIZE(shift)
            IF (weights(op) == 0.0_dp) cycle
            CALL operator_backend%apply(op, block, ab, status)
            IF (status /= 0) RETURN
            IF (continuum(op)) THEN
               CALL operator_backend%adjoint(op, block, second, status)
               IF (status /= 0) RETURN
               image = image + weights(op)**2*((1.0_dp + abs(shift(op))**2)*sb &
                                               - conjg(shift(op))*ab - shift(op)*second)
               cycle
            END IF
            ab(:, :) = ab - shift(op)*sb
            CALL operator_backend%solve(ab, solved, status)
            IF (status /= 0) RETURN
            CALL operator_backend%adjoint(op, solved, second, status)
            IF (status /= 0) RETURN
            CALL operator_backend%apply(0, solved, ss, status)
            IF (status /= 0) RETURN
            image = image + weights(op)**2*(second - conjg(shift(op))*ss)
         END DO
      END SUBROUTINE apply_q
   END SUBROUTINE quadratic_solve_core
 
! **************************************************************************************************
!> \brief Forward an operator application to the public callback interface.
!> \param this ...
!> \param OPERATOR ...
!> \param v ...
!> \param y ...
!> \param status ...
! **************************************************************************************************
   SUBROUTINE callback_apply(this, OPERATOR, v, y, status)
      CLASS(callback_operator_type), INTENT(INOUT) :: this
      INTEGER, INTENT(IN) :: operator
      COMPLEX(dp), INTENT(IN) :: v(:, :)
      COMPLEX(dp), INTENT(OUT) :: y(:, :)
      INTEGER, INTENT(OUT) :: status
 
      CALL this%apply_callback(operator, v, y, status)
   END SUBROUTINE callback_apply
 
! **************************************************************************************************
!> \brief Forward an adjoint application to the public callback interface.
!> \param this ...
!> \param OPERATOR ...
!> \param v ...
!> \param y ...
!> \param status ...
! **************************************************************************************************
   SUBROUTINE callback_adjoint(this, OPERATOR, v, y, status)
      CLASS(callback_operator_type), INTENT(INOUT) :: this
      INTEGER, INTENT(IN) :: operator
      COMPLEX(dp), INTENT(IN) :: v(:, :)
      COMPLEX(dp), INTENT(OUT) :: y(:, :)
      INTEGER, INTENT(OUT) :: status
 
      CALL this%adjoint_callback(operator, v, y, status)
   END SUBROUTINE callback_adjoint
 
! **************************************************************************************************
!> \brief Forward a metric solve to the public callback interface.
!> \param this ...
!> \param v ...
!> \param y ...
!> \param status ...
! **************************************************************************************************
   SUBROUTINE callback_solve(this, v, y, status)
      CLASS(callback_operator_type), INTENT(INOUT) :: this
      COMPLEX(dp), INTENT(IN) :: v(:, :)
      COMPLEX(dp), INTENT(OUT) :: y(:, :)
      INTEGER, INTENT(OUT) :: status
 
      CALL this%solve_callback(v, y, status)
   END SUBROUTINE callback_solve
 
! **************************************************************************************************
!> \brief Hermiticity and finiteness check for the dense reference.
!> \param matrix Input square matrix
!> \return Whether the matrix is finite Hermitian
! **************************************************************************************************
   FUNCTION hermitian(matrix) RESULT(valid)
      COMPLEX(dp), INTENT(IN)                            :: matrix(:, :)
      LOGICAL                                            :: valid
 
      valid = .false.
      IF (.NOT. all(ieee_is_finite(real(matrix, dp))) .OR. .NOT. all(ieee_is_finite(aimag(matrix)))) RETURN
      valid = maxval(abs(matrix - conjg(transpose(matrix)))) <= 1.e-10_dp*max(tiny(1.0_dp), maxval(abs(matrix)))
   END FUNCTION hermitian
 
! **************************************************************************************************
!> \brief Small dense Hermitian eigensolve shared by the oracle and Ritz projection.
!> \param matrix Overwritten with eigenvectors
!> \param values Ascending eigenvalues
!> \param status Return status
! **************************************************************************************************
   SUBROUTINE diagonalize(matrix, values, status)
      COMPLEX(dp), INTENT(INOUT) :: matrix(:, :)
      REAL(dp), INTENT(OUT) :: values(:)
      INTEGER, INTENT(OUT) :: status
      COMPLEX(dp), ALLOCATABLE :: workspace(:)
      REAL(dp), ALLOCATABLE :: rwork(:)
      INTEGER :: n
#if defined(__HAS_IEEE_EXCEPTIONS)
      LOGICAL :: halt(SIZE(IEEE_ALL))
#endif
      n = SIZE(values)
      ALLOCATE (workspace(max(1, 2*n)), rwork(max(1, 3*n - 2)))
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_get_halting_mode(ieee_all, halt)
      CALL ieee_set_halting_mode(ieee_all, .false.)
#endif
      CALL zheev('V', 'L', n, matrix, n, values, workspace, SIZE(workspace), rwork, status)
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_set_halting_mode(ieee_all, halt)
#endif
      IF (status /= 0 .OR. .NOT. all(ieee_is_finite(values))) status = -4
   END SUBROUTINE diagonalize
END MODULE quadratic_pseudospectrum