Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
r""" Sparse matrices over `\ZZ/n\ZZ` for `n` small
This is a compiled implementation of sparse matrices over `\ZZ/n\ZZ` for `n` small.
TODO: - move vectors into a Cython vector class - add _add_ and _mul_ methods.
EXAMPLES::
sage: a = matrix(Integers(37),3,3,range(9),sparse=True); a [0 1 2] [3 4 5] [6 7 8] sage: type(a) <type 'sage.matrix.matrix_modn_sparse.Matrix_modn_sparse'> sage: parent(a) Full MatrixSpace of 3 by 3 sparse matrices over Ring of integers modulo 37 sage: a^2 [15 18 21] [ 5 17 29] [32 16 0] sage: a+a [ 0 2 4] [ 6 8 10] [12 14 16] sage: b = a.new_matrix(2,3,range(6)); b [0 1 2] [3 4 5] sage: a*b Traceback (most recent call last): ... TypeError: unsupported operand parent(s) for *: 'Full MatrixSpace of 3 by 3 sparse matrices over Ring of integers modulo 37' and 'Full MatrixSpace of 2 by 3 sparse matrices over Ring of integers modulo 37' sage: b*a [15 18 21] [ 5 17 29]
::
sage: TestSuite(a).run() sage: TestSuite(b).run()
::
sage: a.echelonize(); a [ 1 0 36] [ 0 1 2] [ 0 0 0] sage: b.echelonize(); b [ 1 0 36] [ 0 1 2] sage: a.pivots() (0, 1) sage: b.pivots() (0, 1) sage: a.rank() 2 sage: b.rank() 2 sage: a[2,2] = 5 sage: a.rank() 3
TESTS::
sage: matrix(Integers(37),0,0,sparse=True).inverse() [] """
#***************************************************************************** # Copyright (C) 2006 William Stein <wstein@gmail.com> # # This program 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. # http://www.gnu.org/licenses/ #*****************************************************************************
from __future__ import absolute_import
from collections import Iterator, Sequence
from cysignals.memory cimport check_calloc, sig_malloc, sig_free from cysignals.signals cimport sig_on, sig_off
from sage.modules.vector_modn_sparse cimport *
from cpython.sequence cimport *
from sage.libs.gmp.mpz cimport mpz_init_set_si cimport sage.matrix.matrix as matrix cimport sage.matrix.matrix_sparse as matrix_sparse cimport sage.matrix.matrix_dense as matrix_dense from sage.rings.finite_rings.integer_mod cimport IntegerMod_int, IntegerMod_abstract
from sage.misc.misc import verbose, get_verbose
import sage.rings.all as rings
from sage.matrix.matrix2 import Matrix as Matrix2 from sage.arith.all import is_prime
from sage.structure.element import is_Vector
cimport sage.structure.element
from sage.data_structures.binary_search cimport * from sage.modules.vector_integer_sparse cimport *
from .matrix_integer_sparse cimport Matrix_integer_sparse from sage.misc.decorators import rename_keyword
################ # TODO: change this to use extern cdef's methods. from sage.rings.fast_arith cimport arith_int cdef arith_int ai ai = arith_int() ################
# The 46341 below is because the mod-n sparse code still uses # int's, even on 64-bit computers. Improving this is # Trac Ticket #12679. MAX_MODULUS = 46341
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): def __cinit__(self, parent, entries, copy, coerce):
# allocate memory cdef Py_ssize_t i, nr, nc cdef int p
def __dealloc__(self): cdef Py_ssize_t i
def __init__(self, parent, entries, copy, coerce): """ Create a sparse matrix over the integers modulo ``n``.
INPUT:
- ``parent`` -- a matrix space
- ``entries`` -- can be one of the following:
* a Python dictionary whose items have the form ``(i, j): x``, where ``0 <= i < nrows``, ``0 <= j < ncols``, and ``x`` is coercible to an element of the integers modulo ``n``. The ``i,j`` entry of ``self`` is set to ``x``. The ``x``'s can be ``0``. * Alternatively, entries can be a list of *all* the entries of the sparse matrix, read row-by-row from top to bottom (so they would be mostly 0).
- ``copy`` -- ignored
- ``coerce`` -- ignored """ cdef int s, z, p cdef Py_ssize_t i, j, k
cdef PyObject** X
# Sparse input format. raise IndexError("invalid entries list") # Dense input format seq = PySequence_Fast(entries,"expected a list") X = PySequence_Fast_ITEMS(seq) k = 0 R = self._base_ring # Get fast access to the entries list. for i from 0 <= i < self._nrows: for j from 0 <= j < self._ncols: z = R(<object>X[k]) if z != 0: set_entry(&self.rows[i], j, z) k = k + 1 else: # scalar?
cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j): cdef IntegerMod_int n
######################################################################## # LEVEL 2 functionality # * def _pickle # * def _unpickle # * cdef _add_ # * cdef _mul_ # * cpdef _cmp_ # * __neg__ # * __invert__ # * __copy__ # * _multiply_classical # * _list -- list of underlying elements (need not be a copy) # * x _dict -- sparse dictionary of underlying elements (need not be a copy) ######################################################################## # def _pickle(self): # def _unpickle(self, data, int version): # use version >= 0 # cpdef _add_(self, right): # cdef _mul_(self, Matrix right): # cpdef int _cmp_(self, Matrix right) except -2: # def __neg__(self): # def __invert__(self): # def __copy__(self): # def _multiply_classical(left, matrix.Matrix _right): # def _list(self):
def _dict(self): """ Unsafe version of the dict method, mainly for internal use. This may return the dict of elements, but as an *unsafe* reference to the underlying dict of the object. It might be dangerous if you change entries of the returned dict.
EXAMPLES::
sage: MS = MatrixSpace(GF(13), 50, 50, sparse=True) sage: m = MS.random_element(density=0.002) sage: m._dict() {(4, 44): 7, (5, 25): 4, (26, 9): 9, (43, 43): 6, (44, 38): 1}
TESTS::
sage: parent(m._dict()[26,9]) Finite Field of size 13 """
cdef Py_ssize_t i, j, k cdef IntegerMod_int n
def _pickle(self): """ TESTS::
sage: M = Matrix( GF(2), [[1,1,1,1,0,0,0,0,0,0]], sparse=True ) sage: loads(dumps(M)) [1 1 1 1 0 0 0 0 0 0] sage: loads(dumps(M)) == M True """
def _unpickle(self, data, version): else: raise ValueError("unknown matrix format")
cdef sage.structure.element.Matrix _matrix_times_matrix_(self, sage.structure.element.Matrix _right): """ This code is implicitly called for multiplying self by another sparse matrix.
EXAMPLES: sage: a = matrix(GF(43), 3, 3, range(9), sparse=True) sage: b = matrix(GF(43), 3, 3, range(10,19), sparse=True) sage: a*b [ 2 5 8] [33 2 14] [21 42 20] sage: a*a [15 18 21] [42 11 23] [26 4 25] sage: c = matrix(GF(43), 3, 8, range(24), sparse=True) sage: a*c [40 0 3 6 9 12 15 18] [26 38 7 19 31 0 12 24] [12 33 11 32 10 31 9 30]
Even though sparse and dense matrices are represented differently, they still compare as equal if they have the same entries: sage: a*b == a._matrix_times_matrix_dense(b) True sage: d = matrix(GF(43), 3, 8, range(24)) sage: a*c == a*d True
TESTS:
The following shows that :trac:`23669` has been addressed::
sage: p = next_prime(2**15) sage: M = Matrix(GF(p), 1,3, lambda i,j: -1, sparse=True); M [32770 32770 32770] sage: M*M.transpose() # previously returned [32738] [3]
""" cdef Matrix_modn_sparse right, ans
cdef c_vector_modint* v
# Build a table that gives the nonzero positions in each column of right cdef Py_ssize_t i, j, k
# Now do the multiplication, getting each row completely before filling it in. cdef int x, y, s
def _matrix_times_matrix_dense(self, sage.structure.element.Matrix _right): """ Multiply self by the sparse matrix _right, and return the result as a dense matrix.
EXAMPLES: sage: a = matrix(GF(10007), 2, [1,2,3,4], sparse=True) sage: b = matrix(GF(10007), 2, 3, [1..6], sparse=True) sage: a * b [ 9 12 15] [19 26 33] sage: c = a._matrix_times_matrix_dense(b); c [ 9 12 15] [19 26 33] sage: type(c) <type 'sage.matrix.matrix_modn_dense_double.Matrix_modn_dense_double'>
sage: a = matrix(GF(2), 20, 20, sparse=True) sage: a*a == a._matrix_times_matrix_dense(a) True sage: type(a._matrix_times_matrix_dense(a)) <type 'sage.matrix.matrix_mod2_dense.Matrix_mod2_dense'> """ cdef Matrix_modn_sparse right cdef matrix_dense.Matrix_dense ans
cdef c_vector_modint* v
# Build a table that gives the nonzero positions in each column of right cdef Py_ssize_t i, j, k
# Now do the multiplication, getting each row completely before filling it in. cdef int x, y, s
#ans._matrix[i][j] = s
######################################################################## # LEVEL 3 functionality (Optional) # * cdef _sub_ # * __deepcopy__ # * __invert__ # * Matrix windows -- only if you need strassen for that base # * Other functions (list them here): # x - echelon form in place # x - nonzero_positions ######################################################################## def swap_rows(self, r1, r2): self.check_bounds_and_mutability(r1,0) self.check_bounds_and_mutability(r2,0) self.swap_rows_c(r1, r2)
cdef swap_rows_c(self, Py_ssize_t n1, Py_ssize_t n2): """ Swap the rows in positions n1 and n2. No bounds checking. """ cdef c_vector_modint tmp
def _echelon_in_place_classical(self): """ Replace self by its reduction to reduced row echelon form.
ALGORITHM: We use Gauss elimination, in a slightly intelligent way, in that we clear each column using a row with the minimum number of nonzero entries.
TODO: Implement switching to a dense method when the matrix gets dense. """
cdef Py_ssize_t i, r, c, min, min_row, start_row cdef int a0, a_inverse, b, do_verb cdef c_vector_modint tmp
tm = verbose('on column %s of %s'%(c, self._ncols), level = 2, caller_name = 'matrix_modn_sparse echelon') #end if # Since there is at least one nonzero entry, the first entry # of the positions list is defined. It is the first position # of a nonzero entry, and it equals c precisely if row r # is a row we could use to clear column c. #endif #endif #endfor # print("min number of entries in a pivoting row = ", min) # Since we can use row r to clear column c, the # entry in position c in row r must be the first nonzero entry. &self.rows[start_row], self.p - b)
def _nonzero_positions_by_row(self, copy=True): """ Returns the list of pairs (i,j) such that self[i,j] != 0.
It is safe to change the resulting list (unless you give the option copy=False).
EXAMPLES:: sage: M = Matrix(GF(7), [[0,0,0,1,0,0,0,0],[0,1,0,0,0,0,1,0]], sparse=True); M [0 0 0 1 0 0 0 0] [0 1 0 0 0 0 1 0] sage: M.nonzero_positions() [(0, 3), (1, 1), (1, 6)]
""" cdef Py_ssize_t i, j
def density(self): """ Return the density of self, i.e., the ratio of the number of nonzero entries of self to the total size of self.
EXAMPLES::
sage: A = matrix(QQ,3,3,[0,1,2,3,0,0,6,7,8],sparse=True) sage: A.density() 2/3
Notice that the density parameter does not ensure the density of a matrix; it is only an upper bound.
::
sage: A = random_matrix(GF(127),200,200,density=0.3, sparse=True) sage: A.density() 2073/8000 """ cdef Py_ssize_t i, nonzero_entries
def transpose(self): """ Return the transpose of self.
EXAMPLES::
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]
``.T`` is a convenient shortcut for the transpose::
sage: A.T [0 2 3] [1 0 0] [0 0 0] """ cdef int i, j cdef c_vector_modint row cdef Matrix_modn_sparse 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
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_rows([2,1]) [6 7 8] [3 4 5] """ cdef int i,k cdef Matrix_modn_sparse A cdef c_vector_modint row
rows = list(rows)
raise IndexError("row %s out of range" % i)
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
cols = list(cols)
cdef _init_linbox(self):
@rename_keyword(deprecation=6094, method="algorithm") def _rank_linbox(self, algorithm): """ See self.rank(). """ return x # the returend pivots list is currently wrong #r, pivots = linbox.rank(1) else: raise TypeError("only GF(p) supported via LinBox")
def rank(self, gauss=False): """ Compute the rank of self.
INPUT:
- ``gauss`` - if True LinBox' Gaussian elimination is used. If False 'Symbolic Reordering' as implemented in LinBox is used. If 'native' the native Sage implementation is used. (default: False)
EXAMPLES::
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: r3 = A.rank(gauss='native') sage: r1 == r2 == r3 True sage: r1 155
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. """
elif gauss is True: return self._rank_linbox(1) elif gauss == "native": return Matrix2.rank(self) else: raise TypeError("parameter 'gauss' not understood") else: return Matrix2.rank(self)
def _solve_right_nonsingular_square(self, B, algorithm=None, check_rank = True): """ If self is a matrix `A`, then this function returns a vector or matrix `X` such that `A X = B`. If `B` is a vector then `X` is a vector and if `B` is a matrix, then `X` is a matrix.
.. NOTE::
In Sage one can also write ``A B`` for ``A.solve_right(B)``, i.e., Sage implements the "the MATLAB/Octave backslash operator".
INPUT:
- ``B`` - a matrix or vector
- ``algorithm`` - one of the following:
- ``'LinBox:BlasElimination'`` - dense elimination
- ``'LinBox:Blackbox'`` - LinBox chooses a Blackbox algorithm
- ``'LinBox:Wiedemann'`` - Wiedemann's algorithm
- ``'generic'`` - use generic implementation (inversion)
- ``None`` - LinBox chooses an algorithm (default)
- ``check_rank`` - check rank before attempting to solve (default: True)
OUTPUT: a matrix or vector
EXAMPLES::
sage: A = matrix(GF(127), 3, [1,2,3,-1,2,5,2,3,1], sparse=True) sage: b = vector(GF(127),[1,2,3]) sage: x = A \ b; x (73, 76, 10) sage: A * x (1, 2, 3) """ cdef Matrix_modn_sparse b cdef Matrix_modn_sparse X cdef c_vector_modint *x
B = B.change_ring(self.base_ring())
return Matrix2.solve_right(self, B)
raise ValueError("self must be of full rank.")
raise ValueError("input matrices must have the same number of rows.")
raise NotImplementedError("input matrix must be square")
matrix = False b = self.matrix_space(1, self.ncols(),sparse=True)(B.list()) else: B = B.sparse_matrix() else: raise TypeError("B must be a matrix or vector over the same base as self")
elif algorithm == "LinBox:BlasElimination": algorithm = 1 elif algorithm == "LinBox:Blackbox": algorithm = 2 elif algorithm == "LinBox:Wiedemann": algorithm = 3 else: raise TypeError("parameter 'algorithm' not understood")
# Convert back to a vector return (X.base_ring() ** X.ncols())(X.list()) else:
def lift(self): """ Return lift of this matrix to a sparse matrix over the integers.
EXAMPLES: sage: a = matrix(GF(7),2,3,[1..6], sparse=True) sage: a.lift() [1 2 3] [4 5 6] sage: a.lift().parent() Full MatrixSpace of 2 by 3 sparse matrices over Integer Ring
Subdivisions are preserved when lifting::
sage: a.subdivide([], [1,1]); a [1||2 3] [4||5 6] sage: a.lift() [1||2 3] [4||5 6] """ cdef Py_ssize_t i, j cdef Matrix_integer_sparse L 0, 0, 0)
cdef mpz_vector* L_row cdef c_vector_modint* A_row raise MemoryError("error allocating space for sparse vector during sparse lift") sig_free(L_row.entries) L_row.entries = NULL raise MemoryError("error allocating space for sparse vector during sparse lift")
|