add LU decomposition

doc/fmprb_mat.txt

@@ -135,6 +135,32 @@ void fmprb_mat_mul(fmprb_mat_t res, const fmprb_mat_t mat1,
     const fmprb_mat_t mat2, long prec)
 
     Sets res to the matrix product of mat1 and mat2. The operands must have
-    compatible dimensions for matrix multiplication. Aliasing of input
+    compatible dimensions for matrix multiplication.
-    and output is currently not supported.
+
+*******************************************************************************
+
+    LU decomposition
+
+*******************************************************************************
+
+int fmprb_mat_lu(long * P, fmprb_mat_t LU, const fmprb_mat_t A, long prec)
+
+    Given an $n \times n$ matrix $A$, computes an LU decomposition $PLU = A$
+    using Gaussian elimination with partial pivoting.
+    The input and output matrices can be the same, performing the
+    decomposition in-place.
+
+    Entry $i$ in the permutation vector is set to the row index in
+    the input matrix corresponding to row $i$ in the output matrix.
+
+    The algorithm succeeds and returns nonzero if it can find $n$ invertible
+    (i.e. not containing zero) pivot entries. This guarantees that the matrix
+    is invertible.
+
+    The algorithm fails and returns zero, leaving the values in $P$ and $LU$
+    undefined, if it cannot find $n$ invertible pivot elements.
+    In this case, either the matrix is singular, the input matrix was
+    computed to insufficient precision, or the LU decomposition was
+    attempted at insufficient precision.
+
 

fmprb.h

@@ -286,9 +286,12 @@ fmprb_print(const fmprb_t x)
 static __inline__ void
 fmprb_printd(const fmprb_t x, long digits)
 {
-    fmpr_printd(fmprb_midref(x), digits);
+    fmpr_printd(fmprb_midref(x), FLINT_ABS(digits));
-    printf(" +/- ");
+    if (digits > 0)
-    fmpr_printd(fmprb_radref(x), 5);
+    {
+        printf(" +/- ");
+        fmpr_printd(fmprb_radref(x), 5);
+    }
 }
 
 static __inline__ void

fmprb_mat.h

@@ -52,6 +52,14 @@ void fmprb_mat_init(fmprb_mat_t mat, long r, long c);
 
 void fmprb_mat_clear(fmprb_mat_t mat);
 
+static __inline__ void
+fmprb_mat_swap(fmprb_mat_t mat1, fmprb_mat_t mat2)
+{
+    fmprb_mat_struct t = *mat1;
+    *mat1 = *mat2;
+    *mat2 = t;
+}
+
 /* Conversions */
 
 void fmprb_mat_set(fmprb_mat_t dest, const fmprb_mat_t src);
@@ -86,5 +94,9 @@ void fmprb_mat_sub(fmprb_mat_t res, const fmprb_mat_t mat1, const fmprb_mat_t ma
 
 void fmprb_mat_mul(fmprb_mat_t res, const fmprb_mat_t mat1, const fmprb_mat_t mat2, long prec);
 
+/* LU decomposition */
+
+int fmprb_mat_lu(long * P, fmprb_mat_t LU, const fmprb_mat_t A, long prec);
+
 #endif
 

fmprb_mat/lu.c

@@ -0,0 +1,136 @@
+/*=============================================================================
+
+    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
+
+******************************************************************************/
+
+#include "fmprb_mat.h"
+
+static __inline__ void
+fmprb_mat_swap_rows(fmprb_mat_t mat, long * perm, long r, long s)
+{
+    if (r != s)
+    {
+        fmprb_struct * u;
+        long t;
+
+        if (perm != NULL)
+        {
+            t = perm[s];
+            perm[s] = perm[r];
+            perm[r] = t;
+        }
+
+        u = mat->rows[s];
+        mat->rows[s] = mat->rows[r];
+        mat->rows[r] = u;
+    }
+}
+
+static __inline__ long
+fmprb_mat_find_pivot_partial(const fmprb_mat_t mat,
+                                    long start_row, long end_row, long c)
+{
+    long best_row, i;
+
+    best_row = -1;
+
+    for (i = start_row; i < end_row; i++)
+    {
+        if (!fmprb_contains_zero(fmprb_mat_entry(mat, i, c)))
+        {
+            if (best_row == -1)
+            {
+                best_row = i;
+            }
+            /* todo: should take the radius into account */
+            else if (fmpr_cmpabs(fmprb_midref(fmprb_mat_entry(mat, i, c)),
+                    fmprb_midref(fmprb_mat_entry(mat, best_row, c))) > 0)
+            {
+                best_row = i;
+            }
+        }
+    }
+
+    return best_row;
+}
+
+int
+fmprb_mat_lu(long * P, fmprb_mat_t LU, const fmprb_mat_t A, long prec)
+{
+    fmprb_t d, e;
+    fmprb_struct ** a;
+    long i, j, m, n, r, length, row, col;
+    int result;
+
+    m = fmprb_mat_nrows(A);
+    n = fmprb_mat_ncols(A);
+
+    result = 1;
+
+    if (m == 0 || n == 0)
+        return result;
+
+    fmprb_mat_set(LU, A);
+
+    a = LU->rows;
+
+    row = col = 0;
+    for (i = 0; i < m; i++)
+        P[i] = i;
+
+    fmprb_init(d);
+    fmprb_init(e);
+
+    while (row < m && col < n)
+    {
+        r = fmprb_mat_find_pivot_partial(LU, row, m, col);
+
+        if (r == -1)
+        {
+            result = 0;
+            break;
+        }
+        else if (r != row)
+            fmprb_mat_swap_rows(LU, P, row, r);
+
+        fmprb_set(d, a[row] + col);
+
+        for (j = row + 1; j < m; j++)
+        {
+            fmprb_div(e, a[j] + col, d, prec);
+            fmprb_neg(e, e);
+            _fmprb_vec_scalar_addmul(a[j] + col,
+                a[row] + col, n - col, e, prec);
+            fmprb_zero(a[j] + col);
+            fmprb_neg(a[j] + row, e);
+        }
+
+        row++;
+        col++;
+    }
+
+    fmprb_clear(d);
+    fmprb_clear(e);
+
+    return result;
+}

fmprb_mat/test/t-lu.c

@@ -0,0 +1,160 @@
+/*=============================================================================
+
+    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
+
+******************************************************************************/
+
+#include "fmprb_mat.h"
+#include "perm.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()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("lu....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        fmpq_mat_t Q;
+        fmprb_mat_t A, LU, P, L, U, T;
+        long i, j, n, qbits, prec, *perm;
+        int result, q_invertible, r_invertible;
+
+        n = n_randint(state, 8);
+        qbits = 1 + n_randint(state, 100);
+        prec = 2 + n_randint(state, 202);
+
+        fmpq_mat_init(Q, n, n);
+        fmprb_mat_init(A, n, n);
+        fmprb_mat_init(LU, n, n);
+        fmprb_mat_init(P, n, n);
+        fmprb_mat_init(L, n, n);
+        fmprb_mat_init(U, n, n);
+        fmprb_mat_init(T, n, n);
+        perm = _perm_init(n);
+
+        fmpq_mat_randtest(Q, state, qbits);
+        q_invertible = fmpq_mat_is_invertible(Q);
+
+        result = 1;
+
+        if (!q_invertible)
+        {
+            fmprb_mat_set_fmpq_mat(A, Q, prec);
+            r_invertible = fmprb_mat_lu(perm, LU, A, prec);
+            if (r_invertible)
+            {
+                printf("FAIL: matrix is singular over Q but not over R\n");
+                printf("n = %ld, prec = %ld\n", n, prec);
+                printf("\n");
+
+                printf("Q = \n"); fmpq_mat_print(Q); printf("\n\n");
+                printf("A = \n"); fmprb_mat_printd(A, 15); printf("\n\n");
+                printf("LU = \n"); fmprb_mat_printd(LU, 15); printf("\n\n");
+            }
+        }
+        else
+        {
+            /* now this must converge */
+            while (1)
+            {
+                fmprb_mat_set_fmpq_mat(A, Q, prec);
+                r_invertible = fmprb_mat_lu(perm, LU, A, prec);
+                if (r_invertible)
+                {
+                    break;
+                }
+                else
+                {
+                    if (prec > 10000)
+                    {
+                        printf("FAIL: failed to converge at 10000 bits\n");
+                        abort();
+                    }
+                    prec *= 2;
+                }
+            }
+
+            fmprb_mat_one(L);
+            for (i = 0; i < n; i++)
+                for (j = 0; j < i; j++)
+                    fmprb_set(fmprb_mat_entry(L, i, j),
+                        fmprb_mat_entry(LU, i, j));
+
+            for (i = 0; i < n; i++)
+                for (j = i; j < n; j++)
+                    fmprb_set(fmprb_mat_entry(U, i, j),
+                        fmprb_mat_entry(LU, i, j));
+
+            for (i = 0; i < n; i++)
+                fmprb_one(fmprb_mat_entry(P, perm[i], i));
+
+            fmprb_mat_mul(T, P, L, prec);
+            fmprb_mat_mul(T, T, U, prec);
+
+            if (!fmprb_mat_contains_fmpq_mat(T, Q))
+            {
+                printf("FAIL (containment, iter = %ld)\n", iter);
+                printf("n = %ld, prec = %ld\n", n, prec);
+                printf("\n");
+
+                printf("Q = \n"); fmpq_mat_print(Q); printf("\n\n");
+                printf("A = \n"); fmprb_mat_printd(A, 15); printf("\n\n");
+                printf("LU = \n"); fmprb_mat_printd(LU, 15); printf("\n\n");
+                printf("L = \n"); fmprb_mat_printd(L, 15); printf("\n\n");
+                printf("U = \n"); fmprb_mat_printd(U, 15); printf("\n\n");
+                printf("P*L*U = \n"); fmprb_mat_printd(T, 15); printf("\n\n");
+
+                abort();
+            }
+        }
+
+        fmpq_mat_clear(Q);
+        fmprb_mat_clear(A);
+        fmprb_mat_clear(LU);
+        fmprb_mat_clear(P);
+        fmprb_mat_clear(L);
+        fmprb_mat_clear(U);
+        fmprb_mat_clear(T);
+        _perm_clear(perm);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}