quadratic_pseudospectrum_unittest.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                                                      !
!--------------------------------------------------------------------------------------------------!
PROGRAM quadratic_pseudospectrum_unittest
   USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
#if defined(__HAS_IEEE_EXCEPTIONS)
   USE ieee_exceptions, ONLY: ieee_all, ieee_get_halting_mode, ieee_set_halting_mode
#endif
   USE kinds, ONLY: dp
   USE libcp2k, ONLY: cp2k_init, cp2k_finalize
   USE quadratic_pseudospectrum, ONLY: quadratic_dense, &
                                       quadratic_result_type
 
   IMPLICIT NONE
   COMPLEX(dp) :: op(2, 2, 4), s(2, 2), change(2, 2), rotated(2, 2, 4), metric(2, 2)
   COMPLEX(dp), ALLOCATABLE :: defect(:, :, :), overlap(:, :), projected(:, :)
   REAL(dp) :: shift(4), weights(4), expected, center
   TYPE(quadratic_result_type) :: dense, iterative, transformed
   INTEGER :: status, i, n, io_unit
 
   CALL cp2k_init()
   io_unit = cp_logger_get_default_io_unit()
   s = cmplx(0.0_dp, 0.0_dp, dp)
   s(1, 1) = 1.0_dp
   s(2, 2) = 1.0_dp
   op = cmplx(0.0_dp, 0.0_dp, dp)
   op(1, 1, 1) = 1.0_dp
   op(2, 2, 1) = -1.0_dp
   op(1, 2, 2) = cmplx(0.0_dp, 1.0_dp, dp)
   op(2, 1, 2) = cmplx(0.0_dp, -1.0_dp, dp)
   shift = [0.3_dp, 0.2_dp, 0.0_dp, 0.0_dp]
   weights = [1.0_dp, 0.7_dp, 0.7_dp, 0.7_dp]
   expected = 1.0_dp + shift(1)**2 + weights(2)**2*(1.0_dp + shift(2)**2) &
              - 2.0_dp*sqrt(shift(1)**2 + weights(2)**4*shift(2)**2)
   CALL quadratic_dense(op, s, shift, weights, 1, .false., 2, 200, 1.e-10_dp, 1.e-10_dp, dense, status)
   IF (status /= 0) error stop 'Quadratic Pauli oracle failed'
   IF (abs(dense%eigenvalues(1) - expected) > 1.e-12_dp) error stop 'Analytic quadratic eigenvalue mismatch'
   CALL quadratic_dense(op, s, shift, weights, 1, .true., 2, 200, 1.e-10_dp, 1.e-10_dp, iterative, status)
   IF (status /= 0) error stop 'Quadratic Pauli iteration failed'
   IF (maxval(abs(dense%eigenvalues - iterative%eigenvalues)) > 1.e-12_dp) error stop 'Pauli Ritz mismatch'
   change = reshape([cmplx(1.2_dp, 0.1_dp, dp), cmplx(0.15_dp, -0.1_dp, dp), &
                     cmplx(0.1_dp, 0.2_dp, dp), cmplx(0.8_dp, 0.0_dp, dp)], [2, 2])
   metric = matmul(conjg(transpose(change)), change)
   DO i = 1, 4
      rotated(:, :, i) = matmul(conjg(transpose(change)), matmul(op(:, :, i), change))
   END DO
   CALL quadratic_dense(rotated, metric, shift, weights, 1, .false., 2, 200, 1.e-10_dp, 1.e-10_dp, &
                        transformed, status)
   IF (status /= 0) error stop 'Nonorthogonal complex basis failed'
   IF (maxval(abs(dense%eigenvalues - transformed%eigenvalues)) > 1.e-12_dp) error stop 'Basis covariance failed'
   IF (maxval(abs(dense%errors - transformed%errors)) > 1.e-12_dp) error stop 'Observable residual covariance failed'
   IF (maxval(abs(dense%means - transformed%means)) > 1.e-12_dp) error stop 'Observable mean covariance failed'
   DO i = 1, 4
      op(:, :, i) = op(:, :, i) + real(i, dp)*s
      shift(i) = shift(i) + real(i, dp)
   END DO
   CALL quadratic_dense(op, s, shift, weights, 1, .false., 2, 200, 1.e-10_dp, 1.e-10_dp, transformed, status)
   IF (status /= 0) error stop 'Shifted quadratic oracle failed'
   IF (maxval(abs(dense%eigenvalues - transformed%eigenvalues)) > 1.e-12_dp) error stop 'Origin invariance failed'
   IF (maxval(abs(dense%errors - transformed%errors)) > 1.e-12_dp) error stop 'Residual origin invariance failed'
 
   n = 40
   center = 20.0_dp
   ALLOCATE (defect(n, n, 4), overlap(n, n), projected(3, 3))
   defect = cmplx(0.0_dp, 0.0_dp, dp)
   overlap = cmplx(0.0_dp, 0.0_dp, dp)
   DO i = 1, n
      overlap(i, i) = 1.0_dp
      defect(i, i, 1) = 5.0_dp
      defect(i, i, 2) = real(i, dp) - center
      IF (i < n) THEN
         defect(i, i + 1, 1) = cmplx(-0.1_dp, 0.02_dp, dp)
         defect(i + 1, i, 1) = conjg(defect(i, i + 1, 1))
      END IF
   END DO
   defect(20, 20, 1) = 0.0_dp
   shift = 0.0_dp
   weights = [1.0_dp, 0.2_dp, 0.2_dp, 0.2_dp]
   CALL quadratic_dense(defect, overlap, shift, weights, 3, .false., n, 500, 1.e-9_dp, 1.e-10_dp, dense, status)
   IF (status /= 0) error stop 'Defect oracle failed'
   IF (abs(dense%states(20, 1))**2 < 0.95_dp) error stop 'Defect state not localized'
   CALL quadratic_dense(defect, overlap, shift, weights, 3, .true., 20, 500, 1.e-9_dp, 1.e-10_dp, iterative, status)
   IF (status /= 0) error stop 'Defect iteration failed'
   IF (maxval(abs(dense%eigenvalues - iterative%eigenvalues)) > 1.e-8_dp) error stop 'Defect Ritz eigenvalue mismatch'
   projected(:, :) = matmul(conjg(transpose(dense%states)), iterative%states)
   IF (sum(abs(projected)**2) < 3.0_dp - 1.e-7_dp) error stop 'Defect invariant subspace mismatch'
   CALL quadratic_dense(defect, overlap, shift, weights, 3, .true., 12, 1, 1.e-12_dp, 1.e-10_dp, iterative, status)
   IF (status /= 1) error stop 'Unconverged quadratic solve accepted'
   ! Exactly degenerate joint states: compare the projector, not individual eigenvectors.
   defect = cmplx(0.0_dp, 0.0_dp, dp)
   DO i = 1, n
      defect(i, i, 1) = 2.0_dp
   END DO
   defect(19, 19, 1) = 0.0_dp
   defect(20, 20, 1) = 0.0_dp
   CALL quadratic_dense(defect, overlap, shift, weights, 2, .true., 8, 200, 1.e-10_dp, 1.e-10_dp, iterative, status)
   IF (status /= 0) error stop 'Degenerate quadratic iteration failed'
   IF (maxval(abs(iterative%eigenvalues)) > 1.e-12_dp) error stop 'Degenerate eigenvalues failed'
   IF (abs(sum(abs(iterative%states(19:20, :))**2) - 2.0_dp) > 1.e-12_dp) error stop 'Degenerate projector failed'
   overlap(1, 1) = -1.0_dp
   CALL quadratic_dense(defect, overlap, shift, weights, 3, .false., n, 200, 1.e-9_dp, 1.e-10_dp, dense, status)
   IF (status /= -2) error stop 'Indefinite metric accepted'
   overlap(1, 1) = 1.e-14_dp
   CALL quadratic_dense(defect, overlap, shift, weights, 3, .false., n, 200, 1.e-9_dp, 1.e-10_dp, dense, status)
   IF (status /= -2) error stop 'Ill-conditioned metric accepted'
   DEALLOCATE (defect, overlap, projected)
   CALL periodic_tests()
   CALL translation_tests()
   CALL folded_translation_tests()
   CALL compressed_translation_tests()
   CALL cp2k_finalize()
CONTAINS
! **************************************************************************************************
!> \brief Compression of a three-state unitary to two AOs retains the missing normal Gram exactly.
! **************************************************************************************************
   SUBROUTINE compressed_translation_tests()
      COMPLEX(dp)                                        :: change(2, 2), changed(2, 2, 2), &
                                                            metric(2, 2), op(2, 2, 2), query(2), &
                                                            transformed(2, 2)
      INTEGER                                            :: i, status
      REAL(dp)                                           :: center, diagonal(2), err, expected(2), &
                                                            split, weights(2)
      TYPE(quadratic_result_type)                        :: full, projected, rotated
 
      metric = 0
      metric(1, 1) = 1
      metric(2, 2) = 1
      op = 0
      op(1, 1, 1) = 0.2_dp
      op(2, 2, 1) = -0.4_dp
      op(1, 2, 2) = 1
      query = [cmplx(0.1_dp, 0.0_dp, dp), exp(cmplx(0.0_dp, 0.73_dp, dp))]
      weights = [1.0_dp, 0.7_dp]
      diagonal = [0.1_dp**2, 0.5_dp**2] + 2*weights(2)**2
      center = sum(diagonal)/2
      split = sqrt(((diagonal(1) - diagonal(2))/2)**2 + weights(2)**4)
      expected = [center - split, center + split]
      CALL quadratic_dense(op, metric, query, weights, 2, .false., 2, 100, 1.e-11_dp, 1.e-10_dp, &
                           full, status, unitary=[.false., .true.])
      IF (status /= 0) error stop 'Compressed unitary solve failed'
      err = maxval(abs(full%eigenvalues - expected))
      IF (err > 1.e-12_dp) error stop 'Full translation normal Gram differs from analytic result'
      IF (maxval(abs(full%errors**2 - full%projected_errors**2 - full%leakage**2)) > 1.e-12_dp) THEN
         error stop 'Continuum/projected residual decomposition failed'
      END IF
      IF (maxval(abs(full%leakage(2, :)**2 - abs(full%states(1, :))**2)) > 1.e-12_dp) THEN
         error stop 'Translation leakage differs from discarded ambient component'
      END IF
      CALL quadratic_dense(op, metric, query, weights, 2, .false., 2, 100, 1.e-11_dp, 1.e-10_dp, projected, status)
      IF (status /= 0) error stop 'Projected unitary comparison failed'
      IF (maxval(abs(projected%eigenvalues - full%eigenvalues)) < 0.1_dp) error stop 'Projection negative control failed'
      change = reshape([cmplx(1.2_dp, 0.1_dp, dp), cmplx(0.15_dp, -0.1_dp, dp), &
                        cmplx(0.1_dp, 0.2_dp, dp), cmplx(0.8_dp, 0.0_dp, dp)], [2, 2])
      transformed = matmul(conjg(transpose(change)), change)
      DO i = 1, 2
         changed(:, :, i) = matmul(conjg(transpose(change)), matmul(op(:, :, i), change))
      END DO
      CALL quadratic_dense(changed, transformed, query, weights, 2, .true., 2, 100, 1.e-11_dp, 1.e-10_dp, &
                           rotated, status, unitary=[.false., .true.])
      IF (status /= 0) error stop 'Compressed unitary metric iteration failed'
      IF (maxval(abs(rotated%eigenvalues - full%eigenvalues)) > 1.e-12_dp) error stop 'Unitary metric covariance failed'
      IF (maxval(abs(rotated%leakage - full%leakage)) > 1.e-12_dp) error stop 'Leakage metric covariance failed'
      op(:, :, 2) = 2*metric
      CALL quadratic_dense(op, metric, query, weights, 2, .false., 2, 100, 1.e-11_dp, 1.e-10_dp, &
                           rotated, status, unitary=[.false., .true.])
      IF (status /= -4) error stop 'Invalid expanding unitary compression accepted'
      IF (io_unit > 0) WRITE (io_unit, '(A,ES15.6)') 'Compressed unitary full-normal analytic error: ', err
   END SUBROUTINE compressed_translation_tests
 
! **************************************************************************************************
!> \brief A trimerized periodic chain must equal its independently assembled three-band Bloch blocks.
! **************************************************************************************************
   SUBROUTINE folded_translation_tests()
      INTEGER, PARAMETER                                 :: cells = 6, n = 3*cells
 
      COMPLEX(dp)                                        :: bloch(3, 3, 2), identity(3, 3), &
                                                            metric(n, n), operators(n, n, 2), &
                                                            phase, query(2)
      INTEGER                                            :: first, i, ie, ik, iq, j, loc, status
      REAL(dp)                                           :: energy(5), expected(n), hopping, &
                                                            joined(n), k, maximum_error, q, &
                                                            weights(2)
      TYPE(quadratic_result_type)                        :: full, sector
 
      operators = 0
      metric = 0
      identity = 0
      DO i = 1, 3
         identity(i, i) = 1
      END DO
      DO i = 1, n
         metric(i, i) = 1
         j = modulo(i, n) + 1
         operators(i, j, 2) = 1
         hopping = 0.8_dp
         IF (modulo(i, 3) == 0) hopping = 1.2_dp
         operators(i, j, 1) = -hopping
         operators(j, i, 1) = -hopping
      END DO
      energy = [-2.2_dp, -1.0_dp, 0.0_dp, 0.7_dp, 2.0_dp]
      weights = [1.0_dp, 0.7_dp]
      maximum_error = 0
      DO ik = 0, 12
         k = -acos(-1.0_dp) + 2*acos(-1.0_dp)*real(ik, dp)/13
         DO ie = 1, SIZE(energy)
            query = [cmplx(energy(ie), 0.0_dp, dp), exp(cmplx(0.0_dp, k, dp))]
            CALL quadratic_dense(operators, metric, query, weights, n, .false., n, 100, &
                                 1.e-10_dp, 1.e-10_dp, full, status)
            IF (status /= 0) error stop 'Full trimerized quadratic scan failed'
            DO iq = 0, cells - 1
               q = 2*acos(-1.0_dp)*real(iq, dp)/cells
               phase = exp(cmplx(0.0_dp, q, dp))
               bloch = 0
               bloch(1, 2, 2) = 1
               bloch(2, 3, 2) = 1
               bloch(3, 1, 2) = phase
               bloch(1, 2, 1) = -0.8_dp
               bloch(2, 3, 1) = -0.8_dp
               bloch(3, 1, 1) = -1.2_dp*phase
               bloch(:, :, 1) = bloch(:, :, 1) + conjg(transpose(bloch(:, :, 1)))
               CALL quadratic_dense(bloch, identity, query, weights, 3, .false., 3, 100, &
                                    1.e-10_dp, 1.e-10_dp, sector, status)
               IF (status /= 0) error stop 'Bloch-block quadratic scan failed'
               first = 3*iq + 1
               joined(first:first + 2) = sector%eigenvalues
            END DO
            DO i = 1, n
               loc = minloc(joined, dim=1)
               expected(i) = joined(loc)
               joined(loc) = huge(1.0_dp)
            END DO
            maximum_error = max(maximum_error, maxval(abs(full%eigenvalues - expected)))
         END DO
      END DO
      IF (maximum_error > 1.e-11_dp) error stop 'Full/Bloch trimerized quadratic spectra differ'
      IF (io_unit > 0) WRITE (io_unit, '(A,I0,A,ES15.6)') 'Trimerized full/Bloch spectra: ', 13*SIZE(energy)*n, &
         ' eigenvalues; maximum error ', maximum_error
   END SUBROUTINE folded_translation_tests
 
! **************************************************************************************************
!> \brief Exact Bloch phases, nonnormal open translations and an independent generalized Gram oracle.
! **************************************************************************************************
   SUBROUTINE translation_tests()
      INTEGER, PARAMETER                                 :: n = 18
 
      COMPLEX(dp) :: changed(n, n, 2), identity(n, n), metric(n, n), offset(2), &
         operators(n, n, 2), projector(n, n), query(2), split(n, n, 3), transform(n, n)
      INTEGER                                            :: i, j, status
      REAL(dp)                                           :: covariance_error, expected(n), &
                                                            gap_error, hopping, mean_error, &
                                                            nonnormal_error, phase, weights(2)
      TYPE(quadratic_result_type)                        :: altered, incorrect, iterative, original
 
      identity = 0
      operators = 0
      DO i = 1, n
         identity(i, i) = 1
         operators(i, modulo(i, n) + 1, 2) = 1
      END DO
      operators(:, :, 1) = -operators(:, :, 2) - conjg(transpose(operators(:, :, 2)))
      phase = 2*acos(-1.0_dp)*5/real(n, dp)
      query = [cmplx(-2*cos(phase) + 0.031_dp, 0.0_dp, dp), exp(cmplx(0.0_dp, phase, dp))]
      weights = [1.0_dp, 1.3_dp]
      CALL quadratic_dense(operators, identity, query, weights, 1, .false., n, 200, &
                           1.e-10_dp, 1.e-10_dp, original, status)
      IF (status /= 0) error stop 'Unitary translation quadratic solve failed'
      IF (abs(original%eigenvalues(1) - 0.031_dp**2) > 1.e-12_dp) error stop 'Exact Bloch energy/phase failed'
      IF (abs(original%complex_means(2, 1) - query(2)) > 1.e-12_dp) error stop 'Bloch translation mean lost phase'
      IF (original%errors(2, 1) > 1.e-12_dp) error stop 'Bloch translation residual failed'
      CALL quadratic_dense(operators, identity, real(query, dp), weights, 1, .false., n, 200, &
                           1.e-10_dp, 1.e-10_dp, altered, status)
      IF (status /= -1) error stop 'Hermitian interface accepted non-Hermitian translation'
 
      ! The truncated translation is nonnormal. Its adjoint cannot be replaced by itself
      ! or by an independent sum of squares of its Hermitian real/imaginary parts.
      operators(n, 1, 2) = 0
      operators(:, :, 1) = 0
      DO i = 1, n - 1
         hopping = 0.8_dp
         IF (modulo(i, 3) == 0) hopping = 1.2_dp
         operators(i, i + 1, 1) = -hopping
         operators(i + 1, i, 1) = -hopping
      END DO
      operators(4, 4, 1) = 0.17_dp
      query = [cmplx(0.33_dp, 0.0_dp, dp), exp(cmplx(0.0_dp, 0.73_dp, dp))]
      CALL quadratic_dense(operators, identity, query, weights, n, .false., n, 200, &
                           1.e-10_dp, 1.e-10_dp, original, status)
      IF (status /= 0) error stop 'Open trimerized translation solve failed'
      CALL gram_reference(operators, identity, query, weights, expected)
      gap_error = maxval(abs(original%eigenvalues - expected))
      IF (gap_error > 1.e-11_dp) error stop 'General quadratic Gram oracle mismatch'
      CALL quadratic_dense(operators, identity, query, weights, 2, .true., 10, 2000, &
                           1.e-10_dp, 1.e-10_dp, iterative, status)
      IF (status /= 0) error stop 'Nonnormal restarted iteration failed'
      IF (maxval(abs(iterative%eigenvalues - expected(:2))) > 1.e-10_dp) error stop 'Nonnormal Ritz spectrum mismatch'
      projector = matmul(iterative%states, conjg(transpose(iterative%states))) - &
                  matmul(original%states(:, :2), conjg(transpose(original%states(:, :2))))
      IF (maxval(abs(projector)) > 1.e-8_dp) error stop 'Nonnormal invariant-subspace mismatch'
      split(:, :, 1) = operators(:, :, 1)
      split(:, :, 2) = 0.5_dp*(operators(:, :, 2) + conjg(transpose(operators(:, :, 2))))
      split(:, :, 3) = cmplx(0.0_dp, -0.5_dp, dp)*(operators(:, :, 2) - conjg(transpose(operators(:, :, 2))))
      CALL quadratic_dense(split, identity, [real(query(1), dp), real(query(2), dp), aimag(query(2))], &
                           [weights(1), weights(2), weights(2)], n, .false., n, 200, &
                           1.e-10_dp, 1.e-10_dp, incorrect, status)
      IF (status /= 0) error stop 'Hermitian-split comparison solve failed'
      nonnormal_error = maxval(abs(original%eigenvalues - incorrect%eigenvalues))
      IF (nonnormal_error < 0.01_dp) error stop 'Open-translation counterexample lost nonnormality'
 
      transform = identity
      DO i = 1, n - 1
         transform(i, i) = 1 + 0.02_dp*i
         transform(i, i + 1) = cmplx(0.21_dp, -0.12_dp, dp)
      END DO
      metric = matmul(conjg(transpose(transform)), transform)
      DO j = 1, 2
         changed(:, :, j) = matmul(conjg(transpose(transform)), matmul(operators(:, :, j), transform))
      END DO
      CALL quadratic_dense(changed, metric, query, weights, n, .false., n, 200, &
                           1.e-10_dp, 1.e-10_dp, altered, status)
      IF (status /= 0) error stop 'Nonnormal complex AO basis failed'
      CALL gram_reference(changed, metric, query, weights, expected)
      covariance_error = maxval(abs(original%eigenvalues - altered%eigenvalues))
      IF (covariance_error > 1.e-11_dp) error stop 'Nonnormal metric covariance failed'
      IF (maxval(abs(altered%eigenvalues - expected)) > 1.e-11_dp) error stop 'Nonorthogonal Gram oracle mismatch'
      mean_error = maxval(abs(original%complex_means - altered%complex_means))
      IF (mean_error > 1.e-10_dp) error stop 'Complex expectation covariance failed'
      IF (maxval(abs(original%errors - altered%errors)) > 1.e-10_dp) error stop 'Complex residual covariance failed'
      offset = [cmplx(0.21_dp, -0.12_dp, dp), cmplx(-0.43_dp, 0.31_dp, dp)]
      DO j = 1, 2
         changed(:, :, j) = changed(:, :, j) + offset(j)*metric
      END DO
      CALL quadratic_dense(changed, metric, query + offset, weights, n, .false., n, 200, &
                           1.e-10_dp, 1.e-10_dp, altered, status)
      IF (status /= 0) error stop 'Complex operator/query origin solve failed'
      IF (maxval(abs(original%eigenvalues - altered%eigenvalues)) > 1.e-11_dp) error stop 'Complex origin covariance failed'
      DO j = 1, 2
         IF (maxval(abs(altered%complex_means(j, :) - original%complex_means(j, :) - offset(j))) > 1.e-10_dp) THEN
            error stop 'Complex expectation origin covariance failed'
         END IF
      END DO
      IF (io_unit > 0) WRITE (io_unit, '(A,4ES15.6)') 'Translation Gram/metric/complex-mean/split-error: ', &
         gap_error, covariance_error, mean_error, nonnormal_error
      changed(:, :, 1) = matmul(conjg(transpose(transform)), matmul(operators(:, :, 1), transform))
      changed(:, :, 2) = metric
      query(2) = 1
      CALL quadratic_dense(changed, metric, query, weights, 2, .true., 10, 2000, &
                           1.e-10_dp, 1.e-10_dp, iterative, status, unitary=[.false., .true.])
      IF (status /= 0) error stop 'Restarted identity translation failed'
      IF (iterative%iterations < 2) error stop 'Identity translation did not exercise restarts'
      IF (maxval(iterative%errors(2, :)) > 5.e-8_dp) error stop 'Identity full residual exceeds roundoff'
      IF (maxval(iterative%projected_errors(2, :)) > 1.e-12_dp) error stop 'Identity projected residual drift'
      IF (maxval(iterative%leakage(2, :)) > 5.e-8_dp) error stop 'Identity projection leakage exceeds roundoff'
   END SUBROUTINE translation_tests
 
! **************************************************************************************************
!> \brief Independent explicit normal-product reference using LU solves and a generalized eigensolve.
!> \param operators Covariant operators
!> \param metric Positive AO overlap
!> \param query Complex operator queries
!> \param weights Nonnegative residual scales
!> \param values Complete ascending generalized spectrum
! **************************************************************************************************
   SUBROUTINE gram_reference(operators, metric, query, weights, values)
      COMPLEX(dp), INTENT(IN) :: operators(:, :, :), metric(:, :), query(:)
      REAL(dp), INTENT(IN) :: weights(:)
      REAL(dp), INTENT(OUT) :: values(:)
      COMPLEX(dp), ALLOCATABLE :: gram(:, :), shifted(:, :), solved(:, :), lu(:, :), work(:)
      REAL(dp), ALLOCATABLE :: rwork(:)
      INTEGER, ALLOCATABLE :: pivots(:)
      INTEGER :: n, i, status
#if defined(__HAS_IEEE_EXCEPTIONS)
      LOGICAL :: halt(SIZE(IEEE_ALL))
#endif
 
      n = SIZE(metric, 1)
      ALLOCATE (gram(n, n), shifted(n, n), solved(n, n), lu(n, n), pivots(n), work(4*n), rwork(4*n))
      gram = 0
      DO i = 1, SIZE(query)
         lu(:, :) = metric
         shifted(:, :) = operators(:, :, i) - query(i)*metric
         solved(:, :) = shifted
#if defined(__HAS_IEEE_EXCEPTIONS)
         CALL ieee_get_halting_mode(ieee_all, halt)
         CALL ieee_set_halting_mode(ieee_all, .false.)
#endif
         CALL zgesv(n, n, lu, n, pivots, solved, n, status)
#if defined(__HAS_IEEE_EXCEPTIONS)
         CALL ieee_set_halting_mode(ieee_all, halt)
#endif
         IF (status /= 0) error stop 'Explicit Gram metric solve failed'
         gram(:, :) = gram + weights(i)**2*matmul(conjg(transpose(shifted)), solved)
      END DO
      IF (maxval(abs(gram - conjg(transpose(gram)))) > 1.e-11_dp) error stop 'Explicit Gram is not Hermitian'
      lu(:, :) = metric
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_get_halting_mode(ieee_all, halt)
      CALL ieee_set_halting_mode(ieee_all, .false.)
#endif
      CALL zhegv(1, 'N', 'L', n, gram, n, lu, n, values, work, SIZE(work), rwork, status)
#if defined(__HAS_IEEE_EXCEPTIONS)
      CALL ieee_set_halting_mode(ieee_all, halt)
#endif
      IF (status /= 0) error stop 'Explicit generalized Gram spectrum failed'
   END SUBROUTINE gram_reference
 
! **************************************************************************************************
!> \brief Three-dimensional trigonometric coordinates, basis covariance and torus-query invariance.
! **************************************************************************************************
   SUBROUTINE periodic_tests()
      INTEGER, PARAMETER                                 :: n = 8
 
      COMPLEX(dp)                                        :: changed(n, n, 7), identity(n, n), &
                                                            metric(n, n), operators(n, n, 7), &
                                                            transform(n, n)
      INTEGER                                            :: i, j, k, status
      REAL(dp)                                           :: diagonal(n), expected(n), move(3), &
                                                            radius(3), shift(7), theta(3), value, &
                                                            weights(7)
      TYPE(quadratic_result_type)                        :: altered, iterative, original
 
      identity = 0
      operators = 0
      theta = [0.21_dp, -0.37_dp, 1.43_dp]
      radius = [0.7_dp, 1.2_dp, 0.9_dp]
      shift(1) = 0.73_dp
      weights(1) = 1
      DO j = 1, 3
         shift(2*j:2*j + 1) = [cos(theta(j)), sin(theta(j))]
         weights(2*j:2*j + 1) = radius(j)
      END DO
      DO i = 1, n
         identity(i, i) = 1
         operators(i, i, 1) = 0.2_dp*i
         DO j = 1, 3
            value = 2*acos(-1.0_dp)*modulo(real(i*(2*j - 1), dp)/real(n + j, dp), 1.0_dp)
            operators(i, i, 2*j) = cos(value)
            operators(i, i, 2*j + 1) = sin(value)
         END DO
         diagonal(i) = sum((weights*abs(operators(i, i, :) - shift))**2)
      END DO
      DO i = 1, n
         k = minloc(diagonal, dim=1)
         expected(i) = diagonal(k)
         diagonal(k) = huge(1.0_dp)
      END DO
      CALL quadratic_dense(operators, identity, shift, weights, n, .false., n, 100, 1.e-10_dp, 1.e-10_dp, original, status)
      IF (status /= 0) error stop 'Periodic quadratic diagonal oracle failed'
      IF (maxval(abs(original%eigenvalues - expected)) > 1.e-12_dp) error stop 'Seven-observable analytic spectrum failed'
      DO i = 1, n - 1
         operators(i, i + 1, 1) = cmplx(-0.12_dp, 0.04_dp, dp)
         operators(i + 1, i, 1) = conjg(operators(i, i + 1, 1))
      END DO
      ! The six initial vectors leave only one or two slots for the three-vector residual block.
      DO k = n - 1, n
         CALL quadratic_dense(operators(:k, :k, :), identity(:k, :k), shift, weights, 3, .false., k, 100, &
                              1.e-10_dp, 1.e-10_dp, original, status)
         IF (status /= 0) error stop 'Partial residual block reference failed'
         CALL quadratic_dense(operators(:k, :k, :), identity(:k, :k), shift, weights, 3, .true., k, 2, &
                              1.e-10_dp, 1.e-10_dp, iterative, status)
         IF (status /= 0) error stop 'Partial residual block did not fill the available space'
         IF (iterative%iterations /= 2) error stop 'Partial residual block iteration count mismatch'
         IF (maxval(abs(original%eigenvalues - iterative%eigenvalues)) > 1.e-12_dp) THEN
            error stop 'Partial residual block eigenvalue mismatch'
         END IF
      END DO
      ! Reuse the full-size reference from the final partial-residual-block check.
      transform = identity
      DO i = 1, n - 1
         transform(i, i + 1) = cmplx(0.13_dp, -0.08_dp, dp)
      END DO
      metric = matmul(conjg(transpose(transform)), transform)
      DO j = 1, 7
         changed(:, :, j) = matmul(conjg(transpose(transform)), matmul(operators(:, :, j), transform))
      END DO
      CALL quadratic_dense(changed, metric, shift, weights, 3, .false., n, 100, 1.e-10_dp, 1.e-10_dp, altered, status)
      IF (status /= 0) error stop 'Periodic complex metric failed'
      IF (maxval(abs(original%eigenvalues - altered%eigenvalues)) > 1.e-12_dp) error stop 'Periodic metric covariance failed'
      CALL quadratic_dense(changed, metric, shift, weights, 3, .true., n, 2, 1.e-10_dp, 1.e-10_dp, iterative, status)
      IF (status /= 0) error stop 'Partial residual block with complex metric failed'
      IF (maxval(abs(altered%eigenvalues - iterative%eigenvalues)) > 1.e-12_dp) THEN
         error stop 'Partial residual block metric covariance failed'
      END IF
      move = [0.43_dp, -0.67_dp, 0.31_dp]
      changed = operators
      DO j = 1, 3
         changed(:, :, 2*j) = cos(move(j))*operators(:, :, 2*j) - sin(move(j))*operators(:, :, 2*j + 1)
         changed(:, :, 2*j + 1) = sin(move(j))*operators(:, :, 2*j) + cos(move(j))*operators(:, :, 2*j + 1)
         shift(2*j:2*j + 1) = [cos(theta(j) + move(j)), sin(theta(j) + move(j))]
      END DO
      CALL quadratic_dense(changed, identity, shift, weights, 3, .false., n, 100, 1.e-10_dp, 1.e-10_dp, altered, status)
      IF (status /= 0) error stop 'Shifted periodic coordinates failed'
  IF (maxval(abs(original%eigenvalues - altered%eigenvalues)) > 1.e-12_dp) error stop 'Periodic coordinate-origin covariance failed'
      DO j = 1, 3
         value = theta(j) + 2*acos(-1.0_dp)*(j - 2)
         shift(2*j:2*j + 1) = [cos(value), sin(value)]
      END DO
      CALL quadratic_dense(operators, identity, shift, weights, 3, .false., n, 100, 1.e-10_dp, 1.e-10_dp, altered, status)
      IF (status /= 0) error stop 'Periodic query wrapping failed'
      IF (maxval(abs(original%eigenvalues - altered%eigenvalues)) > 1.e-12_dp) error stop 'Periodic query invariance failed'
      CALL quadratic_dense(operators, identity, shift, weights(:6), 1, .false., n, 100, 1.e-10_dp, 1.e-10_dp, altered, status)
      IF (status /= -1) error stop 'Mismatched observable weights accepted'
      CALL quadratic_dense(operators(:, :, :0), identity, shift(:0), weights(:0), 1, .false., n, 100, &
                           1.e-10_dp, 1.e-10_dp, altered, status)
      IF (status /= -1) error stop 'Empty observable tuple accepted'
      IF (io_unit > 0) WRITE (io_unit, *) 'Periodic quadratic tuple, metric, iteration and origin checks passed.'
   END SUBROUTINE periodic_tests
END PROGRAM quadratic_pseudospectrum_unittest