1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
/*
Copyright (C) 2010 Fredrik Johansson
Copyright (C) 2020 William Hart
Copyright (C) 2020 Daniel Schultz
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 "nmod.h"
#include "nmod_mat.h"
#include "thread_support.h"
#if FLINT_USES_BLAS
# include <cblas.h>
# include "longlong.h"
#endif
void
nmod_mat_mul(nmod_mat_t C, const nmod_mat_t A, const nmod_mat_t B)
{
slong m = A->r;
slong k = A->c;
slong n = B->c;
slong min_dim = FLINT_MIN(FLINT_MIN(m, k), n);
slong cutoff;
slong flint_num_threads = flint_get_num_threads();
FLINT_ASSERT(C->r == A->r);
FLINT_ASSERT(C->c == B->c);
FLINT_ASSERT(A->c == B->r);
#if FLINT_USES_BLAS
/*
tuning is based on several assumptions:
(1) blas_num_threads >= flint_num_threads.
(2) nmod_mat_mul_blas (with crt) only beats
nmod_mat_mul_classical on square multiplications
of large enough dimension
(3) if nmod_mat_mul_blas beats nmod_mat_mul_classical on
square multiplications of size d, then it beats it on
rectangular multiplications as long as all dimensions are >= d
*/
if (FLINT_BITS == 64 && min_dim > 100)
{
flint_bitcnt_t bits = FLINT_BIT_COUNT(A->mod.n);
if (FLINT_BIT_COUNT(k) + 2*bits < 53 + 5)
{
/* mul_blas definitely avoids the slow crt */
cutoff = 100;
}
else if (flint_num_threads > 1)
{
/* mul_blas with crt is competing against mul_classical_threaded */
bits = FLINT_MAX(bits, 32);
cutoff = 100 + 5*flint_num_threads*bits/2;
}
else
{
/* mul_blas with crt is competing against mul_strassen */
cutoff = 450;
}
if (min_dim > cutoff && nmod_mat_mul_blas(C, A, B))
return;
}
#endif
if (C == A || C == B)
{
nmod_mat_t T;
nmod_mat_init(T, m, n, A->mod.n);
nmod_mat_mul(T, A, B);
nmod_mat_swap_entrywise(C, T);
nmod_mat_clear(T);
return;
}
if (FLINT_BITS == 64 && C->mod.n < 2048)
cutoff = 400;
else
cutoff = 200;
if (flint_num_threads > 1)
nmod_mat_mul_classical_threaded(C, A, B);
else if (min_dim < cutoff)
nmod_mat_mul_classical(C, A, B);
else
nmod_mat_mul_strassen(C, A, B);
}