diff --git a/acb_mat/eig_enclosure_rump.c b/acb_mat/eig_enclosure_rump.c index 13f326e5..ee95999b 100644 --- a/acb_mat/eig_enclosure_rump.c +++ b/acb_mat/eig_enclosure_rump.c @@ -148,7 +148,7 @@ acb_approx_mag(mag_t res, const acb_t x) /* Extract k largest rows to freeze */ static void -partition_X(slong * u, slong * v, const acb_mat_t X) +partition_X_sorted(slong * u, slong * v, const acb_mat_t X, slong prec) { slong i, j, n, k, c; slong * row_idx; @@ -181,7 +181,9 @@ partition_X(slong * u, slong * v, const acb_mat_t X) if (mag_cmp(row_mag + j, row_mag + j + 1) > 0) { mag_swap(row_mag + j, row_mag + j + 1); - c = row_idx[j]; row_idx[j] = row_idx[j + 1]; row_idx[j + 1] = c; + c = row_idx[j]; + row_idx[j] = row_idx[j + 1]; + row_idx[j + 1] = c; } } } @@ -189,6 +191,7 @@ partition_X(slong * u, slong * v, const acb_mat_t X) /* Not frozen rows of the approximation. */ for (i = 0; i < n - k; i++) u[i] = row_idx[i]; + /* Frozen rows of the approximation. */ for (i = 0; i < k; i++) v[i] = row_idx[n - k + i]; @@ -198,6 +201,22 @@ partition_X(slong * u, slong * v, const acb_mat_t X) mag_clear(t); } +static void +partition_X_trivial(slong * u, slong * v, const acb_mat_t X, slong prec) +{ + slong n, k, i; + + n = acb_mat_nrows(X); + k = acb_mat_ncols(X); + + /* Not frozen rows of the approximation. */ + for (i = 0; i < n - k; i++) + u[i] = i; + /* Frozen rows of the approximation. */ + for (i = 0; i < k; i++) + v[i] = n - k + i; +} + void acb_mat_eig_enclosure_rump(acb_t lambda, acb_mat_t J, acb_mat_t X, const acb_mat_t A, const acb_t lambda_approx, const acb_mat_t X_approx, slong prec) @@ -221,7 +240,10 @@ acb_mat_eig_enclosure_rump(acb_t lambda, acb_mat_t J, acb_mat_t X, const acb_mat /* Frozen rows of the approximation. */ v = flint_malloc(sizeof(slong) * k); - partition_X(u, v, X_approx); + if (k == 1) + partition_X_sorted(u, v, X_approx, prec); + else + partition_X_trivial(u, v, X_approx, prec); mag_init(eps); acb_mat_init(R, n, n); diff --git a/acb_mat/test/t-eig_multiple.c b/acb_mat/test/t-eig_multiple.c index c4345b0d..7eba9f72 100644 --- a/acb_mat/test/t-eig_multiple.c +++ b/acb_mat/test/t-eig_multiple.c @@ -142,6 +142,95 @@ int main() acb_clear(b); } + /* Test convergence for DFT matrices */ + for (iter = 0; iter < 50 * arb_test_multiplier(); iter++) + { + acb_mat_t A, R, QC; + acb_ptr E; + acb_t t; + fmpq_mat_t Q, Qinv; + slong i, n, c0, c1, c2, c3; + slong prec; + int algorithm, result; + + n = n_randint(state, 30); + algorithm = n_randint(state, 2); + acb_mat_init(A, n, n); + acb_mat_init(R, n, n); + E = _acb_vec_init(n); + acb_init(t); + acb_mat_init(QC, n, n); + fmpq_mat_init(Q, n, n); + fmpq_mat_init(Qinv, n, n); + + /* The current algorithm is not robust enough. */ +#if 0 + do { + fmpq_mat_randtest(Q, state, 2 + n_randint(state, 10)); + } while (!fmpq_mat_inv(Qinv, Q)); +#else + fmpq_mat_one(Q); + fmpq_mat_one(Qinv); +#endif + + for (prec = 32; ; prec *= 2) + { + if (prec > 10000) + { + flint_printf("FAIL: unsuccessful, prec > 10000\n\n"); + flint_printf("algorithm = %d, iter %wd\n\n", algorithm, iter); + flint_abort(); + } + + acb_mat_dft(A, 0, prec); + +#if 0 + acb_mat_set_fmpq_mat(QC, Q, prec); + acb_mat_mul(A, A, QC, prec); + acb_mat_set_fmpq_mat(QC, Qinv, prec); + acb_mat_mul(A, QC, A, prec); +#endif + + acb_mat_approx_eig_qr(E, NULL, R, A, NULL, 0, prec); + + if (algorithm == 0) + result = acb_mat_eig_multiple_rump(E, A, E, R, prec); + else + result = acb_mat_eig_multiple(E, A, E, R, prec); + + /* Verify the known eigenvalues + multiplicities */ + if (result) + { + c0 = c1 = c2 = c3 = 0; + + for (i = 0; i < n; i++) + { + acb_set_d_d(t, 1.0, 0.0); + c0 += acb_contains(E + i, t); + acb_set_d_d(t, -1.0, 0.0); + c1 += acb_contains(E + i, t); + acb_set_d_d(t, 0.0, 1.0); + c2 += acb_contains(E + i, t); + acb_set_d_d(t, 0.0, -1.0); + c3 += acb_contains(E + i, t); + } + + result = (n == 0 || (c0 == (n+4)/4 && c1 == (n+2)/4 && c2 == (n-1)/4 && c3 == (n+1)/4)); + } + + if (result) + break; + } + + acb_mat_clear(A); + acb_mat_clear(R); + acb_mat_clear(QC); + _acb_vec_clear(E, n); + acb_clear(t); + fmpq_mat_clear(Q); + fmpq_mat_clear(Qinv); + } + flint_randclear(state); flint_cleanup(); flint_printf("PASS\n"); diff --git a/acb_mat/test/t-eig_simple.c b/acb_mat/test/t-eig_simple.c index 095f0d29..94168627 100644 --- a/acb_mat/test/t-eig_simple.c +++ b/acb_mat/test/t-eig_simple.c @@ -156,9 +156,10 @@ int main() /* Test convergence, given companion matrices */ for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++) { - acb_mat_t A, R; + acb_mat_t A, R, QC; acb_ptr E; acb_ptr roots; + fmpq_mat_t Q, Qinv; acb_poly_t f; slong i, j, n, prec, count, count2; int algorithm, success; @@ -170,6 +171,9 @@ int main() acb_poly_init(f); acb_mat_init(A, n, n); acb_mat_init(R, n, n); + fmpq_mat_init(Q, n, n); + fmpq_mat_init(Qinv, n, n); + acb_mat_init(QC, n, n); for (i = 0; i < n; i++) { @@ -182,6 +186,10 @@ int main() goto new_root; } + do { + fmpq_mat_randtest(Q, state, 2 + n_randint(state, 100)); + } while (!fmpq_mat_inv(Qinv, Q)); + success = 0; for (prec = 32; !success; prec *= 2) @@ -189,7 +197,7 @@ int main() if (prec > 10000) { flint_printf("FAIL: unsuccessful, prec > 10000\n\n"); - flint_printf("algorithm = %d\n\n", algorithm); + flint_printf("algorithm = %d, iter %wd\n\n", algorithm, iter); flint_printf("A = \n"); acb_mat_printd(A, 20); flint_printf("\n\n"); flint_printf("R = \n"); acb_mat_printd(R, 20); flint_printf("\n\n"); flint_printf("roots = \n"); @@ -202,7 +210,12 @@ int main() } acb_poly_product_roots(f, roots, n, prec); + acb_mat_companion(A, f, prec); + acb_mat_set_fmpq_mat(QC, Q, prec); + acb_mat_mul(A, A, QC, prec); + acb_mat_set_fmpq_mat(QC, Qinv, prec); + acb_mat_mul(A, QC, A, prec); acb_mat_approx_eig_qr(E, NULL, R, A, NULL, 0, prec); @@ -250,6 +263,9 @@ int main() } } + fmpq_mat_clear(Q); + fmpq_mat_clear(Qinv); + acb_mat_clear(QC); acb_mat_clear(A); acb_mat_clear(R); acb_poly_clear(f); diff --git a/doc/source/acb_mat.rst b/doc/source/acb_mat.rst index 4ef2e654..4f584ae1 100644 --- a/doc/source/acb_mat.rst +++ b/doc/source/acb_mat.rst @@ -581,10 +581,17 @@ Component and error operations Eigenvalues and eigenvectors ------------------------------------------------------------------------------- -The functions in this section are experimental. There may be classes +The functions in this section are experimental. There are classes of matrices where the algorithms fail to converge even as *prec* is increased, or for which the error bounds are much worse -than necessary. Manually balancing badly scaled matrices may help. +than necessary. In some cases, it can help to manually precondition +the matrix *A* by applying a similarity transformation `T^{-1} A T`. + +* If *A* is badly scaled, take `T` to be a matrix such that the entries + of `T^{-1} A T` are more uniform (this is known as balancing). +* Simply taking `T` to be a random invertible matrix can help if an + algorithm fails to converge despite `A` being well-scaled. (This + can be the case when dealing with multiple eigenvalues.) .. function:: 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)