# HG changeset patch # User Martin Albrecht # Date 1189866063 -3600 # Node ID faf211fd9c6707064e1649a7eb1fb7e374d9654e # Parent 2cce131bce328a2ca283d5fb590b7aa8ddab5200 first steps towards faster sparse linear algerba over GF(p) diff -r 2cce131bce32 -r faf211fd9c67 sage/libs/linbox/linbox.pxd --- a/sage/libs/linbox/linbox.pxd Sat Sep 15 14:32:28 2007 +0100 +++ b/sage/libs/linbox/linbox.pxd Sat Sep 15 15:21:03 2007 +0100 @@ -1,7 +1,6 @@ include "../../ext/cdefs.pxi" include "../../ext/cdefs.pxi" -cdef extern from "../libs/m4ri/packedmatrix.h": - ctypedef struct packedmatrix +include '../../modules/vector_modn_sparse_h.pxi' ctypedef size_t mod_int @@ -18,14 +17,16 @@ cdef class Linbox_modn_dense: size_t B_nr, size_t B_nc) cdef unsigned long rank(self) except -1 +cdef class Linbox_modn_sparse: + cdef c_vector_modint *rows + cdef size_t nrows + cdef size_t ncols + cdef unsigned int modulus -## cdef class Linbox_mod2_dense: -## cdef packedmatrix *matrix - -## cdef set(self, packedmatrix *matrix) -## cdef int echelonize(self) -## cdef matrix_matrix_multiply(self, packedmatrix *ans, packedmatrix *B) -## cdef unsigned long rank(self) except -1 + cdef set(self, int modulus, size_t nrows, size_t ncols, c_vector_modint *rows) + cdef object rank(self, int reorder) + #cdef void solve(self, void *x, void *b) + cdef class Linbox_integer_dense: cdef mpz_t** matrix diff -r 2cce131bce32 -r faf211fd9c67 sage/libs/linbox/linbox.pyx --- a/sage/libs/linbox/linbox.pyx Sat Sep 15 14:32:28 2007 +0100 +++ b/sage/libs/linbox/linbox.pyx Sat Sep 15 15:21:03 2007 +0100 @@ -87,82 +87,35 @@ cdef class Linbox_modn_dense: r = linbox_modn_dense_rank(self.n, self.matrix, self.nrows, self.ncols) return r - - - - - -########################################################################## -## Dense matrices GF(2) -########################################################################## - -## cdef extern from "linbox_wrap.h": -## ctypedef struct packedmatrix - -## cdef int linbox_mod2_dense_echelonize(packedmatrix *m) - -## cdef void linbox_mod2_dense_minpoly(mod_int **mp, size_t* degree, packedmatrix *matrix, int do_minpoly) - -## cdef int linbox_mod2_dense_matrix_matrix_multiply(packedmatrix *ans, packedmatrix *A, packedmatrix *B) - -## cdef int linbox_mod2_dense_rank(packedmatrix *m) - -## cdef class Linbox_mod2_dense: -## cdef set(self, packedmatrix *matrix): -## self.matrix = matrix - -## cdef int echelonize(self): -## cdef int r -## r = linbox_mod2_dense_echelonize(self.matrix) -## return r - -## def minpoly(self): -## return self._poly(True) - -## def charpoly(self): -## return self._poly(False) - -## def _poly(self, minpoly): +########################################################################## +## Sparse matices modulo p. +########################################################################## + +cdef extern from "linbox_wrap.h": + int linbox_modn_sparse_matrix_rank(unsigned long modulus, size_t nrows, size_t ncols, void* rows, int reorder) #, int **pivots) + void linbox_modn_sparse_matrix_solve(unsigned long modulus, size_t numrows, size_t numcols, void **x, void *a, void *b) + +cdef class Linbox_modn_sparse: + cdef set(self, int modulus, size_t nrows, size_t ncols, c_vector_modint *rows): + self.rows = rows + self.nrows = nrows + self.ncols = ncols + self.modulus = modulus + + cdef object rank(self, int reorder): + #cdef int *_pivots + cdef int r + r = linbox_modn_sparse_matrix_rank(self.modulus, self.nrows, self.ncols, self.rows, reorder) + + #pivots = [_pivots[i] for i in range(r)] + #free(_pivots) + return r#, pivots + +## cdef void solve(self, void *x, void *b): ## """ -## INPUT: -## as given - -## OUTPUT: -## coefficients of charpoly or minpoly as a Python list ## """ -## cdef mod_int *f -## cdef size_t degree -## linbox_mod2_dense_minpoly(&f, °ree, self.matrix, minpoly) - -## v = [] -## cdef Py_ssize_t i -## for i from 0 <= i <= degree: -## v.append(f[i]) -## linbox_modn_dense_delete_array(f) -## return v - -## cdef matrix_matrix_multiply(self, -## packedmatrix *ans, -## packedmatrix *B): -## cdef int e - -## e = linbox_mod2_dense_matrix_matrix_multiply(ans, self.matrix, B) -## if e: -## raise RuntimeError, "error doing matrix matrix multiply mod2 using linbox" - -## cdef unsigned long rank(self) except -1: -## cdef unsigned long r -## r = linbox_mod2_dense_rank(self.matrix) -## return r - - -########################################################################## -## Sparse matices modulo p. -########################################################################## -cdef class Linbox_modn_sparse: - pass - - +## linbox_modn_sparse_matrix_solve(self.modulus, self.nrows, self.ncols, &x, self.rows, b) + ########################################################################## ## Sparse matrices over ZZ ########################################################################## diff -r 2cce131bce32 -r faf211fd9c67 sage/matrix/matrix_modn_sparse.pxd --- a/sage/matrix/matrix_modn_sparse.pxd Sat Sep 15 14:32:28 2007 +0100 +++ b/sage/matrix/matrix_modn_sparse.pxd Sat Sep 15 15:21:03 2007 +0100 @@ -6,3 +6,5 @@ cdef class Matrix_modn_sparse(matrix_spa cdef c_vector_modint* rows cdef public int p cdef swap_rows_c(self, Py_ssize_t n1, Py_ssize_t n2) + + cdef _init_linbox(self) diff -r 2cce131bce32 -r faf211fd9c67 sage/matrix/matrix_modn_sparse.pyx --- a/sage/matrix/matrix_modn_sparse.pyx Sat Sep 15 14:32:28 2007 +0100 +++ b/sage/matrix/matrix_modn_sparse.pyx Sat Sep 15 15:21:03 2007 +0100 @@ -79,6 +79,13 @@ from sage.rings.integer_mod cimport Inte from sage.misc.misc import verbose, get_verbose, graphics_filename +from sage.rings.integer import Integer + +from sage.matrix.matrix0 import Matrix as Matrix0 +from sage.rings.arith import is_prime + +from sage.structure.element import is_Vector + ################ # TODO: change this to use extern cdef's methods. from sage.ext.arith cimport arith_int @@ -88,6 +95,10 @@ ai = arith_int() import sage.ext.multi_modular MAX_MODULUS = sage.ext.multi_modular.MAX_MODULUS + +from sage.libs.linbox.linbox cimport Linbox_modn_sparse +cdef Linbox_modn_sparse linbox +linbox = Linbox_modn_sparse() cdef class Matrix_modn_sparse(matrix_sparse.Matrix_sparse): @@ -462,3 +473,197 @@ cdef class Matrix_modn_sparse(matrix_spa for j from 0 <= j < self.rows[i].num_nonzero: k+=1 return QQ(k)/QQ(self.nrows()*self.ncols()) + + cdef _init_linbox(self): + _sig_on + linbox.set(self.p, self._nrows, self._ncols, self.rows) + _sig_off + + def _rank_linbox(self): + """ + See self.rank(). + """ + if is_prime(self.p): + x = self.fetch('rank') + if not x is None: + return x + self._init_linbox() + _sig_on + # the returend pivots list is currently wrong + #r, pivots = linbox.rank(1) + r = linbox.rank(1) + r = Integer(r) + _sig_off + self.cache('rank', r) + return r + else: + raise TypeError, "only GF(p) supported via LinBox" + + def rank(self, gauss=False): + """ + Compute the rank of self. + + INPUT: + gauss -- if True Gaussian elimination is used. If False + 'Symbolic Reordering' as implemented in LinBox + is used. (default: False) + + EXAMPLE: + sage: A = random_matrix(GF(127),200,200,density=0.01,sparse=True) + sage: r1 = A.rank(gauss=False) + sage: r2 = A.rank(gauss=True) + sage: r1 == r2 + True + + ALGORITHM: Uses LinBox or native implementation. + + REFERENCES: Jean-Guillaume Dumas and Gilles Villars. 'Computing the + Rank of Large Sparse Matrices over Finite Fields'. Proc. CASC'2002, + The Fifth International Workshop on Computer Algebra in Scientific Computing, + Big Yalta, Crimea, Ukraine, 22-27 sept. 2002, Springer-Verlag, + http://perso.ens-lyon.fr/gilles.villard/BIBLIOGRAPHIE/POSTSCRIPT/rankjgd.ps + + NOTE: + For very sparse matrices Gaussian elimination is faster because + it barly has anything to do. If the fill in needs to be considered, + 'Symbolic Reordering' is usually much faster. + """ + x = self.fetch('rank') + if not x is None: return x + + if is_prime(self.p) and not gauss: + return self._rank_linbox() + else: + return Matrix0.rank(self) + + def transpose(self): + """ + Return the transpose of self. + + EXAMPLE: + sage: A = matrix(GF(127),3,3,[0,1,0,2,0,0,3,0,0],sparse=True) + sage: A + [0 1 0] + [2 0 0] + [3 0 0] + sage: A.transpose() + [0 2 3] + [1 0 0] + [0 0 0] + """ + cdef int i, j + cdef c_vector_modint row + cdef Matrix_modn_sparse B + + B = self.new_matrix(nrows = self.ncols(), ncols = self.nrows()) + for i from 0 <= i < self._nrows: + row = self.rows[i] + for j from 0 <= j < row.num_nonzero: + set_entry(&B.rows[row.positions[j]], i, row.entries[j]) + return B + + def matrix_from_rows(self, rows): + """ + Return the matrix constructed from self using rows with indices + in the rows list. + + INPUT: + rows -- list or tuple of row indices + + EXAMPLE: + sage: M = MatrixSpace(GF(127),3,3,sparse=True) + sage: A = M(range(9)); A + [0 1 2] + [3 4 5] + [6 7 8] + sage: A.matrix_from_rows([2,1]) + [6 7 8] + [3 4 5] + """ + cdef int i,k + cdef Matrix_modn_sparse A + cdef c_vector_modint row + + if not isinstance(rows, (list, tuple)): + raise TypeError, "rows must be a list of integers" + + + A = self.new_matrix(nrows = len(rows)) + + k = 0 + for ii in rows: + i = ii + if i < 0 or i >= self.nrows(): + raise IndexError, "row %s out of range"%i + + row = self.rows[i] + for j from 0 <= j < row.num_nonzero: + set_entry(&A.rows[k], row.positions[j], row.entries[j]) + k += 1 + return A + + + def matrix_from_columns(self, cols): + """ + Return the matrix constructed from self using columns with + indices in the columns list. + + EXAMPLES: + sage: M = MatrixSpace(GF(127),3,3,sparse=True) + sage: A = M(range(9)); A + [0 1 2] + [3 4 5] + [6 7 8] + sage: A.matrix_from_columns([2,1]) + [2 1] + [5 4] + [8 7] + """ + cdef int i,j + cdef Matrix_modn_sparse A + cdef c_vector_modint row + + if not isinstance(cols, (list, tuple)): + raise TypeError, "rows must be a list of integers" + + A = self.new_matrix(ncols = len(cols)) + + cols = dict(zip([int(e) for e in cols],range(len(cols)))) + + for i from 0 <= i < self.nrows(): + row = self.rows[i] + for j from 0 <= j < row.num_nonzero: + if int(row.positions[j]) in cols: + set_entry(&A.rows[i], cols[int(row.positions[j])], row.entries[j]) + return A + +## def _solve(self, Matrix_modn_sparse B): +## """ +## """ +## cdef c_vector_modint *X + +## cdef Matrix_modn_sparse C + +## if not self.is_square(): +## raise NotImplementedError, "input matrix must be square" + +## if self.rank() != self.nrows(): +## raise ValueError, "input matrix must have full rank but it doesn't" + +## self._init_linbox() + +## matrix = True +## if is_Vector(B): +## matrix = False +## C = self.matrix_space(1,self.nrows())(B.list()) + +## #_sig_on +## linbox.solve(&X, C.rows) +## #_sig_off + +## ## X = None +## ## if not matrix: +## ## # Convert back to a vector +## ## return (X.base_ring() ** X.nrows())(X.list()) +## ## else: +## ## return X