mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
approximate eigenvectors + better test code
This commit is contained in:
fredrik
parent
5535cba2b2
commit
0e2867460c
4 changed files with 490 additions and 127 deletions
|
|
@ -21,6 +21,17 @@ Todo items present in the original code:
|
|||
|
||||
#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 void
|
||||
acb_get_mid(acb_t res, const acb_t x)
|
||||
{
|
||||
|
|
@ -64,6 +75,31 @@ acb_approx_div_arb(acb_t res, const acb_t x, const arb_t y, slong prec)
|
|||
arf_div(arb_midref(acb_imagref(res)), arb_midref(acb_imagref(x)), arb_midref(y), prec, ARF_RND_DOWN);
|
||||
}
|
||||
|
||||
static void
|
||||
acb_approx_inv(acb_t z, const acb_t x, slong prec)
|
||||
{
|
||||
arf_set(arb_midref(acb_realref(z)), arb_midref(acb_realref(x)));
|
||||
arf_set(arb_midref(acb_imagref(z)), arb_midref(acb_imagref(x)));
|
||||
|
||||
mag_zero(arb_radref(acb_realref(z)));
|
||||
mag_zero(arb_radref(acb_imagref(z)));
|
||||
|
||||
acb_inv(z, z, prec);
|
||||
|
||||
mag_zero(arb_radref(acb_realref(z)));
|
||||
mag_zero(arb_radref(acb_imagref(z)));
|
||||
}
|
||||
|
||||
static void
|
||||
acb_approx_div(acb_t z, const acb_t x, const acb_t y, slong prec)
|
||||
{
|
||||
acb_t t;
|
||||
acb_init(t);
|
||||
acb_approx_inv(t, y, prec);
|
||||
acb_approx_mul(z, x, t, prec);
|
||||
acb_clear(t);
|
||||
}
|
||||
|
||||
void
|
||||
acb_mat_approx_qr_step(acb_mat_t A, acb_mat_t Q, slong n0, slong n1, const acb_t shift, slong prec)
|
||||
{
|
||||
|
|
@ -316,7 +352,7 @@ acb_mat_approx_hessenberg_reduce_0(acb_mat_t A, acb_ptr T, slong prec)
|
|||
|
||||
arf_ui_div(scale_inv, 1, scale, prec, ARF_RND_DOWN);
|
||||
|
||||
if (arf_is_zero(scale_inv))
|
||||
if (arf_is_zero(scale))
|
||||
{
|
||||
acb_zero(T + i);
|
||||
acb_zero(acb_mat_entry(A, i, i - 1));
|
||||
|
|
@ -515,6 +551,207 @@ acb_mat_approx_hessenberg_reduce_1(acb_mat_t A, acb_srcptr T, slong prec)
|
|||
acb_clear(t);
|
||||
}
|
||||
|
||||
/* Right eigenvectors of a triu matrix. No aliasing. */
|
||||
void
|
||||
acb_mat_approx_eig_triu_r(acb_mat_t ER, const acb_mat_t A, slong prec)
|
||||
{
|
||||
slong i, j, k, n;
|
||||
mag_t tm, smin, unfl, simin, smlnum, rmax;
|
||||
acb_t r, s, t;
|
||||
|
||||
n = acb_mat_nrows(A);
|
||||
|
||||
acb_mat_one(ER);
|
||||
|
||||
acb_init(r);
|
||||
acb_init(s);
|
||||
acb_init(t);
|
||||
mag_init(tm);
|
||||
mag_init(smin);
|
||||
mag_init(smlnum);
|
||||
mag_init(unfl);
|
||||
mag_init(simin);
|
||||
mag_init(rmax);
|
||||
|
||||
mag_set_ui_2exp_si(unfl, 1, -30 * prec);
|
||||
mag_mul_ui(smlnum, unfl, n);
|
||||
mag_mul_2exp_si(smlnum, smlnum, prec);
|
||||
mag_set_ui_2exp_si(simin, 1, prec / 2);
|
||||
mag_one(rmax);
|
||||
|
||||
for (i = 1; i < n; i++)
|
||||
{
|
||||
acb_set(s, acb_mat_entry(A, i, i));
|
||||
|
||||
/* smin = max(eps * abs(s), smlnum) */
|
||||
acb_approx_mag(smin, s);
|
||||
mag_mul_2exp_si(smin, smin, -prec);
|
||||
mag_max(smin, smin, smlnum);
|
||||
|
||||
for (j = i - 1; j >= 0; j--)
|
||||
{
|
||||
acb_approx_dot(r, NULL, 0, A->rows[j] + j + 1, 1, ER->rows[i] + j + 1, 1, i - j, prec);
|
||||
acb_approx_sub(t, acb_mat_entry(A, j, j), s, prec);
|
||||
|
||||
/* if abs(t) < smin: t = smin */
|
||||
acb_approx_mag(tm, t);
|
||||
if (mag_cmp(tm, smin) < 0)
|
||||
{
|
||||
acb_zero(t);
|
||||
arf_set_mag(arb_midref(acb_realref(t)), smin);
|
||||
}
|
||||
|
||||
acb_approx_div(acb_mat_entry(ER, i, j), r, t, prec);
|
||||
acb_neg(acb_mat_entry(ER, i, j), acb_mat_entry(ER, i, j));
|
||||
|
||||
acb_approx_mag(tm, r);
|
||||
mag_max(rmax, rmax, tm);
|
||||
if (mag_cmp(rmax, simin) > 0)
|
||||
{
|
||||
arb_t b;
|
||||
arb_init(b);
|
||||
arf_set_mag(arb_midref(b), rmax);
|
||||
|
||||
for (k = j; k < i + 1; k++)
|
||||
{
|
||||
acb_approx_div_arb(acb_mat_entry(ER, i, k),
|
||||
acb_mat_entry(ER, i, k), b, prec);
|
||||
}
|
||||
|
||||
mag_one(rmax);
|
||||
}
|
||||
}
|
||||
|
||||
if (mag_cmp_2exp_si(rmax, 0) != 0)
|
||||
{
|
||||
arb_t b;
|
||||
arb_init(b);
|
||||
arf_set_mag(arb_midref(b), rmax);
|
||||
|
||||
for (k = 0; k < i + 1; k++)
|
||||
{
|
||||
acb_approx_div_arb(acb_mat_entry(ER, i, k),
|
||||
acb_mat_entry(ER, i, k), b, prec);
|
||||
}
|
||||
|
||||
arb_clear(b);
|
||||
}
|
||||
}
|
||||
|
||||
acb_mat_transpose(ER, ER);
|
||||
|
||||
acb_clear(r);
|
||||
acb_clear(s);
|
||||
acb_clear(t);
|
||||
mag_clear(tm);
|
||||
mag_clear(smin);
|
||||
mag_clear(smlnum);
|
||||
mag_clear(unfl);
|
||||
mag_clear(simin);
|
||||
mag_clear(rmax);
|
||||
}
|
||||
|
||||
/* Left eigenvectors of a triu matrix. No aliasing. */
|
||||
void
|
||||
acb_mat_approx_eig_triu_l(acb_mat_t EL, const acb_mat_t A, slong prec)
|
||||
{
|
||||
slong i, j, k, n;
|
||||
mag_t tm, smin, unfl, simin, smlnum, rmax;
|
||||
acb_t r, s, t;
|
||||
acb_mat_t AT;
|
||||
|
||||
n = acb_mat_nrows(A);
|
||||
acb_mat_init(AT, n, n);
|
||||
|
||||
acb_mat_one(EL);
|
||||
acb_mat_transpose(AT, A);
|
||||
|
||||
acb_init(r);
|
||||
acb_init(s);
|
||||
acb_init(t);
|
||||
mag_init(tm);
|
||||
mag_init(smin);
|
||||
mag_init(smlnum);
|
||||
mag_init(unfl);
|
||||
mag_init(simin);
|
||||
mag_init(rmax);
|
||||
|
||||
mag_set_ui_2exp_si(unfl, 1, -30 * prec);
|
||||
mag_mul_ui(smlnum, unfl, n);
|
||||
mag_mul_2exp_si(smlnum, smlnum, prec);
|
||||
mag_set_ui_2exp_si(simin, 1, prec / 2);
|
||||
mag_one(rmax);
|
||||
|
||||
for (i = 0; i < n - 1; i++)
|
||||
{
|
||||
acb_set(s, acb_mat_entry(AT, i, i));
|
||||
|
||||
/* smin = max(eps * abs(s), smlnum) */
|
||||
acb_approx_mag(smin, s);
|
||||
mag_mul_2exp_si(smin, smin, -prec);
|
||||
mag_max(smin, smin, smlnum);
|
||||
|
||||
for (j = i + 1; j < n; j++)
|
||||
{
|
||||
acb_approx_dot(r, NULL, 0, EL->rows[i] + i, 1, AT->rows[j] + i, 1, j - i, prec);
|
||||
acb_approx_sub(t, acb_mat_entry(AT, j, j), s, prec);
|
||||
|
||||
/* if abs(t) < smin: t = smin */
|
||||
acb_approx_mag(tm, t);
|
||||
if (mag_cmp(tm, smin) < 0)
|
||||
{
|
||||
acb_zero(t);
|
||||
arf_set_mag(arb_midref(acb_realref(t)), smin);
|
||||
}
|
||||
|
||||
acb_approx_div(acb_mat_entry(EL, i, j), r, t, prec);
|
||||
acb_neg(acb_mat_entry(EL, i, j), acb_mat_entry(EL, i, j));
|
||||
|
||||
acb_approx_mag(tm, r);
|
||||
mag_max(rmax, rmax, tm);
|
||||
if (mag_cmp(rmax, simin) > 0)
|
||||
{
|
||||
arb_t b;
|
||||
arb_init(b);
|
||||
arf_set_mag(arb_midref(b), rmax);
|
||||
|
||||
for (k = i; k < j + 1; k++)
|
||||
{
|
||||
acb_approx_div_arb(acb_mat_entry(EL, i, k),
|
||||
acb_mat_entry(EL, i, k), b, prec);
|
||||
}
|
||||
|
||||
mag_one(rmax);
|
||||
}
|
||||
}
|
||||
|
||||
if (mag_cmp_2exp_si(rmax, 0) != 0)
|
||||
{
|
||||
arb_t b;
|
||||
arb_init(b);
|
||||
arf_set_mag(arb_midref(b), rmax);
|
||||
|
||||
for (k = i; k < n; k++)
|
||||
{
|
||||
acb_approx_div_arb(acb_mat_entry(EL, i, k),
|
||||
acb_mat_entry(EL, i, k), b, prec);
|
||||
}
|
||||
|
||||
arb_clear(b);
|
||||
}
|
||||
}
|
||||
|
||||
acb_clear(r);
|
||||
acb_clear(s);
|
||||
acb_clear(t);
|
||||
mag_clear(tm);
|
||||
mag_clear(smin);
|
||||
mag_clear(smlnum);
|
||||
mag_clear(unfl);
|
||||
mag_clear(simin);
|
||||
mag_clear(rmax);
|
||||
}
|
||||
|
||||
int
|
||||
acb_mat_approx_hessenberg_qr(acb_mat_t A, acb_mat_t Q, const mag_t tol, slong maxiter, slong prec)
|
||||
{
|
||||
|
|
@ -727,15 +964,9 @@ acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, co
|
|||
{
|
||||
slong n, i, j;
|
||||
acb_ptr T;
|
||||
acb_mat_t Acopy;
|
||||
acb_mat_t Acopy, Q;
|
||||
int result;
|
||||
|
||||
if (L != NULL || R != NULL)
|
||||
{
|
||||
flint_printf("computation of eigenvectors not implemented\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
n = acb_mat_nrows(A);
|
||||
|
||||
T = _acb_vec_init(n);
|
||||
|
|
@ -744,15 +975,45 @@ acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, co
|
|||
|
||||
acb_mat_approx_hessenberg_reduce_0(Acopy, T, prec);
|
||||
|
||||
if (L != NULL || R != NULL)
|
||||
{
|
||||
acb_mat_init(Q, n, n);
|
||||
acb_mat_set(Q, Acopy);
|
||||
acb_mat_approx_hessenberg_reduce_1(Q, T, prec);
|
||||
}
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
for (j = i + 2; j < n; j++)
|
||||
acb_zero(acb_mat_entry(Acopy, j, i));
|
||||
|
||||
result = acb_mat_approx_hessenberg_qr(Acopy, NULL, tol, maxiter, prec);
|
||||
result = acb_mat_approx_hessenberg_qr(Acopy,
|
||||
(L != NULL || R != NULL) ? Q : NULL, tol, maxiter, prec);
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
acb_set(E + i, acb_mat_entry(Acopy, i, i));
|
||||
|
||||
if (R != NULL)
|
||||
{
|
||||
acb_mat_t ER;
|
||||
acb_mat_init(ER, n, n);
|
||||
acb_mat_approx_eig_triu_r(ER, Acopy, prec);
|
||||
acb_mat_approx_mul(R, Q, ER, prec);
|
||||
acb_mat_clear(ER);
|
||||
}
|
||||
|
||||
if (L != NULL)
|
||||
{
|
||||
acb_mat_t EL;
|
||||
acb_mat_init(EL, n, n);
|
||||
acb_mat_approx_eig_triu_l(EL, Acopy, prec);
|
||||
acb_mat_conjugate_transpose(Q, Q);
|
||||
acb_mat_approx_mul(L, EL, Q, prec);
|
||||
acb_mat_clear(EL);
|
||||
}
|
||||
|
||||
if (L != NULL || R != NULL)
|
||||
acb_mat_clear(Q);
|
||||
|
||||
_acb_vec_clear(T, n);
|
||||
acb_mat_clear(Acopy);
|
||||
|
||||
|
|
|
|||
|
|
@ -21,70 +21,167 @@ int main()
|
|||
|
||||
flint_randinit(state);
|
||||
|
||||
/* Test DFT matrices */
|
||||
for (iter = 0; iter < 50 * arb_test_multiplier(); iter++)
|
||||
/* Test random & DFT matrices */
|
||||
for (iter = 0; iter < 200 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
acb_mat_t A;
|
||||
acb_mat_t A, L, R;
|
||||
acb_ptr E;
|
||||
acb_t t;
|
||||
mag_t b;
|
||||
slong i, n, prec, goal, c0, c1, c2, c3;
|
||||
slong i, j, n, prec, goal, c0, c1, c2, c3;
|
||||
int wantL, wantR, result, dft;
|
||||
|
||||
n = n_randint(state, 30);
|
||||
dft = n_randint(state, 2);
|
||||
if (dft)
|
||||
n = n_randint(state, 30);
|
||||
else
|
||||
n = n_randint(state, 15);
|
||||
goal = 2 + n_randint(state, 100);
|
||||
wantL = n_randint(state, 2);
|
||||
wantR = n_randint(state, 2);
|
||||
|
||||
acb_mat_init(A, n, n);
|
||||
acb_mat_init(L, n, n);
|
||||
acb_mat_init(R, n, n);
|
||||
acb_init(t);
|
||||
mag_init(b);
|
||||
E = _acb_vec_init(n);
|
||||
|
||||
for (prec = 32; ; prec *= 2)
|
||||
{
|
||||
acb_mat_dft(A, 0, prec);
|
||||
acb_mat_approx_eig_qr(E, NULL, NULL, A, NULL, 0, prec);
|
||||
|
||||
c0 = c1 = c2 = c3 = 0;
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
if (dft)
|
||||
{
|
||||
acb_set_d_d(t, 1.0, 0.0);
|
||||
acb_sub(t, t, E + i, prec);
|
||||
acb_get_mag(b, t);
|
||||
c0 += (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
|
||||
acb_set_d_d(t, -1.0, 0.0);
|
||||
acb_sub(t, t, E + i, prec);
|
||||
acb_get_mag(b, t);
|
||||
c1 += (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
|
||||
acb_set_d_d(t, 0.0, 1.0);
|
||||
acb_sub(t, t, E + i, prec);
|
||||
acb_get_mag(b, t);
|
||||
c2 += (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
|
||||
acb_set_d_d(t, 0.0, -1.0);
|
||||
acb_sub(t, t, E + i, prec);
|
||||
acb_get_mag(b, t);
|
||||
c3 += (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
acb_mat_dft(A, 0, prec);
|
||||
}
|
||||
else
|
||||
{
|
||||
acb_mat_randtest(A, state, 2 + n_randint(state, 200), 5);
|
||||
acb_mat_get_mid(A, A);
|
||||
}
|
||||
|
||||
if (n == 0 || (c0 == (n+4)/4 && c1 == (n+2)/4 && c2 == (n-1)/4 && c3 == (n+1)/4))
|
||||
acb_mat_approx_eig_qr(E, wantL ? L : NULL, wantR ? R : NULL, A, NULL, 0, prec);
|
||||
|
||||
if (dft)
|
||||
{
|
||||
/* Verify the known eigenvalues + multiplicities */
|
||||
c0 = c1 = c2 = c3 = 0;
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
acb_set_d_d(t, 1.0, 0.0);
|
||||
acb_sub(t, t, E + i, prec);
|
||||
acb_get_mag(b, t);
|
||||
c0 += (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
|
||||
acb_set_d_d(t, -1.0, 0.0);
|
||||
acb_sub(t, t, E + i, prec);
|
||||
acb_get_mag(b, t);
|
||||
c1 += (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
|
||||
acb_set_d_d(t, 0.0, 1.0);
|
||||
acb_sub(t, t, E + i, prec);
|
||||
acb_get_mag(b, t);
|
||||
c2 += (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
|
||||
acb_set_d_d(t, 0.0, -1.0);
|
||||
acb_sub(t, t, E + i, prec);
|
||||
acb_get_mag(b, t);
|
||||
c3 += (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
}
|
||||
|
||||
result = (n == 0 || (c0 == (n+4)/4 && c1 == (n+2)/4 && c2 == (n-1)/4 && c3 == (n+1)/4));
|
||||
}
|
||||
else
|
||||
{
|
||||
result = 1;
|
||||
}
|
||||
|
||||
if (result && wantL)
|
||||
{
|
||||
acb_mat_t LA, D;
|
||||
acb_mat_init(LA, n, n);
|
||||
acb_mat_init(D, n, n);
|
||||
|
||||
/* Check LA - lambda L = 0 */
|
||||
acb_mat_approx_mul(LA, L, A, prec);
|
||||
for (i = 0; i < n; i++)
|
||||
acb_set(acb_mat_entry(D, i, i), E + i);
|
||||
acb_mat_approx_mul(D, D, L, prec);
|
||||
acb_mat_sub(LA, LA, D, prec);
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
for (j = 0; j < n; j++)
|
||||
{
|
||||
acb_get_mag(b, acb_mat_entry(LA, i, j));
|
||||
result = result && (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
}
|
||||
}
|
||||
|
||||
acb_mat_clear(LA);
|
||||
acb_mat_clear(D);
|
||||
}
|
||||
|
||||
if (result && wantR)
|
||||
{
|
||||
acb_mat_t AR, D;
|
||||
acb_mat_init(AR, n, n);
|
||||
acb_mat_init(D, n, n);
|
||||
|
||||
/* Check AR - R lambda = 0 */
|
||||
acb_mat_approx_mul(AR, A, R, prec);
|
||||
for (i = 0; i < n; i++)
|
||||
acb_set(acb_mat_entry(D, i, i), E + i);
|
||||
acb_mat_approx_mul(D, R, D, prec);
|
||||
acb_mat_sub(AR, AR, D, prec);
|
||||
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
for (j = 0; j < n; j++)
|
||||
{
|
||||
acb_get_mag(b, acb_mat_entry(AR, i, j));
|
||||
result = result && (mag_cmp_2exp_si(b, -goal) < 0);
|
||||
}
|
||||
}
|
||||
|
||||
acb_mat_clear(AR);
|
||||
acb_mat_clear(D);
|
||||
}
|
||||
|
||||
if (result)
|
||||
break;
|
||||
|
||||
if (prec > 300)
|
||||
if (prec > 2000)
|
||||
{
|
||||
flint_printf("FAIL (convergence, DFT matrix)\n\n");
|
||||
flint_printf("FAIL (convergence, dft = %d)\n\n", dft);
|
||||
flint_printf("n = %wd\n\n", n);
|
||||
acb_mat_printd(A, 10);
|
||||
flint_printf("\n\n");
|
||||
for (i = 0; i < n; i++)
|
||||
{
|
||||
acb_printn(E + i, 50, 0);
|
||||
flint_printf("\n");
|
||||
}
|
||||
flint_printf("\n");
|
||||
if (wantL)
|
||||
{
|
||||
flint_printf("L = \n");
|
||||
acb_mat_printd(L, 10);
|
||||
flint_printf("\n\n");
|
||||
}
|
||||
if (wantR)
|
||||
{
|
||||
flint_printf("R = \n");
|
||||
acb_mat_printd(R, 10);
|
||||
flint_printf("\n\n");
|
||||
}
|
||||
flint_abort();
|
||||
}
|
||||
}
|
||||
|
||||
acb_mat_clear(A);
|
||||
acb_mat_clear(L);
|
||||
acb_mat_clear(R);
|
||||
_acb_vec_clear(E, n);
|
||||
acb_clear(t);
|
||||
mag_clear(b);
|
||||
|
|
|
|||
|
|
@ -5,12 +5,23 @@
|
|||
|
||||
An :type:`acb_mat_t` represents a dense matrix over the complex numbers,
|
||||
implemented as an array of entries of type :type:`acb_struct`.
|
||||
|
||||
The dimension (number of rows and columns) of a matrix is fixed at
|
||||
initialization, and the user must ensure that inputs and outputs to
|
||||
an operation have compatible dimensions. The number of rows or columns
|
||||
in a matrix can be zero.
|
||||
|
||||
.. note::
|
||||
|
||||
Methods prefixed with *acb_mat_approx* treat all input entries
|
||||
as floating-point numbers (ignoring the radii of the balls) and
|
||||
compute floating-point output (balls with zero radius) representing
|
||||
approximate solutions *without error bounds*.
|
||||
All other methods compute rigorous error bounds.
|
||||
The *approx* methods are typically useful for computing initial
|
||||
values or preconditioners for rigorous solvers.
|
||||
Some users may also find *approx* methods useful
|
||||
for doing ordinary numerical linear algebra in applications where
|
||||
error bounds are not needed.
|
||||
|
||||
Types, macros and constants
|
||||
-------------------------------------------------------------------------------
|
||||
|
|
@ -268,6 +279,13 @@ Arithmetic
|
|||
Sets *res* to *mat* raised to the power *exp*. Requires that *mat*
|
||||
is a square matrix.
|
||||
|
||||
.. function:: void acb_mat_approx_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
|
||||
|
||||
Approximate matrix multiplication. The input radii are ignored and
|
||||
the output matrix is set to an approximate floating-point result.
|
||||
For performance reasons, the radii in the output matrix will *not*
|
||||
necessarily be written (zeroed), but will remain zero if they
|
||||
are already zeroed in *res* before calling this function.
|
||||
|
||||
Scalar arithmetic
|
||||
-------------------------------------------------------------------------------
|
||||
|
|
@ -435,6 +453,24 @@ Gaussian elimination and solving
|
|||
versions and additionally handles small or triangular matrices
|
||||
by direct formulas.
|
||||
|
||||
.. function:: void acb_mat_approx_solve_triu(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec)
|
||||
|
||||
.. function:: void acb_mat_approx_solve_tril(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec)
|
||||
|
||||
.. function:: int acb_mat_approx_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)
|
||||
|
||||
.. function:: void acb_mat_approx_solve_lu_precomp(acb_mat_t X, const slong * perm, const acb_mat_t A, const acb_mat_t B, slong prec)
|
||||
|
||||
.. function:: int acb_mat_approx_solve(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec)
|
||||
|
||||
These methods perform approximate solving *without any error control*.
|
||||
The radii in the input matrices are ignored, the computations are done
|
||||
numerically with floating-point arithmetic (using ordinary
|
||||
Gaussian elimination and triangular solving, accelerated through
|
||||
the use of block recursive strategies for large matrices), and the
|
||||
output matrices are set to the approximate floating-point results with
|
||||
zeroed error bounds.
|
||||
|
||||
Characteristic polynomial
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -484,66 +520,27 @@ Component and error operations
|
|||
|
||||
Adds *err* in-place to the radii of the entries of *mat*.
|
||||
|
||||
Approximate linear algebra
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
Matrix multiplication
|
||||
.....................
|
||||
|
||||
.. function:: void acb_mat_approx_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
|
||||
|
||||
Approximate matrix multiplication. The input radii are ignored and
|
||||
the output matrix is set to an approximate floating-point result.
|
||||
The radii in the output matrix will *not* necessarily be zeroed.
|
||||
|
||||
Solving
|
||||
.....................
|
||||
|
||||
.. function:: void acb_mat_approx_solve_triu(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec)
|
||||
|
||||
.. function:: void acb_mat_approx_solve_tril(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec)
|
||||
|
||||
.. function:: int acb_mat_approx_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)
|
||||
|
||||
.. function:: void acb_mat_approx_solve_lu_precomp(acb_mat_t X, const slong * perm, const acb_mat_t A, const acb_mat_t B, slong prec)
|
||||
|
||||
.. function:: int acb_mat_approx_solve(acb_mat_t X, const acb_mat_t A, const acb_mat_t B, slong prec)
|
||||
|
||||
These methods perform approximate solving *without any error control*.
|
||||
The radii in the input matrices are ignored, the computations are done
|
||||
numerically with floating-point arithmetic (using ordinary
|
||||
Gaussian elimination and triangular solving, accelerated through
|
||||
the use of block recursive strategies for large matrices), and the
|
||||
output matrices are set to the approximate floating-point results with
|
||||
zeroed error bounds.
|
||||
|
||||
Approximate solutions are useful for computing preconditioning matrices
|
||||
for certified solutions. Some users may also find these methods useful
|
||||
for doing ordinary numerical linear algebra in applications where
|
||||
error bounds are not needed.
|
||||
|
||||
Eigenvalues and eigenvectors
|
||||
............................
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
.. 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)
|
||||
|
||||
Computes floating-point approximations of the eigenvalues of the
|
||||
given *n* by *n* matrix *A*. Approximations of all the *n*
|
||||
eigenvalues (counting repetitions) of *A* are
|
||||
written to the vector *E*, in no particular order.
|
||||
Uses the implicitly shifted QR algorithm with reduction
|
||||
to Hessenberg form.
|
||||
|
||||
Not implemented: if *L* and/or *R* are non-*NULL*, approximations of the left
|
||||
eigenvectors are written to the rows of *L* and approximations
|
||||
of the right eigenvectors are written to the columns of *R*,
|
||||
respectively.
|
||||
Computes floating-point approximations of all the *n* eigenvalues
|
||||
(and optionally eigenvectors) of the
|
||||
given *n* by *n* matrix *A*. The approximations of the
|
||||
eigenvalues are written to the vector *E*, in no particular order.
|
||||
If *L* is not *NULL*, approximations of the corresponding left
|
||||
eigenvectors are written to the rows of *L*. If *R* is not *NULL*,
|
||||
approximations of the corresponding
|
||||
right eigenvectors are written to the columns of *R*.
|
||||
|
||||
The parameters *tol* and *maxiter* can be used to control the
|
||||
target numerical error and the maximum allowed number of iterations
|
||||
target numerical error and the maximum number of iterations
|
||||
allowed before giving up. Passing *NULL* and 0 respectively results
|
||||
in default values being used.
|
||||
|
||||
Uses the implicitly shifted QR algorithm with reduction
|
||||
to Hessenberg form.
|
||||
No guarantees are made about the accuracy of the output. A nonzero
|
||||
return value indicates that the QR iteration converged numerically,
|
||||
but this is only a heuristic termination test and does not imply
|
||||
|
|
|
|||
|
|
@ -5,12 +5,23 @@
|
|||
|
||||
An :type:`arb_mat_t` represents a dense matrix over the real numbers,
|
||||
implemented as an array of entries of type :type:`arb_struct`.
|
||||
|
||||
The dimension (number of rows and columns) of a matrix is fixed at
|
||||
initialization, and the user must ensure that inputs and outputs to
|
||||
an operation have compatible dimensions. The number of rows or columns
|
||||
in a matrix can be zero.
|
||||
|
||||
.. note::
|
||||
|
||||
Methods prefixed with *arb_mat_approx* treat all input entries
|
||||
as floating-point numbers (ignoring the radii of the balls) and
|
||||
compute floating-point output (balls with zero radius) representing
|
||||
approximate solutions *without error bounds*.
|
||||
All other methods compute rigorous error bounds.
|
||||
The *approx* methods are typically useful for computing initial
|
||||
values or preconditioners for rigorous solvers.
|
||||
Some users may also find *approx* methods useful
|
||||
for doing ordinary numerical linear algebra in applications where
|
||||
error bounds are not needed.
|
||||
|
||||
Types, macros and constants
|
||||
-------------------------------------------------------------------------------
|
||||
|
|
@ -294,6 +305,12 @@ Arithmetic
|
|||
This function assumes that all exponents are small and is unsafe
|
||||
for general use.
|
||||
|
||||
.. function:: void arb_mat_approx_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)
|
||||
|
||||
Approximate matrix multiplication. The input radii are ignored and
|
||||
the output matrix is set to an approximate floating-point result.
|
||||
The radii in the output matrix will *not* necessarily be zeroed.
|
||||
|
||||
Scalar arithmetic
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -468,6 +485,29 @@ Gaussian elimination and solving
|
|||
versions and additionally handles small or triangular matrices
|
||||
by direct formulas.
|
||||
|
||||
.. function:: void arb_mat_approx_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)
|
||||
|
||||
.. function:: void arb_mat_approx_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)
|
||||
|
||||
.. function:: int arb_mat_approx_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec)
|
||||
|
||||
.. function:: void arb_mat_approx_solve_lu_precomp(arb_mat_t X, const slong * perm, const arb_mat_t A, const arb_mat_t B, slong prec)
|
||||
|
||||
.. function:: int arb_mat_approx_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)
|
||||
|
||||
These methods perform approximate solving *without any error control*.
|
||||
The radii in the input matrices are ignored, the computations are done
|
||||
numerically with floating-point arithmetic (using ordinary
|
||||
Gaussian elimination and triangular solving, accelerated through
|
||||
the use of block recursive strategies for large matrices), and the
|
||||
output matrices are set to the approximate floating-point results with
|
||||
zeroed error bounds.
|
||||
|
||||
Approximate solutions are useful for computing preconditioning matrices
|
||||
for certified solutions. Some users may also find these methods useful
|
||||
for doing ordinary numerical linear algebra in applications where
|
||||
error bounds are not needed.
|
||||
|
||||
Cholesky decomposition and solving
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -659,35 +699,3 @@ Component and error operations
|
|||
.. function:: void arb_mat_add_error_mag(arb_mat_t mat, const mag_t err)
|
||||
|
||||
Adds *err* in-place to the radii of the entries of *mat*.
|
||||
|
||||
Approximate solving
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
.. function:: void arb_mat_approx_mul(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec)
|
||||
|
||||
Approximate matrix multiplication. The input radii are ignored and
|
||||
the output matrix is set to an approximate floating-point result.
|
||||
The radii in the output matrix will *not* necessarily be zeroed.
|
||||
|
||||
.. function:: void arb_mat_approx_solve_triu(arb_mat_t X, const arb_mat_t U, const arb_mat_t B, int unit, slong prec)
|
||||
|
||||
.. function:: void arb_mat_approx_solve_tril(arb_mat_t X, const arb_mat_t L, const arb_mat_t B, int unit, slong prec)
|
||||
|
||||
.. function:: int arb_mat_approx_lu(slong * P, arb_mat_t LU, const arb_mat_t A, slong prec)
|
||||
|
||||
.. function:: void arb_mat_approx_solve_lu_precomp(arb_mat_t X, const slong * perm, const arb_mat_t A, const arb_mat_t B, slong prec)
|
||||
|
||||
.. function:: int arb_mat_approx_solve(arb_mat_t X, const arb_mat_t A, const arb_mat_t B, slong prec)
|
||||
|
||||
These methods perform approximate solving *without any error control*.
|
||||
The radii in the input matrices are ignored, the computations are done
|
||||
numerically with floating-point arithmetic (using ordinary
|
||||
Gaussian elimination and triangular solving, accelerated through
|
||||
the use of block recursive strategies for large matrices), and the
|
||||
output matrices are set to the approximate floating-point results with
|
||||
zeroed error bounds.
|
||||
|
||||
Approximate solutions are useful for computing preconditioning matrices
|
||||
for certified solutions. Some users may also find these methods useful
|
||||
for doing ordinary numerical linear algebra in applications where
|
||||
error bounds are not needed.
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue