From 3e3127d4360fe1fa0583a65cd912f032b21fa97b Mon Sep 17 00:00:00 2001 From: alex Date: Wed, 6 Apr 2016 13:39:59 -0400 Subject: [PATCH] TST: add frobenius norm test files --- acb_mat/test/t-frobenius_norm.c | 278 ++++++++++++++++++++++++++++++++ arb_mat/test/t-frobenius_norm.c | 235 +++++++++++++++++++++++++++ 2 files changed, 513 insertions(+) create mode 100644 acb_mat/test/t-frobenius_norm.c create mode 100644 arb_mat/test/t-frobenius_norm.c diff --git a/acb_mat/test/t-frobenius_norm.c b/acb_mat/test/t-frobenius_norm.c new file mode 100644 index 00000000..197722fa --- /dev/null +++ b/acb_mat/test/t-frobenius_norm.c @@ -0,0 +1,278 @@ +/*============================================================================= + + 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 "acb_mat.h" + +static void +_acb_mat_set_fmpq_mat_fmpq_mat(acb_mat_t C, + const fmpq_mat_t rsrc, const fmpq_mat_t isrc, slong prec) +{ + slong i, j; + for (i = 0; i < acb_mat_nrows(C); i++) + { + for (j = 0; j < acb_mat_ncols(C); j++) + { + arb_set_fmpq(acb_realref(acb_mat_entry(C, i, j)), + fmpq_mat_entry(rsrc, i, j), prec); + arb_set_fmpq(acb_imagref(acb_mat_entry(C, i, j)), + fmpq_mat_entry(isrc, i, j), prec); + } + } +} + +static void +_acb_mat_conjugate_transpose(acb_mat_t B, const acb_mat_t A) +{ + slong i, j; + acb_mat_transpose(B, A); + for (i = 0; i < acb_mat_nrows(B); i++) + for (j = 0; j < acb_mat_ncols(B); j++) + acb_conj(acb_mat_entry(B, i, j), acb_mat_entry(B, i, j)); +} + +static void +_fmpq_mat_sum_of_squares(fmpq_t res, const fmpq_mat_t Q) +{ + slong i, j; + fmpq_zero(res); + for (i = 0; i < fmpq_mat_nrows(Q); i++) + { + for (j = 0; j < fmpq_mat_ncols(Q); j++) + { + fmpq_addmul(res, fmpq_mat_entry(Q, i, j), fmpq_mat_entry(Q, i, j)); + } + } +} + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("frobenius_norm...."); + fflush(stdout); + + flint_randinit(state); + + /* compare to the exact rational norm */ + for (iter = 0; iter < 10000; iter++) + { + fmpq_mat_t Qr, Qi; + fmpq_t q; + acb_mat_t A; + slong n, qbits, prec; + + n = n_randint(state, 8); + qbits = 1 + n_randint(state, 100); + prec = 2 + n_randint(state, 200); + + fmpq_mat_init(Qr, n, n); + fmpq_mat_init(Qi, n, n); + fmpq_init(q); + + acb_mat_init(A, n, n); + + fmpq_mat_randtest(Qr, state, qbits); + fmpq_mat_randtest(Qi, state, qbits); + + /* compute the square of the exact rational norm */ + { + fmpq_t qr, qi; + fmpq_init(qr); + fmpq_init(qi); + _fmpq_mat_sum_of_squares(qr, Qr); + _fmpq_mat_sum_of_squares(qi, Qi); + fmpq_add(q, qr, qi); + fmpq_clear(qr); + fmpq_clear(qi); + } + + _acb_mat_set_fmpq_mat_fmpq_mat(A, Qr, Qi, prec); + + /* check that the arb interval contains the exact value */ + { + arb_t a; + arb_init(a); + + acb_mat_frobenius_norm(a, A, prec); + arb_mul(a, a, a, prec); + + if (!arb_contains_fmpq(a, q)) + { + flint_printf("FAIL (containment, iter = %wd)\n", iter); + flint_printf("n = %wd, prec = %wd\n", n, prec); + flint_printf("\n"); + + flint_printf("Qr = \n"); fmpq_mat_print(Qr); flint_printf("\n\n"); + flint_printf("Qi = \n"); fmpq_mat_print(Qi); flint_printf("\n\n"); + flint_printf("frobenius_norm(Q)^2 = \n"); + fmpq_print(q); flint_printf("\n\n"); + + flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n"); + flint_printf("frobenius_norm(A)^2 = \n"); + arb_printd(a, 15); flint_printf("\n\n"); + flint_printf("frobenius_norm(A)^2 = \n"); + arb_print(a); flint_printf("\n\n"); + + abort(); + } + + arb_clear(a); + } + + /* check that the upper bound is not less than the exact value */ + { + mag_t b; + fmpq_t y; + + mag_init(b); + fmpq_init(y); + + acb_mat_bound_frobenius_norm(b, A); + mag_mul(b, b, b); + mag_get_fmpq(y, b); + + if (fmpq_cmp(q, y) > 0) + { + flint_printf("FAIL (bound, iter = %wd)\n", iter); + flint_printf("n = %wd, prec = %wd\n", n, prec); + flint_printf("\n"); + + flint_printf("Qr = \n"); fmpq_mat_print(Qr); flint_printf("\n\n"); + flint_printf("Qi = \n"); fmpq_mat_print(Qi); flint_printf("\n\n"); + flint_printf("frobenius_norm(Q)^2 = \n"); + fmpq_print(q); flint_printf("\n\n"); + + flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n"); + flint_printf("bound_frobenius_norm(A)^2 = \n"); + mag_printd(b, 15); flint_printf("\n\n"); + flint_printf("bound_frobenius_norm(A)^2 = \n"); + mag_print(b); flint_printf("\n\n"); + + abort(); + } + + mag_clear(b); + fmpq_clear(y); + } + + fmpq_mat_clear(Qr); + fmpq_mat_clear(Qi); + fmpq_clear(q); + acb_mat_clear(A); + } + + /* check trace(A^H A) = frobenius_norm(A)^2 */ + for (iter = 0; iter < 10000; iter++) + { + slong m, n, prec; + acb_mat_t A, AH, AHA; + acb_t t; + + prec = 2 + n_randint(state, 200); + + m = n_randint(state, 10); + n = n_randint(state, 10); + + acb_mat_init(A, m, n); + acb_mat_init(AH, n, m); + acb_mat_init(AHA, n, n); + acb_init(t); + + acb_mat_randtest(A, state, 2 + n_randint(state, 100), 10); + _acb_mat_conjugate_transpose(AH, A); + acb_mat_mul(AHA, AH, A, prec); + acb_mat_trace(t, AHA, prec); + acb_sqrt(t, t, prec); + + /* check the norm bound */ + { + mag_t low, frobenius; + + mag_init(low); + acb_get_mag_lower(low, t); + + mag_init(frobenius); + acb_mat_bound_frobenius_norm(frobenius, A); + + if (mag_cmp(low, frobenius) > 0) + { + flint_printf("FAIL (frobenius norm bound)\n", iter); + flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec); + flint_printf("\n"); + + flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n"); + + flint_printf("lower(sqrt(trace(A^H A))) = \n"); + mag_printd(low, 15); flint_printf("\n\n"); + + flint_printf("bound_frobenius_norm(A) = \n"); + mag_printd(frobenius, 15); flint_printf("\n\n"); + + abort(); + } + + mag_clear(low); + mag_clear(frobenius); + } + + /* check the norm interval */ + { + arb_t frobenius; + + arb_init(frobenius); + acb_mat_frobenius_norm(frobenius, A, prec); + + if (!arb_overlaps(acb_realref(t), frobenius)) + { + flint_printf("FAIL (frobenius norm overlap)\n", iter); + flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec); + flint_printf("\n"); + + flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n"); + + flint_printf("sqrt(trace(A^H A)) = \n"); + acb_printd(t, 15); flint_printf("\n\n"); + + flint_printf("frobenius_norm(A) = \n"); + arb_printd(frobenius, 15); flint_printf("\n\n"); + + abort(); + } + + arb_clear(frobenius); + } + + acb_mat_clear(A); + acb_mat_clear(AH); + acb_mat_clear(AHA); + acb_clear(t); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arb_mat/test/t-frobenius_norm.c b/arb_mat/test/t-frobenius_norm.c new file mode 100644 index 00000000..0110f761 --- /dev/null +++ b/arb_mat/test/t-frobenius_norm.c @@ -0,0 +1,235 @@ +/*============================================================================= + + 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" + +static void +_fmpq_mat_sum_of_squares(fmpq_t res, const fmpq_mat_t Q) +{ + slong i, j; + fmpq_zero(res); + for (i = 0; i < fmpq_mat_nrows(Q); i++) + { + for (j = 0; j < fmpq_mat_ncols(Q); j++) + { + fmpq_addmul(res, fmpq_mat_entry(Q, i, j), fmpq_mat_entry(Q, i, j)); + } + } +} + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("frobenius_norm...."); + fflush(stdout); + + flint_randinit(state); + + /* compare to the exact rational norm */ + for (iter = 0; iter < 10000; iter++) + { + fmpq_mat_t Q; + fmpq_t q; + arb_mat_t A; + slong n, qbits, prec; + + n = n_randint(state, 8); + qbits = 1 + n_randint(state, 100); + prec = 2 + n_randint(state, 200); + + fmpq_mat_init(Q, n, n); + fmpq_init(q); + + arb_mat_init(A, n, n); + + fmpq_mat_randtest(Q, state, qbits); + _fmpq_mat_sum_of_squares(q, Q); + + arb_mat_set_fmpq_mat(A, Q, prec); + + /* check that the arb interval contains the exact value */ + { + arb_t a; + arb_init(a); + + arb_mat_frobenius_norm(a, A, prec); + arb_mul(a, a, a, prec); + + if (!arb_contains_fmpq(a, q)) + { + flint_printf("FAIL (containment, iter = %wd)\n", iter); + flint_printf("n = %wd, prec = %wd\n", n, prec); + flint_printf("\n"); + + flint_printf("Q = \n"); fmpq_mat_print(Q); flint_printf("\n\n"); + flint_printf("frobenius_norm(Q)^2 = \n"); + fmpq_print(q); flint_printf("\n\n"); + + flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n"); + flint_printf("frobenius_norm(A)^2 = \n"); + arb_printd(a, 15); flint_printf("\n\n"); + flint_printf("frobenius_norm(A)^2 = \n"); + arb_print(a); flint_printf("\n\n"); + + abort(); + } + + arb_clear(a); + } + + /* check that the upper bound is not less than the exact value */ + { + mag_t b; + fmpq_t y; + + mag_init(b); + fmpq_init(y); + + arb_mat_bound_frobenius_norm(b, A); + mag_mul(b, b, b); + mag_get_fmpq(y, b); + + if (fmpq_cmp(q, y) > 0) + { + flint_printf("FAIL (bound, iter = %wd)\n", iter); + flint_printf("n = %wd, prec = %wd\n", n, prec); + flint_printf("\n"); + + flint_printf("Q = \n"); fmpq_mat_print(Q); flint_printf("\n\n"); + flint_printf("frobenius_norm(Q)^2 = \n"); + fmpq_print(q); flint_printf("\n\n"); + + flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n"); + flint_printf("bound_frobenius_norm(A)^2 = \n"); + mag_printd(b, 15); flint_printf("\n\n"); + flint_printf("bound_frobenius_norm(A)^2 = \n"); + mag_print(b); flint_printf("\n\n"); + + abort(); + } + + mag_clear(b); + fmpq_clear(y); + } + + fmpq_mat_clear(Q); + fmpq_clear(q); + arb_mat_clear(A); + } + + /* check trace(A^T A) = frobenius_norm(A)^2 */ + for (iter = 0; iter < 10000; iter++) + { + slong m, n, prec; + arb_mat_t A, AT, ATA; + arb_t t; + + prec = 2 + n_randint(state, 200); + + m = n_randint(state, 10); + n = n_randint(state, 10); + + arb_mat_init(A, m, n); + arb_mat_init(AT, n, m); + arb_mat_init(ATA, n, n); + arb_init(t); + + arb_mat_randtest(A, state, 2 + n_randint(state, 100), 10); + arb_mat_transpose(AT, A); + arb_mat_mul(ATA, AT, A, prec); + arb_mat_trace(t, ATA, prec); + arb_sqrt(t, t, prec); + + /* check the norm bound */ + { + mag_t low, frobenius; + + mag_init(low); + arb_get_mag_lower(low, t); + + mag_init(frobenius); + arb_mat_bound_frobenius_norm(frobenius, A); + + if (mag_cmp(low, frobenius) > 0) + { + flint_printf("FAIL (bound)\n", iter); + flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec); + flint_printf("\n"); + + flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n"); + + flint_printf("lower(sqrt(trace(A^T A))) = \n"); + mag_printd(low, 15); flint_printf("\n\n"); + + flint_printf("bound_frobenius_norm(A) = \n"); + mag_printd(frobenius, 15); flint_printf("\n\n"); + + abort(); + } + + mag_clear(low); + mag_clear(frobenius); + } + + /* check the norm interval */ + { + arb_t frobenius; + + arb_init(frobenius); + arb_mat_frobenius_norm(frobenius, A, prec); + + if (!arb_overlaps(t, frobenius)) + { + flint_printf("FAIL (overlap)\n", iter); + flint_printf("m = %wd, n = %wd, prec = %wd\n", m, n, prec); + flint_printf("\n"); + + flint_printf("A = \n"); arb_mat_printd(A, 15); flint_printf("\n\n"); + + flint_printf("sqrt(trace(A^T A)) = \n"); + arb_printd(t, 15); flint_printf("\n\n"); + + flint_printf("frobenius_norm(A) = \n"); + arb_printd(frobenius, 15); flint_printf("\n\n"); + + abort(); + } + + arb_clear(frobenius); + } + + arb_mat_clear(A); + arb_mat_clear(AT); + arb_mat_clear(ATA); + arb_clear(t); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +}