From 260c530197a599b92a5ab3306ab9e9f365432ce5 Mon Sep 17 00:00:00 2001 From: alex Date: Tue, 23 Feb 2016 13:46:29 -0500 Subject: [PATCH] MAINT: document sparsity-related arb_mat functions and move them to a separate file --- arb_mat.h | 18 ++++++++++++ arb_mat/adjacency.c | 62 ++++++++++++++++++++++++++++++++++++++++++ arb_mat/exp.c | 51 +++++++++------------------------- doc/source/arb_mat.rst | 16 +++++++++++ 4 files changed, 109 insertions(+), 38 deletions(-) create mode 100644 arb_mat/adjacency.c diff --git a/arb_mat.h b/arb_mat.h index 9690adfd..833f6ed4 100644 --- a/arb_mat.h +++ b/arb_mat.h @@ -299,6 +299,24 @@ void arb_mat_charpoly(arb_poly_t cp, const arb_mat_t mat, slong prec); void arb_mat_trace(arb_t trace, const arb_mat_t mat, slong prec); +/* Sparsity structure */ + +void _arb_mat_entrywise_nonzero_round_up(fmpz_mat_t A, const arb_mat_t src); + +ARB_MAT_INLINE void +arb_mat_adjacency(fmpz_mat_t A, const arb_mat_t src) +{ + _arb_mat_entrywise_nonzero_round_up(A, src); +} + +slong _arb_mat_count_nonzero_round_up(const arb_mat_t src); + +ARB_MAT_INLINE slong +arb_mat_count_nonzero(const arb_mat_t src) +{ + return _arb_mat_count_nonzero_round_up(src); +} + #ifdef __cplusplus } #endif diff --git a/arb_mat/adjacency.c b/arb_mat/adjacency.c new file mode 100644 index 00000000..ef65d695 --- /dev/null +++ b/arb_mat/adjacency.c @@ -0,0 +1,62 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with FLINT; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2016 Arb authors + +******************************************************************************/ + +#include "arb_mat.h" + +void +_arb_mat_entrywise_nonzero_round_up(fmpz_mat_t A, const arb_mat_t src) +{ + slong i, j; + fmpz_mat_zero(A); + for (i = 0; i < arb_mat_nrows(src); i++) + { + for (j = 0; j < arb_mat_ncols(src); j++) + { + if (!arb_is_zero(arb_mat_entry(src, i, j))) + { + fmpz_one(fmpz_mat_entry(A, i, j)); + } + } + } +} + +slong +_arb_mat_count_nonzero_round_up(const arb_mat_t src) +{ + slong nnz; + slong i, j; + nnz = 0; + for (i = 0; i < arb_mat_nrows(src); i++) + { + for (j = 0; j < arb_mat_ncols(src); j++) + { + if (!arb_is_zero(arb_mat_entry(src, i, j))) + { + nnz++; + } + } + } + return nnz; +} diff --git a/arb_mat/exp.c b/arb_mat/exp.c index 0d80f60f..aa3a33ef 100644 --- a/arb_mat/exp.c +++ b/arb_mat/exp.c @@ -32,7 +32,7 @@ /* Warshall's algorithm */ void -_fmpz_mat_transitive_closure(fmpz_mat_t A) +_fmpz_mat_transitive_closure(fmpz_mat_t B, fmpz_mat_t A) { slong k, i, j, dim; dim = fmpz_mat_nrows(A); @@ -43,17 +43,22 @@ _fmpz_mat_transitive_closure(fmpz_mat_t A) abort(); } + if (A != B) + { + fmpz_mat_set(B, A); + } + for (k = 0; k < dim; k++) { for (i = 0; i < dim; i++) { for (j = 0; j < dim; j++) { - if (fmpz_is_zero(fmpz_mat_entry(A, i, j)) && - !fmpz_is_zero(fmpz_mat_entry(A, i, k)) && - !fmpz_is_zero(fmpz_mat_entry(A, k, j))) + if (fmpz_is_zero(fmpz_mat_entry(B, i, j)) && + !fmpz_is_zero(fmpz_mat_entry(B, i, k)) && + !fmpz_is_zero(fmpz_mat_entry(B, k, j))) { - fmpz_one(fmpz_mat_entry(A, i, j)); + fmpz_one(fmpz_mat_entry(B, i, j)); } } } @@ -71,37 +76,6 @@ _arb_mat_is_diagonal(const arb_mat_t A) return 1; } -int -_arb_mat_any_is_zero(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 (arb_is_zero(arb_mat_entry(A, i, j))) - return 1; - return 0; -} - -void -_arb_mat_exp_get_structure(fmpz_mat_t C, const arb_mat_t A) -{ - slong i, j, dim; - - dim = arb_mat_nrows(A); - fmpz_mat_zero(C); - for (i = 0; i < dim; i++) - { - for (j = 0; j < dim; j++) - { - if (!arb_is_zero(arb_mat_entry(A, i, j))) - { - fmpz_one(fmpz_mat_entry(C, i, j)); - } - } - } - _fmpz_mat_transitive_closure(C); -} - void _arb_mat_exp_set_structure(arb_mat_t B, const fmpz_mat_t C) { @@ -296,11 +270,12 @@ arb_mat_exp(arb_mat_t B, const arb_mat_t A, slong prec) fmpz_mat_t S; int using_structure; - using_structure = _arb_mat_any_is_zero(A); + using_structure = arb_mat_count_nonzero(A) < dim * dim; if (using_structure) { fmpz_mat_init(S, dim, dim); - _arb_mat_exp_get_structure(S, A); + arb_mat_adjacency(S, A); + _fmpz_mat_transitive_closure(S, S); } q = pow(wp, 0.25); /* wanted magnitude */ diff --git a/doc/source/arb_mat.rst b/doc/source/arb_mat.rst index e9f1875c..f0787fa0 100644 --- a/doc/source/arb_mat.rst +++ b/doc/source/arb_mat.rst @@ -324,3 +324,19 @@ Special functions Sets *trace* to the trace of the matrix, i.e. the sum of entries on the main diagonal of *mat*. The matrix is required to be square. + +Sparsity structure +------------------------------------------------------------------------------- + +.. function:: void _arb_mat_entrywise_nonzero_round_up(fmpz_mat_t A, const arb_mat_t src) + +.. function:: void arb_mat_adjacency(fmpz_mat_t A, const arb_mat_t src) + + Each entry of *A* is set to 0 or 1 depending on whether or not + the corresponding entry of *src* is zero. + +.. function:: slong _arb_mat_count_nonzero_round_up(const arb_mat_t src) + +.. function:: slong arb_mat_count_nonzero(const arb_mat_t src) + + Returns the number of nonzero entries of the matrix.