Use entries + stride instead of entries + array of row pointers in matrix types

src/acb_mat.h

@@ -25,7 +25,7 @@
 extern "C" {
 #endif
 
-#define acb_mat_entry(mat,i,j) ((mat)->rows[i] + (j))
+#define acb_mat_entry(mat,i,j) ((mat)->entries + (i) * (mat)->stride + (j))
 #define acb_mat_nrows(mat) ((mat)->r)
 #define acb_mat_ncols(mat) ((mat)->c)
 
@@ -54,9 +54,8 @@ void acb_mat_swap_entrywise(acb_mat_t mat1, acb_mat_t mat2);
 void acb_mat_window_init(acb_mat_t window, const acb_mat_t mat, slong r1, slong c1, slong r2, slong c2);
 
 ACB_MAT_INLINE void
-acb_mat_window_clear(acb_mat_t window)
+acb_mat_window_clear(acb_mat_t FLINT_UNUSED(window))
 {
-    flint_free(window->rows);
 }
 
 /* Conversions */
@@ -226,17 +225,7 @@ void acb_mat_vector_mul_col(acb_ptr res, const acb_mat_t A, acb_srcptr v, slong
 
 /* Solving */
 
-ACB_MAT_INLINE void
-acb_mat_swap_rows(acb_mat_t mat, slong * perm, slong r, slong s)
-{
-    if (r != s)
-    {
-        if (perm != NULL)
-            FLINT_SWAP(slong, perm[r], perm[s]);
-
-        FLINT_SWAP(acb_ptr, mat->rows[r], mat->rows[s]);
-    }
-}
+void 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);

src/acb_mat/approx_eig_qr.c

@@ -584,7 +584,7 @@ acb_mat_approx_eig_triu_r(acb_mat_t ER, const acb_mat_t A, slong prec)
 
         for (j = i - 1; j >= 0; j--)
         {
-            acb_approx_dot(r, NULL, 0, A->rows[j] + j + 1, 1, ER->rows[i] + j + 1, 1, i - j, prec);
+            acb_approx_dot(r, NULL, 0, acb_mat_entry(A, j, j + 1), 1, acb_mat_entry(ER, i, j + 1), 1, i - j, prec);
             acb_approx_sub(t, acb_mat_entry(A, j, j), s, prec);
 
             /* if abs(t) < smin: t = smin */
@@ -688,7 +688,7 @@ acb_mat_approx_eig_triu_l(acb_mat_t EL, const acb_mat_t A, slong prec)
 
         for (j = i + 1; j < n; j++)
         {
-            acb_approx_dot(r, NULL, 0, EL->rows[i] + i, 1, AT->rows[j] + i, 1, j - i, prec);
+            acb_approx_dot(r, NULL, 0, acb_mat_entry(EL, i, i), 1, acb_mat_entry(AT, j, i), 1, j - i, prec);
             acb_approx_sub(t, acb_mat_entry(AT, j, j), s, prec);
 
             /* if abs(t) < smin: t = smin */

src/acb_mat/approx_lu.c

@@ -9,27 +9,31 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include <string.h>
 #include "acb.h"
 #include "acb_mat.h"
 
 static void
-_apply_permutation(slong * AP, acb_mat_t A, slong * P,
-    slong n, slong offset)
+_apply_permutation(slong * AP, acb_mat_t A, const slong * P,
+    slong num_rows, slong row_offset, slong num_cols, slong col_offset)
 {
-    if (n != 0)
+    if (num_rows != 0)
     {
-        acb_ptr * Atmp;
+        acb_ptr Atmp;
         slong * APtmp;
         slong i;
 
-        Atmp = flint_malloc(sizeof(acb_ptr) * n);
-        APtmp = flint_malloc(sizeof(slong) * n);
+        /* todo: reduce memory allocation */
+        Atmp = flint_malloc(sizeof(acb_struct) * num_rows * num_cols);
+        APtmp = flint_malloc(sizeof(slong) * num_rows);
 
-        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 < num_rows; i++)
+            memcpy(Atmp + i * num_cols, acb_mat_entry(A, P[i] + row_offset, col_offset), sizeof(acb_struct) * num_cols);
+        for (i = 0; i < num_rows; i++)
+            memcpy(acb_mat_entry(A, i + row_offset, col_offset), Atmp + i * num_cols, sizeof(acb_struct) * num_cols);
 
-        for (i = 0; i < n; i++) APtmp[i] = AP[P[i] + offset];
-        for (i = 0; i < n; i++) AP[i + offset] = APtmp[i];
+        for (i = 0; i < num_rows; i++) APtmp[i] = AP[P[i] + row_offset];
+        for (i = 0; i < num_rows; i++) AP[i + row_offset] = APtmp[i];
 
         flint_free(Atmp);
         flint_free(APtmp);
@@ -86,7 +90,6 @@ int
 acb_mat_approx_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;
 
@@ -98,8 +101,6 @@ acb_mat_approx_lu_classical(slong * P, acb_mat_t LU, const acb_mat_t A, slong pr
 
     acb_mat_get_mid(LU, A);
 
-    a = LU->rows;
-
     row = col = 0;
     for (i = 0; i < m; i++)
         P[i] = i;
@@ -121,16 +122,16 @@ acb_mat_approx_lu_classical(slong * P, acb_mat_t LU, const acb_mat_t A, slong pr
         else if (r != row)
             acb_mat_swap_rows(LU, P, row, r);
 
-        _acb_approx_inv(d, a[row] + col, prec);
+        _acb_approx_inv(d, acb_mat_entry(LU, row, col), prec);
 
         for (j = row + 1; j < m; j++)
         {
-            _acb_approx_mul(e, a[j] + col, d, prec);
+            _acb_approx_mul(e, acb_mat_entry(LU, j, col), d, prec);
             acb_neg(e, e);
-            _acb_vec_approx_scalar_addmul(a[j] + col,
-                a[row] + col, n - col, e, prec);
-            acb_zero(a[j] + col);
-            acb_neg(a[j] + row, e);
+            _acb_vec_approx_scalar_addmul(acb_mat_entry(LU, j, col),
+                acb_mat_entry(LU, row, col), n - col, e, prec);
+            acb_zero(acb_mat_entry(LU, j, col));
+            acb_neg(acb_mat_entry(LU, j, row), e);
         }
 
         row++;
@@ -182,7 +183,7 @@ acb_mat_approx_lu_recursive(slong * P, acb_mat_t LU, const acb_mat_t A, slong pr
     /* r1 = rank of A0 */
     r1 = FLINT_MIN(m, n1);
 
-    _apply_permutation(P, LU, P1, m, 0);
+    _apply_permutation(P, LU, P1, m, 0, n - n1, n1);
 
     acb_mat_window_init(A00, LU, 0, 0, r1, r1);
     acb_mat_window_init(A10, LU, r1, 0, m, r1);
@@ -206,8 +207,7 @@ acb_mat_approx_lu_recursive(slong * P, acb_mat_t LU, const acb_mat_t A, slong pr
     if (!r2)
         r1 = r2 = 0;
     else
-        _apply_permutation(P, LU, P1, m - r1, r1);
-
+        _apply_permutation(P, LU, P1, m - r1, r1, n1, 0);
     flint_free(P1);
     acb_mat_window_clear(A00);
     acb_mat_window_clear(A01);

src/acb_mat/approx_mul.c

@@ -62,26 +62,11 @@ acb_mat_approx_mul_classical(acb_mat_t C, const acb_mat_t A, const acb_mat_t B,
     }
     else
     {
-        acb_ptr tmp;
-        TMP_INIT;
-
-        TMP_START;
-        tmp = TMP_ALLOC(sizeof(acb_struct) * br * bc);
-
-        for (i = 0; i < br; i++)
-            for (j = 0; j < bc; j++)
-                tmp[j * br + i] = *acb_mat_entry(B, i, j);
-
         for (i = 0; i < ar; i++)
-        {
             for (j = 0; j < bc; j++)
-            {
                 acb_approx_dot(acb_mat_entry(C, i, j), NULL, 0,
-                    A->rows[i], 1, tmp + j * br, 1, br, prec);
-            }
-        }
-
-        TMP_END;
+                    acb_mat_entry(A, i, 0), 1,
+                    acb_mat_entry(B, 0, j), B->stride, br, prec);
     }
 }
 

src/acb_mat/approx_solve_lu_precomp.c

@@ -28,9 +28,9 @@
         for (c = 0; c < m; c++)
         {
             for (i = 0; i < n; i++)
-                tmp[i] = B->rows[perm[i]][c];
+                tmp[i] = *acb_mat_entry(B, perm[i], c);
             for (i = 0; i < n; i++)
-                X->rows[i][c] = tmp[i];
+                *acb_mat_entry(X, i, c) = tmp[i];
         }
 
         flint_free(tmp);

src/acb_mat/approx_solve_tril.c

@@ -57,7 +57,7 @@ acb_mat_approx_solve_tril_classical(acb_mat_t X,
 
         for (j = 0; j < n; j++)
         {
-            acb_approx_dot(s, acb_mat_entry(B, j, i), 1, L->rows[j], 1, tmp, 1, j, prec);
+            acb_approx_dot(s, acb_mat_entry(B, j, i), 1, acb_mat_entry(L, j, 0), 1, tmp, 1, j, prec);
 
             if (!unit)
                 acb_approx_div(tmp + j, s, acb_mat_entry(L, j, j), t, prec);

src/acb_mat/approx_solve_triu.c

@@ -57,7 +57,7 @@ acb_mat_approx_solve_triu_classical(acb_mat_t X, const acb_mat_t U,
 
         for (j = n - 1; j >= 0; j--)
         {
-            acb_approx_dot(s, acb_mat_entry(B, j, i), 1, U->rows[j] + j + 1, 1, tmp + j + 1, 1, n - j - 1, prec);
+            acb_approx_dot(s, acb_mat_entry(B, j, i), 1, acb_mat_entry(U, j, j + 1), 1, tmp + j + 1, 1, n - j - 1, prec);
 
             if (!unit)
                 acb_approx_div(tmp + j, s, acb_mat_entry(U, j, j), t, prec);

src/acb_mat/charpoly.c

@@ -9,69 +9,23 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include "gr.h"
+#include "gr_mat.h"
 #include "acb_poly.h"
 #include "acb_mat.h"
 
 void _acb_mat_charpoly(acb_ptr cp, const acb_mat_t mat, slong prec)
 {
-    const slong n = mat->r;
-
-    if (n == 0)
-    {
-        acb_one(cp);
-    }
-    else if (n == 1)
+    if (!acb_mat_is_finite(mat))
     {
-        acb_neg(cp + 0, acb_mat_entry(mat, 0, 0));
-        acb_one(cp + 1);
+        _acb_vec_indeterminate(cp, mat->r + 1);
     }
     else
     {
-        slong i, k, t;
-        acb_ptr a, A, s;
-
-        a = _acb_vec_init(n * n);
-        A = a + (n - 1) * n;
-
-        _acb_vec_zero(cp, n + 1);
-        acb_neg(cp + 0, acb_mat_entry(mat, 0, 0));
-
-        for (t = 1; t < n; t++)
-        {
-            for (i = 0; i <= t; i++)
-            {
-                acb_set(a + 0 * n + i, acb_mat_entry(mat, i, t));
-            }
-
-            acb_set(A + 0, acb_mat_entry(mat, t, t));
-
-            for (k = 1; k < t; k++)
-            {
-                for (i = 0; i <= t; i++)
-                {
-                    s = a + k * n + i;
-                    acb_dot(s, NULL, 0, mat->rows[i], 1, a + (k - 1) * n, 1, t + 1, prec);
-                }
-
-                acb_set(A + k, a + k * n + t);
-            }
-
-            acb_dot(A + t, NULL, 0, mat->rows[t], 1, a + (t - 1) * n, 1, t + 1, prec);
-
-            for (k = 0; k <= t; k++)
-            {
-                acb_dot(cp + k, cp + k, 1, A, 1, cp + k - 1, -1, k, prec);
-                acb_sub(cp + k, cp + k, A + k, prec);
-            }
-        }
-
-        /* Shift all coefficients up by one */
-        for (i = n; i > 0; i--)
-            acb_swap(cp + i, cp + (i - 1));
-
-        acb_one(cp + 0);
-        _acb_poly_reverse(cp, cp, n + 1, n + 1);
-        _acb_vec_clear(a, n * n);
+        gr_ctx_t ctx;
+        gr_ctx_init_complex_acb(ctx, prec);
+        if (_gr_mat_charpoly_berkowitz(cp, (const gr_mat_struct *) mat, ctx) != GR_SUCCESS)
+            _acb_vec_indeterminate(cp, mat->r + 1);
     }
 }
 

src/acb_mat/clear.c

@@ -16,8 +16,5 @@ void
 acb_mat_clear(acb_mat_t mat)
 {
     if (mat->entries != NULL)
-    {
         _acb_vec_clear(mat->entries, mat->r * mat->c);
-        flint_free(mat->rows);
-    }
 }

src/acb_mat/det_lu.c

@@ -16,12 +16,10 @@ slong
 acb_mat_gauss_partial(acb_mat_t A, slong prec)
 {
     acb_t e;
-    acb_ptr * a;
     slong j, m, n, r, rank, row, col, sign;
 
     m = A->r;
     n = A->c;
-    a = A->rows;
     rank = row = col = 0;
     sign = 1;
 
@@ -45,9 +43,9 @@ acb_mat_gauss_partial(acb_mat_t A, slong prec)
 
         for (j = row + 1; j < m; j++)
         {
-            acb_div(e, a[j] + col, a[row] + col, prec);
+            acb_div(e, acb_mat_entry(A, j, col), acb_mat_entry(A, row, col), prec);
             acb_neg(e, e);
-            _acb_vec_scalar_addmul(a[j] + col + 1, a[row] + col + 1, n - col - 1, e, prec);
+            _acb_vec_scalar_addmul(acb_mat_entry(A, j, col + 1), acb_mat_entry(A, row, col + 1), n - col - 1, e, prec);
         }
 
         row++;
@@ -110,7 +108,7 @@ acb_mat_det_lu_inplace(acb_t det, acb_mat_t A, slong prec)
 
         for (i = rank; i < n; i++)
         {
-            acb_vec_get_arf_2norm_squared_bound(t, A->rows[i] + rank,  n - rank, MAG_BITS);
+            acb_vec_get_arf_2norm_squared_bound(t, acb_mat_entry(A, i, rank),  n - rank, MAG_BITS);
             arf_mul(d, d, t, MAG_BITS, ARF_RND_UP);
         }
 

src/acb_mat/det_precond.c

@@ -9,19 +9,32 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include <string.h>
 #include "perm.h"
 #include "acb.h"
 #include "acb_mat.h"
 
 static void
-_apply_permutation(acb_mat_t A, slong * P, slong n)
+_apply_permutation(acb_mat_t A, const slong * P, slong num_rows)
 {
-    acb_ptr * Atmp;
-    slong i;
-    Atmp = flint_malloc(sizeof(acb_ptr) * n);
-    for (i = 0; i < n; i++) Atmp[i] = A->rows[P[i]];
-    for (i = 0; i < n; i++) A->rows[i] = Atmp[i];
-    flint_free(Atmp);
+    if (num_rows != 0)
+    {
+        acb_ptr Atmp;
+        slong i;
+        slong row_offset = 0;
+        slong col_offset = 0;
+        slong num_cols = A->c;
+
+        /* todo: reduce memory allocation */
+        Atmp = flint_malloc(sizeof(acb_struct) * num_rows * num_cols);
+
+        for (i = 0; i < num_rows; i++)
+            memcpy(Atmp + i * num_cols, acb_mat_entry(A, P[i] + row_offset, col_offset), sizeof(acb_struct) * num_cols);
+        for (i = 0; i < num_rows; i++)
+            memcpy(acb_mat_entry(A, i + row_offset, col_offset), Atmp + i * num_cols, sizeof(acb_struct) * num_cols);
+
+        flint_free(Atmp);
+    }
 }
 
 /* Enclosure of det(I + eps) using Gershgorin circles.