Merge pull request #2034 from fredrik-johansson/matmul

Matrix multiplication for nfloat

doc/source/gr.rst

@@ -844,6 +844,23 @@ Ordering methods
     of *x* is less than, equal or greater than the absolute value of *y*.
     This may return ``GR_DOMAIN`` if the ring is not an ordered ring.
 
+.. function:: truth_t gr_le(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+              truth_t gr_lt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+              truth_t gr_ge(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+              truth_t gr_gt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+              truth_t gr_abs_le(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+              truth_t gr_abs_lt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+              truth_t gr_abs_ge(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+              truth_t gr_abs_gt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+
+    Wrappers of ``gr_cmp`` and ``gr_cmpabs`` returning truth values
+    for the comparison operations ``<=``, ``<``, ``>=``, ``>``.
+
+.. function:: int gr_min(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+              int gr_max(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx)
+
+    Minimum and maximum value.
+
 Enclosure and interval methods
 ........................................................................
 

doc/source/gr_mat.rst

@@ -154,6 +154,7 @@ Assignment and special values
 .. function:: int gr_mat_set(gr_mat_t res, const gr_mat_t mat, gr_ctx_t ctx)
               int gr_mat_set_fmpz_mat(gr_mat_t res, const fmpz_mat_t mat, gr_ctx_t ctx)
               int gr_mat_set_fmpq_mat(gr_mat_t res, const fmpq_mat_t mat, gr_ctx_t ctx)
+              int gr_mat_set_gr_mat_other(gr_mat_t res, const gr_mat_t mat, gr_ctx_t mat_ctx, gr_ctx_t ctx)
 
     Sets *res* to the value of *mat*.
 
@@ -213,6 +214,58 @@ Basic row, column and entry operations
     This predicate is always decidable (even if the underlying ring
     is not computable), returning ``T_TRUE`` or ``T_FALSE``.
 
+Entrywise operations
+-------------------------------------------------------------------------------
+
+.. function:: int gr_mat_entrywise_unary_op(gr_mat_t res, gr_method_unary_op f, const gr_mat_t mat, gr_ctx_t ctx)
+
+    Sets *res* to the application of the function *f* to the
+    entries of matrix *mat*. Returns ``GR_DOMAIN`` if the matrix dimensions do not match.
+
+.. function:: int gr_mat_entrywise_binary_op(gr_mat_t res, gr_method_binary_op f, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+
+    Sets *res* to the application of the function *f*
+    to the entries of *mat1* as first argument and the entries of *mat2*
+    as second argument.
+    Returns ``GR_DOMAIN`` if the matrix dimensions do not match.
+
+.. function:: int gr_mat_entrywise_binary_op_scalar(gr_mat_t res, gr_method_binary_op f, const gr_mat_t mat, gr_srcptr c, gr_ctx_t ctx)
+
+    Sets *res* to the application of the function *f*
+    to the entries of *mat* as first argument and the scalar *c*
+    as second argument.
+    Returns ``GR_DOMAIN`` if the matrix dimensions do not match.
+
+.. function:: truth_t gr_mat_entrywise_unary_predicate_all(gr_method_unary_predicate f, const gr_mat_t mat, gr_ctx_t ctx)
+              truth_t gr_mat_entrywise_unary_predicate_any(gr_method_unary_predicate f, const gr_mat_t mat, gr_ctx_t ctx)
+
+    Returns whether the predicate *f* is true for all entries,
+    respectively for any entry, in the matrix *mat*.
+
+.. function:: truth_t gr_mat_entrywise_binary_predicate_all(gr_method_binary_predicate f, const gr_mat_t mat1, const gr_mat_t mat2, gr_ctx_t ctx)
+
+    Returns whether the binary predicate *f* is true for all entries
+    in *mat1* paired with the corresponding entries in *mat2*.
+    Returns ``T_FALSE`` if the matrix dimensions are not compatible.
+
+Norms
+-------------------------------------------------------------------------------
+
+.. function:: int gr_mat_norm_max(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+
+    Max norm: `\max_{i,j} |a_{i,j}|`.
+
+.. function:: int gr_mat_norm_1(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+
+    1-norm (largest absolute column sum): `\max_{1\le j \le n} \sum_{i=1}^m |a_{i,j}|`.
+
+.. function:: int gr_mat_norm_inf(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+
+    Infinity-norm (largest absolute row sum): `\max_{1\le i \le m} \sum_{j=1}^n |a_{i,j}|`.
+
+.. function:: int gr_mat_norm_frobenius(gr_ptr res, const gr_mat_t mat, gr_ctx_t ctx)
+
+    Frobenius norm: `\sqrt{\sum_{i,j} |a_{i,j}|^2}`.
 
 Arithmetic
 -------------------------------------------------------------------------------
@@ -802,6 +855,19 @@ on each test iteration, otherwise the given ring is tested.
     Tests the given function ``solve_impl`` for correctness as an implementation
     of :func:`gr_mat_nonsingular_solve_tril` / :func:`gr_mat_nonsingular_solve_triu`.
 
+.. function:: void gr_mat_test_approx_mul_max_norm(gr_method_mat_binary_op mul_impl, gr_srcptr rel_tol, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+
+    Tests the given implementation of matrix multiplication for accuracy
+    over an approximate numerical ring by checking that
+    `|C-AB| \le |A||B| rel\_tol` holds in the max norm,
+    using classical multiplication for reference.
+
+.. function:: void gr_mat_test_approx_mul_pos_entrywise_accurate(gr_method_mat_binary_op mul_impl, gr_srcptr rel_tol, flint_rand_t state, slong iters, slong maxn, gr_ctx_t ctx)
+
+    Tests the given implementation of matrix multiplication for accuracy
+    over an approximate numerical ring by generating nonnegative matrices
+    and checking that the entrywise relative error compared to
+    classical multiplication does not exceed *rel_tol*.
 
 .. raw:: latex
 

doc/source/nfloat.rst

@@ -314,6 +314,16 @@ code for reduced overhead.
 .. function:: int _nfloat_vec_dot(nfloat_ptr res, nfloat_srcptr initial, int subtract, nfloat_srcptr x, nfloat_srcptr y, slong len, gr_ctx_t ctx)
               int _nfloat_vec_dot_rev(nfloat_ptr res, nfloat_srcptr initial, int subtract, nfloat_srcptr x, nfloat_srcptr y, slong len, gr_ctx_t ctx)
 
+Matrix functions
+-------------------------------------------------------------------------------
+
+.. function:: int nfloat_mat_mul_fixed_classical(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+              int nfloat_mat_mul_waksman(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+              int nfloat_mat_mul_block(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, slong min_block_size, gr_ctx_t ctx)
+              int nfloat_mat_mul(gr_mat_t C, const gr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
+
+    Different implementations of matrix multiplication.
+
 Internal functions
 -------------------------------------------------------------------------------
 

src/gr.h

@@ -1109,6 +1109,25 @@ GR_INLINE WARN_UNUSED_RESULT int gr_cmpabs(int * res, gr_srcptr x, gr_srcptr y,
 GR_INLINE WARN_UNUSED_RESULT int gr_cmp_other(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx) { return GR_BINARY_OP_OTHER_GET_INT(ctx, CMP_OTHER)(res, x, y, y_ctx, ctx); }
 GR_INLINE WARN_UNUSED_RESULT int gr_cmpabs_other(int * res, gr_srcptr x, gr_srcptr y, gr_ctx_t y_ctx, gr_ctx_t ctx) { return GR_BINARY_OP_OTHER_GET_INT(ctx, CMPABS_OTHER)(res, x, y, y_ctx, ctx); }
 
+#define __GR_CMP(cfun, expr) \
+    int cmp; \
+    if ((cfun)(&cmp, x, y, ctx) != GR_SUCCESS) \
+        return T_UNKNOWN; \
+    return (expr) ? T_TRUE : T_FALSE; \
+
+GR_INLINE WARN_UNUSED_RESULT truth_t gr_le(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { __GR_CMP(gr_cmp, cmp <= 0) }
+GR_INLINE WARN_UNUSED_RESULT truth_t gr_lt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { __GR_CMP(gr_cmp, cmp < 0) }
+GR_INLINE WARN_UNUSED_RESULT truth_t gr_ge(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { __GR_CMP(gr_cmp, cmp >= 0) }
+GR_INLINE WARN_UNUSED_RESULT truth_t gr_gt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { __GR_CMP(gr_cmp, cmp > 0) }
+
+GR_INLINE WARN_UNUSED_RESULT truth_t gr_abs_le(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { __GR_CMP(gr_cmpabs, cmp <= 0) }
+GR_INLINE WARN_UNUSED_RESULT truth_t gr_abs_lt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { __GR_CMP(gr_cmpabs, cmp < 0) }
+GR_INLINE WARN_UNUSED_RESULT truth_t gr_abs_ge(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { __GR_CMP(gr_cmpabs, cmp >= 0) }
+GR_INLINE WARN_UNUSED_RESULT truth_t gr_abs_gt(gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { __GR_CMP(gr_cmpabs, cmp > 0) }
+
+GR_INLINE WARN_UNUSED_RESULT int gr_min(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_OP(ctx, MIN)(res, x, y, ctx); }
+GR_INLINE WARN_UNUSED_RESULT int gr_max(gr_ptr res, gr_srcptr x, gr_srcptr y, gr_ctx_t ctx) { return GR_BINARY_OP(ctx, MAX)(res, x, y, ctx); }
+
 GR_INLINE WARN_UNUSED_RESULT int gr_gen(gr_ptr res, gr_ctx_t ctx) { return GR_CONSTANT_OP(ctx, GEN)(res, ctx); }
 GR_INLINE WARN_UNUSED_RESULT int gr_gens(gr_vec_t res, gr_ctx_t ctx) { return GR_VEC_CTX_OP(ctx, GENS)(res, ctx); }
 GR_INLINE WARN_UNUSED_RESULT int gr_gens_recursive(gr_vec_t res, gr_ctx_t ctx) { return GR_VEC_CTX_OP(ctx, GENS_RECURSIVE)(res, ctx); }

src/gr/acb.c

@@ -1004,6 +1004,32 @@ _gr_acb_cmpabs(int * res, const acb_t x, const acb_t y, const gr_ctx_t ctx)
     return _gr_acb_cmp(res, t, u, ctx);
 }
 
+int
+_gr_acb_min(acb_t res, const acb_t x, const acb_t y, const gr_ctx_t ctx)
+{
+    if (arb_is_zero(acb_imagref(x)) && arb_is_zero(acb_imagref(y)))
+    {
+        arb_min(acb_realref(res), acb_realref(x), acb_realref(y), ACB_CTX_PREC(ctx));
+        arb_zero(acb_imagref(res));
+        return GR_SUCCESS;
+    }
+    else
+        return GR_UNABLE;
+}
+
+int
+_gr_acb_max(acb_t res, const acb_t x, const acb_t y, const gr_ctx_t ctx)
+{
+    if (arb_is_zero(acb_imagref(x)) && arb_is_zero(acb_imagref(y)))
+    {
+        arb_max(acb_realref(res), acb_realref(x), acb_realref(y), ACB_CTX_PREC(ctx));
+        arb_zero(acb_imagref(res));
+        return GR_SUCCESS;
+    }
+    else
+        return GR_UNABLE;
+}
+
 int
 _gr_acb_pi(acb_t res, const gr_ctx_t ctx)
 {
@@ -1330,13 +1356,7 @@ _gr_acb_gamma_fmpq(acb_t res, const fmpq_t x, const gr_ctx_t ctx)
     }
 }
 
-
-int
-_gr_acb_rgamma(acb_t res, const acb_t x, const gr_ctx_t ctx)
-{
-    acb_rgamma(res, x, ACB_CTX_PREC(ctx));
-    return GR_SUCCESS;
-}
+DEF_FUNC(rgamma)
 
 int
 _gr_acb_lgamma(acb_t res, const acb_t x, const gr_ctx_t ctx)
@@ -1642,12 +1662,7 @@ int _gr_acb_stieltjes(acb_t res, const fmpz_t n, const acb_t a, const gr_ctx_t c
     return acb_is_finite(res) ? GR_SUCCESS : GR_UNABLE;
 }
 
-int
-_gr_acb_dirichlet_eta(acb_t res, const acb_t x, const gr_ctx_t ctx)
-{
-    acb_dirichlet_eta(res, x, ACB_CTX_PREC(ctx));
-    return GR_SUCCESS;
-}
+DEF_FUNC(dirichlet_eta)
 
 /* todo
 int
@@ -2220,6 +2235,8 @@ gr_method_tab_input _acb_methods_input[] =
     {GR_METHOD_ARG,             (gr_funcptr) _gr_acb_arg},
     {GR_METHOD_CMP,             (gr_funcptr) _gr_acb_cmp},
     {GR_METHOD_CMPABS,          (gr_funcptr) _gr_acb_cmpabs},
+    {GR_METHOD_MIN,             (gr_funcptr) _gr_acb_min},
+    {GR_METHOD_MAX,             (gr_funcptr) _gr_acb_max},
     {GR_METHOD_PI,              (gr_funcptr) _gr_acb_pi},
     {GR_METHOD_EXP,             (gr_funcptr) _gr_acb_exp},
     {GR_METHOD_EXPM1,           (gr_funcptr) _gr_acb_expm1},