diff --git a/acb_mat.h b/acb_mat.h
index dce02deb..af3c6db0 100644
--- a/acb_mat.h
+++ b/acb_mat.h
@@ -94,6 +94,8 @@ void acb_mat_set_round_arb_mat(acb_mat_t dest, const arb_mat_t src, slong prec);
void acb_mat_randtest(acb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits);
+void acb_mat_randtest_eig(acb_mat_t A, flint_rand_t state, acb_srcptr E, slong prec);
+
/* I/O */
void acb_mat_fprintd(FILE * file, const acb_mat_t mat, slong digits);
@@ -406,6 +408,8 @@ void acb_mat_det(acb_t det, const acb_mat_t A, slong prec);
int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec);
+void acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec);
+
/* Special functions */
void acb_mat_exp_taylor_sum(acb_mat_t S, const acb_mat_t A, slong N, slong prec);
diff --git a/acb_mat/eig_global_enclosure.c b/acb_mat/eig_global_enclosure.c
new file mode 100644
index 00000000..c61c8308
--- /dev/null
+++ b/acb_mat/eig_global_enclosure.c
@@ -0,0 +1,63 @@
+/*
+ Copyright 2018 Fredrik Johansson
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "acb_mat.h"
+
+void
+acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec)
+{
+ acb_mat_t Y, R1, R2;
+ slong i, j, n;
+ mag_t r1, r2;
+
+ n = acb_mat_nrows(A);
+
+ acb_mat_init(Y, n, n);
+ acb_mat_init(R1, n, n);
+ acb_mat_init(R2, n, n);
+ mag_init(r1);
+ mag_init(r2);
+
+ /* Y ~= inv(R) */
+ acb_mat_one(R1);
+ acb_mat_approx_solve(Y, R, R1, prec);
+
+ /* R2 = Y*R - I */
+ acb_mat_mul(R2, Y, R, prec);
+ for (i = 0; i < n; i++)
+ acb_sub_ui(acb_mat_entry(R2, i, i), acb_mat_entry(R2, i, i), 1, prec);
+
+ acb_mat_bound_inf_norm(r2, R2);
+
+ if (mag_cmp_2exp_si(r2, 0) < 0)
+ {
+ /* R1 = Y*(AR - RD) */
+ acb_mat_mul(R2, A, R, prec);
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
+ acb_submul(acb_mat_entry(R2, i, j), acb_mat_entry(R, i, j), E + j, prec);
+ acb_mat_mul(R1, Y, R2, prec);
+
+ acb_mat_bound_inf_norm(r1, R1);
+ mag_geom_series(r2, r2, 0);
+ mag_mul(eps, r1, r2);
+ }
+ else
+ {
+ mag_inf(eps);
+ }
+
+ acb_mat_clear(R1);
+ acb_mat_clear(R2);
+ acb_mat_clear(Y);
+ mag_clear(r1);
+ mag_clear(r2);
+}
diff --git a/acb_mat/randtest_eig.c b/acb_mat/randtest_eig.c
new file mode 100644
index 00000000..c87d1ab5
--- /dev/null
+++ b/acb_mat/randtest_eig.c
@@ -0,0 +1,58 @@
+/*
+ Copyright (C) 2018 Fredrik Johansson
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "acb_mat.h"
+
+void
+acb_mat_randtest_eig(acb_mat_t A, flint_rand_t state, acb_srcptr E, slong prec)
+{
+ slong n, i, j, ebits;
+ acb_mat_t U, Q;
+
+ n = acb_mat_nrows(A);
+ ebits = 1 + n_randint(state, 5);
+
+ acb_mat_init(U, n, n);
+ acb_mat_init(Q, n, n);
+
+ /* Skew-Hermitian matrix */
+ acb_mat_randtest(Q, state, prec, 1);
+ if (n_randint(state, 2))
+ acb_mat_get_mid(Q, Q);
+ for (i = 0; i < n; i++)
+ {
+ for (j = i + 1; j < n; j++)
+ {
+ acb_neg(acb_mat_entry(Q, i, j), acb_mat_entry(Q, j, i));
+ acb_conj(acb_mat_entry(Q, i, j), acb_mat_entry(Q, i, j));
+ }
+ arb_zero(acb_realref(acb_mat_entry(Q, i, i)));
+ }
+
+ acb_mat_exp(Q, Q, prec);
+
+ acb_mat_randtest(U, state, prec, ebits);
+ if (n_randint(state, 2))
+ acb_mat_get_mid(U, U);
+ for (i = 0; i < n; i++)
+ for (j = 0; j < i; j++)
+ acb_zero(acb_mat_entry(U, i, j));
+
+ for (i = 0; i < n; i++)
+ acb_set(acb_mat_entry(U, i, i), E + i);
+
+ acb_mat_mul(U, Q, U, prec);
+ acb_mat_conjugate_transpose(Q, Q);
+ acb_mat_mul(A, U, Q, prec);
+
+ acb_mat_clear(U);
+ acb_mat_clear(Q);
+}
diff --git a/acb_mat/test/t-eig_global_enclosure.c b/acb_mat/test/t-eig_global_enclosure.c
new file mode 100644
index 00000000..dfe85007
--- /dev/null
+++ b/acb_mat/test/t-eig_global_enclosure.c
@@ -0,0 +1,107 @@
+/*
+ Copyright (C) 2018 Fredrik Johansson
+
+ This file is part of Arb.
+
+ Arb is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 2.1 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "acb_mat.h"
+
+int main()
+{
+ slong iter;
+ flint_rand_t state;
+
+ flint_printf("eig_global_enclosure....");
+ fflush(stdout);
+
+ flint_randinit(state);
+
+ for (iter = 0; iter < 250 * arb_test_multiplier(); iter++)
+ {
+ acb_mat_t A, X;
+ acb_ptr E, F, B;
+ mag_t eps;
+ slong i, j, n, prec;
+ int success;
+
+ n = n_randint(state, 8);
+ prec = 2 + n_randint(state, 100);
+
+ acb_mat_init(A, n, n);
+ acb_mat_init(X, n, n);
+ E = _acb_vec_init(n);
+ F = _acb_vec_init(n);
+ B = _acb_vec_init(n);
+ mag_init(eps);
+
+ acb_mat_randtest_eig(A, state, E, prec);
+ acb_mat_approx_eig_qr(F, NULL, X, A, NULL, 0, prec);
+
+ /* perturb F further */
+ if (n_randint(state, 2))
+ {
+ for (i = 0; i < n; i++)
+ {
+ acb_randtest(B, state, prec, 1);
+ acb_mul_2exp_si(B, B, -n_randint(state, prec));
+ acb_add(F + i, F + i, B, prec);
+ }
+ }
+
+ acb_mat_eig_global_enclosure(eps, A, F, X, prec);
+
+ _acb_vec_set(B, F, n);
+ for (i = 0; i < n; i++)
+ acb_add_error_mag(B + i, eps);
+
+ for (i = 0; i < n; i++)
+ {
+ success = 0;
+
+ for (j = 0; j < n; j++)
+ {
+ if (acb_contains(B + j, E + i))
+ {
+ success = 1;
+ break;
+ }
+ }
+
+ if (!success)
+ {
+ flint_printf("FAIL\n\n");
+ flint_printf("A =\n");
+ acb_mat_printd(A, 20);
+ flint_printf("\n\ni = %wd\n\n", i);
+ flint_printf("\nE =\n");
+ for (j = 0; j < n; j++)
+ {
+ acb_printn(E + j, 20, 0); flint_printf("\n");
+ }
+ flint_printf("\nB =\n");
+ for (j = 0; j < n; j++)
+ {
+ acb_printn(B + j, 20, 0); flint_printf("\n");
+ }
+ flint_abort();
+ }
+ }
+
+ acb_mat_clear(A);
+ acb_mat_clear(X);
+ _acb_vec_clear(E, n);
+ _acb_vec_clear(F, n);
+ _acb_vec_clear(B, n);
+ mag_clear(eps);
+ }
+
+ flint_randclear(state);
+ flint_cleanup();
+ flint_printf("PASS\n");
+ return EXIT_SUCCESS;
+}
diff --git a/doc/source/acb_mat.rst b/doc/source/acb_mat.rst
index 4108e1a3..9db75a21 100644
--- a/doc/source/acb_mat.rst
+++ b/doc/source/acb_mat.rst
@@ -104,6 +104,13 @@ Random generation
Sets *mat* to a random matrix with up to *prec* bits of precision
and with exponents of width up to *mag_bits*.
+.. function:: void acb_mat_randtest_eig(acb_mat_t mat, flint_rand_t state, acb_srcptr E, slong prec)
+
+ Sets *mat* to a random matrix with the prescribed eigenvalues
+ supplied as the vector *E*. The output matrix is required to be
+ square. We generate a random unitary matrix via a matrix
+ exponential, and then evaluate an inverse Schur decomposition.
+
Input and output
-------------------------------------------------------------------------------
@@ -545,3 +552,46 @@ Eigenvalues and eigenvectors
return value indicates that the QR iteration converged numerically,
but this is only a heuristic termination test and does not imply
any statement whatsoever about error bounds.
+
+.. function:: void acb_mat_eig_global_enclosure(mag_t eps, const acb_mat_t A, acb_srcptr E, const acb_mat_t R, slong prec)
+
+ Given an *n* by *n* matrix *A*, a length-*n* vector *E*
+ containing approximations of the eigenvalues of *A*,
+ and an *n* by *n* matrix *R* containing approximations of
+ the corresponding eigenvectors, computes a rigorous bound
+ `\varepsilon` such that every eigenvalue `\lambda` of *A* satisfies
+ `|\lambda - \hat \lambda_k| \le \varepsilon`
+ for some `\hat \lambda_k` in *E*.
+
+ In other words, the union of the balls
+ `B_k = \{z : |z - \hat \lambda_k| \le \varepsilon\}` is guaranteed to
+ be an enclosure of all eigenvalues of *A*.
+
+ Note that there is no guarantee that each ball `B_k` can be
+ identified with a single eigenvalue: it is possible that some
+ balls contain several eigenvalues while other balls contain
+ no eigenvalues. In other words, this method is not powerful enough
+ to compute isolating balls for the individual eigenvalues (or even
+ for clusters of eigenvalues other than the whole spectrum).
+ Nevertheless, in practice the balls `B_k` will represent
+ eigenvalues one-to-one with high probability if the
+ given approximations are good.
+
+ The output can be used to certify
+ that all eigenvalues of *A* lie in some region of the complex
+ plane (such as a specific half-plane, strip, disk, or annulus)
+ without the need to certify the individual eigenvalues.
+ The output is easily converted into lower or upper bounds
+ for the absolute values or real or complex parts
+ of the spectrum, and with high probability these bounds will be tight.
+ Using :func:`acb_add_error_mag` and :func:`acb_union`, the output
+ can also be converted to a single :type:`acb_t` enclosing
+ the whole spectrum of *A* in a rectangle, but note that to
+ test whether a condition holds for all eigenvalues of *A*, it
+ is typically better to iterate over the individual balls `B_k`.
+
+ This function implements the fast algorithm in Theorem 1 in
+ [Miy2010]_ which extends the Bauer-Fike theorem. Approximations
+ can, for instance, be computed using :func:`acb_mat_approx_eig_qr`.
+ No assumptions are made about the structure of *A* or the
+ quality of the given approximations.
diff --git a/doc/source/credits.rst b/doc/source/credits.rst
index 11de7fa3..52a4d995 100644
--- a/doc/source/credits.rst
+++ b/doc/source/credits.rst
@@ -225,6 +225,8 @@ Bibliography
.. [Kri2013] \A. Krishnamoorthy and D. Menon, "Matrix Inversion Using Cholesky Decomposition" Proc. of the International Conference on Signal Processing Algorithms, Architectures, Arrangements, and Applications (SPA-2013), pp. 70-72, 2013.
+.. [Miy2010] \S. Miyajima, "Fast enclosure for all eigenvalues in generalized eigenvalue problems", Journal of Computational and Applied Mathematics, 233 (2010), 2994-3004, https://dx.doi.org/10.1016/j.cam.2009.11.048
+
.. [MPFR2012] The MPFR team, "MPFR Algorithms" (2012), http://www.mpfr.org/algo.html
.. [NIST2012] National Institute of Standards and Technology, *Digital Library of Mathematical Functions* (2012), http://dlmf.nist.gov/
@@ -246,4 +248,3 @@ Bibliography
.. [Tak2000] \D. Takahashi, "A fast algorithm for computing large Fibonacci numbers", Information Processing Letters 75 (2000) 243-246, http://www.ii.uni.wroc.pl/~lorys/IPL/article75-6-1.pdf
.. [Tre2008] \L. N. Trefethen, "Is Gauss Quadrature Better than Clenshaw-Curtis?", SIAM Review, 50:1 (2008), 67-87, https://doi.org/10.1137/060659831
-