ENH: arb_mat_mul_entrywise

This commit is contained in:
alex 2016-02-23 17:59:42 -05:00
parent f29fae67f6
commit 1e3bb266d0
4 changed files with 198 additions and 0 deletions

View file

@ -144,6 +144,8 @@ void arb_mat_mul_classical(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, sl
void arb_mat_mul_threaded(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec); void arb_mat_mul_threaded(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec);
void arb_mat_mul_entrywise(arb_mat_t res, const arb_mat_t mat1, const arb_mat_t mat2, slong prec);
void arb_mat_sqr(arb_mat_t B, const arb_mat_t A, slong prec); void arb_mat_sqr(arb_mat_t B, const arb_mat_t A, slong prec);
void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec); void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec);

49
arb_mat/mul_entrywise.c Normal file
View file

@ -0,0 +1,49 @@
/*=============================================================================
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 ARB; 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_mul_entrywise(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)
{
slong i, j;
if (arb_mat_nrows(A) != arb_mat_nrows(B) ||
arb_mat_ncols(A) != arb_mat_ncols(B))
{
flint_printf("arb_mat_mul_entrywise: incompatible dimensions\n");
abort();
}
for (i = 0; i < arb_mat_nrows(A); i++)
{
for (j = 0; j < arb_mat_ncols(A); j++)
{
arb_mul(arb_mat_entry(C, i, j),
arb_mat_entry(A, i, j),
arb_mat_entry(B, i, j), prec);
}
}
}

View file

@ -0,0 +1,142 @@
/*=============================================================================
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 ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2016 Arb authors
******************************************************************************/
#include "arb_mat.h"
void
_fmpq_mat_mul_entrywise(fmpq_mat_t C, const fmpq_mat_t A, const fmpq_mat_t B)
{
slong i, j;
for (i = 0; i < fmpq_mat_nrows(A); i++)
{
for (j = 0; j < fmpq_mat_ncols(A); j++)
{
fmpq_mul(fmpq_mat_entry(C, i, j),
fmpq_mat_entry(A, i, j),
fmpq_mat_entry(B, i, j));
}
}
}
int main()
{
slong iter;
flint_rand_t state;
flint_printf("mul_entrywise....");
fflush(stdout);
flint_randinit(state);
for (iter = 0; iter < 10000; iter++)
{
slong m, n, qbits1, qbits2, rbits1, rbits2, rbits3;
fmpq_mat_t A, B, C;
arb_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, 10);
n = n_randint(state, 10);
fmpq_mat_init(A, m, n);
fmpq_mat_init(B, m, n);
fmpq_mat_init(C, m, n);
arb_mat_init(a, m, n);
arb_mat_init(b, m, n);
arb_mat_init(c, m, n);
arb_mat_init(d, m, n);
fmpq_mat_randtest(A, state, qbits1);
fmpq_mat_randtest(B, state, qbits2);
_fmpq_mat_mul_entrywise(C, A, B);
arb_mat_set_fmpq_mat(a, A, rbits1);
arb_mat_set_fmpq_mat(b, B, rbits2);
arb_mat_mul_entrywise(c, a, b, rbits3);
if (!arb_mat_contains_fmpq_mat(c, C))
{
flint_printf("FAIL\n\n");
flint_printf("m = %wd, n = %wd, bits3 = %wd\n", m, n, 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 = "); arb_mat_printd(a, 15); flint_printf("\n\n");
flint_printf("b = "); arb_mat_printd(b, 15); flint_printf("\n\n");
flint_printf("c = "); arb_mat_printd(c, 15); flint_printf("\n\n");
abort();
}
/* test aliasing with a */
if (arb_mat_nrows(a) == arb_mat_nrows(c) &&
arb_mat_ncols(a) == arb_mat_ncols(c))
{
arb_mat_set(d, a);
arb_mat_mul_entrywise(d, d, b, rbits3);
if (!arb_mat_equal(d, c))
{
flint_printf("FAIL (aliasing 1)\n\n");
abort();
}
}
/* test aliasing with b */
if (arb_mat_nrows(b) == arb_mat_nrows(c) &&
arb_mat_ncols(b) == arb_mat_ncols(c))
{
arb_mat_set(d, b);
arb_mat_mul_entrywise(d, a, d, rbits3);
if (!arb_mat_equal(d, c))
{
flint_printf("FAIL (aliasing 2)\n\n");
abort();
}
}
fmpq_mat_clear(A);
fmpq_mat_clear(B);
fmpq_mat_clear(C);
arb_mat_clear(a);
arb_mat_clear(b);
arb_mat_clear(c);
arb_mat_clear(d);
}
flint_randclear(state);
flint_cleanup();
flint_printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -176,6 +176,11 @@ Arithmetic
if the matrices are sufficiently large and more than one thread if the matrices are sufficiently large and more than one thread
can be used. can be used.
.. function:: void arb_mat_mul_entrywise(arb_mat_t C, const arb_mat_t A, const arb_mat_t B, slong prec)
Sets *C* to the entrywise product of *A* and *B*.
The operands must have the same dimensions.
.. function:: void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec) .. function:: void arb_mat_sqr_classical(arb_mat_t B, const arb_mat_t A, slong prec)
.. function:: void arb_mat_sqr(arb_mat_t res, const arb_mat_t mat, slong prec) .. function:: void arb_mat_sqr(arb_mat_t res, const arb_mat_t mat, slong prec)