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

Merge pull request #105 from argriffing/bool-mat-for-exp-sparsity-structure

Browse source 
MAINT: use bool_mat for structure detection in matrix exp
This commit is contained in:
Fredrik Johansson 2016-03-03 00:44:02 +01:00
commit 2e6683c34a
2 changed files with 139 additions and 192 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

174
acb_mat/exp.c
View file

@ -25,77 +25,10 @@
#include "double_extras.h" #include "double_extras.h"
#include "acb_mat.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); 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) */ /* evaluates the truncated Taylor series (assumes no aliasing) */
void void
_acb_mat_exp_taylor(acb_mat_t S, const acb_mat_t A, slong N, slong prec) _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 void
acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec) acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
{ {
slong i, j, dim, wp, N, q, r; slong i, j, dim, nz;
mag_t norm, err; bool_mat_t C;
acb_mat_t T;
int is_real;
if (!acb_mat_is_square(A)) if (!acb_mat_is_square(A))
{ {
@ -210,20 +166,47 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
return; return;
} }
/* todo: generalize to (possibly permuted) block diagonal structure */ nz = _acb_mat_count_is_zero(A);
if (_acb_mat_is_diagonal(A))
if (nz == dim * dim)
{ {
if (B != A) acb_mat_one(B);
{
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);
}
return; return;
} }
if (nz == 0)
{
bool_mat_init(C, dim, dim);
bool_mat_complement(C, C);
}
else
{
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))
{
_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); is_real = acb_mat_is_real(A);
wp = prec + 3 * FLINT_BIT_COUNT(prec); wp = prec + 3 * FLINT_BIT_COUNT(prec);
@ -234,22 +217,6 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
acb_mat_bound_inf_norm(norm, A); acb_mat_bound_inf_norm(norm, A);
if (mag_is_zero(norm))
{
acb_mat_one(B);
}
else
{
fmpz_mat_t S;
int using_structure;
using_structure = _acb_mat_any_is_zero(A);
if (using_structure)
{
fmpz_mat_init(S, dim, dim);
_acb_mat_exp_get_structure(S, A);
}
q = pow(wp, 0.25); /* wanted magnitude */ q = pow(wp, 0.25); /* wanted magnitude */
if (mag_cmp_2exp_si(norm, 2 * wp) > 0) /* too big */ if (mag_cmp_2exp_si(norm, 2 * wp) > 0) /* too big */
@ -271,21 +238,17 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
{ {
for (i = 0; i < dim; i++) for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++) for (j = 0; j < dim; j++)
if (bool_mat_get_entry(C, i, j))
arb_add_error_mag(acb_realref(acb_mat_entry(B, i, j)), err); arb_add_error_mag(acb_realref(acb_mat_entry(B, i, j)), err);
} }
else else
{ {
for (i = 0; i < dim; i++) for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++) for (j = 0; j < dim; j++)
if (bool_mat_get_entry(C, i, j))
acb_add_error_mag(acb_mat_entry(B, i, j), err); 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);
}
for (i = 0; i < r; i++) for (i = 0; i < r; i++)
{ {
acb_mat_sqr(T, B, wp); acb_mat_sqr(T, B, wp);
@ -296,10 +259,11 @@ acb_mat_exp(acb_mat_t B, const acb_mat_t A, slong prec)
for (j = 0; j < dim; j++) for (j = 0; j < dim; j++)
acb_set_round(acb_mat_entry(B, i, j), acb_set_round(acb_mat_entry(B, i, j),
acb_mat_entry(B, i, j), prec); acb_mat_entry(B, i, j), prec);
}
mag_clear(norm); mag_clear(norm);
mag_clear(err); mag_clear(err);
acb_mat_clear(T); acb_mat_clear(T);
} }
bool_mat_clear(C);
}

131
arb_mat/exp.c
View file

@ -26,46 +26,10 @@
#include "fmpz_mat.h" #include "fmpz_mat.h"
#include "double_extras.h" #include "double_extras.h"
#include "arb_mat.h" #include "arb_mat.h"
#include "fmpz_mat_extras.h" #include "bool_mat.h"
#define LOG2_OVER_E 0.25499459743395350926 #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 slong
_arb_mat_exp_choose_N(const mag_t norm, slong prec) _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 void
arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec) arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
{ {
slong i, j, dim, wp, N, q, r; slong i, j, dim, nz;
mag_t norm, err; bool_mat_t C;
arb_mat_t T;
if (!arb_mat_is_square(A)) if (!arb_mat_is_square(A))
{ {
@ -203,20 +181,46 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
return; return;
} }
/* todo: generalize to (possibly permuted) block diagonal structure */ nz = arb_mat_count_is_zero(A);
if (_arb_mat_is_diagonal(A))
if (nz == dim * dim)
{ {
if (B != A) arb_mat_one(B);
{
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);
}
return; return;
} }
if (nz == 0)
{
bool_mat_init(C, dim, dim);
bool_mat_complement(C, C);
}
else
{
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))
{
_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); wp = prec + 3 * FLINT_BIT_COUNT(prec);
mag_init(norm); mag_init(norm);
@ -225,23 +229,6 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
arb_mat_bound_inf_norm(norm, A); arb_mat_bound_inf_norm(norm, A);
if (mag_is_zero(norm))
{
arb_mat_one(B);
}
else
{
fmpz_mat_t S;
int using_structure;
using_structure = arb_mat_count_is_zero(A) > 0;
if (using_structure)
{
fmpz_mat_init(S, dim, dim);
arb_mat_entrywise_not_is_zero(S, A);
fmpz_mat_transitive_closure(S, S);
}
q = pow(wp, 0.25); /* wanted magnitude */ q = pow(wp, 0.25); /* wanted magnitude */
if (mag_cmp_2exp_si(norm, 2 * wp) > 0) /* too big */ if (mag_cmp_2exp_si(norm, 2 * wp) > 0) /* too big */
@ -261,14 +248,9 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
for (i = 0; i < dim; i++) for (i = 0; i < dim; i++)
for (j = 0; j < dim; j++) for (j = 0; j < dim; j++)
if (bool_mat_get_entry(C, i, j))
arb_add_error_mag(arb_mat_entry(B, i, j), err); 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);
}
for (i = 0; i < r; i++) for (i = 0; i < r; i++)
{ {
arb_mat_sqr(T, B, wp); arb_mat_sqr(T, B, wp);
@ -279,10 +261,11 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec)
for (j = 0; j < dim; j++) for (j = 0; j < dim; j++)
arb_set_round(arb_mat_entry(B, i, j), arb_set_round(arb_mat_entry(B, i, j),
arb_mat_entry(B, i, j), prec); arb_mat_entry(B, i, j), prec);
}
mag_clear(norm); mag_clear(norm);
mag_clear(err); mag_clear(err);
arb_mat_clear(T); arb_mat_clear(T);
} }
bool_mat_clear(C);
}