Ticket #13177: trac_13177.patch
File trac_13177.patch, 39.4 KB (added by , 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 977 977 978 978 Extension('sage.matrix.matrix_modn_sparse', 979 979 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']), 981 983 982 984 Extension('sage.matrix.matrix_mpolynomial_dense', 983 985 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 1 1 include "../../ext/cdefs.pxi" 2 2 3 include '../../modules/vector_modn_sparse_h.pxi' 3 from sage.matrix.matrix_modn_sparse cimport c_vector_modint 4 4 5 5 from sage.matrix.matrix_integer_dense cimport mod_int 6 6 
sage/matrix/matrix2.pyx
diff git a/sage/matrix/matrix2.pyx b/sage/matrix/matrix2.pyx
a b 5714 5714 sage: m == transformation_matrix * m_original 5715 5715 True 5716 5716 """ 5717 if self.fetch('in_echelon_form') and not kwds.get('transformation',False): 5718 return 5717 5719 self.check_mutability() 5718 5720 5719 5721 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 1 1 r""" 2 Sparse matrices over `\ZZ/n\ZZ` for `n` small2 Sparse matrices over `\ZZ/n\ZZ` for small `n`, currently at most 46341. 3 3 4 4 This is a compiled implementation of sparse matrices over 5 `\ZZ/n\ZZ` for `n` small.5 `\ZZ/n\ZZ` for small `n`. 6 6 7 TODO:  move vectors into a Cython vector class  add _add_ and 8 _mul_ methods. 7 TODO: 8  Move vectors into a Cython vector class 9  Add _add_ and _mul_ methods. 9 10 10 11 EXAMPLES:: 11 12 … … 102 103 from matrix_integer_sparse cimport Matrix_integer_sparse 103 104 from sage.misc.decorators import rename_keyword 104 105 105 ################ 106 # TODO: change this to use extern cdef's methods.106 ############################################ 107 # TODO: Change this to use extern cdef's methods. 107 108 from sage.rings.fast_arith cimport arith_int 108 109 cdef arith_int ai 109 110 ai = arith_int() 110 ################ 111 ############################################ 111 112 112 113 # The 46341 below is because the modn sparse code still uses 113 114 # int's, even on 64bit computers. Improving this is … … 119 120 cdef Linbox_modn_sparse linbox 120 121 linbox = Linbox_modn_sparse() 121 122 123 from sage.libs.linbox.modular cimport ModDoubleField, ModDoubleFieldElement 124 from sage.libs.linbox.echelonform cimport EchelonFormDomainDouble 125 from sage.libs.linbox.echelonform cimport BlasMatrix 126 122 127 cdef class Matrix_modn_sparse(matrix_sparse.Matrix_sparse): 123 128 124 129 ######################################################################## … … 130 135 # x * get_unsafe 131 136 # x * __richcmp__  always the same 132 137 ######################################################################## 133 def __cinit__(self, parent, entries, copy, coerce ):138 def __cinit__(self, parent, entries, copy, coerce, bint init_rows=True): 134 139 matrix.Matrix.__init__(self, parent) 135 140 136 141 # allocate memory … … 142 147 p = parent.base_ring().order() 143 148 self.p = p 144 149 145 146 150 self.rows = <c_vector_modint*> sage_malloc(nr*sizeof(c_vector_modint)) 147 151 if self.rows == NULL: 148 152 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)152 153 154 if init_rows: 155 for i in range(nr): 156 init_c_vector_modint(&self.rows[i], p, nc, 0) 153 157 154 158 def __dealloc__(self): 155 159 cdef int i 156 for i from 0 <= i < self._nrows:160 for i in range(self._nrows): 157 161 clear_c_vector_modint(&self.rows[i]) 158 162 sage_free(self.rows) 159 163 160 def __init__(self, parent, entries, copy, coerce ):164 def __init__(self, parent, entries, copy, coerce, bint init_rows=True): 161 165 """ 162 166 Create a sparse matrix modulo n. 163 167 … … 181 185 """ 182 186 cdef int s, z, p 183 187 cdef Py_ssize_t i, j, k 188 cdef list entries1 184 189 185 cdef void** X186 187 190 if entries is None: 188 191 return 189 192 … … 199 202 set_entry(&self.rows[i], j, z) 200 203 elif isinstance(entries, list): 201 204 # Dense input format 205 entries1 = entries 202 206 if len(entries) != self._nrows * self._ncols: 203 207 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)206 208 k = 0 207 209 R = self._base_ring 208 210 # 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]) 212 214 if z != 0: 213 215 set_entry(&self.rows[i], j, z) 214 216 k = k + 1 … … 219 221 return 220 222 if self._nrows != self._ncols: 221 223 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): 223 225 set_entry(&self.rows[i], i, s) 224 225 226 226 227 cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, value): 227 228 set_entry(&self.rows[i], j, (<IntegerMod_int> value).ivalue) … … 290 291 cdef Py_ssize_t i, j, k 291 292 d = {} 292 293 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): 295 296 n = IntegerMod_int.__new__(IntegerMod_int) 296 297 IntegerMod_abstract.__init__(n, self._base_ring) 297 298 n.ivalue = self.rows[i].entries[j] … … 348 349 sage: a*c == a*d 349 350 True 350 351 """ 351 cdef Matrix_modn_sparse right, ans 352 right = _right 352 cdef Matrix_modn_sparse ans, right = _right 353 353 354 354 cdef c_vector_modint* v 355 cdef Py_ssize_t i, j, k 356 cdef int x, y, s 355 357 356 # Build a table that givesthe nonzero positions in each column of right358 # Build a table with the nonzero positions in each column of right 357 359 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): 360 361 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): 362 363 nonzero_positions_in_columns[v.positions[j]].add(i) 363 364 364 365 ans = self.new_matrix(self._nrows, right._ncols) 365 366 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): 370 370 v = &self.rows[i] 371 for j from 0 <= j < right._ncols:371 for j in range(right._ncols): 372 372 s = 0 373 373 c = nonzero_positions_in_columns[j] 374 for k from 0 <= k < v.num_nonzero:374 for k in range(v.num_nonzero): 375 375 if v.positions[k] in c: 376 376 y = get_entry(&right.rows[v.positions[k]], j) 377 377 x = v.entries[k] * y … … 379 379 set_entry(&ans.rows[i], j, s) 380 380 return ans 381 381 382 383 382 def _matrix_times_matrix_dense(self, sage.structure.element.Matrix _right): 384 383 """ 385 Multiply self by the sparse matrix _right, and return the384 Multiply ``self`` by the sparse matrix ``_right``, and return the 386 385 result as a dense matrix. 387 386 388 EXAMPLES: 387 EXAMPLES:: 388 389 389 sage: a = matrix(GF(10007), 2, [1,2,3,4], sparse=True) 390 390 sage: b = matrix(GF(10007), 2, 3, [1..6], sparse=True) 391 391 sage: a * b … … 403 403 sage: type(a._matrix_times_matrix_dense(a)) 404 404 <type 'sage.matrix.matrix_mod2_dense.Matrix_mod2_dense'> 405 405 """ 406 cdef Matrix_modn_sparse right 406 cdef Matrix_modn_sparse right = _right 407 407 cdef matrix_dense.Matrix_dense ans 408 right = _right409 408 410 409 cdef c_vector_modint* v 410 cdef Py_ssize_t i, j, k 411 cdef int x, y, s 411 412 412 # Build a table that givesthe nonzero positions in each column of right413 # Build a table for the nonzero positions in each column of right 413 414 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): 416 416 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): 418 418 nonzero_positions_in_columns[v.positions[j]].add(i) 419 419 420 420 ans = self.new_matrix(self._nrows, right._ncols, sparse=False) 421 421 422 422 # 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): 426 424 v = &self.rows[i] 427 for j from 0 <= j < right._ncols:425 for j in range(right._ncols): 428 426 s = 0 429 427 c = nonzero_positions_in_columns[j] 430 for k from 0 <= k < v.num_nonzero:428 for k in range(v.num_nonzero): 431 429 if v.positions[k] in c: 432 430 y = get_entry(&right.rows[v.positions[k]], j) 433 431 x = v.entries[k] * y 434 432 s = (s + x)%self.p 435 433 ans.set_unsafe_int(i, j, s) 436 #ans._matrix[i][j] = s437 434 return ans 438 435 439 436 ######################################################################## … … 460 457 self.rows[n1] = self.rows[n2] 461 458 self.rows[n2] = tmp 462 459 463 def _echelon_in_place_classical(self ):460 def _echelon_in_place_classical(self, density_cutoff=0.05, bint finish_dense_part=True): 464 461 """ 465 Replace selfby its reduction to reduced row echelon form.462 Replace ``self`` by its reduction to reduced row echelon form. 466 463 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. 473 473 """ 474 474 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 476 480 self.check_mutability() 477 481 478 482 cdef Py_ssize_t i, r, c, min, min_row, start_row … … 484 488 tm = verbose(caller_name = 'sparse_matrix_pyx matrix_modint echelon') 485 489 do_verb = (get_verbose() >= 2) 486 490 487 for c from 0 <= c < self._ncols: 491 cdef Matrix_modn_sparse E 492 for c in range(self._ncols): 488 493 if do_verb and (c % fifth == 0 and c>0): 489 494 tm = verbose('on column %s of %s'%(c, self._ncols), 490 495 level = 2, … … 506 511 #endfor 507 512 if min_row != 1: 508 513 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 510 529 pivots.append(c) 511 530 # Since we can use row r to clear column c, the 512 531 # entry in position c in row r must be the first nonzero entry. … … 516 535 scale_c_vector_modint(&self.rows[r], a_inverse) 517 536 self.swap_rows_c(r, start_row) 518 537 sig_on() 519 for i from 0 <= i < self._nrows:538 for i in range(self._nrows): 520 539 if i != start_row: 521 540 b = get_entry(&self.rows[i], c) 522 541 if b != 0: … … 530 549 531 550 self.cache('pivots',tuple(pivots)) 532 551 self.cache('in_echelon_form',True) 552 if finish_dense_part: 553 return self.nrows() 533 554 534 555 def _nonzero_positions_by_row(self, copy=True): 535 556 """ … … 552 573 return x 553 574 nzp = [] 554 575 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): 557 578 nzp.append((i,self.rows[i].positions[j])) 558 579 self.cache('nonzero_positions', nzp) 559 580 if copy: … … 574 595 the four values are nonzero. 575 596 576 597 INPUT: 577 578 598 579 599  ``filename``  either a path or None in which case a 580 600 filename in the current directory is chosen automatically 581 601 (default:None) … … 628 648 setPixel = im.setPixel 629 649 colorExact = im.colorExact 630 650 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): 633 653 x = <int>(invblk * self.rows[i].positions[j]) 634 654 y = <int>(invblk * i) 635 655 r,g,b = colorComponents( getPixel((x,y))) … … 642 662 643 663 def density(self): 644 664 """ 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``. 648 667 649 668 EXAMPLES:: 650 669 651 670 sage: A = matrix(QQ,3,3,[0,1,2,3,0,0,6,7,8],sparse=True) 652 671 sage: A.density() 653 672 2/3 654 673 655 674 Notice that the density parameter does not ensure the density 656 675 of a matrix; it is only an upper bound. 657 676 658 677 :: 659 678 660 679 sage: A = random_matrix(GF(127),200,200,density=0.3, sparse=True) 661 680 sage: A.density() 662 681 2073/8000 663 682 """ 664 cdef Py_ssize_t i, n onzero_entries683 cdef Py_ssize_t i, nz = 0 665 684 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 669 687 670 return rings.ZZ(n onzero_entries)/rings.ZZ(self._nrows*self._ncols)688 return rings.ZZ(nz) / rings.ZZ(self._nrows*self._ncols) 671 689 672 690 def transpose(self): 673 691 """ 674 Return the transpose of self.675 692 Return the transpose of ``self``. 693 676 694 EXAMPLE:: 677 695 678 696 sage: A = matrix(GF(127),3,3,[0,1,0,2,0,0,3,0,0],sparse=True) 679 697 sage: A 680 698 [0 1 0] … … 695 713 cdef int i, j 696 714 cdef c_vector_modint row 697 715 cdef Matrix_modn_sparse B 698 716 699 717 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): 701 719 row = self.rows[i] 702 for j from 0 <= j < row.num_nonzero:720 for j in range(row.num_nonzero): 703 721 set_entry(&B.rows[row.positions[j]], i, row.entries[j]) 704 722 if self._subdivisions is not None: 705 723 row_divs, col_divs = self.subdivisions() 706 724 B.subdivide(col_divs, row_divs) 707 725 return B 708 726 709 727 def matrix_from_rows(self, rows): 710 728 """ 711 Return the matrix constructed from selfusing rows with indices in729 Return the matrix constructed from ``self`` using rows with indices in 712 730 the rows list. 713 731 714 732 INPUT: 715 716 733 717 734  ``rows``  list or tuple of row indices 718 719 735 720 736 EXAMPLE:: 721 737 722 738 sage: M = MatrixSpace(GF(127),3,3,sparse=True) 723 739 sage: A = M(range(9)); A 724 740 [0 1 2] … … 728 744 [6 7 8] 729 745 [3 4 5] 730 746 """ 731 cdef int i,k747 cdef Py_ssize_t i, k 732 748 cdef Matrix_modn_sparse A 733 cdef c_vector_modint row734 749 735 750 if not isinstance(rows, (list, tuple)): 736 751 raise TypeError, "rows must be a list of integers" 737 752 738 739 753 A = self.new_matrix(nrows = len(rows)) 740 754 741 k = 0 742 for ii in rows: 743 i = ii 755 for k, i in enumerate(rows): 744 756 if i < 0 or i >= self.nrows(): 745 757 raise IndexError, "row %s out of range"%i 746 758 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 751 761 return A 752 753 762 754 763 def matrix_from_columns(self, cols): 755 764 """ … … 768 777 [5 4] 769 778 [8 7] 770 779 """ 771 cdef int i,j780 cdef Py_ssize_t i, j 772 781 cdef Matrix_modn_sparse A 773 782 cdef c_vector_modint row 774 783 … … 779 788 780 789 cols = dict(zip([int(e) for e in cols],range(len(cols)))) 781 790 782 for i from 0 <= i < self.nrows():791 for i in range(self._nrows): 783 792 row = self.rows[i] 784 for j from 0 <= j < row.num_nonzero:793 for j in range(row.num_nonzero): 785 794 if int(row.positions[j]) in cols: 786 795 set_entry(&A.rows[i], cols[int(row.positions[j])], row.entries[j]) 787 796 return A … … 814 823 815 824 def rank(self, gauss=False): 816 825 """ 817 Compute the rank of self.818 826 Compute the rank of ``self``. 827 819 828 INPUT: 820 821 829 822 830  ``gauss``  if True LinBox' Gaussian elimination is 823 831 used. If False 'Symbolic Reordering' as implemented in LinBox is 824 832 used. If 'native' the native Sage implementation is used. (default: 825 833 False) 826 827 834 828 835 EXAMPLE:: 829 836 830 837 sage: A = random_matrix(GF(127),200,200,density=0.01,sparse=True) 831 838 sage: r1 = A.rank(gauss=False) 832 839 sage: r2 = A.rank(gauss=True) … … 835 842 True 836 843 sage: r1 837 844 155 838 845 839 846 ALGORITHM: Uses LinBox or native implementation. 840 847 841 848 REFERENCES: 842 849 843 850  JeanGuillaume Dumas and Gilles Villars. 'Computing the Rank … … 846 853 on Computer Algebra in Scientific Computing, Big Yalta, 847 854 Crimea, Ukraine, 2227 sept. 2002, SpringerVerlag, 848 855 http://perso.enslyon.fr/gilles.villard/BIBLIOGRAPHIE/POSTSCRIPT/rankjgd.ps 849 856 850 857 .. note:: 851 858 852 859 For very sparse matrices Gaussian elimination is faster … … 873 880 874 881 def _solve_right_nonsingular_square(self, B, algorithm=None, check_rank = True): 875 882 """ 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 881 887 .. note:: 882 888 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 887 892 INPUT: 888 889 893 890 894  ``B``  a matrix or vector 891 895 892 896  ``algorithm``  one of the following: … … 905 909 906 910  ``check_rank``  check rank before attempting to 907 911 solve (default: True) 908 909 912 910 913 OUTPUT: a matrix or vector 911 914 912 915 EXAMPLES:: 913 916 914 917 sage: A = matrix(GF(127), 3, [1,2,3,1,2,5,2,3,1], sparse=True) 915 918 sage: b = vector(GF(127),[1,2,3]) 916 919 sage: x = A \ b; x … … 982 985 """ 983 986 Return lift of this matrix to a sparse matrix over the integers. 984 987 985 EXAMPLES: 988 EXAMPLES:: 989 986 990 sage: a = matrix(GF(7),2,3,[1..6], sparse=True) 987 991 sage: a.lift() 988 992 [1 2 3] … … 1007 1011 1008 1012 cdef mpz_vector* L_row 1009 1013 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])) 1012 1018 A_row = &(self.rows[i]) 1013 1019 sage_free(L_row.entries) 1014 1020 L_row.entries = <mpz_t*> sage_malloc(sizeof(mpz_t)*A_row.num_nonzero) … … 1021 1027 sage_free(L_row.entries) 1022 1028 L_row.entries = NULL 1023 1029 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): 1025 1031 L_row.positions[j] = A_row.positions[j] 1026 1032 mpz_init_set_si(L_row.entries[j], A_row.entries[j]) 1027 1033 L.subdivide(self.subdivisions()) 1028 1034 return L 1029 1035 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  Partiallyechelonized 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[r1].num_nonzero == 0: 1105 raise ValueError, "Row r1 has no nonzero entries" 1106 c = self.rows[r1].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(ir, 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 nonpivots 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 nonpivot columns: %s"%len(non_pivots) 1153 assert set(pivots + non_pivots) == set(range(n)) 1154 1155 if verb: 1156 print "Copying nonpivot 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(ir, 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(ir, j): 1236 Z.rows[i].positions[nz] = j+c 1237 Z.rows[i].entries[nz] = <long>E.getEntry(ir, 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 offdiagonal 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`` inplace 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 13 13 """ 14 14 v.entries = <int*>sage_malloc(num_nonzero*sizeof(int)) 15 15 if v.entries == NULL: 16 raise MemoryError, "Error allocating memory "16 raise MemoryError, "Error allocating memory for sparse vector." 17 17 v.positions = <Py_ssize_t*>sage_malloc(num_nonzero*sizeof(Py_ssize_t)) 18 18 if v.positions == NULL: 19 19 sage_free(v.entries) 20 raise MemoryError, "Error allocating memory "20 raise MemoryError, "Error allocating memory for sparse vector." 21 21 return 0 22 22 23 23 cdef int init_c_vector_modint(c_vector_modint* v, int p, Py_ssize_t degree, … … 35 35 v.p = p 36 36 return 0 37 37 38 cdef 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 38 64 cdef void clear_c_vector_modint(c_vector_modint* v): 39 65 sage_free(v.entries) 40 66 sage_free(v.positions) … … 51 77 i = 0 52 78 j = n1 53 79 while i<=j: 54 if i == j:55 if v[i] == x:56 return i57 return 158 80 k = (i+j)/2 59 81 if v[k] > x: 60 82 j = k1 61 83 elif v[k] < x: 62 84 i = k+1 63 else: # only possibility is thatv[k] == x85 else: # v[k] == x 64 86 return k 65 87 return 1 66 88 … … 93 115 j = k1 94 116 elif v[k] < x: 95 117 i = k+1 96 else: # only possibility is thatv[k] == x118 else: # v[k] == x 97 119 ins[0] = k 98 120 return k 99 # end while100 121 ins[0] = j+1 101 122 return 1 102 123 … … 107 128 """ 108 129 if n >= v.degree or n < 0: 109 130 raise IndexError, "Index must be between 0 and the degree minus 1." 110 return 1111 131 cdef Py_ssize_t m 112 132 m = binary_search0_modn(v.positions, v.num_nonzero, n) 113 133 if m == 1: … … 133 153 """ 134 154 if n < 0 or n >= v.degree: 135 155 raise IndexError, "Index (=%s) must be between 0 and %s."%(n, v.degree1) 136 return 1137 156 cdef Py_ssize_t i, m, ins 138 cdef Py_ssize_t m2, ins2139 cdef Py_ssize_t *pos140 cdef int *e141 157 142 158 x = x % v.p 143 159 if x<0: x = x + v.p 144 160 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 mth 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 nonzero x 154 164 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[i1] = e[i] 164 v.positions[i1] = pos[i] 165 sage_free(e) 166 sage_free(pos) 165 else: # Remove entry as x == 0 167 166 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) 168 173 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[i1] 190 v.positions[i] = pos[i1] 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[i1] 182 v.positions[i] = v.positions[i1] 183 184 v.entries[ins] = x 185 v.positions[ins] = n 186 187 cdef 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 193 213 194 214 cdef int add_c_vector_modint_init(c_vector_modint* sum, c_vector_modint* v, 195 215 c_vector_modint* w, int multiple) except 1: … … 198 218 """ 199 219 if v.p != w.p: 200 220 raise ArithmeticError, "The vectors must be modulo the same prime." 201 return 1202 221 if v.degree != w.degree: 203 222 raise ArithmeticError, "The vectors must have the same degree." 204 return 1205 223 206 224 cdef int s 207 225 cdef Py_ssize_t nz, i, j, k … … 276 294 if scalar < 0: 277 295 scalar = scalar + v.p 278 296 cdef Py_ssize_t i 279 for i from 0 <= i < v.num_nonzero:297 for i in range(v.num_nonzero): 280 298 v.entries[i] = (v.entries[i] * scalar) % v.p 281 299 return 0 282 300 283