hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 17:31:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

MAINT: use bool_mat for structure detection for matrix exp

Browse source
This commit is contained in:
alex 2016-03-02 14:56:06 -05:00
parent 298b670941
commit 247964e4d2
2 changed files with 139 additions and 192 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

190
acb_mat/exp.c
View file

@ -25,77 +25,10 @@
#include "double_extras.h"
#include "acb_mat.h"
#include "fmpz_mat_extras.h"
#include "bool_mat.h"
slong _arb_mat_exp_choose_N(const mag_t norm, slong prec);
int
_acb_mat_is_diagonal(const acb_mat_t A)
{
slong i, j;
for (i = 0; i < acb_mat_nrows(A); i++)
for (j = 0; j < acb_mat_ncols(A); j++)
if (i != j && !acb_is_zero(acb_mat_entry(A, i, j)))
return 0;
return 1;
}
int
_acb_mat_any_is_zero(const acb_mat_t A)
{
slong i, j;
for (i = 0; i < acb_mat_nrows(A); i++)
for (j = 0; j < acb_mat_ncols(A); j++)
if (acb_is_zero(acb_mat_entry(A, i, j)))
return 1;
return 0;
}
void
_acb_mat_exp_get_structure(fmpz_mat_t C, const acb_mat_t A)
{
slong i, j, dim;
dim = acb_mat_nrows(A);
fmpz_mat_zero(C);
for (i = 0; i < dim; i++)
{
for (j = 0; j < dim; j++)
{
if (!acb_is_zero(acb_mat_entry(A, i, j)))
{
fmpz_one(fmpz_mat_entry(C, i, j));
}
}
}
fmpz_mat_transitive_closure(C, C);
}
void
_acb_mat_exp_set_structure(acb_mat_t B, const fmpz_mat_t C)
{
slong i, j, dim;
dim = acb_mat_nrows(B);
for (i = 0; i < dim; i++)
{
for (j = 0; j < dim; j++)
{
if (fmpz_is_zero(fmpz_mat_entry(C, i, j)))
{
if (i == j)
{
acb_one(acb_mat_entry(B, i, j));
}
else
{
acb_zero(acb_mat_entry(B, i, j));
}
}
}
}
}
/* evaluates the truncated Taylor series (assumes no aliasing) */
void
_acb_mat_exp_taylor(acb_mat_t S, const acb_mat_t A, slong N, slong prec)
@ -185,13 +118,36 @@ _acb_mat_exp_taylor(acb_mat_t S, const acb_mat_t A, slong N, slong prec)
}
}
static slong
_acb_mat_count_is_zero(const acb_mat_t A)
{
slong nz, i, j;
for (nz = 0, i = 0; i < acb_mat_nrows(A); i++)
for (j = 0; j < acb_mat_ncols(A); j++)
nz += acb_is_zero(acb_mat_entry(A, i, j));
return nz;
}
static void
_acb_mat_exp_diagonal(acb_mat_t B, const acb_mat_t A, slong prec)
{
slong n, i;
n = acb_mat_nrows(A);
if (B != A)
{
acb_mat_zero(B);
}
for (i = 0; i < n; i++)
{
acb_exp(acb_mat_entry(B, i, i), acb_mat_entry(A, i, i), prec);
}
}
void
acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
{
slong i, j, dim, wp, N, q, r;
mag_t norm, err;
acb_mat_t T;
int is_real;
slong i, j, dim, nz;
bool_mat_t C;
if (!acb_mat_is_square(A))
{
@ -210,45 +166,56 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
return;
}
/* todo: generalize to (possibly permuted) block diagonal structure */
if (_acb_mat_is_diagonal(A))
nz = _acb_mat_count_is_zero(A);
if (nz == dim * dim)
{
if (B != A)
{
acb_mat_zero(B);
}
for (i = 0; i < dim; i++)
{
acb_exp(acb_mat_entry(B, i, i), acb_mat_entry(A, i, i), prec);
}
acb_mat_one(B);
return;
}
is_real = acb_mat_is_real(A);
wp = prec + 3 * FLINT_BIT_COUNT(prec);
mag_init(norm);
mag_init(err);
acb_mat_init(T, dim, dim);
acb_mat_bound_inf_norm(norm, A);
if (mag_is_zero(norm))
if (nz == 0)
{
acb_mat_one(B);
bool_mat_init(C, dim, dim);
bool_mat_complement(C, C);
}
else
{
fmpz_mat_t S;
int using_structure;
using_structure = _acb_mat_any_is_zero(A);
if (using_structure)
bool_mat_t S;
bool_mat_init(S, dim, dim);
for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++)
bool_mat_set_entry(S, i, j, !acb_is_zero(acb_mat_entry(A, i, j)));
if (bool_mat_is_diagonal(S))
{
fmpz_mat_init(S, dim, dim);
_acb_mat_exp_get_structure(S, A);
_acb_mat_exp_diagonal(B, A, prec);
bool_mat_clear(S);
return;
}
else
{
bool_mat_init(C, dim, dim);
bool_mat_transitive_closure(C, S);
bool_mat_clear(S);
}
}
/* evaluate using scaling and squaring of truncated taylor series */
{
slong wp, N, q, r;
mag_t norm, err;
acb_mat_t T;
int is_real;
is_real = acb_mat_is_real(A);
wp = prec + 3 * FLINT_BIT_COUNT(prec);
mag_init(norm);
mag_init(err);
acb_mat_init(T, dim, dim);
acb_mat_bound_inf_norm(norm, A);
q = pow(wp, 0.25); /* wanted magnitude */
@ -271,19 +238,15 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
{
for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++)
arb_add_error_mag(acb_realref(acb_mat_entry(B, i, j)), err);
if (bool_mat_get_entry(C, i, j))
arb_add_error_mag(acb_realref(acb_mat_entry(B, i, j)), err);
}
else
{
for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++)
acb_add_error_mag(acb_mat_entry(B, i, j), err);
}
if (using_structure)
{
_acb_mat_exp_set_structure(B, S);
fmpz_mat_clear(S);
if (bool_mat_get_entry(C, i, j))
acb_add_error_mag(acb_mat_entry(B, i, j), err);
}
for (i = 0; i < r; i++)
@ -296,10 +259,11 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
for (j = 0; j < dim; j++)
acb_set_round(acb_mat_entry(B, i, j),
acb_mat_entry(B, i, j), prec);
mag_clear(norm);
mag_clear(err);
acb_mat_clear(T);
}
mag_clear(norm);
mag_clear(err);
acb_mat_clear(T);
bool_mat_clear(C);
}

141
arb_mat/exp.c
View file

@ -26,46 +26,10 @@
#include "fmpz_mat.h"
#include "double_extras.h"
#include "arb_mat.h"
#include "fmpz_mat_extras.h"
#include "bool_mat.h"
#define LOG2_OVER_E 0.25499459743395350926
int
_arb_mat_is_diagonal(const arb_mat_t A)
{
slong i, j;
for (i = 0; i < arb_mat_nrows(A); i++)
for (j = 0; j < arb_mat_ncols(A); j++)
if (i != j && !arb_is_zero(arb_mat_entry(A, i, j)))
return 0;
return 1;
}
void
_arb_mat_exp_set_structure(arb_mat_t B, const fmpz_mat_t C)
{
slong i, j, dim;
dim = arb_mat_nrows(B);
for (i = 0; i < dim; i++)
{
for (j = 0; j < dim; j++)
{
if (fmpz_is_zero(fmpz_mat_entry(C, i, j)))
{
if (i == j)
{
arb_one(arb_mat_entry(B, i, j));
}
else
{
arb_zero(arb_mat_entry(B, i, j));
}
}
}
}
}
slong
_arb_mat_exp_choose_N(const mag_t norm, slong prec)
{
@ -179,12 +143,26 @@ _arb_mat_exp_taylor(arb_mat_t S, const arb_mat_t A, slong N, slong prec)
}
}
static void
_arb_mat_exp_diagonal(arb_mat_t B, const arb_mat_t A, slong prec)
{
slong n, i;
n = arb_mat_nrows(A);
if (B != A)
{
arb_mat_zero(B);
}
for (i = 0; i < n; i++)
{
arb_exp(arb_mat_entry(B, i, i), arb_mat_entry(A, i, i), prec);
}
}
void
arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
{
slong i, j, dim, wp, N, q, r;
mag_t norm, err;
arb_mat_t T;
slong i, j, dim, nz;
bool_mat_t C;
if (!arb_mat_is_square(A))
{
@ -203,44 +181,53 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
return;
}
/* todo: generalize to (possibly permuted) block diagonal structure */
if (_arb_mat_is_diagonal(A))
nz = arb_mat_count_is_zero(A);
if (nz == dim * dim)
{
if (B != A)
{
arb_mat_zero(B);
}
for (i = 0; i < dim; i++)
{
arb_exp(arb_mat_entry(B, i, i), arb_mat_entry(A, i, i), prec);
}
arb_mat_one(B);
return;
}
wp = prec + 3 * FLINT_BIT_COUNT(prec);
mag_init(norm);
mag_init(err);
arb_mat_init(T, dim, dim);
arb_mat_bound_inf_norm(norm, A);
if (mag_is_zero(norm))
if (nz == 0)
{
arb_mat_one(B);
bool_mat_init(C, dim, dim);
bool_mat_complement(C, C);
}
else
{
fmpz_mat_t S;
int using_structure;
using_structure = arb_mat_count_is_zero(A) > 0;
if (using_structure)
bool_mat_t S;
bool_mat_init(S, dim, dim);
for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++)
bool_mat_set_entry(S, i, j, !arb_is_zero(arb_mat_entry(A, i, j)));
if (bool_mat_is_diagonal(S))
{
fmpz_mat_init(S, dim, dim);
arb_mat_entrywise_not_is_zero(S, A);
fmpz_mat_transitive_closure(S, S);
_arb_mat_exp_diagonal(B, A, prec);
bool_mat_clear(S);
return;
}
else
{
bool_mat_init(C, dim, dim);
bool_mat_transitive_closure(C, S);
bool_mat_clear(S);
}
}
/* evaluate using scaling and squaring of truncated taylor series */
{
slong wp, N, q, r;
mag_t norm, err;
arb_mat_t T;
wp = prec + 3 * FLINT_BIT_COUNT(prec);
mag_init(norm);
mag_init(err);
arb_mat_init(T, dim, dim);
arb_mat_bound_inf_norm(norm, A);
q = pow(wp, 0.25); /* wanted magnitude */
@ -261,13 +248,8 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++)
arb_add_error_mag(arb_mat_entry(B, i, j), err);
if (using_structure)
{
_arb_mat_exp_set_structure(B, S);
fmpz_mat_clear(S);
}
if (bool_mat_get_entry(C, i, j))
arb_add_error_mag(arb_mat_entry(B, i, j), err);
for (i = 0; i < r; i++)
{
@ -279,10 +261,11 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
for (j = 0; j < dim; j++)
arb_set_round(arb_mat_entry(B, i, j),
arb_mat_entry(B, i, j), prec);
mag_clear(norm);
mag_clear(err);
arb_mat_clear(T);
}
mag_clear(norm);
mag_clear(err);
arb_mat_clear(T);
bool_mat_clear(C);
}