MAINT: document sparsity-related arb_mat functions and move them to a separate file

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

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;
+}

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)) &&
+                if (fmpz_is_zero(fmpz_mat_entry(B, i, j)) &&
-                    !fmpz_is_zero(fmpz_mat_entry(A, i, k)) &&
+                    !fmpz_is_zero(fmpz_mat_entry(B, i, k)) &&
-                    !fmpz_is_zero(fmpz_mat_entry(A, k, j)))
+                    !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 */

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.