implement acb_mat_mul_reorder

acb_mat.h

@@ -175,11 +175,10 @@ void acb_mat_add(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slon
 
 void acb_mat_sub(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec);
 
-void acb_mat_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec);
-
 void acb_mat_mul_classical(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec);
-
 void acb_mat_mul_threaded(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec);
+void acb_mat_mul_reorder(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec);
+void acb_mat_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec);
 
 void acb_mat_mul_entrywise(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec);
 

acb_mat/mul.c

@@ -11,14 +11,79 @@
 
 #include "acb_mat.h"
 
+static slong
+acb_mat_bits(const acb_mat_t A)
+{
+    slong b, t, i, ar, ac;
+
+    ar = acb_mat_nrows(A);
+    ac = acb_mat_ncols(A);
+
+    b = 0;
+    for (i = 0; i < ar; i++)
+    {
+        t = _arb_vec_bits((arb_srcptr) A->rows[i], 2 * ac);
+        b = FLINT_MAX(b, t);
+    }
+
+    return b;
+}
+
+int acb_mat_is_lagom(const acb_mat_t A)
+{
+    slong i, j, M, N;
+
+    M = acb_mat_nrows(A);
+    N = acb_mat_ncols(A);
+
+    for (i = 0; i < M; i++)
+    {
+        for (j = 0; j < N; j++)
+        {
+            if (!ARB_IS_LAGOM(acb_realref(acb_mat_entry(A, i, j))) ||
+                !ARB_IS_LAGOM(acb_imagref(acb_mat_entry(A, i, j))))
+                return 0;
+        }
+    }
+
+    return 1;
+}
+
 void
 acb_mat_mul(acb_mat_t C, const acb_mat_t A, const acb_mat_t B, slong prec)
 {
+    slong ar, ac, br, bc, n, abits, bbits, bits;
+
+    ar = acb_mat_nrows(A);
+    ac = acb_mat_ncols(A);
+    br = acb_mat_nrows(B);
+    bc = acb_mat_ncols(B);
+
+    if (ac != br || ar != acb_mat_nrows(C) || bc != acb_mat_ncols(C))
+    {
+        flint_printf("acb_mat_mul: incompatible dimensions\n");
+        flint_abort();
+    }
+
+    n = FLINT_MIN(ar, ac);
+    n = FLINT_MIN(ac, bc);
+
+    if (n >= 5)
+    {
+        abits = acb_mat_bits(A);
+        bbits = acb_mat_bits(B);
+
+        bits = FLINT_MIN(prec, FLINT_MAX(abits, bbits));
+
+        if (bits < 8000 && n >= 5 + bits / 64)
+        {
+            acb_mat_mul_reorder(C, A, B, prec);
+            return;
+        }
+    }
+
     if (flint_get_num_threads() > 1 &&
-        ((double) acb_mat_nrows(A) *
-         (double) acb_mat_nrows(B) *
-         (double) acb_mat_ncols(B) *
-         (double) prec > 100000))
+        ((double) ar * (double) ac * (double) bc * (double) prec > 100000))
     {
         acb_mat_mul_threaded(C, A, B, prec);
     }

acb_mat/mul_reorder.c

@@ -0,0 +1,286 @@
+/*
+    Copyright (C) 2018 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb 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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_mat.h"
+
+static void
+copy_re_shallow(arb_mat_t X, const acb_mat_t A)
+{
+    slong M, N, i, j;
+    M = arb_mat_nrows(X);
+    N = arb_mat_ncols(X);
+    for (i = 0; i < M; i++)
+        for (j = 0; j < N; j++)
+            *arb_mat_entry(X, i, j) = *acb_realref(acb_mat_entry(A, i, j));
+}
+
+static void
+copy_im_shallow(arb_mat_t X, const acb_mat_t A)
+{
+    slong M, N, i, j;
+    M = arb_mat_nrows(X);
+    N = arb_mat_ncols(X);
+    for (i = 0; i < M; i++)
+        for (j = 0; j < N; j++)
+            *arb_mat_entry(X, i, j) = *acb_imagref(acb_mat_entry(A, i, j));
+}
+
+static void
+clear_shallow(arb_mat_t X)
+{
+    slong M, N, i, j;
+    M = arb_mat_nrows(X);
+    N = arb_mat_ncols(X);
+    for (i = 0; i < M; i++)
+        for (j = 0; j < N; j++)
+            arb_init(arb_mat_entry(X, i, j));
+}
+
+/* todo: only use if prec (also look at bits!) < 3200 (except of course: pure real/...) */
+
+void
+acb_mat_mul_reorder(acb_mat_t C, const acb_mat_t A, const acb_mat_t B, slong prec)
+{
+    arb_mat_t X, Y, Z, W;
+    slong M, N, P;
+    slong i, j;
+
+    M = acb_mat_nrows(A);
+    N = acb_mat_ncols(A);
+    P = acb_mat_ncols(B);
+
+    if (acb_mat_is_real(A))
+    {
+        if (acb_mat_is_real(B))
+        {
+            arb_mat_init(X, M, N);
+            arb_mat_init(Y, N, P);
+            arb_mat_init(Z, M, P);
+
+            copy_re_shallow(X, A);
+            copy_re_shallow(Y, B);
+
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    arb_zero(acb_imagref(acb_mat_entry(C, i, j)));
+
+            if (A == C || B == C)
+            {
+                arb_mat_mul(Z, X, Y, prec);
+
+                for (i = 0; i < M; i++)
+                    for (j = 0; j < P; j++)
+                        arb_swap(acb_realref(acb_mat_entry(C, i, j)), acb_mat_entry(Z, i, j));
+            }
+            else
+            {
+                copy_re_shallow(Z, C);
+                arb_mat_mul(Z, X, Y, prec);
+
+                for (i = 0; i < M; i++)
+                    for (j = 0; j < P; j++)
+                        *acb_realref(acb_mat_entry(C, i, j)) = *arb_mat_entry(Z, i, j);
+
+                clear_shallow(Z);
+            }
+
+            clear_shallow(X);
+            clear_shallow(Y);
+
+            arb_mat_clear(X);
+            arb_mat_clear(Y);
+            arb_mat_clear(Z);
+        }
+        else
+        {
+            arb_mat_init(X, M, N);
+            arb_mat_init(Y, N, P);
+            arb_mat_init(Z, M, P);
+
+            /* (reA*reB, reA*imB) */
+
+            copy_re_shallow(X, A);
+            copy_re_shallow(Y, B);
+
+            if (A == C || B == C)
+            {
+                arb_mat_t T;
+                arb_mat_init(T, M, P);
+
+                arb_mat_mul(T, X, Y, prec);
+                copy_im_shallow(Y, B);
+                arb_mat_mul(Z, X, Y, prec);
+
+                for (i = 0; i < M; i++)
+                    for (j = 0; j < P; j++)
+                        arb_swap(acb_realref(acb_mat_entry(C, i, j)), acb_mat_entry(T, i, j));
+
+                for (i = 0; i < M; i++)
+                    for (j = 0; j < P; j++)
+                        arb_swap(acb_imagref(acb_mat_entry(C, i, j)), acb_mat_entry(Z, i, j));
+
+                arb_mat_clear(T);
+            }
+            else
+            {
+                copy_re_shallow(Z, C);
+                arb_mat_mul(Z, X, Y, prec);
+
+                for (i = 0; i < M; i++)
+                    for (j = 0; j < P; j++)
+                        *acb_realref(acb_mat_entry(C, i, j)) = *arb_mat_entry(Z, i, j);
+
+                copy_im_shallow(Z, C);
+                copy_im_shallow(Y, B);
+                arb_mat_mul(Z, X, Y, prec);
+
+                for (i = 0; i < M; i++)
+                    for (j = 0; j < P; j++)
+                        *acb_imagref(acb_mat_entry(C, i, j)) = *arb_mat_entry(Z, i, j);
+
+                clear_shallow(Z);
+            }
+
+            clear_shallow(X);
+            clear_shallow(Y);
+
+            arb_mat_clear(X);
+            arb_mat_clear(Y);
+            arb_mat_clear(Z);
+        }
+    }
+    else if (acb_mat_is_real(B))
+    {
+        arb_mat_init(X, M, N);
+        arb_mat_init(Y, N, P);
+        arb_mat_init(Z, M, P);
+
+        /* (reA*reB, imA*reB) */
+
+        copy_re_shallow(X, A);
+        copy_re_shallow(Y, B);
+
+        if (A == C || B == C)
+        {
+            arb_mat_t T;
+            arb_mat_init(T, M, P);
+
+            arb_mat_mul(T, X, Y, prec);
+            copy_im_shallow(X, A);
+            arb_mat_mul(Z, X, Y, prec);
+
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    arb_swap(acb_realref(acb_mat_entry(C, i, j)), acb_mat_entry(T, i, j));
+
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    arb_swap(acb_imagref(acb_mat_entry(C, i, j)), acb_mat_entry(Z, i, j));
+
+            arb_mat_clear(T);
+        }
+        else
+        {
+            copy_re_shallow(Z, C);
+            arb_mat_mul(Z, X, Y, prec);
+
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    *acb_realref(acb_mat_entry(C, i, j)) = *arb_mat_entry(Z, i, j);
+
+            copy_im_shallow(Z, C);
+            copy_im_shallow(X, A);
+            arb_mat_mul(Z, X, Y, prec);
+
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    *acb_imagref(acb_mat_entry(C, i, j)) = *arb_mat_entry(Z, i, j);
+
+            clear_shallow(Z);
+        }
+
+        clear_shallow(X);
+        clear_shallow(Y);
+
+        arb_mat_clear(X);
+        arb_mat_clear(Y);
+        arb_mat_clear(Z);
+    }
+    else
+    {
+        arb_mat_init(X, M, N);
+        arb_mat_init(Y, N, P);
+        arb_mat_init(Z, M, P);
+        arb_mat_init(W, M, P);
+
+        copy_re_shallow(X, A);
+        copy_re_shallow(Y, B);
+        arb_mat_mul(Z, X, Y, prec);
+
+        copy_im_shallow(X, A);
+        copy_im_shallow(Y, B);
+        arb_mat_mul(W, X, Y, prec);
+
+        if (A == C || B == C)
+        {
+            arb_mat_t T;
+            arb_mat_init(T, M, P);
+            arb_mat_sub(T, Z, W, prec);
+
+            copy_re_shallow(X, A);
+            arb_mat_mul(Z, X, Y, prec);
+
+            copy_im_shallow(X, A);
+            copy_re_shallow(Y, B);
+            arb_mat_mul(W, X, Y, prec);
+
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    arb_swap(acb_realref(acb_mat_entry(C, i, j)),
+                        arb_mat_entry(T, i, j));
+
+            arb_mat_clear(T);
+
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    arb_add(acb_imagref(acb_mat_entry(C, i, j)),
+                        arb_mat_entry(Z, i, j), arb_mat_entry(W, i, j), prec);
+        }
+        else
+        {
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    arb_sub(acb_realref(acb_mat_entry(C, i, j)),
+                        arb_mat_entry(Z, i, j), arb_mat_entry(W, i, j), prec);
+
+            copy_re_shallow(X, A);
+            arb_mat_mul(Z, X, Y, prec);
+
+            copy_im_shallow(X, A);
+            copy_re_shallow(Y, B);
+            arb_mat_mul(W, X, Y, prec);
+
+            for (i = 0; i < M; i++)
+                for (j = 0; j < P; j++)
+                    arb_add(acb_imagref(acb_mat_entry(C, i, j)),
+                        arb_mat_entry(Z, i, j), arb_mat_entry(W, i, j), prec);
+        }
+
+        clear_shallow(X);
+        clear_shallow(Y);
+
+        arb_mat_clear(X);
+        arb_mat_clear(Y);
+        arb_mat_clear(Z);
+        arb_mat_clear(W);
+    }
+}
+

acb_mat/randtest.c

@@ -14,15 +14,23 @@
 void
 acb_mat_randtest(acb_mat_t mat, flint_rand_t state, slong prec, slong mag_bits)
 {
-    slong i, j;
+    slong i, j, density;
+
+    density = n_randint(state, 100);
 
     if (n_randint(state, 2))
         for (i = 0; i < acb_mat_nrows(mat); i++)
             for (j = 0; j < acb_mat_ncols(mat); j++)
-                acb_randtest(acb_mat_entry(mat, i, j), state, prec, mag_bits);
+                if (n_randint(state, 100) < density)
+                    acb_randtest(acb_mat_entry(mat, i, j), state, prec, mag_bits);
+                else
+                    acb_zero(acb_mat_entry(mat, i, j));
     else
         for (i = 0; i < acb_mat_nrows(mat); i++)
             for (j = 0; j < acb_mat_ncols(mat); j++)
-                acb_randtest_precise(acb_mat_entry(mat, i, j), state, prec, mag_bits);
+                if (n_randint(state, 100) < density)
+                    acb_randtest_precise(acb_mat_entry(mat, i, j), state, prec, mag_bits);
+                else
+                    acb_zero(acb_mat_entry(mat, i, j));
 }
 

acb_mat/test/t-mul_reorder.c

@@ -0,0 +1,374 @@
+/*
+    Copyright (C) 2012, 2018 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb 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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_mat.h"
+
+void
+_acb_mat_init_randtest(acb_mat_t mat, slong r, slong c, flint_rand_t state)
+{
+    acb_mat_init(mat, r, c);
+    acb_mat_randtest(mat, state, 2 + n_randint(state, 200), 10);
+}
+
+void
+_acb_mat_nprintd(const char * name, acb_mat_t mat)
+{
+    flint_printf("%s = ", name);
+    acb_mat_printd(mat, 15);
+    flint_printf("\n\n");
+}
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("mul_reorder....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        slong m, n, k, qbits1, qbits2, rbits1, rbits2, rbits3;
+        fmpq_mat_t A, B, C;
+        acb_mat_t a, b, c, d;
+
+        qbits1 = 2 + n_randint(state, 200);
+        qbits2 = 2 + n_randint(state, 200);
+        rbits1 = 2 + n_randint(state, 200);
+        rbits2 = 2 + n_randint(state, 200);
+        rbits3 = 2 + n_randint(state, 200);
+
+        m = n_randint(state, 8);
+        n = n_randint(state, 8);
+        k = n_randint(state, 8);
+
+        fmpq_mat_init(A, m, n);
+        fmpq_mat_init(B, n, k);
+        fmpq_mat_init(C, m, k);
+
+        acb_mat_init(a, m, n);
+        acb_mat_init(b, n, k);
+        _acb_mat_init_randtest(c, m, k, state);
+        _acb_mat_init_randtest(d, m, k, state);
+
+        fmpq_mat_randtest(A, state, qbits1);
+        fmpq_mat_randtest(B, state, qbits2);
+        fmpq_mat_mul(C, A, B);
+
+        acb_mat_set_fmpq_mat(a, A, rbits1);
+        acb_mat_set_fmpq_mat(b, B, rbits2);
+        acb_mat_mul_reorder(c, a, b, rbits3);
+
+        if (!acb_mat_contains_fmpq_mat(c, C))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("m = %wd, n = %wd, k = %wd, bits3 = %wd\n", m, n, k, rbits3);
+
+            flint_printf("A = "); fmpq_mat_print(A); flint_printf("\n\n");
+            flint_printf("B = "); fmpq_mat_print(B); flint_printf("\n\n");
+            flint_printf("C = "); fmpq_mat_print(C); flint_printf("\n\n");
+
+            flint_printf("a = "); acb_mat_printd(a, 15); flint_printf("\n\n");
+            flint_printf("b = "); acb_mat_printd(b, 15); flint_printf("\n\n");
+            flint_printf("c = "); acb_mat_printd(c, 15); flint_printf("\n\n");
+
+            flint_abort();
+        }
+
+        /* test aliasing with a */
+        if (acb_mat_nrows(a) == acb_mat_nrows(c) &&
+            acb_mat_ncols(a) == acb_mat_ncols(c))
+        {
+            acb_mat_set(d, a);
+            acb_mat_mul_reorder(d, d, b, rbits3);
+            if (!acb_mat_equal(d, c))
+            {
+                flint_printf("FAIL (aliasing 1)\n\n");
+                flint_abort();
+            }
+        }
+
+        /* test aliasing with b */
+        if (acb_mat_nrows(b) == acb_mat_nrows(c) &&
+            acb_mat_ncols(b) == acb_mat_ncols(c))
+        {
+            acb_mat_set(d, b);
+            acb_mat_mul_reorder(d, a, d, rbits3);
+            if (!acb_mat_equal(d, c))
+            {
+                flint_printf("FAIL (aliasing 2)\n\n");
+                flint_abort();
+            }
+        }
+
+        fmpq_mat_clear(A);
+        fmpq_mat_clear(B);
+        fmpq_mat_clear(C);
+
+        acb_mat_clear(a);
+        acb_mat_clear(b);
+        acb_mat_clear(c);
+        acb_mat_clear(d);
+    }
+
+    /* general aliasing test */
+    for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+    {
+        slong m, n;
+        slong rbits;
+        acb_mat_t a, b, c, d;
+
+        rbits = 2 + n_randint(state, 200);
+
+        m = n_randint(state, 8);
+        n = n_randint(state, 8);
+
+        _acb_mat_init_randtest(a, m, n, state);
+        _acb_mat_init_randtest(b, n, n, state);
+        _acb_mat_init_randtest(c, m, n, state);
+        _acb_mat_init_randtest(d, m, n, state);
+
+        acb_mat_mul_reorder(c, a, b, rbits);
+        acb_mat_set(d, a);
+        acb_mat_mul_reorder(d, d, b, rbits);
+
+        if (!acb_mat_equal(c, d))
+        {
+            flint_printf("FAIL (aliasing 3)\n\n");
+            _acb_mat_nprintd("a", a);
+            _acb_mat_nprintd("b", b);
+            _acb_mat_nprintd("c", d);
+            _acb_mat_nprintd("d", d);
+            flint_abort();
+        }
+
+        acb_mat_clear(a);
+        acb_mat_clear(b);
+        acb_mat_clear(c);
+        acb_mat_clear(d);
+
+        _acb_mat_init_randtest(a, m, m, state);
+        _acb_mat_init_randtest(b, m, n, state);
+        _acb_mat_init_randtest(c, m, n, state);
+        _acb_mat_init_randtest(d, m, n, state);
+
+        acb_mat_mul_reorder(c, a, b, rbits);
+        acb_mat_set(d, b);
+        acb_mat_mul_reorder(d, a, d, rbits);
+
+        if (!acb_mat_equal(c, d))
+        {
+            flint_printf("FAIL (aliasing 4)\n\n");
+            _acb_mat_nprintd("a", a);
+            _acb_mat_nprintd("b", b);
+            _acb_mat_nprintd("c", d);
+            _acb_mat_nprintd("d", d);
+            flint_abort();
+        }
+
+        acb_mat_clear(a);
+        acb_mat_clear(b);
+        acb_mat_clear(c);
+        acb_mat_clear(d);
+    }
+
+    /* check algebraic properties like associativity and distributivity */
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        slong m, n, k, l;
+        slong rbits;
+        acb_mat_t a, b, c, d, ab, ac, bd, cd, s;
+
+        rbits = 2 + n_randint(state, 200);
+
+        m = n_randint(state, 8);
+        n = n_randint(state, 8);
+        k = n_randint(state, 8);
+        l = n_randint(state, 8);
+
+        _acb_mat_init_randtest(a, m, n, state);
+        _acb_mat_init_randtest(b, n, k, state);
+        _acb_mat_init_randtest(c, n, k, state);
+        _acb_mat_init_randtest(d, k, l, state);
+
+        _acb_mat_init_randtest(ab, m, k, state);
+        _acb_mat_init_randtest(ac, m, k, state);
+        _acb_mat_init_randtest(bd, n, l, state);
+        _acb_mat_init_randtest(cd, n, l, state);
+        _acb_mat_init_randtest(s, n, k, state);
+
+        acb_mat_mul_reorder(ab, a, b, rbits);
+        acb_mat_mul_reorder(ac, a, c, rbits);
+        acb_mat_mul_reorder(bd, b, d, rbits);
+        acb_mat_mul_reorder(cd, c, d, rbits);
+        acb_mat_add(s, b, c, rbits);
+
+        /* check associativity of multiplication */
+        /* (A*B)*D = A*(B*D) */
+        {
+            acb_mat_t lhs, rhs;
+
+            _acb_mat_init_randtest(lhs, m, l, state);
+            _acb_mat_init_randtest(rhs, m, l, state);
+
+            acb_mat_mul_reorder(lhs, ab, d, rbits);
+            acb_mat_mul_reorder(rhs, a, bd, rbits);
+
+            if (!acb_mat_overlaps(lhs, rhs))
+            {
+                flint_printf("FAIL\n\n");
+                flint_printf("m, n, k, l = %wd, %wd, %wd, %wd\n", m, n, k, l);
+                flint_printf("rbits = %wd\n", rbits);
+
+                _acb_mat_nprintd("a", a);
+                _acb_mat_nprintd("b", b);
+                _acb_mat_nprintd("d", d);
+                _acb_mat_nprintd("(a*b)*d", lhs);
+                _acb_mat_nprintd("a*(b*d)", rhs);
+
+                flint_abort();
+            }
+
+            acb_mat_clear(lhs);
+            acb_mat_clear(rhs);
+        }
+
+        /* check left distributivity of multiplication over addition */
+        /* A*(B + C) = A*B + A*C */
+        {
+            acb_mat_t lhs, rhs;
+
+            _acb_mat_init_randtest(lhs, m, k, state);
+            _acb_mat_init_randtest(rhs, m, k, state);
+
+            acb_mat_mul_reorder(lhs, a, s, rbits);
+            acb_mat_add(rhs, ab, ac, rbits);
+
+            if (!acb_mat_overlaps(lhs, rhs))
+            {
+                flint_printf("FAIL\n\n");
+                flint_printf("m, n, k, l = %wd, %wd, %wd, %wd\n", m, n, k, l);
+                flint_printf("rbits = %wd\n", rbits);
+
+                _acb_mat_nprintd("a", a);
+                _acb_mat_nprintd("b", b);
+                _acb_mat_nprintd("c", c);
+                _acb_mat_nprintd("a*(b + c)", lhs);
+                _acb_mat_nprintd("a*b + b*c", rhs);
+
+                flint_abort();
+            }
+
+            acb_mat_clear(lhs);
+            acb_mat_clear(rhs);
+        }
+
+        /* check right distributivity of multiplication over addition */
+        /* (B + C)*D = B*D + C*D */
+        {
+            acb_mat_t lhs, rhs;
+
+            _acb_mat_init_randtest(lhs, n, l, state);
+            _acb_mat_init_randtest(rhs, n, l, state);
+
+            acb_mat_mul_reorder(lhs, s, d, rbits);
+            acb_mat_add(rhs, bd, cd, rbits);
+
+            if (!acb_mat_overlaps(lhs, rhs))
+            {
+                flint_printf("FAIL\n\n");
+                flint_printf("m, n, k, l = %wd, %wd, %wd, %wd\n", m, n, k, l);
+                flint_printf("rbits = %wd\n", rbits);
+
+                _acb_mat_nprintd("b", b);
+                _acb_mat_nprintd("c", c);
+                _acb_mat_nprintd("d", d);
+                _acb_mat_nprintd("(b + c)*d", lhs);
+                _acb_mat_nprintd("b*d + c*d", rhs);
+
+                flint_abort();
+            }
+
+            acb_mat_clear(lhs);
+            acb_mat_clear(rhs);
+        }
+
+        /* check left multiplicative identity I*D = D */
+        {
+            acb_mat_t one, lhs;
+
+            _acb_mat_init_randtest(one, k, k, state);
+            _acb_mat_init_randtest(lhs, k, l, state);
+
+            acb_mat_one(one);
+            acb_mat_mul_reorder(lhs, one, d, rbits);
+
+            if (!acb_mat_contains(lhs, d))
+            {
+                flint_printf("FAIL\n\n");
+                flint_printf("k = %wd, l = %wd\n", k, l);
+                flint_printf("rbits = %wd\n", rbits);
+
+                _acb_mat_nprintd("identity * d", lhs);
+                _acb_mat_nprintd("d", d);
+
+                flint_abort();
+            }
+
+            acb_mat_clear(one);
+            acb_mat_clear(lhs);
+        }
+
+        /* check right multiplicative identity A*I = A */
+        {
+            acb_mat_t one, lhs;
+
+            _acb_mat_init_randtest(one, n, n, state);
+            _acb_mat_init_randtest(lhs, m, n, state);
+
+            acb_mat_one(one);
+            acb_mat_mul_reorder(lhs, a, one, rbits);
+
+            if (!acb_mat_contains(lhs, a))
+            {
+                flint_printf("FAIL\n\n");
+                flint_printf("m = %wd, n = %wd\n", m, n);
+                flint_printf("rbits = %wd\n", rbits);
+
+                _acb_mat_nprintd("a * identity", lhs);
+                _acb_mat_nprintd("a", a);
+
+                flint_abort();
+            }
+
+            acb_mat_clear(one);
+            acb_mat_clear(lhs);
+        }
+
+        acb_mat_clear(a);
+        acb_mat_clear(b);
+        acb_mat_clear(c);
+        acb_mat_clear(d);
+        acb_mat_clear(ab);
+        acb_mat_clear(ac);
+        acb_mat_clear(bd);
+        acb_mat_clear(cd);
+        acb_mat_clear(s);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/acb_mat.rst

@@ -193,20 +193,26 @@ Arithmetic
     Sets *res* to the difference of *mat1* and *mat2*. The operands must have
     the same dimensions.
 
-.. function:: void acb_mat_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
-
 .. function:: void acb_mat_mul_classical(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
 
 .. function:: void acb_mat_mul_threaded(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
 
+.. function:: void acb_mat_mul_reorder(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
+
+.. function:: void acb_mat_mul(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)
+
     Sets *res* to the matrix product of *mat1* and *mat2*. The operands must have
     compatible dimensions for matrix multiplication.
 
-    The *threaded* version splits the computation
-    over the number of threads returned by *flint_get_num_threads()*.
-    The default version automatically calls the *threaded* version
-    if the matrices are sufficiently large and more than one thread
-    can be used.
+    The *classical* version performs matrix multiplication in the trivial way.
+
+    The *threaded* version performs classical multiplication but splits the
+    computation over the number of threads returned by *flint_get_num_threads()*.
+
+    The *reorder* version reorders the data and performs one to four real
+    matrix multiplications via :func:`arb_mat_mul`.
+
+    The default version chooses an algorithm automatically.
 
 .. function:: void acb_mat_mul_entrywise(acb_mat_t res, const acb_mat_t mat1, const acb_mat_t mat2, slong prec)