block LU algorithms for acb_mat

acb_mat.h

@@ -66,6 +66,16 @@ acb_mat_swap(acb_mat_t mat1, acb_mat_t mat2)
     *mat2 = t;
 }
 
+/* Window matrices */
+
+void acb_mat_window_init(acb_mat_t window, const acb_mat_t mat, slong r1, slong c1, slong r2, slong c2);
+
+ARB_MAT_INLINE void
+acb_mat_window_clear(acb_mat_t window)
+{
+    flint_free(window->rows);
+}
+
 /* Conversions */
 
 void acb_mat_set(acb_mat_t dest, const acb_mat_t src);
@@ -358,6 +368,16 @@ acb_mat_swap_rows(acb_mat_t mat, slong * perm, slong r, slong s)
 slong acb_mat_find_pivot_partial(const acb_mat_t mat,
                                     slong start_row, slong end_row, slong c);
 
+void acb_mat_solve_tril_classical(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec);
+void acb_mat_solve_tril_recursive(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec);
+void acb_mat_solve_tril(acb_mat_t X, const acb_mat_t L, const acb_mat_t B, int unit, slong prec);
+
+void acb_mat_solve_triu_classical(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec);
+void acb_mat_solve_triu_recursive(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec);
+void acb_mat_solve_triu(acb_mat_t X, const acb_mat_t U, const acb_mat_t B, int unit, slong prec);
+
+int acb_mat_lu_classical(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec);
+int acb_mat_lu_recursive(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec);
 int acb_mat_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec);
 
 void acb_mat_solve_lu_precomp(acb_mat_t X, const slong * perm,

acb_mat/lu.c

@@ -14,60 +14,9 @@
 int
 acb_mat_lu(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)
 {
-    acb_t d, e;
+    if (acb_mat_nrows(A) < 8 || acb_mat_ncols(A) < 8)
-    acb_ptr * a;
+        return acb_mat_lu_classical(P, LU, A, prec);
-    slong i, j, m, n, r, row, col;
+    else
-    int result;
+        return acb_mat_lu_recursive(P, LU, A, prec);
-
-    if (acb_mat_is_empty(A))
-        return 1;
-
-    m = acb_mat_nrows(A);
-    n = acb_mat_ncols(A);
-
-    acb_mat_set(LU, A);
-
-    a = LU->rows;
-
-    row = col = 0;
-    for (i = 0; i < m; i++)
-        P[i] = i;
-
-    acb_init(d);
-    acb_init(e);
-
-    result = 1;
-
-    while (row < m && col < n)
-    {
-        r = acb_mat_find_pivot_partial(LU, row, m, col);
-
-        if (r == -1)
-        {
-            result = 0;
-            break;
-        }
-        else if (r != row)
-            acb_mat_swap_rows(LU, P, row, r);
-
-        acb_set(d, a[row] + col);
-
-        for (j = row + 1; j < m; j++)
-        {
-            acb_div(e, a[j] + col, d, prec);
-            acb_neg(e, e);
-            _acb_vec_scalar_addmul(a[j] + col,
-                a[row] + col, n - col, e, prec);
-            acb_zero(a[j] + col);
-            acb_neg(a[j] + row, e);
-        }
-
-        row++;
-        col++;
-    }
-
-    acb_clear(d);
-    acb_clear(e);
-
-    return result;
 }
+

acb_mat/lu_classical.c

@@ -0,0 +1,74 @@
+/*
+    Copyright (C) 2012 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"
+
+int
+acb_mat_lu_classical(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)
+{
+    acb_t d, e;
+    acb_ptr * a;
+    slong i, j, m, n, r, row, col;
+    int result;
+
+    if (acb_mat_is_empty(A))
+        return 1;
+
+    m = acb_mat_nrows(A);
+    n = acb_mat_ncols(A);
+
+    acb_mat_set(LU, A);
+
+    a = LU->rows;
+
+    row = col = 0;
+    for (i = 0; i < m; i++)
+        P[i] = i;
+
+    acb_init(d);
+    acb_init(e);
+
+    result = 1;
+
+    while (row < m && col < n)
+    {
+        r = acb_mat_find_pivot_partial(LU, row, m, col);
+
+        if (r == -1)
+        {
+            result = 0;
+            break;
+        }
+        else if (r != row)
+            acb_mat_swap_rows(LU, P, row, r);
+
+        acb_set(d, a[row] + col);
+
+        for (j = row + 1; j < m; j++)
+        {
+            acb_div(e, a[j] + col, d, prec);
+            acb_neg(e, e);
+            _acb_vec_scalar_addmul(a[j] + col,
+                a[row] + col, n - col, e, prec);
+            acb_zero(a[j] + col);
+            acb_neg(a[j] + row, e);
+        }
+
+        row++;
+        col++;
+    }
+
+    acb_clear(d);
+    acb_clear(e);
+
+    return result;
+}
+

acb_mat/lu_recursive.c

@@ -0,0 +1,113 @@
+/*
+    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
+_apply_permutation(slong * AP, acb_mat_t A, slong * P,
+    slong n, slong offset)
+{
+    if (n != 0)
+    {
+        acb_ptr * Atmp;
+        slong * APtmp;
+        slong i;
+
+        Atmp = flint_malloc(sizeof(acb_ptr) * n);
+        APtmp = flint_malloc(sizeof(slong) * n);
+
+        for (i = 0; i < n; i++) Atmp[i] = A->rows[P[i] + offset];
+        for (i = 0; i < n; i++) A->rows[i + offset] = Atmp[i];
+
+        for (i = 0; i < n; i++) APtmp[i] = AP[P[i] + offset];
+        for (i = 0; i < n; i++) AP[i + offset] = APtmp[i];
+
+        flint_free(Atmp);
+        flint_free(APtmp);
+    }
+}
+
+int
+acb_mat_lu_recursive(slong * P, acb_mat_t LU, const acb_mat_t A, slong prec)
+{
+    slong i, m, n, r1, r2, n1;
+    acb_mat_t A0, A1, A00, A01, A10, A11;
+    slong * P1;
+
+    m = A->r;
+    n = A->c;
+
+    if (m <= 1 || n <= 1)
+    {
+        return acb_mat_lu_classical(P, LU, A, prec);
+    }
+
+    if (LU != A)
+        acb_mat_set(LU, A);
+
+    n1 = n / 2;
+
+    for (i = 0; i < m; i++)
+        P[i] = i;
+
+    P1 = flint_malloc(sizeof(slong) * m);
+    acb_mat_window_init(A0, LU, 0, 0, m, n1);
+    acb_mat_window_init(A1, LU, 0, n1, m, n);
+
+    r1 = acb_mat_lu(P1, A0, A0, prec);
+
+    if (!r1)
+    {
+        flint_free(P1);
+        acb_mat_window_clear(A0);
+        acb_mat_window_clear(A1);
+        return 0;
+    }
+
+    /* r1 = rank of A0 */
+    r1 = FLINT_MIN(m, n1);
+
+    _apply_permutation(P, LU, P1, m, 0);
+
+    acb_mat_window_init(A00, LU, 0, 0, r1, r1);
+    acb_mat_window_init(A10, LU, r1, 0, m, r1);
+    acb_mat_window_init(A01, LU, 0, n1, r1, n);
+    acb_mat_window_init(A11, LU, r1, n1, m, n);
+
+    acb_mat_solve_tril(A01, A00, A01, 1, prec);
+
+    {
+        /* acb_mat_submul(A11, A11, A10, A01, prec); */
+        acb_mat_t T;
+        acb_mat_init(T, A10->r, A01->c);
+        acb_mat_mul(T, A10, A01, prec);
+        acb_mat_sub(A11, A11, T, prec);
+        acb_mat_clear(T);
+    }
+
+    r2 = acb_mat_lu(P1, A11, A11, prec);
+
+    if (!r2)
+        r1 = r2 = 0;
+    else
+        _apply_permutation(P, LU, P1, m - r1, r1);
+
+    flint_free(P1);
+    acb_mat_window_clear(A00);
+    acb_mat_window_clear(A01);
+    acb_mat_window_clear(A10);
+    acb_mat_window_clear(A11);
+    acb_mat_window_clear(A0);
+    acb_mat_window_clear(A1);
+
+    return r1 && r2;
+}
+

acb_mat/solve_lu_precomp.c

@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2012 Fredrik Johansson
+    Copyright (C) 2012,2018 Fredrik Johansson
 
     This file is part of Arb.
 
@@ -46,6 +46,15 @@ acb_mat_solve_lu_precomp(acb_mat_t X, const slong * perm,
         }
     }
 
+    /* todo: solve_tril and solve_triu have some overhead; should be
+       able to eliminate the basecase code below */
+    if (n >= 8 && m >= 8)
+    {
+        acb_mat_solve_tril(X, A, X, 1, prec);
+        acb_mat_solve_triu(X, A, X, 0, prec);
+        return;
+    }
+
     for (c = 0; c < m; c++)
     {
         /* solve Ly = b */

acb_mat/solve_tril.c

@@ -0,0 +1,108 @@
+/*
+    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"
+
+void
+acb_mat_solve_tril_classical(acb_mat_t X,
+        const acb_mat_t L, const acb_mat_t B, int unit, slong prec)
+{
+    slong i, j, k, n, m;
+    acb_ptr tmp;
+    acb_t s;
+
+    n = L->r;
+    m = B->c;
+
+    acb_init(s);
+    tmp = _acb_vec_init(n);
+
+    for (i = 0; i < m; i++)
+    {
+        for (j = 0; j < n; j++)
+            acb_set(tmp + j, acb_mat_entry(X, j, i));
+
+        for (j = 0; j < n; j++)
+        {
+            acb_zero(s);
+            for (k = 0; k < j; k++)
+                acb_addmul(s, L->rows[j] + k, tmp + k, prec);
+            acb_sub(s, acb_mat_entry(B, j, i), s, prec);
+            if (!unit)
+                acb_div(s, s, acb_mat_entry(L, j, j), prec);
+            acb_set(tmp + j, s);
+        }
+
+        for (j = 0; j < n; j++)
+            acb_set(acb_mat_entry(X, j, i), tmp + j);
+    }
+
+    _acb_vec_clear(tmp, n);
+    acb_clear(s);
+}
+
+void
+acb_mat_solve_tril_recursive(acb_mat_t X,
+        const acb_mat_t L, const acb_mat_t B, int unit, slong prec)
+{
+    acb_mat_t LA, LC, LD, XX, XY, BX, BY, T;
+    slong r, n, m;
+
+    n = L->r;
+    m = B->c;
+    r = n / 2;
+
+    if (n == 0 || m == 0)
+        return;
+
+    /*
+    Denoting inv(M) by M^, we have:
+
+    [A 0]^ [X]  ==  [A^          0 ] [X]  ==  [A^ X]
+    [C D]  [Y]  ==  [-D^ C A^    D^] [Y]  ==  [D^ (Y - C A^ X)]
+    */
+    acb_mat_window_init(LA, L, 0, 0, r, r);
+    acb_mat_window_init(LC, L, r, 0, n, r);
+    acb_mat_window_init(LD, L, r, r, n, n);
+    acb_mat_window_init(BX, B, 0, 0, r, m);
+    acb_mat_window_init(BY, B, r, 0, n, m);
+    acb_mat_window_init(XX, X, 0, 0, r, m);
+    acb_mat_window_init(XY, X, r, 0, n, m);
+
+    acb_mat_solve_tril(XX, LA, BX, unit, prec);
+
+    /* acb_mat_submul(XY, BY, LC, XX); */
+    acb_mat_init(T, LC->r, BX->c);
+    acb_mat_mul(T, LC, XX, prec);
+    acb_mat_sub(XY, BY, T, prec);
+    acb_mat_clear(T);
+
+    acb_mat_solve_tril(XY, LD, XY, unit, prec);
+
+    acb_mat_window_clear(LA);
+    acb_mat_window_clear(LC);
+    acb_mat_window_clear(LD);
+    acb_mat_window_clear(BX);
+    acb_mat_window_clear(BY);
+    acb_mat_window_clear(XX);
+    acb_mat_window_clear(XY);
+}
+
+void
+acb_mat_solve_tril(acb_mat_t X, const acb_mat_t L,
+                                    const acb_mat_t B, int unit, slong prec)
+{
+    if (B->r < 8 || B->c < 8)
+        acb_mat_solve_tril_classical(X, L, B, unit, prec);
+    else
+        acb_mat_solve_tril_recursive(X, L, B, unit, prec);
+}
+

acb_mat/solve_triu.c

@@ -0,0 +1,110 @@
+/*
+    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"
+
+void
+acb_mat_solve_triu_classical(acb_mat_t X, const acb_mat_t U,
+    const acb_mat_t B, int unit, slong prec)
+{
+    slong i, j, k, n, m;
+    acb_ptr tmp;
+    acb_t s;
+
+    n = U->r;
+    m = B->c;
+
+    tmp = _acb_vec_init(n);
+    acb_init(s);
+
+    for (i = 0; i < m; i++)
+    {
+        for (j = 0; j < n; j++)
+            acb_set(tmp + j, acb_mat_entry(X, j, i));
+
+        for (j = n - 1; j >= 0; j--)
+        {
+            acb_zero(s);
+            for (k = 0; k < n - j - 1; k++)
+                acb_addmul(s, U->rows[j] + j + 1 + k, tmp + j + 1 + k, prec);
+
+            acb_sub(s, acb_mat_entry(B, j, i), s, prec);
+            if (!unit)
+                acb_div(s, s, acb_mat_entry(U, j, j), prec);
+
+            acb_set(tmp + j, s);
+        }
+
+        for (j = 0; j < n; j++)
+            acb_set(acb_mat_entry(X, j, i), tmp + j);
+    }
+
+    _acb_vec_clear(tmp, n);
+    acb_clear(s);
+}
+
+void
+acb_mat_solve_triu_recursive(acb_mat_t X,
+        const acb_mat_t U, const acb_mat_t B, int unit, slong prec)
+{
+    acb_mat_t UA, UB, UD, XX, XY, BX, BY, T;
+    slong r, n, m;
+
+    n = U->r;
+    m = B->c;
+    r = n / 2;
+
+    if (n == 0 || m == 0)
+        return;
+
+    /*
+    Denoting inv(M) by M^, we have:
+    [A B]^ [X]  ==  [A^ (X - B D^ Y)]
+    [0 D]  [Y]  ==  [    D^ Y       ]
+    */
+
+    acb_mat_window_init(UA, U, 0, 0, r, r);
+    acb_mat_window_init(UB, U, 0, r, r, n);
+    acb_mat_window_init(UD, U, r, r, n, n);
+    acb_mat_window_init(BX, B, 0, 0, r, m);
+    acb_mat_window_init(BY, B, r, 0, n, m);
+    acb_mat_window_init(XX, X, 0, 0, r, m);
+    acb_mat_window_init(XY, X, r, 0, n, m);
+
+    acb_mat_solve_triu(XY, UD, BY, unit, prec);
+
+    /* acb_mat_submul(XX, BX, UB, XY); */
+    acb_mat_init(T, UB->r, XY->c);
+    acb_mat_mul(T, UB, XY, prec);
+    acb_mat_sub(XX, BX, T, prec);
+    acb_mat_clear(T);
+
+    acb_mat_solve_triu(XX, UA, XX, unit, prec);
+
+    acb_mat_window_clear(UA);
+    acb_mat_window_clear(UB);
+    acb_mat_window_clear(UD);
+    acb_mat_window_clear(BX);
+    acb_mat_window_clear(BY);
+    acb_mat_window_clear(XX);
+    acb_mat_window_clear(XY);
+}
+
+void
+acb_mat_solve_triu(acb_mat_t X, const acb_mat_t U,
+                                    const acb_mat_t B, int unit, slong prec)
+{
+    if (B->r < 8 || B->c < 8)
+        acb_mat_solve_triu_classical(X, U, B, unit, prec);
+    else
+        acb_mat_solve_triu_recursive(X, U, B, unit, prec);
+}
+

acb_mat/test/t-lu.c

@@ -39,7 +39,7 @@ int main()
         slong i, j, n, qbits, prec, *perm;
         int q_invertible, r_invertible;
 
-        n = n_randint(state, 8);
+        n = n_randint(state, 10);
         qbits = 1 + n_randint(state, 100);
         prec = 2 + n_randint(state, 202);
 

acb_mat/test/t-lu_recursive.c

@@ -0,0 +1,177 @@
+/*
+    Copyright (C) 2012 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"
+
+int fmpq_mat_is_invertible(const fmpq_mat_t A)
+{
+    int r;
+    fmpq_t t;
+    fmpq_init(t);
+    fmpq_mat_det(t, A);
+    r = !fmpq_is_zero(t);
+    fmpq_clear(t);
+    return r;
+}
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("lu_recursive....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* Dummy test with rectangular matrices. Rectangular matrices are
+       not actually supported (the output may be bogus), but the algorithm
+       should at least not crash. */
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        slong m, n, prec;
+        slong *perm;
+        acb_mat_t A, LU;
+
+        n = n_randint(state, 20);
+        m = n_randint(state, 20);
+        prec = 2 + n_randint(state, 200);
+
+        acb_mat_init(A, n, m);
+        acb_mat_init(LU, n, m);
+        perm = _perm_init(n);
+
+        acb_mat_randtest(A, state, prec, 10);
+
+        if (n_randint(state, 2))
+        {
+            acb_mat_lu_recursive(perm, LU, A, prec);
+        }
+        else
+        {
+            acb_mat_set(LU, A);
+            acb_mat_lu_recursive(perm, LU, LU, prec);
+        }
+
+        acb_mat_clear(A);
+        acb_mat_clear(LU);
+        _perm_clear(perm);
+    }
+
+    for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
+    {
+        fmpq_mat_t Q;
+        acb_mat_t A, LU, P, L, U, T;
+        slong i, j, n, qbits, prec, *perm;
+        int q_invertible, r_invertible;
+
+        n = n_randint(state, 20);
+        qbits = 1 + n_randint(state, 100);
+        prec = 2 + n_randint(state, 202);
+
+        fmpq_mat_init(Q, n, n);
+        acb_mat_init(A, n, n);
+        acb_mat_init(LU, n, n);
+        acb_mat_init(P, n, n);
+        acb_mat_init(L, n, n);
+        acb_mat_init(U, n, n);
+        acb_mat_init(T, n, n);
+        perm = _perm_init(n);
+
+        fmpq_mat_randtest(Q, state, qbits);
+        q_invertible = fmpq_mat_is_invertible(Q);
+
+        if (!q_invertible)
+        {
+            acb_mat_set_fmpq_mat(A, Q, prec);
+            r_invertible = acb_mat_lu_recursive(perm, LU, A, prec);
+            if (r_invertible)
+            {
+                flint_printf("FAIL: matrix is singular over Q but not over R\n");
+                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("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
+                flint_printf("LU = \n"); acb_mat_printd(LU, 15); flint_printf("\n\n");
+            }
+        }
+        else
+        {
+            /* now this must converge */
+            while (1)
+            {
+                acb_mat_set_fmpq_mat(A, Q, prec);
+                r_invertible = acb_mat_lu_recursive(perm, LU, A, prec);
+                if (r_invertible)
+                {
+                    break;
+                }
+                else
+                {
+                    if (prec > 10000)
+                    {
+                        flint_printf("FAIL: failed to converge at 10000 bits\n");
+                        flint_abort();
+                    }
+                    prec *= 2;
+                }
+            }
+
+            acb_mat_one(L);
+            for (i = 0; i < n; i++)
+                for (j = 0; j < i; j++)
+                    acb_set(acb_mat_entry(L, i, j),
+                        acb_mat_entry(LU, i, j));
+
+            for (i = 0; i < n; i++)
+                for (j = i; j < n; j++)
+                    acb_set(acb_mat_entry(U, i, j),
+                        acb_mat_entry(LU, i, j));
+
+            for (i = 0; i < n; i++)
+                acb_one(acb_mat_entry(P, perm[i], i));
+
+            acb_mat_mul(T, P, L, prec);
+            acb_mat_mul(T, T, U, prec);
+
+            if (!acb_mat_contains_fmpq_mat(T, 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("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
+                flint_printf("LU = \n"); acb_mat_printd(LU, 15); flint_printf("\n\n");
+                flint_printf("L = \n"); acb_mat_printd(L, 15); flint_printf("\n\n");
+                flint_printf("U = \n"); acb_mat_printd(U, 15); flint_printf("\n\n");
+                flint_printf("P*L*U = \n"); acb_mat_printd(T, 15); flint_printf("\n\n");
+
+                flint_abort();
+            }
+        }
+
+        fmpq_mat_clear(Q);
+        acb_mat_clear(A);
+        acb_mat_clear(LU);
+        acb_mat_clear(P);
+        acb_mat_clear(L);
+        acb_mat_clear(U);
+        acb_mat_clear(T);
+        _perm_clear(perm);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_mat/test/t-solve.c

@@ -28,8 +28,8 @@ int main()
         slong n, m, qbits, prec;
         int q_invertible, r_invertible, r_invertible2;
 
-        n = n_randint(state, 8);
+        n = n_randint(state, 10);
-        m = n_randint(state, 8);
+        m = n_randint(state, 10);
         qbits = 1 + n_randint(state, 30);
         prec = 2 + n_randint(state, 200);
 

acb_mat/test/t-solve_lu.c

@@ -0,0 +1,135 @@
+/*
+    Copyright (C) 2012 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"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("solve_lu....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+    {
+        fmpq_mat_t Q, QX, QB;
+        acb_mat_t A, X, B;
+        slong n, m, qbits, prec;
+        int q_invertible, r_invertible, r_invertible2;
+
+        n = n_randint(state, 10);
+        m = n_randint(state, 10);
+        qbits = 1 + n_randint(state, 30);
+        prec = 2 + n_randint(state, 200);
+
+        fmpq_mat_init(Q, n, n);
+        fmpq_mat_init(QX, n, m);
+        fmpq_mat_init(QB, n, m);
+
+        acb_mat_init(A, n, n);
+        acb_mat_init(X, n, m);
+        acb_mat_init(B, n, m);
+
+        fmpq_mat_randtest(Q, state, qbits);
+        fmpq_mat_randtest(QB, state, qbits);
+
+        q_invertible = fmpq_mat_solve_fraction_free(QX, Q, QB);
+
+        if (!q_invertible)
+        {
+            acb_mat_set_fmpq_mat(A, Q, prec);
+            r_invertible = acb_mat_solve_lu(X, A, B, prec);
+            if (r_invertible)
+            {
+                flint_printf("FAIL: matrix is singular over Q but not over R\n");
+                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("QX = \n"); fmpq_mat_print(QX); flint_printf("\n\n");
+                flint_printf("QB = \n"); fmpq_mat_print(QB); flint_printf("\n\n");
+                flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
+                flint_abort();
+            }
+        }
+        else
+        {
+            /* now this must converge */
+            while (1)
+            {
+                acb_mat_set_fmpq_mat(A, Q, prec);
+                acb_mat_set_fmpq_mat(B, QB, prec);
+
+                r_invertible = acb_mat_solve_lu(X, A, B, prec);
+                if (r_invertible)
+                {
+                    break;
+                }
+                else
+                {
+                    if (prec > 10000)
+                    {
+                        flint_printf("FAIL: failed to converge at 10000 bits\n");
+                        flint_printf("Q = \n"); fmpq_mat_print(Q); flint_printf("\n\n");
+                        flint_printf("QX = \n"); fmpq_mat_print(QX); flint_printf("\n\n");
+                        flint_printf("QB = \n"); fmpq_mat_print(QB); flint_printf("\n\n");
+                        flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
+                        flint_abort();
+                    }
+                    prec *= 2;
+                }
+            }
+
+            if (!acb_mat_contains_fmpq_mat(X, QX))
+            {
+                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("QB = \n"); fmpq_mat_print(QB); flint_printf("\n\n");
+                flint_printf("QX = \n"); fmpq_mat_print(QX); flint_printf("\n\n");
+
+                flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
+                flint_printf("B = \n"); acb_mat_printd(B, 15); flint_printf("\n\n");
+                flint_printf("X = \n"); acb_mat_printd(X, 15); flint_printf("\n\n");
+
+                flint_abort();
+            }
+
+            /* test aliasing */
+            r_invertible2 = acb_mat_solve_lu(B, A, B, prec);
+            if (!acb_mat_equal(X, B) || r_invertible != r_invertible2)
+            {
+                flint_printf("FAIL (aliasing)\n");
+                flint_printf("A = \n"); acb_mat_printd(A, 15); flint_printf("\n\n");
+                flint_printf("B = \n"); acb_mat_printd(B, 15); flint_printf("\n\n");
+                flint_printf("X = \n"); acb_mat_printd(X, 15); flint_printf("\n\n");
+                flint_abort();
+            }
+        }
+
+        fmpq_mat_clear(Q);
+        fmpq_mat_clear(QB);
+        fmpq_mat_clear(QX);
+        acb_mat_clear(A);
+        acb_mat_clear(B);
+        acb_mat_clear(X);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_mat/test/t-solve_tril.c

@@ -0,0 +1,99 @@
+/*
+    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"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("solve_tril....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
+    {
+        acb_mat_t A, X, B, Y;
+        slong rows, cols, prec, i, j;
+        int unit;
+
+        prec = 2 + n_randint(state, 200);
+        rows = n_randint(state, 30);
+        cols = n_randint(state, 30);
+        unit = n_randint(state, 2);
+
+        acb_mat_init(A, rows, rows);
+        acb_mat_init(B, rows, cols);
+        acb_mat_init(X, rows, cols);
+        acb_mat_init(Y, rows, cols);
+
+        acb_mat_randtest(A, state, prec, 10);
+        acb_mat_randtest(X, state, prec, 10);
+        acb_mat_randtest(Y, state, prec, 10);
+
+        for (i = 0; i < rows; i++)
+        {
+            if (unit)
+                acb_one(acb_mat_entry(A, i, i));
+            else
+                acb_set_ui(acb_mat_entry(A, i, i), 1 + n_randint(state, 100));
+
+            for (j = i + 1; j < rows; j++)
+                acb_zero(acb_mat_entry(A, i, j));
+        }
+
+        acb_mat_mul(B, A, X, prec);
+
+        if (unit)  /* check that diagonal entries are ignored */
+        {
+            for (i = 0; i < rows; i++)
+                acb_set_ui(acb_mat_entry(A, i, i), 1 + n_randint(state, 100));
+        }
+
+        /* Check Y = A^(-1) * (A * X) = X */
+        acb_mat_solve_tril(Y, A, B, unit, prec);
+
+        if (!acb_mat_overlaps(Y, X))
+        {
+            flint_printf("FAIL\n");
+            flint_printf("A = \n"); acb_mat_printd(A, 10); flint_printf("\n\n");
+            flint_printf("B = \n"); acb_mat_printd(B, 10); flint_printf("\n\n");
+            flint_printf("X = \n"); acb_mat_printd(X, 10); flint_printf("\n\n");
+            flint_printf("Y = \n"); acb_mat_printd(Y, 10); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        /* Check aliasing */
+        acb_mat_solve_tril(B, A, B, unit, prec);
+        if (!acb_mat_equal(B, Y))
+        {
+            flint_printf("FAIL (aliasing)\n");
+            flint_printf("A = \n"); acb_mat_printd(A, 10); flint_printf("\n\n");
+            flint_printf("B = \n"); acb_mat_printd(B, 10); flint_printf("\n\n");
+            flint_printf("X = \n"); acb_mat_printd(X, 10); flint_printf("\n\n");
+            flint_printf("Y = \n"); acb_mat_printd(Y, 10); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        acb_mat_clear(A);
+        acb_mat_clear(B);
+        acb_mat_clear(X);
+        acb_mat_clear(Y);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_mat/test/t-solve_triu.c

@@ -0,0 +1,99 @@
+/*
+    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"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("solve_triu....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
+    {
+        acb_mat_t A, X, B, Y;
+        slong rows, cols, prec, i, j;
+        int unit;
+
+        prec = 2 + n_randint(state, 200);
+        rows = n_randint(state, 30);
+        cols = n_randint(state, 30);
+        unit = n_randint(state, 2);
+
+        acb_mat_init(A, rows, rows);
+        acb_mat_init(B, rows, cols);
+        acb_mat_init(X, rows, cols);
+        acb_mat_init(Y, rows, cols);
+
+        acb_mat_randtest(A, state, prec, 10);
+        acb_mat_randtest(X, state, prec, 10);
+        acb_mat_randtest(Y, state, prec, 10);
+
+        for (i = 0; i < rows; i++)
+        {
+            if (unit)
+                acb_one(acb_mat_entry(A, i, i));
+            else
+                acb_set_ui(acb_mat_entry(A, i, i), 1 + n_randint(state, 100));
+
+            for (j = 0; j < i; j++)
+                acb_zero(acb_mat_entry(A, i, j));
+        }
+
+        acb_mat_mul(B, A, X, prec);
+
+        if (unit)  /* check that diagonal entries are ignored */
+        {
+            for (i = 0; i < rows; i++)
+                acb_set_ui(acb_mat_entry(A, i, i), 1 + n_randint(state, 100));
+        }
+
+        /* Check Y = A^(-1) * (A * X) = X */
+        acb_mat_solve_triu(Y, A, B, unit, prec);
+
+        if (!acb_mat_overlaps(Y, X))
+        {
+            flint_printf("FAIL\n");
+            flint_printf("A = \n"); acb_mat_printd(A, 10); flint_printf("\n\n");
+            flint_printf("B = \n"); acb_mat_printd(B, 10); flint_printf("\n\n");
+            flint_printf("X = \n"); acb_mat_printd(X, 10); flint_printf("\n\n");
+            flint_printf("Y = \n"); acb_mat_printd(Y, 10); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        /* Check aliasing */
+        acb_mat_solve_triu(B, A, B, unit, prec);
+        if (!acb_mat_equal(B, Y))
+        {
+            flint_printf("FAIL (aliasing)\n");
+            flint_printf("A = \n"); acb_mat_printd(A, 10); flint_printf("\n\n");
+            flint_printf("B = \n"); acb_mat_printd(B, 10); flint_printf("\n\n");
+            flint_printf("X = \n"); acb_mat_printd(X, 10); flint_printf("\n\n");
+            flint_printf("Y = \n"); acb_mat_printd(Y, 10); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        acb_mat_clear(A);
+        acb_mat_clear(B);
+        acb_mat_clear(X);
+        acb_mat_clear(Y);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_mat/window_init.c

@@ -0,0 +1,29 @@
+/*
+    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"
+
+void
+acb_mat_window_init(acb_mat_t window, const acb_mat_t mat,
+    slong r1, slong c1, slong r2, slong c2)
+{
+    slong i;
+    window->entries = NULL;
+
+    window->rows = flint_malloc((r2 - r1) * sizeof(acb_ptr));
+
+    for (i = 0; i < r2 - r1; i++)
+        window->rows[i] = mat->rows[r1 + i] + c1;
+
+    window->r = r2 - r1;
+    window->c = c2 - c1;
+}
+