Sparse3

CMakeLists.txt

@@ -259,7 +259,7 @@ set(_BUILD_DIRS
     acb_theta       dirichlet       bernoulli       hypgeom
 
     gr              gr_generic      gr_vec          gr_mat
-    gr_poly         gr_mpoly        gr_special
+    gr_poly         gr_mpoly        gr_sparse_vec   gr_sparse_mat        gr_special
 
     calcium
     fmpz_mpoly_q

Makefile.in

@@ -199,7 +199,8 @@ HEADER_DIRS :=                                                              \
         acb_theta       dirichlet       bernoulli       hypgeom             \
                                                                             \
         gr              gr_generic      gr_vec          gr_mat              \
-        gr_poly         gr_mpoly        gr_special                          \
+        gr_poly         gr_mpoly        gr_sparse_vec   gr_sparse_mat       \
+        gr_special                                                          \
                                                                             \
         calcium                                                             \
         fmpz_mpoly_q                                                        \

doc/source/gr_mat.rst

@@ -199,20 +199,26 @@ Basic row, column and entry operations
     ``c`` is the number of columns of ``mat``. If ``perm`` is non-``NULL``, the
     permutation of the columns will also be applied to ``perm``.
 
-.. function:: truth_t gr_mat_is_empty(const gr_mat_t mat, gr_ctx_t ctx)
+.. macro:: gr_mat_is_empty(mat, ctx)
 
     Returns whether *mat* is an empty matrix, having either zero
     rows or zero column. This predicate is always decidable (even if
     the underlying ring is not computable), returning
     ``T_TRUE`` or ``T_FALSE``.
 
-.. function:: truth_t gr_mat_is_square(const gr_mat_t mat, gr_ctx_t ctx)
+.. macro:: gr_mat_is_square(mat, ctx)
 
     Returns whether *mat* is a square matrix, having the same number
     of rows as columns (not the same thing as being a perfect square!).
     This predicate is always decidable (even if the underlying ring
     is not computable), returning ``T_TRUE`` or ``T_FALSE``.
 
+.. macro:: gr_mat_is_compatible(mat1, mat2, ctx)
+
+    Returns whether *mat1* and *mat2* have the same dimensions.
+    This predicate is always decidable (even if the underlying ring
+    is not computable), returning ``T_TRUE`` or ``T_FALSE``.
+
 
 Arithmetic
 -------------------------------------------------------------------------------

doc/source/index.rst

@@ -57,6 +57,7 @@ Generic rings
    gr_generic.rst
    gr_special.rst
    gr_vec.rst
+   gr_sparse_vec.rst
    gr_mat.rst
    gr_poly.rst
    gr_mpoly.rst

src/fmpz_mod_mat/det.c

@@ -39,4 +39,4 @@ void fmpz_mod_mat_det(fmpz_t res, const fmpz_mod_mat_t mat, const fmpz_mod_ctx_t
             GR_MUST_SUCCEED(gr_mat_det_berkowitz(res, (const gr_mat_struct *) mat, gr_ctx));
         }
     }
-}
+}

src/gr_mat.h

@@ -99,6 +99,8 @@ WARN_UNUSED_RESULT int gr_mat_randops(gr_mat_t mat, flint_rand_t state, slong co
 WARN_UNUSED_RESULT int gr_mat_randpermdiag(int * parity, gr_mat_t mat, flint_rand_t state, gr_ptr diag, slong n, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_mat_randrank(gr_mat_t mat, flint_rand_t state, slong rank, gr_ctx_t ctx);
 
+#define gr_mat_is_compatible(mat1, mat2, ctx) (((mat1)->r == (mat2)->r && (mat1)->c == (mat2)->c) ? T_TRUE : T_FALSE)
+
 GR_MAT_INLINE truth_t
 gr_mat_is_empty(const gr_mat_t mat, gr_ctx_t FLINT_UNUSED(ctx))
 {
@@ -135,17 +137,50 @@ WARN_UNUSED_RESULT int gr_mat_sub(gr_mat_t res, const gr_mat_t mat1, const gr_ma
 
 /* todo: test, wrap; div; more conversions */
 WARN_UNUSED_RESULT int gr_mat_add_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_add_scalar_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_add_scalar_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_add_scalar_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_add_scalar_fmpq(gr_mat_t res, const gr_mat_t mat, fmpq_t x, gr_ctx_t ctx);
+
 WARN_UNUSED_RESULT int gr_mat_sub_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_sub_scalar_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_sub_scalar_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_sub_scalar_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_sub_scalar_fmpq(gr_mat_t res, const gr_mat_t mat, fmpq_t x, gr_ctx_t ctx);
+
 WARN_UNUSED_RESULT int gr_mat_mul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_mul_scalar_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_mul_scalar_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_mul_scalar_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_mul_scalar_fmpq(gr_mat_t res, const gr_mat_t mat, fmpq_t x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_mul_scalar_2exp_si(gr_mat_t res, const gr_mat_t mat, slong e, gr_ctx_t ctx);
+
+WARN_UNUSED_RESULT int gr_mat_div_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_div_scalar_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_div_scalar_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_div_scalar_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_div_scalar_fmpq(gr_mat_t res, const gr_mat_t mat, fmpq_t x, gr_ctx_t ctx);
+
+WARN_UNUSED_RESULT int gr_mat_divexact_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_divexact_scalar_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_divexact_scalar_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_divexact_scalar_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_divexact_scalar_fmpq(gr_mat_t res, const gr_mat_t mat, fmpq_t x, gr_ctx_t ctx);
+
 WARN_UNUSED_RESULT int gr_mat_addmul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_mat_submul_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx);
-WARN_UNUSED_RESULT int gr_mat_div_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx);
+
+WARN_UNUSED_RESULT int gr_mat_mul_vec(gr_ptr v, const gr_mat_t A, gr_srcptr u, gr_ctx_t ctx);
 
 WARN_UNUSED_RESULT int gr_mat_mul_classical(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_mat_mul_strassen(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_mat_mul_generic(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_mat_mul(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
 
+WARN_UNUSED_RESULT int gr_mat_addmul_classical(gr_mat_t D, const gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_addmul_generic(gr_mat_t D, const gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_mat_addmul(gr_mat_t D, const gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx);
+
 /* todo */
 GR_MAT_INLINE WARN_UNUSED_RESULT int
 gr_mat_sqr(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
@@ -214,6 +249,27 @@ WARN_UNUSED_RESULT int gr_mat_hilbert(gr_mat_t mat, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_mat_hadamard(gr_mat_t mat, gr_ctx_t ctx);
 /* todo: dft, dct */
 
+#define _GR_MAT_SHALLOW_TRANSPOSE(AT, A, sz) { \
+    for (j = 0; j < A->c; j++) \
+    { \
+        for (i = 0; i < A->r; i++) \
+        { \
+            set_shallow(GR_MAT_ENTRY(AT, j, i, sz), GR_MAT_ENTRY(A, i, j, sz), ctx); \
+        } \
+    } \
+}
+
+#define _GR_MAT_INIT_SHALLOW_TRANSPOSE(AT, A, sz) { \
+    AT->r = A->c; \
+    AT->c = A->r; \
+    AT->rows = flint_malloc(AT->r * sizeof(gr_ptr)); \
+    AT->entries = TMP_ALLOC(sz * A->c * A->r); \
+    for (j = 0; j < A->c; j++) \
+        AT->rows[j] = GR_ENTRY(AT->entries, j * A->r, sz); \
+    _GR_MAT_SHALLOW_TRANSPOSE(AT, A, sz) \
+}
+
+
 WARN_UNUSED_RESULT int gr_mat_transpose(gr_mat_t B, const gr_mat_t A, gr_ctx_t ctx);
 
 WARN_UNUSED_RESULT int gr_mat_nonsingular_solve_tril_classical(gr_mat_t X, const gr_mat_t L, const gr_mat_t B, int unit, gr_ctx_t ctx);

src/gr_mat/add_scalar.c

@@ -11,35 +11,56 @@
 
 #include "gr_mat.h"
 
-int
-gr_mat_add_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx)
-{
-    slong i, j, r, c;
-    slong sz = ctx->sizeof_elem;
-    int status = GR_SUCCESS;
-
-    r = gr_mat_nrows(res, ctx);
-    c = gr_mat_ncols(res, ctx);
-
-    if (res == mat)
-    {
-        for (i = 0; i < FLINT_MIN(r, c); i++)
-            status |= gr_add(GR_MAT_ENTRY(res, i, i, sz), GR_MAT_ENTRY(res, i, i, sz), x, ctx);
-    }
-    else
-    {
-        for (i = 0; i < r; i++)
-        {
-            for (j = 0; j < c; j++)
-            {
-                /* todo: vectorize */
-                if (i == j)
-                    status |= gr_add(GR_MAT_ENTRY(res, i, j, sz), GR_MAT_ENTRY(mat, i, j, sz), x, ctx);
-                else
-                    status |= gr_set(GR_MAT_ENTRY(res, i, j, sz), GR_MAT_ENTRY(mat, i, j, sz), ctx);
-            }
-        }
-    }
-
-    return status;
-}
+#define GR_MAT_ADD_SCALAR(FUNC, res, mat, x, ctx)        \
+    slong i, j, r, c, sz = (ctx)->sizeof_elem;           \
+    int status = GR_SUCCESS;                             \
+    r = gr_mat_nrows(res, ctx);                          \
+    c = gr_mat_ncols(res, ctx);                          \
+    if (res == mat)                                      \
+    {                                                    \
+        for (i = 0; i < FLINT_MIN(r, c); i++)            \
+            status |= (FUNC)(                            \
+                GR_MAT_ENTRY(res, i, i, sz),             \
+                GR_MAT_ENTRY(res, i, i, sz),             \
+                x, ctx                                   \
+            );                                           \
+    }                                                    \
+    else                                                 \
+    {                                                    \
+        for (i = 0; i < r; i++)                          \
+        {                                                \
+            for (j = 0; j < c; j++)                      \
+            {                                            \
+                /* todo: vectorize */                    \
+                if (i == j)                              \
+                    status |= (FUNC)(                    \
+                        GR_MAT_ENTRY(res, i, j, sz),     \
+                        GR_MAT_ENTRY(mat, i, j, sz),     \
+                        x, ctx                           \
+                    );                                   \
+                else                                     \
+                    status |= gr_set(                    \
+                        GR_MAT_ENTRY(res, i, j, sz),     \
+                        GR_MAT_ENTRY(mat, i, j, sz),     \
+                        ctx                              \
+                    );                                   \
+            }                                            \
+        }                                                \
+    }                                                    \
+    return status;                                       \
+
+
+int gr_mat_add_scalar(gr_mat_t res, const gr_mat_t mat, gr_srcptr x, gr_ctx_t ctx)
+{ GR_MAT_ADD_SCALAR(gr_add, res, mat, x, ctx) }
+
+int gr_mat_add_scalar_si(gr_mat_t res, const gr_mat_t mat, slong x, gr_ctx_t ctx)
+{ GR_MAT_ADD_SCALAR(gr_add_si, res, mat, x, ctx) }
+
+int gr_mat_add_scalar_ui(gr_mat_t res, const gr_mat_t mat, ulong x, gr_ctx_t ctx)
+{ GR_MAT_ADD_SCALAR(gr_add_ui, res, mat, x, ctx) }
+
+int gr_mat_add_scalar_fmpz(gr_mat_t res, const gr_mat_t mat, fmpz_t x, gr_ctx_t ctx)
+{ GR_MAT_ADD_SCALAR(gr_add_fmpz, res, mat, x, ctx) }
+
+int gr_mat_add_scalar_fmpq(gr_mat_t res, const gr_mat_t mat, fmpq_t x, gr_ctx_t ctx)
+{ GR_MAT_ADD_SCALAR(gr_add_fmpq, res, mat, x, ctx) }

src/gr_mat/addmul.c

@@ -0,0 +1,24 @@
+/*
+    Copyright (C) 2022 Fredrik Johansson
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include "gr_mat.h"
+
+int
+gr_mat_addmul_generic(gr_mat_t D, const gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+{
+    return gr_mat_mul_classical(C, A, B, ctx);
+}
+
+int
+gr_mat_addmul(gr_mat_t D, const gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+{
+    return gr_mat_addmul_generic(D, C, A, B, ctx);
+}

src/gr_mat/addmul_classical.c

@@ -0,0 +1,117 @@
+/*
+    Copyright (C) 2022 Fredrik Johansson
+
+    This file is part of FLINT.
+
+    FLINT is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 3 of the License, or
+    (at your option) any later version.  See <https://www.gnu.org/licenses/>.
+*/
+
+#include <stdint.h>
+#include "gr_vec.h"
+#include "gr_mat.h"
+
+int
+gr_mat_addmul_classical(gr_mat_t D, const gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+{
+    slong ar, ac, br, bc, i, j, sz;
+    int status;
+
+    ar = gr_mat_nrows(A, ctx);
+    ac = gr_mat_ncols(A, ctx);
+    br = gr_mat_nrows(B, ctx);
+    bc = gr_mat_ncols(B, ctx);
+
+    if (gr_mat_is_compatible(C, D, ctx) == T_FALSE || ac != br || ar != gr_mat_nrows(C, ctx) || bc != gr_mat_ncols(C, ctx))
+        return GR_DOMAIN;
+
+    if (br == 0)
+    {
+        return GR_SUCCESS;
+    }
+
+    status = GR_SUCCESS;
+
+    if (A == D || B == D)
+    {
+        gr_mat_t T;
+        gr_mat_init(T, ar, bc, ctx);
+        status |= gr_mat_addmul_classical(T, C, A, B, ctx);
+        status |= gr_mat_swap_entrywise(T, D, ctx);
+        gr_mat_clear(T, ctx);
+        return status;
+    }
+
+    sz = ctx->sizeof_elem;
+
+    if (br == 1)
+    {
+        for (i = 0; i < ar; i++)
+        {
+            for (j = 0; j < bc; j++)
+            {
+                if (C != D)
+                    status |= gr_set(GR_MAT_ENTRY(D, i, j, sz), GR_MAT_ENTRY(C, i, j, sz), ctx);
+                status |= gr_addmul(GR_MAT_ENTRY(D, i, j, sz),
+                                    GR_MAT_ENTRY(A, i, 0, sz),
+                                    GR_MAT_ENTRY(B, 0, j, sz), ctx);
+            }
+        }
+    }
+    else
+    {
+        gr_ptr tmp;
+        gr_method_void_unary_op set_shallow = GR_VOID_UNARY_OP(ctx, SET_SHALLOW);
+        TMP_INIT;
+
+        TMP_START;
+        tmp = TMP_ALLOC(sz * br * bc);
+
+        /* Make a shallow transpose so that we can use dot products.
+           Inline common sizes. (Caution: are we sure about the alignment?
+           Some asserts would be nice here.)
+           Todo: we may want inlining in nonsingular_solve etc. as well. */
+        for (i = 0; i < br; i++)
+        {
+            for (j = 0; j < bc; j++)
+            {
+                switch (sz)
+                {
+#if 0
+                    case 1:
+                        ((int8_t *) GR_ENTRY(tmp, j * br + i, 1))[0] = ((int8_t *) GR_MAT_ENTRY(B, i, j, 1))[0];
+                        break;
+                    case 2:
+                        ((int16_t *) GR_ENTRY(tmp, j * br + i, 2))[0] = ((int16_t *) GR_MAT_ENTRY(B, i, j, 2))[0];
+                        break;
+                    case 4:
+                        ((int32_t *) GR_ENTRY(tmp, j * br + i, 4))[0] = ((int32_t *) GR_MAT_ENTRY(B, i, j, 4))[0];
+                        break;
+#if FLINT_BITS == 64
+                    case 8:
+                        ((int64_t *) GR_ENTRY(tmp, j * br + i, 8))[0] = ((int64_t *) GR_MAT_ENTRY(B, i, j, 8))[0];
+                        break;
+#endif
+#endif
+                    default:
+                        set_shallow(GR_ENTRY(tmp, j * br + i, sz), GR_MAT_ENTRY(B, i, j, sz), ctx);
+                }
+            }
+        }
+
+        for (i = 0; i < ar; i++)
+        {
+            for (j = 0; j < bc; j++)
+            {
+                status |= _gr_vec_dot(GR_MAT_ENTRY(D, i, j, sz), GR_MAT_ENTRY(C, i, j, sz), 0,
+                    GR_MAT_ENTRY(A, i, 0, sz), GR_ENTRY(tmp, j * br, sz), br, ctx);
+            }
+        }
+
+        TMP_END;
+    }
+
+    return status;
+}