mirror of
https://github.com/vale981/arb
synced 2025-03-04 17:01:40 -05:00
add acb_mat_eig_multiple_rump, handling overlapping eigenvalues
This commit is contained in:
fredrik
parent
eecc160028
commit
bbf6860121
4 changed files with 365 additions and 0 deletions
|
|
@ -423,6 +423,8 @@ int acb_mat_eig_simple_vdhoeven_mourrain(acb_ptr E, acb_mat_t L, acb_mat_t R,
|
|||
int acb_mat_eig_simple(acb_ptr E, acb_mat_t L, acb_mat_t R,
|
||||
const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec);
|
||||
|
||||
int acb_mat_eig_multiple_rump(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec);
|
||||
|
||||
/* Special functions */
|
||||
|
||||
void acb_mat_exp_taylor_sum(acb_mat_t S, const acb_mat_t A, slong N, slong prec);
|
||||
|
|
|
|||
181
acb_mat/eig_multiple_rump.c
Normal file
181
acb_mat/eig_multiple_rump.c
Normal file
|
|
@ -0,0 +1,181 @@
|
|||
/*
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_mat.h"
|
||||
|
||||
static void
|
||||
acb_approx_mag(mag_t res, const acb_t x)
|
||||
{
|
||||
mag_t t;
|
||||
mag_init(t);
|
||||
arf_get_mag(res, arb_midref(acb_realref(x)));
|
||||
arf_get_mag(t, arb_midref(acb_imagref(x)));
|
||||
mag_hypot(res, res, t);
|
||||
mag_clear(t);
|
||||
}
|
||||
|
||||
static int
|
||||
close(const acb_t x, const acb_t y, const mag_t eps)
|
||||
{
|
||||
arf_t t;
|
||||
mag_t a, b;
|
||||
int result;
|
||||
|
||||
mag_init(a);
|
||||
mag_init(b);
|
||||
arf_init(t);
|
||||
arf_sub(t, arb_midref(acb_realref(x)), arb_midref(acb_realref(y)), MAG_BITS, ARF_RND_UP);
|
||||
arf_get_mag(a, t);
|
||||
arf_sub(t, arb_midref(acb_imagref(x)), arb_midref(acb_imagref(y)), MAG_BITS, ARF_RND_UP);
|
||||
arf_get_mag(b, t);
|
||||
|
||||
mag_hypot(a, a, b);
|
||||
result = (mag_cmp(a, eps) <= 0);
|
||||
|
||||
mag_clear(a);
|
||||
mag_clear(b);
|
||||
arf_clear(t);
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
int
|
||||
acb_mat_eig_multiple_rump(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
|
||||
{
|
||||
slong c, i, j, k, n;
|
||||
acb_mat_t X;
|
||||
acb_ptr F;
|
||||
int result;
|
||||
slong iter;
|
||||
mag_t escale, eps, tm, um;
|
||||
slong ** cluster;
|
||||
slong * cluster_size;
|
||||
slong num_clusters;
|
||||
|
||||
n = acb_mat_nrows(A);
|
||||
|
||||
if (n == 0)
|
||||
return 1;
|
||||
|
||||
cluster = flint_malloc(sizeof(slong *) * n);
|
||||
for (i = 0; i < n; i++)
|
||||
cluster[i] = flint_malloc(sizeof(slong) * n);
|
||||
cluster_size = flint_malloc(sizeof(slong) * n);
|
||||
|
||||
mag_init(eps);
|
||||
mag_init(escale);
|
||||
mag_init(tm);
|
||||
mag_init(um);
|
||||
|
||||
mag_zero(escale);
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
acb_approx_mag(tm, E_approx + i);
|
||||
mag_max(escale, escale, tm);
|
||||
}
|
||||
|
||||
/* todo: when num_clusters = 1, could use fallback global enclosure */
|
||||
|
||||
/* todo: initial clustering could be allowed to be zero */
|
||||
/* M 2^(-0.75p), M 2^(-0.5p), M 2^(-0.25p), ... */
|
||||
for (iter = 0; iter < 2; iter++)
|
||||
{
|
||||
mag_mul_2exp_si(eps, escale, -prec + (iter + 1) * prec/4);
|
||||
|
||||
/* Group the eigenvalue approximations. */
|
||||
num_clusters = 0;
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
int new_cluster = 1;
|
||||
|
||||
for (j = 0; j < num_clusters && new_cluster; j++)
|
||||
{
|
||||
if (close(E_approx + i, E_approx + cluster[j][0], eps))
|
||||
{
|
||||
cluster[j][cluster_size[j]] = i;
|
||||
cluster_size[j]++;
|
||||
new_cluster = 0;
|
||||
}
|
||||
}
|
||||
|
||||
if (new_cluster)
|
||||
{
|
||||
cluster[num_clusters][0] = i;
|
||||
cluster_size[num_clusters] = 1;
|
||||
num_clusters++;
|
||||
}
|
||||
}
|
||||
|
||||
result = 1;
|
||||
|
||||
F = _acb_vec_init(num_clusters);
|
||||
|
||||
for (c = 0; c < num_clusters && result; c++)
|
||||
{
|
||||
k = cluster_size[c];
|
||||
|
||||
acb_mat_init(X, n, k);
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
for (j = 0; j < k; j++)
|
||||
acb_set(acb_mat_entry(X, i, j), acb_mat_entry(R_approx, i, cluster[c][j]));
|
||||
|
||||
acb_mat_eig_enclosure_rump(F + c, NULL, X, A, E_approx + cluster[c][0], X, prec);
|
||||
|
||||
if (!acb_is_finite(F + c))
|
||||
result = 0;
|
||||
|
||||
acb_mat_clear(X);
|
||||
}
|
||||
|
||||
for (i = 0; i < num_clusters; i++)
|
||||
{
|
||||
for (j = i + 1; j < num_clusters; j++)
|
||||
{
|
||||
if (acb_overlaps(F + i, F + j))
|
||||
result = 0;
|
||||
}
|
||||
}
|
||||
|
||||
if (result)
|
||||
{
|
||||
i = 0;
|
||||
for (c = 0; c < num_clusters; c++)
|
||||
{
|
||||
for (j = 0; j < cluster_size[c]; j++)
|
||||
{
|
||||
acb_set(E + i, F + c);
|
||||
i++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
_acb_vec_clear(F, num_clusters);
|
||||
|
||||
if (result)
|
||||
break;
|
||||
}
|
||||
|
||||
if (!result)
|
||||
_acb_vec_indeterminate(E, n);
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
flint_free(cluster[i]);
|
||||
flint_free(cluster);
|
||||
flint_free(cluster_size);
|
||||
|
||||
mag_clear(eps);
|
||||
mag_clear(escale);
|
||||
mag_clear(tm);
|
||||
mag_clear(um);
|
||||
|
||||
return result;
|
||||
}
|
||||
153
acb_mat/test/t-eig_multiple_rump.c
Normal file
153
acb_mat/test/t-eig_multiple_rump.c
Normal file
|
|
@ -0,0 +1,153 @@
|
|||
/*
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "acb_mat.h"
|
||||
|
||||
static void
|
||||
acb_get_mid(acb_t res, const acb_t x)
|
||||
{
|
||||
arb_get_mid_arb(acb_realref(res), acb_realref(x));
|
||||
arb_get_mid_arb(acb_imagref(res), acb_imagref(x));
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
flint_rand_t state;
|
||||
|
||||
flint_printf("eig_multiple_rump....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
acb_mat_t A, R;
|
||||
acb_ptr E, F;
|
||||
acb_t b;
|
||||
slong i, j, n, prec, count, count2;
|
||||
int result;
|
||||
|
||||
n = n_randint(state, 8);
|
||||
prec = 2 + n_randint(state, 200);
|
||||
|
||||
acb_init(b);
|
||||
acb_mat_init(A, n, n);
|
||||
acb_mat_init(R, n, n);
|
||||
E = _acb_vec_init(n);
|
||||
F = _acb_vec_init(n);
|
||||
|
||||
if (n_randint(state, 10) != 0)
|
||||
{
|
||||
for (i = 0; i < n; i++)
|
||||
acb_randtest(E + i, state, prec, 2);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* Randomly repeat eigenvalues. */
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
if (i == 0 || n_randint(state, 2))
|
||||
acb_randtest(E + i, state, prec, 2);
|
||||
else
|
||||
acb_set(E + i, E + n_randint(state, i));
|
||||
}
|
||||
}
|
||||
|
||||
if (n_randint(state, 2))
|
||||
{
|
||||
for (i = 0; i < n; i++)
|
||||
acb_get_mid(E + i, E + i);
|
||||
}
|
||||
|
||||
acb_mat_randtest_eig(A, state, E, prec);
|
||||
acb_mat_approx_eig_qr(F, NULL, R, A, NULL, 0, prec);
|
||||
|
||||
/* Perturb F further. */
|
||||
if (n_randint(state, 10) == 0)
|
||||
{
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
/* Perturb R further. */
|
||||
if (n_randint(state, 10) == 0)
|
||||
{
|
||||
j = n_randint(state, n);
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
acb_randtest(b, state, prec, 1);
|
||||
acb_mul_2exp_si(b, b, -10 - n_randint(state, prec));
|
||||
acb_add(acb_mat_entry(R, i, j), acb_mat_entry(R, i, j), b, prec);
|
||||
}
|
||||
}
|
||||
|
||||
result = acb_mat_eig_multiple_rump(F, A, E, R, prec);
|
||||
|
||||
if (result)
|
||||
{
|
||||
count = 0;
|
||||
count2 = 0;
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
for (j = 0; j < n; j++)
|
||||
{
|
||||
if (j == 0 || !acb_equal(F + j, F + j - 1))
|
||||
count += acb_contains(F + j, E + i);
|
||||
}
|
||||
|
||||
for (j = 0; j < n; j++)
|
||||
{
|
||||
if (j == 0 || !acb_equal(F + j, F + j - 1))
|
||||
count2 += acb_overlaps(F + j, E + i);
|
||||
}
|
||||
}
|
||||
|
||||
if (count != n || count2 != n)
|
||||
{
|
||||
flint_printf("FAIL: count\n\n");
|
||||
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("count = %wd, count2 = %wd\n\n", count, count2);
|
||||
flint_printf("E = \n");
|
||||
for (j = 0; j < n; j++)
|
||||
{
|
||||
acb_printd(E + j, 20);
|
||||
flint_printf("\n");
|
||||
}
|
||||
flint_printf("F = \n");
|
||||
for (j = 0; j < n; j++)
|
||||
{
|
||||
acb_printd(F + j, 20);
|
||||
flint_printf("\n");
|
||||
}
|
||||
flint_abort();
|
||||
}
|
||||
}
|
||||
|
||||
acb_mat_clear(A);
|
||||
acb_mat_clear(R);
|
||||
_acb_vec_clear(E, n);
|
||||
_acb_vec_clear(F, n);
|
||||
acb_clear(b);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -692,3 +692,32 @@ than necessary. Manually balancing badly scaled matrices may help.
|
|||
complexity `O(n^3)`.
|
||||
|
||||
The default version currently uses *vdhoeven_mourrain*.
|
||||
|
||||
By design, these functions terminate instead of attempting to
|
||||
compute eigenvalue clusters if some eigenvalues cannot be isolated.
|
||||
To compute all eigenvalues of a matrix allowing for overlap,
|
||||
:func:`acb_mat_eig_multiple_rump` may be used as a fallback.
|
||||
|
||||
.. function:: int acb_mat_eig_multiple_rump(acb_ptr E, const acb_mat_t A, acb_srcptr E_approx, const acb_mat_t R_approx, slong prec)
|
||||
|
||||
Computes all the eigenvalues of the given *n* by *n* matrix *A*.
|
||||
On success, the output vector *E* contains *n* complex intervals,
|
||||
each representing one eigenvalue of *A* with the correct
|
||||
multiplicities in case of overlap.
|
||||
The output intervals are either disjoint or identical, and
|
||||
identical intervals are guaranteed to be grouped consecutively.
|
||||
Each complete run of *k* identical intervals thus represents a cluster of
|
||||
exactly *k* eigenvalues which could not be separated from each
|
||||
other at the current precision, but which could be isolated
|
||||
from the other `n - k` eigenvalues of the matrix.
|
||||
|
||||
The user supplies approximations *E_approx* and *R_approx*
|
||||
of the eigenvalues and the right eigenvectors.
|
||||
The initial 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.
|
||||
|
||||
This function groups approximate eigenvalues that are close and
|
||||
calls :func:`acb_mat_eig_enclosure_rump` repeatedly to validate
|
||||
each cluster. The complexity is `O(m n^3)` for *m* clusters.
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue