Ticket #13177: trac_13177.patch

File trac_13177.patch, 39.4 KB (added by was, 10 years ago)
  • module_list.py

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1340936331 21600
    # Node ID 2c8f5d2a39d13e77b512699e7e48a0d464a1782d
    # Parent  1f1edda93473594deb3368c79867a19b4beff53e
    Trac #13177: optimize matrix_modn_dense echelon form by implementing switchover to dense linbox echelon
    
    diff --git a/module_list.py b/module_list.py
    a b  
    977977
    978978    Extension('sage.matrix.matrix_modn_sparse',
    979979              sources = ['sage/matrix/matrix_modn_sparse.pyx'],
    980               libraries = ['gmp']),
     980              language="c++",
     981              libraries = ['gmp', 'linbox', 'givaro', BLAS, BLAS2],
     982              extra_compile_args = ['-DDISABLE_COMMENTATOR']),
    981983
    982984    Extension('sage.matrix.matrix_mpolynomial_dense',
    983985              sources = ['sage/matrix/matrix_mpolynomial_dense.pyx'],
  • sage/libs/linbox/linbox.pxd

    diff --git a/sage/libs/linbox/linbox.pxd b/sage/libs/linbox/linbox.pxd
    a b  
    11include "../../ext/cdefs.pxi"
    22
    3 include '../../modules/vector_modn_sparse_h.pxi'
     3from sage.matrix.matrix_modn_sparse cimport c_vector_modint
    44
    55from sage.matrix.matrix_integer_dense cimport mod_int
    66
  • sage/matrix/matrix2.pyx

    diff --git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
    a b  
    57145714            sage: m == transformation_matrix * m_original
    57155715            True
    57165716        """
     5717        if self.fetch('in_echelon_form') and not kwds.get('transformation',False):
     5718            return
    57175719        self.check_mutability()
    57185720
    57195721        if algorithm == 'default':
  • sage/matrix/matrix_modn_sparse.pyx

    diff --git a/sage/matrix/matrix_modn_sparse.pyx b/sage/matrix/matrix_modn_sparse.pyx
    a b  
    11r"""
    2 Sparse matrices over `\ZZ/n\ZZ` for `n` small
     2Sparse matrices over `\ZZ/n\ZZ` for small `n`, currently at most 46341.
    33
    44This is a compiled implementation of sparse matrices over
    5 `\ZZ/n\ZZ` for `n` small.
     5`\ZZ/n\ZZ` for small `n`.
    66
    7 TODO: - move vectors into a Cython vector class - add _add_ and
    8 _mul_ methods.
     7TODO:
     8    - Move vectors into a Cython vector class
     9    - Add _add_ and _mul_ methods.
    910
    1011EXAMPLES::
    1112
     
    102103from matrix_integer_sparse cimport Matrix_integer_sparse
    103104from sage.misc.decorators import rename_keyword
    104105
    105 ################
    106 # TODO: change this to use extern cdef's methods.
     106############################################
     107# TODO: Change this to use extern cdef's methods.
    107108from sage.rings.fast_arith cimport arith_int
    108109cdef arith_int ai
    109110ai = arith_int()
    110 ################
     111############################################
    111112
    112113# The 46341 below is because the mod-n sparse code still uses
    113114# int's, even on 64-bit computers.  Improving this is
     
    119120cdef Linbox_modn_sparse linbox
    120121linbox = Linbox_modn_sparse()
    121122
     123from sage.libs.linbox.modular cimport ModDoubleField, ModDoubleFieldElement
     124from sage.libs.linbox.echelonform cimport EchelonFormDomainDouble
     125from sage.libs.linbox.echelonform cimport  BlasMatrix
     126
    122127cdef class Matrix_modn_sparse(matrix_sparse.Matrix_sparse):   
    123128
    124129    ########################################################################
     
    130135    # x * get_unsafe
    131136    # x * __richcmp__    -- always the same
    132137    ########################################################################
    133     def __cinit__(self, parent, entries, copy, coerce):
     138    def __cinit__(self, parent, entries, copy, coerce, bint init_rows=True):
    134139        matrix.Matrix.__init__(self, parent)
    135140
    136141        # allocate memory
     
    142147        p = parent.base_ring().order()
    143148        self.p = p
    144149
    145 
    146150        self.rows = <c_vector_modint*> sage_malloc(nr*sizeof(c_vector_modint))
    147151        if self.rows == NULL:
    148152            raise MemoryError, "error allocating memory for sparse matrix"
    149        
    150         for i from 0 <= i < nr:
    151             init_c_vector_modint(&self.rows[i], p, nc, 0)
    152153
     154        if init_rows:
     155            for i in range(nr):
     156                init_c_vector_modint(&self.rows[i], p, nc, 0)
    153157
    154158    def __dealloc__(self):
    155159        cdef int i
    156         for i from 0 <= i < self._nrows:
     160        for i in range(self._nrows):
    157161            clear_c_vector_modint(&self.rows[i])
    158162        sage_free(self.rows)
    159163
    160     def __init__(self, parent, entries, copy, coerce):
     164    def __init__(self, parent, entries, copy, coerce, bint init_rows=True):
    161165        """
    162166        Create a sparse matrix modulo n.
    163167       
     
    181185        """
    182186        cdef int s, z, p
    183187        cdef Py_ssize_t i, j, k
     188        cdef list entries1
    184189       
    185         cdef void** X
    186 
    187190        if entries is None:
    188191            return
    189192
     
    199202                    set_entry(&self.rows[i], j, z)
    200203        elif isinstance(entries, list):
    201204            # Dense input format
     205            entries1 = entries
    202206            if len(entries) != self._nrows * self._ncols:
    203207                raise TypeError, "list of entries must be a dictionary of (i,j):x or a dense list of n * m elements"
    204             seq = PySequence_Fast(entries,"expected a list")
    205             X = PySequence_Fast_ITEMS(seq)
    206208            k = 0
    207209            R = self._base_ring
    208210            # Get fast access to the entries list.
    209             for i from 0 <= i < self._nrows:
    210                 for  j from 0 <= j < self._ncols:
    211                     z = R(<object>X[k])
     211            for i in range(self._nrows):
     212                for  j in range(self._ncols):
     213                    z = R(entries[k])
    212214                    if z != 0:
    213215                        set_entry(&self.rows[i], j, z)
    214216                    k = k + 1
     
    219221                return
    220222            if self._nrows != self._ncols:
    221223                raise TypeError, "matrix must be square to initialize with a scalar."
    222             for i from 0 <= i < self._nrows:
     224            for i in range(self._nrows):
    223225                set_entry(&self.rows[i], i, s)
    224                
    225226
    226227    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value):
    227228        set_entry(&self.rows[i], j, (<IntegerMod_int> value).ivalue)
     
    290291        cdef Py_ssize_t i, j, k
    291292        d = {}
    292293        cdef IntegerMod_int n
    293         for i from 0 <= i < self._nrows:
    294             for j from 0 <= j < self.rows[i].num_nonzero:
     294        for i in range(self._nrows):
     295            for j in range(self.rows[i].num_nonzero):
    295296                n = IntegerMod_int.__new__(IntegerMod_int)
    296297                IntegerMod_abstract.__init__(n, self._base_ring)
    297298                n.ivalue = self.rows[i].entries[j]
     
    348349            sage: a*c == a*d
    349350            True
    350351        """
    351         cdef Matrix_modn_sparse right, ans
    352         right = _right
     352        cdef Matrix_modn_sparse ans, right = _right
    353353
    354354        cdef c_vector_modint* v
     355        cdef Py_ssize_t i, j, k
     356        cdef int x, y, s
    355357
    356         # Build a table that gives the nonzero positions in each column of right
     358        # Build a table with the nonzero positions in each column of right
    357359        nonzero_positions_in_columns = [set([]) for _ in range(right._ncols)]
    358         cdef Py_ssize_t i, j, k
    359         for i from 0 <= i < right._nrows:
     360        for i in range(right._nrows):
    360361            v = &(right.rows[i])
    361             for j from 0 <= j < right.rows[i].num_nonzero:
     362            for j in range(right.rows[i].num_nonzero):
    362363                nonzero_positions_in_columns[v.positions[j]].add(i)
    363364
    364365        ans = self.new_matrix(self._nrows, right._ncols)
    365366       
    366         # Now do the multiplication, getting each row completely before filling it in.
    367         cdef int x, y, s
    368 
    369         for i from 0 <= i < self._nrows:
     367        # Now do the multiplication, getting each row completely before
     368        # filling it in
     369        for i in range(self._nrows):
    370370            v = &self.rows[i]
    371             for j from 0 <= j < right._ncols:
     371            for j in range(right._ncols):
    372372                s = 0
    373373                c = nonzero_positions_in_columns[j]
    374                 for k from 0 <= k < v.num_nonzero:
     374                for k in range(v.num_nonzero):
    375375                    if v.positions[k] in c:
    376376                        y = get_entry(&right.rows[v.positions[k]], j)
    377377                        x = v.entries[k] * y
     
    379379                set_entry(&ans.rows[i], j, s)
    380380        return ans
    381381
    382 
    383382    def _matrix_times_matrix_dense(self, sage.structure.element.Matrix _right):
    384383        """
    385         Multiply self by the sparse matrix _right, and return the
     384        Multiply ``self`` by the sparse matrix ``_right``, and return the
    386385        result as a dense matrix.
    387386
    388         EXAMPLES:
     387        EXAMPLES::
     388
    389389            sage: a = matrix(GF(10007), 2, [1,2,3,4], sparse=True)
    390390            sage: b = matrix(GF(10007), 2, 3, [1..6], sparse=True)
    391391            sage: a * b
     
    403403            sage: type(a._matrix_times_matrix_dense(a))
    404404            <type 'sage.matrix.matrix_mod2_dense.Matrix_mod2_dense'>
    405405        """
    406         cdef Matrix_modn_sparse right
     406        cdef Matrix_modn_sparse right = _right
    407407        cdef matrix_dense.Matrix_dense ans
    408         right = _right
    409408
    410409        cdef c_vector_modint* v
     410        cdef Py_ssize_t i, j, k
     411        cdef int x, y, s
    411412
    412         # Build a table that gives the nonzero positions in each column of right
     413        # Build a table for the nonzero positions in each column of right
    413414        nonzero_positions_in_columns = [set([]) for _ in range(right._ncols)]
    414         cdef Py_ssize_t i, j, k
    415         for i from 0 <= i < right._nrows:
     415        for i in range(right._nrows):
    416416            v = &(right.rows[i])
    417             for j from 0 <= j < right.rows[i].num_nonzero:
     417            for j in range(right.rows[i].num_nonzero):
    418418                nonzero_positions_in_columns[v.positions[j]].add(i)
    419419
    420420        ans = self.new_matrix(self._nrows, right._ncols, sparse=False)
    421        
     421
    422422        # Now do the multiplication, getting each row completely before filling it in.
    423         cdef int x, y, s
    424        
    425         for i from 0 <= i < self._nrows:
     423        for i in range(self._nrows):
    426424            v = &self.rows[i]
    427             for j from 0 <= j < right._ncols:
     425            for j in range(right._ncols):
    428426                s = 0
    429427                c = nonzero_positions_in_columns[j]
    430                 for k from 0 <= k < v.num_nonzero:
     428                for k in range(v.num_nonzero):
    431429                    if v.positions[k] in c:
    432430                        y = get_entry(&right.rows[v.positions[k]], j)
    433431                        x = v.entries[k] * y
    434432                        s = (s + x)%self.p
    435433                ans.set_unsafe_int(i, j, s)
    436                 #ans._matrix[i][j] = s
    437434        return ans
    438435
    439436    ########################################################################
     
    460457        self.rows[n1] = self.rows[n2]
    461458        self.rows[n2] = tmp
    462459
    463     def _echelon_in_place_classical(self):
     460    def _echelon_in_place_classical(self, density_cutoff=0.05, bint finish_dense_part=True):
    464461        """
    465         Replace self by its reduction to reduced row echelon form.
     462        Replace ``self`` by its reduction to reduced row echelon form.
    466463       
    467         ALGORITHM: We use Gauss elimination, in a slightly intelligent way,
    468         in that we clear each column using a row with the minimum number of
    469         nonzero entries.
    470        
    471         TODO: Implement switching to a dense method when the matrix gets
    472         dense.
     464        ALGORITHM:
     465
     466        We use Gauss elimination, in a slightly intelligent way, in that
     467        we clear each column using a row with the minimum number of nonzero
     468        entries.
     469
     470        If at some step the minimum row density in the remaining bottom right
     471        matrix drops below ``density_cutoff``, we switch to a dense echelon
     472        computation.
    473473        """
    474474        x = self.fetch('in_echelon_form')
    475         if not x is None and x: return  # already known to be in echelon form
     475        if x:
     476            if finish_dense_part:
     477                return self.nrows() # already known to be in echelon form
     478            else:
     479                return
    476480        self.check_mutability()
    477481
    478482        cdef Py_ssize_t i, r, c, min, min_row,  start_row
     
    484488        tm = verbose(caller_name = 'sparse_matrix_pyx matrix_modint echelon')
    485489        do_verb = (get_verbose() >= 2)
    486490
    487         for c from 0 <= c < self._ncols:
     491        cdef Matrix_modn_sparse E
     492        for c in range(self._ncols):
    488493            if do_verb and (c % fifth == 0 and c>0):
    489494                tm = verbose('on column %s of %s'%(c, self._ncols),
    490495                             level = 2,
     
    506511            #endfor
    507512            if min_row != -1:
    508513                r = min_row
    509                 #print "min number of entries in a pivoting row = ", min
     514                row_density = float(min) / self._ncols
     515                if do_verb and (c % fifth == 0 and c>0):
     516                    print "lower row density = ", row_density
     517                if row_density >= density_cutoff:
     518                    if not finish_dense_part:
     519                        return start_row
     520                    E, new_pivots = self._dense_elimination(start_row, do_verb)
     521                    pivots.extend(new_pivots)
     522                    # Now swap underlying matrix entries.
     523                    (E.rows, self.rows) = (self.rows, E.rows)
     524                    # Store/cache result
     525                    self.cache('in_echelon_form',True)
     526                    self.cache('pivots',tuple(pivots))
     527                    return
     528                   
    510529                pivots.append(c)
    511530                # Since we can use row r to clear column c, the
    512531                # entry in position c in row r must be the first nonzero entry.
     
    516535                    scale_c_vector_modint(&self.rows[r], a_inverse)
    517536                self.swap_rows_c(r, start_row)
    518537                sig_on()
    519                 for i from 0 <= i < self._nrows:
     538                for i in range(self._nrows):
    520539                    if i != start_row:
    521540                        b = get_entry(&self.rows[i], c)
    522541                        if b != 0:
     
    530549           
    531550        self.cache('pivots',tuple(pivots))
    532551        self.cache('in_echelon_form',True)
     552        if finish_dense_part:
     553            return self.nrows()
    533554
    534555    def _nonzero_positions_by_row(self, copy=True):
    535556        """
     
    552573            return x
    553574        nzp = []
    554575        cdef Py_ssize_t i, j
    555         for i from 0 <= i < self._nrows:
    556             for j from 0 <= j < self.rows[i].num_nonzero:
     576        for i in range(self._nrows):
     577            for j in range(self.rows[i].num_nonzero):
    557578                nzp.append((i,self.rows[i].positions[j]))
    558579        self.cache('nonzero_positions', nzp)
    559580        if copy:
     
    574595        the four values are nonzero.
    575596       
    576597        INPUT:
    577        
    578        
     598
    579599        -  ``filename`` - either a path or None in which case a
    580600           filename in the current directory is chosen automatically
    581601           (default:None)
     
    628648        setPixel = im.setPixel
    629649        colorExact = im.colorExact
    630650
    631         for i from 0 <= i < self._nrows:
    632             for j from 0 <= j < self.rows[i].num_nonzero:
     651        for i in range(self._nrows):
     652            for j in range(self.rows[i].num_nonzero):
    633653                x = <int>(invblk * self.rows[i].positions[j])
    634654                y = <int>(invblk * i)
    635655                r,g,b = colorComponents( getPixel((x,y)))
     
    642662
    643663    def density(self):
    644664        """
    645         Return the density of self, i.e., the ratio of the number of
    646         nonzero entries of self to the total size of self.
    647        
     665        Return the density of ``self``, i.e., the ratio of the number of
     666        nonzero entries of ``self`` to the total size of ``self``.
    648667
    649668        EXAMPLES::
    650        
     669
    651670            sage: A = matrix(QQ,3,3,[0,1,2,3,0,0,6,7,8],sparse=True)
    652671            sage: A.density()
    653672            2/3
    654            
     673
    655674        Notice that the density parameter does not ensure the density
    656675        of a matrix; it is only an upper bound.
    657        
     676
    658677        ::
    659        
     678
    660679            sage: A = random_matrix(GF(127),200,200,density=0.3, sparse=True)
    661680            sage: A.density()
    662681            2073/8000
    663682        """
    664         cdef Py_ssize_t i, nonzero_entries
     683        cdef Py_ssize_t i, nz = 0
    665684
    666         nonzero_entries = 0
    667         for i from 0 <= i < self._nrows:
    668             nonzero_entries += self.rows[i].num_nonzero
     685        for i in range(self._nrows):
     686            nz += self.rows[i].num_nonzero
    669687           
    670         return rings.ZZ(nonzero_entries)/rings.ZZ(self._nrows*self._ncols)
     688        return rings.ZZ(nz) / rings.ZZ(self._nrows*self._ncols)
    671689
    672690    def transpose(self):
    673691        """
    674         Return the transpose of self.
    675        
     692        Return the transpose of ``self``.
     693
    676694        EXAMPLE::
    677        
     695
    678696            sage: A = matrix(GF(127),3,3,[0,1,0,2,0,0,3,0,0],sparse=True)
    679697            sage: A
    680698            [0 1 0]
     
    695713        cdef int i, j
    696714        cdef c_vector_modint row
    697715        cdef Matrix_modn_sparse B
    698        
     716
    699717        B = self.new_matrix(nrows = self.ncols(), ncols = self.nrows())
    700         for i from 0 <= i < self._nrows:
     718        for i in range(self._nrows):
    701719            row = self.rows[i]
    702             for j from 0 <= j < row.num_nonzero:
     720            for j in range(row.num_nonzero):
    703721                set_entry(&B.rows[row.positions[j]], i, row.entries[j])
    704722        if self._subdivisions is not None:
    705723            row_divs, col_divs = self.subdivisions()
    706724            B.subdivide(col_divs, row_divs)
    707725        return B
    708                                
     726
    709727    def matrix_from_rows(self, rows):
    710728        """
    711         Return the matrix constructed from self using rows with indices in
     729        Return the matrix constructed from ``self`` using rows with indices in
    712730        the rows list.
    713        
     731
    714732        INPUT:
    715        
    716        
     733
    717734        -  ``rows`` - list or tuple of row indices
    718        
    719        
     735
    720736        EXAMPLE::
    721        
     737
    722738            sage: M = MatrixSpace(GF(127),3,3,sparse=True)
    723739            sage: A = M(range(9)); A
    724740            [0 1 2]
     
    728744            [6 7 8]
    729745            [3 4 5]
    730746        """
    731         cdef int i,k
     747        cdef Py_ssize_t i, k
    732748        cdef Matrix_modn_sparse A
    733         cdef c_vector_modint row
    734749
    735750        if not isinstance(rows, (list, tuple)):
    736751            raise TypeError, "rows must be a list of integers"
    737752
    738 
    739753        A = self.new_matrix(nrows = len(rows))
    740754
    741         k = 0
    742         for ii in rows:
    743             i = ii
     755        for k, i in enumerate(rows):
    744756            if i < 0 or i >= self.nrows():
    745757                raise IndexError, "row %s out of range"%i
    746758
    747             row = self.rows[i]
    748             for j from 0 <= j < row.num_nonzero:
    749                 set_entry(&A.rows[k], row.positions[j], row.entries[j])
    750             k += 1
     759            set_entries(&A.rows[k], &self.rows[i])
     760
    751761        return A
    752            
    753762
    754763    def matrix_from_columns(self, cols):
    755764        """
     
    768777            [5 4]
    769778            [8 7]
    770779        """
    771         cdef int i,j
     780        cdef Py_ssize_t i, j
    772781        cdef Matrix_modn_sparse A
    773782        cdef c_vector_modint row
    774783
     
    779788
    780789        cols = dict(zip([int(e) for e in cols],range(len(cols))))
    781790
    782         for i from 0 <= i < self.nrows():
     791        for i in range(self._nrows):
    783792            row = self.rows[i]
    784             for j from 0 <= j < row.num_nonzero:
     793            for j in range(row.num_nonzero):
    785794                if int(row.positions[j]) in cols:
    786795                    set_entry(&A.rows[i], cols[int(row.positions[j])], row.entries[j])
    787796        return A
     
    814823
    815824    def rank(self, gauss=False):
    816825        """
    817         Compute the rank of self.
    818        
     826        Compute the rank of ``self``.
     827
    819828        INPUT:
    820        
    821        
     829
    822830        -  ``gauss`` - if True LinBox' Gaussian elimination is
    823831           used. If False 'Symbolic Reordering' as implemented in LinBox is
    824832           used. If 'native' the native Sage implementation is used. (default:
    825833           False)
    826        
    827        
     834
    828835        EXAMPLE::
    829        
     836
    830837            sage: A = random_matrix(GF(127),200,200,density=0.01,sparse=True)
    831838            sage: r1 = A.rank(gauss=False)
    832839            sage: r2 = A.rank(gauss=True)
     
    835842            True
    836843            sage: r1
    837844            155
    838        
     845
    839846        ALGORITHM: Uses LinBox or native implementation.
    840        
     847
    841848        REFERENCES:
    842849
    843850        - Jean-Guillaume Dumas and Gilles Villars. 'Computing the Rank
     
    846853          on Computer Algebra in Scientific Computing, Big Yalta,
    847854          Crimea, Ukraine, 22-27 sept. 2002, Springer-Verlag,
    848855          http://perso.ens-lyon.fr/gilles.villard/BIBLIOGRAPHIE/POSTSCRIPT/rankjgd.ps
    849        
     856
    850857        .. note::
    851858
    852859           For very sparse matrices Gaussian elimination is faster
     
    873880
    874881    def _solve_right_nonsingular_square(self, B, algorithm=None, check_rank = True):
    875882        """
    876         If self is a matrix `A`, then this function returns a
    877         vector or matrix `X` such that `A X = B`. If
    878         `B` is a vector then `X` is a vector and if
    879         `B` is a matrix, then `X` is a matrix.
    880        
     883        If ``self`` is a matrix `A`, then this function returns a vector
     884        or matrix `X` such that `A X = B`. If `B` is a vector then `X` is
     885        a vector and if `B` is a matrix, then `X` is a matrix.
     886
    881887        .. note::
    882888
    883            In Sage one can also write ``A  B`` for
    884            ``A.solve_right(B)``, i.e., Sage implements the "the
    885            MATLAB/Octave backslash operator".
    886        
     889           In Sage one can also write ``A  B`` for ``A.solve_right(B)``,
     890           i.e., Sage implements the "the MATLAB/Octave backslash operator".
     891
    887892        INPUT:
    888        
    889        
     893
    890894        -  ``B`` - a matrix or vector
    891895       
    892896        -  ``algorithm`` - one of the following:
     
    905909       
    906910        -  ``check_rank`` - check rank before attempting to
    907911           solve (default: True)
    908        
    909        
     912
    910913        OUTPUT: a matrix or vector
    911        
     914
    912915        EXAMPLES::
    913        
     916
    914917            sage: A = matrix(GF(127), 3, [1,2,3,-1,2,5,2,3,1], sparse=True)
    915918            sage: b = vector(GF(127),[1,2,3])
    916919            sage: x = A \ b; x
     
    982985        """
    983986        Return lift of this matrix to a sparse matrix over the integers.
    984987
    985         EXAMPLES:
     988        EXAMPLES::
     989
    986990            sage: a = matrix(GF(7),2,3,[1..6], sparse=True)
    987991            sage: a.lift()
    988992            [1 2 3]
     
    10071011
    10081012        cdef mpz_vector* L_row
    10091013        cdef c_vector_modint* A_row
    1010         for i from 0 <= i < self._nrows:
    1011             L_row = &(L._matrix[i])
     1014        for i in range(self._nrows):
     1015            # The evil cast here is to get around #included pxi's and C++ rules.  In C mode
     1016            # this compiles fine:
     1017            L_row = <mpz_vector*> (&(L._matrix[i])) 
    10121018            A_row = &(self.rows[i])
    10131019            sage_free(L_row.entries)
    10141020            L_row.entries = <mpz_t*> sage_malloc(sizeof(mpz_t)*A_row.num_nonzero)
     
    10211027                sage_free(L_row.entries)
    10221028                L_row.entries = NULL
    10231029                raise MemoryError, "error allocating space for sparse vector during sparse lift"
    1024             for j from 0 <= j < A_row.num_nonzero:
     1030            for j in range(A_row.num_nonzero):
    10251031                L_row.positions[j] = A_row.positions[j]
    10261032                mpz_init_set_si(L_row.entries[j], A_row.entries[j])
    10271033        L.subdivide(self.subdivisions())
    10281034        return L
    10291035
     1036    def _dense_elimination(self, Py_ssize_t r, bint verb):
     1037        """
     1038        Dense part of the elimination, using the obvious algorithm, which took
     1039        Sebastian Pancratz and William Stein an hour to work out :-).
     1040
     1041        INPUT:
     1042       
     1043            - Partially-echelonized matrix ``self``
     1044            - Index `r` of first dense row
     1045            - Boolean ``verb`` to trigger verbose output
     1046
     1047        OUTPUT:
     1048
     1049            - Sparse matrix in reduced row echelon form
     1050            - List of pivot columns with index at least `c`,
     1051              where `c` is the index of the first dense column
     1052
     1053        TODO:
     1054
     1055            - Include ``sig_on``, ``sig_off``
     1056
     1057        ALGORITHM:
     1058
     1059                                   c
     1060            -------------------------------------
     1061            |                     |             |
     1062            | Already echelonized |             |
     1063            | part of the sparse  |     D       |
     1064            | matrix              |             |
     1065            |                     |             |
     1066            -------------------------------------
     1067          r |                     |             |
     1068            |         0           |    A, E     |
     1069            |                     |             |
     1070            -------------------------------------
     1071
     1072            `A` is set to ``self[r:,c:]``, and `E` is the echelon form of
     1073            the matrix `A`.
     1074
     1075            `D` is at first set to a dense copy of ``self[0:r,c:]`` but
     1076            then corrected so that it is the appropriate part of the
     1077            echelon form of ``self`` using the matrix `E`.
     1078
     1079        EXAMPLES::
     1080
     1081            sage: p = previous_prime(2**15)
     1082            sage: A = matrix(GF(p), 3, 3, [1,1,2,0,3,15,0,1,5], sparse=True); A
     1083            [ 1  1  2]
     1084            [ 0  3 15]
     1085            [ 0  1  5]
     1086            sage: E,P = A._dense_elimination(1,0);
     1087            sage: print E, P
     1088            [    1     0 32746]
     1089            [    0     1     5]
     1090            [    0     0     0] [1]
     1091        """
     1092        if self._nrows == 0 or self._ncols == 0:
     1093            return self, []
     1094
     1095        if r < 0 or r >= self._nrows:
     1096            raise IndexError
     1097
     1098        # First dense column
     1099        cdef Py_ssize_t c
     1100
     1101        if r == 0:
     1102            c = 0
     1103        else:
     1104            if self.rows[r-1].num_nonzero == 0:
     1105                raise ValueError, "Row r-1 has no nonzero entries"
     1106            c = self.rows[r-1].positions[0] + 1
     1107
     1108        if verb:
     1109            print "r = %s, c = %s"%(r,c)
     1110
     1111        if (r == self._nrows or c == self._ncols):
     1112            return self, []
     1113
     1114        # Dimensions of the bottom right matrix
     1115        cdef Py_ssize_t m = self._nrows - r
     1116        cdef Py_ssize_t n = self._ncols - c
     1117
     1118        cdef ModDoubleField *F = new ModDoubleField(self.p)
     1119        cdef BlasMatrix[ModDoubleFieldElement] *A = \
     1120            new BlasMatrix[ModDoubleFieldElement](m, n)
     1121
     1122        cdef Py_ssize_t i, j, k, nz
     1123
     1124        if verb:
     1125            print "Copying dense matrix over to Linbox..."
     1126        for i in range(r, self._nrows):
     1127            nz = self.rows[i].num_nonzero
     1128            for k in range(nz):
     1129               A.setEntry(i-r, self.rows[i].positions[k]-c,
     1130                          <ModDoubleFieldElement> self.rows[i].entries[k])
     1131
     1132        if verb:
     1133            print "Echelonizing %sx%s matrix in linbox..."%(m,n)
     1134        cdef EchelonFormDomainDouble *EF = new EchelonFormDomainDouble(F[0])
     1135        cdef BlasMatrix[ModDoubleFieldElement] *E = new BlasMatrix[ModDoubleFieldElement](m, n)
     1136        cdef int rk = EF.rowReducedEchelon(E[0], A[0])
     1137        if verb:
     1138            print "Rank =", rk
     1139
     1140        # Compute pivots and non-pivots
     1141        cdef Py_ssize_t j0 = 0
     1142        cdef list pivots = [], non_pivots = []
     1143        for i in range(rk):
     1144            for j in range(j0, n):
     1145                if not F.isZero(E.getEntry(i, j)):
     1146                    non_pivots.extend(range(pivots[-1]+1 if pivots else 0, j))
     1147                    pivots.append(j)
     1148                    j0 = j+1
     1149                    break
     1150        non_pivots.extend(range(pivots[-1]+1 if pivots else 0, n))
     1151        if verb:
     1152            print "Number of non-pivot columns: %s"%len(non_pivots)
     1153        assert set(pivots + non_pivots) == set(range(n))
     1154
     1155        if verb:
     1156            print "Copying non-pivot columns into a dense matrix"
     1157        cdef long* D = <long*> sage_calloc(n*r, sizeof(long))
     1158
     1159        if D == NULL:
     1160            del F, A, E, EF
     1161            raise MemoryError, "Failed to allocate the dense matrix D."
     1162
     1163        cdef Py_ssize_t k0
     1164        for i in range(r):
     1165            k0 = 0
     1166            while k0 < self.rows[i].num_nonzero and self.rows[i].positions[k0] < c:
     1167                k0 += 1
     1168            for k in range(k0, self.rows[i].num_nonzero):
     1169                D[(self.rows[i].positions[k]-c)*r + i] = self.rows[i].entries[k]
     1170
     1171        if verb:
     1172            print "Doing the elimination step"
     1173        cdef Py_ssize_t p
     1174        cdef long t, a
     1175
     1176        for i in range(rk):
     1177            p = pivots[i]
     1178            for j in range(p+1, n):
     1179                t = <long> E.getEntry(i,j)
     1180                if t:
     1181                    for k in range(r):
     1182                        a = (D[j*r + k] - D[p*r + k]*t) % self.p
     1183                        if a < 0:
     1184                            a += self.p
     1185                        D[j*r + k] = a
     1186
     1187        if verb:
     1188            print "Constructing final matrix"
     1189
     1190        cdef Matrix_modn_sparse Z = Matrix_modn_sparse(self._parent,
     1191                                    None, False, False, init_rows=False)
     1192
     1193        if verb:
     1194            print "Auxilary vector"
     1195
     1196        cdef c_vector_modint v
     1197        init_c_vector_modint(&v, self.p, self._ncols, self._ncols)
     1198       
     1199        if verb:
     1200            print "Copying over top part of matrix"
     1201
     1202        for i in range(r):
     1203            nz = 0
     1204            for k in range(self.rows[i].num_nonzero):
     1205                if self.rows[i].positions[k] < c:
     1206                    v.positions[nz] = self.rows[i].positions[k]
     1207                    v.entries[nz] = self.rows[i].entries[k]
     1208                    nz += 1
     1209                else:
     1210                    break
     1211            for j in range(len(non_pivots)):
     1212                if D[non_pivots[j]*r + i]:
     1213                    v.positions[nz] = c+non_pivots[j]
     1214                    v.entries[nz] = D[non_pivots[j]*r + i]
     1215                    nz += 1
     1216
     1217            init_c_vector_modint(&Z.rows[i], self.p, self._ncols, nz)
     1218            nz -= 1
     1219            while nz >= 0:
     1220                Z.rows[i].positions[nz] = v.positions[nz]
     1221                Z.rows[i].entries[nz] = v.entries[nz]
     1222                nz -= 1
     1223
     1224        if verb:
     1225            print "Copying over bottom right matrix"
     1226
     1227        for i in range(r, self._nrows):
     1228            nz = 0
     1229            for j in range(n):
     1230                if <long>E.getEntry(i-r, j):
     1231                    nz += 1
     1232            init_c_vector_modint(&Z.rows[i], self.p, self._ncols, nz)
     1233            nz = 0
     1234            for j in range(n):
     1235                if <long>E.getEntry(i-r, j):
     1236                    Z.rows[i].positions[nz] = j+c
     1237                    Z.rows[i].entries[nz] = <long>E.getEntry(i-r, j)
     1238                    nz += 1
     1239
     1240        if verb:
     1241            print "Free memory allocated above"
     1242
     1243        del F, A, E, EF
     1244        sage_free(D)
     1245        clear_c_vector_modint(&v)
     1246
     1247        return Z, [i+c for i in pivots]
     1248
     1249    def _subtract_scalar(self, a, bint inplace=False):
     1250        """
     1251        Returns the matrix ``self`` minus `a`.
     1252
     1253        This has the same entries as ``self`` for off-diagonal entries
     1254        but entries on the diagonal have `a` subtracted from them.
     1255
     1256        INPUT:
     1257
     1258            -  ``a`` -- Scalar;  an integer, int, long, or something that
     1259                        can be coerced into the base ring
     1260
     1261            -  ``inplace`` --  If ``True``, modifies ``self`` in-place
     1262
     1263        OUTPUT:
     1264
     1265            -  ``self`` minus `a`
     1266
     1267        EXAMPLES::
     1268
     1269            sage: p = previous_prime(2**15)
     1270            sage: a = matrix(GF(p), 3, 3, [1,1,2,0,3,15,0,1,5], sparse=True)
     1271            sage: a._subtract_scalar(5)
     1272            [32745     1     2]
     1273            [    0 32747    15]
     1274            [    0     1     0]
     1275        """
     1276        cdef Py_ssize_t m = self._nrows
     1277        cdef Py_ssize_t n = self._ncols
     1278
     1279        cdef Py_ssize_t i
     1280        cdef int b, p = self.p, x
     1281
     1282        cdef Integer y
     1283        cdef IntegerMod_int z
     1284
     1285        if PY_TYPE_CHECK(a, Integer):
     1286            y = a % Integer(p)
     1287            b = mpz_get_si(y.value)
     1288        elif PyInt_Check(a):
     1289            b = a % p
     1290        elif PyLong_Check(a):
     1291            b = a % <long> p
     1292        else:
     1293            # Coerce into the base ring
     1294            R = self._base_ring
     1295            z = R(a)
     1296            b = z.ivalue
     1297
     1298        if b < 0:
     1299            b = b + p
     1300
     1301        if inplace:
     1302
     1303            self.check_mutability()
     1304
     1305            for i in range(m):
     1306                if i >= n:
     1307                    break
     1308                x = get_entry(&self.rows[i], i) - b
     1309                set_entry(&self.rows[i], i, x)
     1310            return self
     1311
     1312        else:
     1313            A = self.__copy__()
     1314            return A._subtract_scalar(a, inplace=True)
     1315
  • sage/modules/vector_modn_sparse_c.pxi

    diff --git a/sage/modules/vector_modn_sparse_c.pxi b/sage/modules/vector_modn_sparse_c.pxi
    a b  
    1313    """
    1414    v.entries = <int*>sage_malloc(num_nonzero*sizeof(int))
    1515    if v.entries == NULL:
    16         raise MemoryError, "Error allocating memory"
     16        raise MemoryError, "Error allocating memory for sparse vector."
    1717    v.positions = <Py_ssize_t*>sage_malloc(num_nonzero*sizeof(Py_ssize_t))
    1818    if v.positions == NULL:
    1919        sage_free(v.entries)
    20         raise MemoryError, "Error allocating memory"
     20        raise MemoryError, "Error allocating memory for sparse vector."
    2121    return 0
    2222
    2323cdef int init_c_vector_modint(c_vector_modint* v, int p, Py_ssize_t degree,
     
    3535    v.p = p
    3636    return 0
    3737
     38cdef int reallocate_c_vector_modint(c_vector_modint* v, Py_ssize_t num_nonzero) except -1:
     39    """
     40    Reallocate memory for a c_vector_modint.
     41    """
     42    if num_nonzero == 0:
     43        sage_free(v.entries)
     44        sage_free(v.positions)
     45        v.entries   = NULL
     46        v.positions = NULL
     47        return 0
     48
     49    cdef int *e
     50    cdef Py_ssize_t *p
     51
     52    e = <int*>         sage_realloc(v.entries,   num_nonzero*sizeof(int))
     53    p = <Py_ssize_t *> sage_realloc(v.positions, num_nonzero*sizeof(Py_ssize_t))
     54
     55    if e == NULL or p == NULL:
     56        sage_free(v.entries   if e == NULL else e)
     57        sage_free(v.positions if p == NULL else p)
     58        raise MemoryError, "Error allocating memory for sparse vector."
     59
     60    v.entries   = e
     61    v.positions = p
     62    return 0
     63
    3864cdef void clear_c_vector_modint(c_vector_modint* v):
    3965    sage_free(v.entries)
    4066    sage_free(v.positions)
     
    5177    i = 0
    5278    j = n-1
    5379    while i<=j:
    54         if i == j:
    55             if v[i] == x:
    56                 return i
    57             return -1
    5880        k = (i+j)/2
    5981        if v[k] > x:
    6082            j = k-1
    6183        elif v[k] < x:
    6284            i = k+1
    63         else:   # only possibility is that v[k] == x
     85        else:  # v[k] == x
    6486            return k
    6587    return -1
    6688
     
    93115            j = k-1
    94116        elif v[k] < x:
    95117            i = k+1
    96         else:   # only possibility is that v[k] == x
     118        else:  # v[k] == x
    97119            ins[0] = k
    98120            return k
    99     # end while
    100121    ins[0] = j+1
    101122    return -1
    102123
     
    107128    """
    108129    if n >= v.degree or n < 0:
    109130        raise IndexError, "Index must be between 0 and the degree minus 1."
    110         return -1
    111131    cdef Py_ssize_t m
    112132    m = binary_search0_modn(v.positions, v.num_nonzero, n)
    113133    if m == -1:
     
    133153    """
    134154    if n < 0 or n >= v.degree:
    135155        raise IndexError, "Index (=%s) must be between 0 and %s."%(n, v.degree-1)
    136         return -1   
    137156    cdef Py_ssize_t i, m, ins
    138     cdef Py_ssize_t m2, ins2
    139     cdef Py_ssize_t *pos
    140     cdef int *e
    141157
    142158    x = x % v.p
    143159    if x<0: x = x + v.p
    144160    m = binary_search_modn(v.positions, v.num_nonzero, n, &ins)
    145    
    146     if m != -1:
    147         # The position n was found in the array of positions.
    148         # Now there are two cases:
    149         #   1. x =/= 0, which is easy, and
    150         #   2. x = 0, in which case we have to recopy
    151         #      positions and entries, without the m-th
    152         #      element, and change num_nonzero.
    153         if x != 0:   # case 1
     161
     162    if m != -1:      # Found position n in the vector
     163        if x != 0:   # Replace with non-zero x
    154164            v.entries[m] = x
    155         else:        # case 2
    156             e = v.entries
    157             pos = v.positions
    158             allocate_c_vector_modint(v, v.num_nonzero - 1)
    159             for i from 0 <= i < m:
    160                 v.entries[i] = e[i]
    161                 v.positions[i] = pos[i]
    162             for i from m < i < v.num_nonzero:
    163                 v.entries[i-1] = e[i]
    164                 v.positions[i-1] = pos[i]
    165             sage_free(e)
    166             sage_free(pos)
     165        else:        # Remove entry as x == 0
    167166            v.num_nonzero = v.num_nonzero - 1
     167
     168            for i in range(m, v.num_nonzero):
     169                v.entries[i]   = v.entries[i+1]
     170                v.positions[i] = v.positions[i+1]
     171
     172            reallocate_c_vector_modint(v, v.num_nonzero)
    168173    else:
    169         # Allocate new memory and copy over elements from the
    170         # old array.  This is similar to case 2 above,
    171         # except we are inserting a new entry rather than
    172         # deleting an old one.  The new entry should be inserted
    173         # at position ins, which was computed using binary search.
    174         #
    175         # There is one exception -- if the new entry is 0, we
    176         # do nothing and return.
    177         if x == 0:
    178             return 0
    179         v.num_nonzero = v.num_nonzero + 1
    180         e = v.entries
    181         pos = v.positions
    182         allocate_c_vector_modint(v, v.num_nonzero)
    183         for i from 0 <= i < ins:
    184             v.entries[i] = e[i]
    185             v.positions[i] = pos[i]
    186         v.entries[ins] = x
    187         v.positions[ins] = n
    188         for i from ins < i < v.num_nonzero:
    189             v.entries[i] = e[i-1]
    190             v.positions[i] = pos[i-1]
    191         sage_free(e)
    192         sage_free(pos)
     174        # Did not find position n, push higher entries up by one
     175        # and insert new value x
     176        if x != 0:
     177            v.num_nonzero = v.num_nonzero + 1
     178            reallocate_c_vector_modint(v, v.num_nonzero)
     179
     180            for i in range(v.num_nonzero - 1, ins, -1):
     181                v.entries[i]   = v.entries[i-1]
     182                v.positions[i] = v.positions[i-1]
     183
     184            v.entries[ins]   = x
     185            v.positions[ins] = n
     186
     187cdef int set_entries(c_vector_modint* v1, c_vector_modint* v2) except -1:
     188    """
     189    Sets the vector ``v1`` to the value of the vector ``v2``.
     190    """
     191    if v1.p != v2.p:
     192        raise ArithmeticError, "The vectors must be modulo the same prime."
     193    if v1.degree != v2.degree:
     194        raise ArithmeticError, "The vectors must have the same degree."
     195
     196    if v2.num_nonzero == 0:
     197        v1.num_nonzero = 0
     198        sage_free(v1.positions)
     199        sage_free(v1.entries)
     200        v1.positions = NULL
     201        v1.entries   = NULL
     202        return 0
     203
     204    cdef Py_ssize_t i
     205
     206    reallocate_c_vector_modint(v1, v2.num_nonzero)
     207
     208    for i in range(v2.num_nonzero):
     209        v1.entries[i]   = v2.entries[i]
     210        v1.positions[i] = v2.positions[i]
     211    v1.num_nonzero = v2.num_nonzero
     212    return 0
    193213
    194214cdef int add_c_vector_modint_init(c_vector_modint* sum, c_vector_modint* v,
    195215                                  c_vector_modint* w, int multiple) except -1:
     
    198218    """
    199219    if v.p != w.p:
    200220        raise ArithmeticError, "The vectors must be modulo the same prime."
    201         return -1
    202221    if v.degree != w.degree:
    203222        raise ArithmeticError, "The vectors must have the same degree."
    204         return -1
    205223
    206224    cdef int s
    207225    cdef Py_ssize_t nz, i, j, k
     
    276294    if scalar < 0:
    277295        scalar = scalar + v.p
    278296    cdef Py_ssize_t i
    279     for i from 0 <= i < v.num_nonzero:
     297    for i in range(v.num_nonzero):
    280298        v.entries[i] = (v.entries[i] * scalar) % v.p
    281299    return 0
    282300
    283