# Changeset 3111:b790f0e707f0

Ignore:
Timestamp:
02/21/07 02:04:30 (6 years ago)
Branch:
default
Message:

snapshot -- first version of sparse rational matrices.

Files:
6 edited
2 moved

Unmodified
Removed
• ## sage/matrix/matrix_integer_dense.pyx

 r3110 cdef Linbox_integer_dense linbox linbox = Linbox_integer_dense() #import sage.misc.misc #USE_LINBOX_POLY = not sage.misc.misc.is_64bit() USE_LINBOX_POLY = True cdef sage.structure.element.Matrix _matrix_times_matrix_c_impl(self, sage.structure.element.Matrix right): return self._multiply_classical(right) # NOTE -- the multimodular matrix multiply implementation # breaks on 64-bit machines; e..g, the following doctests # *all* fail if multimodular matrix multiply is enabled # on sage.math.washington.edu: #sage -t  devel/sage-main/sage/modular/modsym/modsym.py #sage -t  devel/sage-main/sage/modular/modsym/space.py #sage -t  devel/sage-main/sage/modular/modsym/subspace.py #sage -t  devel/sage-main/sage/modular/hecke/hecke_operator.py #sage -t  devel/sage-main/sage/modular/hecke/module.py ############# # see the tune_multiplication function below. if n <= 20: return self._multiply_classical(right) return self._multiply_multi_modular(right) ##         a = self.height(); b = right.height() ##         # waiting for multiply_multi_modular to get fixed, and not assume all matrix entries ##         # are between 0 and prod - 1. ##         if float(max(a,b)) / float(n) >= 0.70: ##             return self._multiply_classical(right) ##         else: ##             return self._multiply_multi_modular(right) a = self.height(); b = right.height() if float(max(a,b)) / float(n) >= 0.70: return self._multiply_classical(right) else: return self._multiply_multi_modular(right) cdef ModuleElement _lmul_c_impl(self, RingElement right):
• ## sage/matrix/matrix_modn_sparse.pyx

 r3109 coerce -- ignored """ cdef int s, y, z, p cdef int s, z, p cdef Py_ssize_t i, j, k cdef object seq cdef void** X R = self._base_ring for ij, x in entries.iteritems(): y = x z = R(y) z = R(x) if z != 0: i, j = ij  # nothing better to do since this is user input, which may be bogus. for i from 0 <= i < self._nrows: for  j from 0 <= j < self._ncols: set_entry(&self.rows[i], j, R(
• ## sage/matrix/matrix_rational_sparse.pyx.orig

 r3092 START_PRIME = 20011  # used for multi-modular algorithms cdef class Vector_mpq cdef void Vector_mpq_rescale(Vector_mpq w, mpq_t x): scale_mpq_vector(&w.v, x) cdef class Vector_mpq: """ Vector_mpq -- a sparse vector of GMP rationals.  This is a Python extension type that wraps the C implementation of sparse vectors modulo a small prime. """ cdef mpq_vector v def __init__(self, int degree, int num_nonzero=0, entries=[], sort=True): cdef int i init_mpq_vector(&self.v, degree, num_nonzero) if entries != []: if len(entries) != num_nonzero: raise ValueError, "length of entries (=%s) must equal num_nonzero (=%s)"%(len(entries), num_nonzero) if sort: entries = list(entries) # copy so as not to modify passed in list entries.sort() for i from 0 <= i < num_nonzero: s = str(entries[i][1]) mpq_set_str(self.v.entries[i], s, 0) self.v.positions[i] = entries[i][0] def __dealloc__(self): clear_mpq_vector(&self.v) def __getitem__(self, int n): cdef mpq_t x cdef sage.rings.rational.Rational a mpq_init(x) mpq_vector_get_entry(&x, &self.v, n) a = sage.rings.rational.Rational() a.set_from_mpq(x) mpq_clear(x) return a def cmp(self, Vector_mpq other): return mpq_vector_cmp(&self.v, &other.v) def __richcmp__(Vector_mpq self, x, int op): if not isinstance(x, Vector_mpq): return -1 cdef int n n = self.cmp(x) if op == 0: return bool(n < 0) elif op == 1: return bool(n <= 0) elif op == 2: return bool(n == 0) elif op == 3: return bool(n != 0) elif op == 4: return bool(n > 0) elif op == 5: return bool(n >= 0) def __setitem__(self, int n, x): cdef object s s = str(x) mpq_vector_set_entry_str(&self.v, n, s) def __repr__(self): return str(list(self)) def degree(self): return self.v.degree def num_nonzero(self): return self.v.num_nonzero def list(self): return mpq_vector_to_list(&self.v) cdef void rescale(self, mpq_t x): scale_mpq_vector(&self.v, x) def __add__(Vector_mpq self, Vector_mpq other): cdef mpq_vector z1, *z2 cdef Vector_mpq w cdef mpq_t ONE mpq_init(ONE) mpq_set_si(ONE,1,1) add_mpq_vector_init(&z1, &self.v, &other.v, ONE) mpq_clear(ONE) w = Vector_mpq(self.v.degree) z2 = &(w.v) clear_mpq_vector(z2)   # free memory wasted on allocated w z2.entries = z1.entries z2.positions = z1.positions z2.num_nonzero = z1.num_nonzero # at this point we do *not* free z1, since it is referenced by w. return w def __sub__(Vector_mpq self, Vector_mpq other): return self + other*(-1) def copy(self): cdef int i cdef Vector_mpq w w = Vector_mpq(self.v.degree, self.v.num_nonzero) for i from 0 <= i < self.v.num_nonzero: mpq_set(w.v.entries[i], self.v.entries[i]) w.v.positions[i] = self.v.positions[i] return w def __mul__(x, y): if isinstance(x, Vector_mpq): self = x other = y elif isinstance(y, Vector_mpq): self = y other = x else: raise TypeError, "Invalid types." cdef object s, z cdef mpq_t t z = self.copy() mpq_init(t) s = str(other) mpq_set_str(t, s, 0) Vector_mpq_rescale(z, t) mpq_clear(t) return z def randomize(self, int sparcity, bound=3): """ randomize(self, int sparcity, exact=False): The sparcity is a bound on the number of nonzeros per row. """ cdef int i for i from 0 <= i < sparcity: self[random.randrange(self.v.degree)] = random.randrange(1,bound) ############################################################# cdef int i self.rows = sage_malloc(nrows*sizeof(mpq_vector)) # TODO: check worked  -- bug self.is_init = init if self.is_init:
• ## sage/matrix/matrix_space.py

 r3109 import matrix_rational_dense ##import matrix_rational_sparse import matrix_rational_sparse ## import matrix_cyclo_dense return matrix_modn_sparse.Matrix_modn_sparse # the default elif sage.rings.rational_field.is_RationalField(R): return matrix_rational_sparse.Matrix_rational_sparse return matrix_generic_sparse.Matrix_generic_sparse
• ## sage/modules/vector_rational_sparse.pyx

 r3092 include 'vector_rational_sparse_c.pxi' cdef class Vector_mpq cdef void Vector_mpq_rescale(Vector_mpq w, mpq_t x): scale_mpq_vector(&w.v, x) cdef class Vector_mpq: """ Vector_mpq -- a sparse vector of GMP rationals.  This is a Python extension type that wraps the C implementation of sparse vectors modulo a small prime. """ cdef mpq_vector v def __init__(self, int degree, int num_nonzero=0, entries=[], sort=True): cdef int i init_mpq_vector(&self.v, degree, num_nonzero) if entries != []: if len(entries) != num_nonzero: raise ValueError, "length of entries (=%s) must equal num_nonzero (=%s)"%(len(entries), num_nonzero) if sort: entries = list(entries) # copy so as not to modify passed in list entries.sort() for i from 0 <= i < num_nonzero: s = str(entries[i][1]) mpq_set_str(self.v.entries[i], s, 0) self.v.positions[i] = entries[i][0] def __dealloc__(self): clear_mpq_vector(&self.v) def __getitem__(self, int n): cdef mpq_t x cdef sage.rings.rational.Rational a mpq_init(x) mpq_vector_get_entry(&x, &self.v, n) a = sage.rings.rational.Rational() a.set_from_mpq(x) mpq_clear(x) return a def cmp(self, Vector_mpq other): return mpq_vector_cmp(&self.v, &other.v) def __richcmp__(Vector_mpq self, x, int op): if not isinstance(x, Vector_mpq): return -1 cdef int n n = self.cmp(x) if op == 0: return bool(n < 0) elif op == 1: return bool(n <= 0) elif op == 2: return bool(n == 0) elif op == 3: return bool(n != 0) elif op == 4: return bool(n > 0) elif op == 5: return bool(n >= 0) def __setitem__(self, int n, x): cdef object s s = str(x) mpq_vector_set_entry_str(&self.v, n, s) def __repr__(self): return str(list(self)) def degree(self): return self.v.degree def num_nonzero(self): return self.v.num_nonzero def list(self): return mpq_vector_to_list(&self.v) cdef void rescale(self, mpq_t x): scale_mpq_vector(&self.v, x) def __add__(Vector_mpq self, Vector_mpq other): cdef mpq_vector z1, *z2 cdef Vector_mpq w cdef mpq_t ONE mpq_init(ONE) mpq_set_si(ONE,1,1) add_mpq_vector_init(&z1, &self.v, &other.v, ONE) mpq_clear(ONE) w = Vector_mpq(self.v.degree) z2 = &(w.v) clear_mpq_vector(z2)   # free memory wasted on allocated w z2.entries = z1.entries z2.positions = z1.positions z2.num_nonzero = z1.num_nonzero # at this point we do *not* free z1, since it is referenced by w. return w def __sub__(Vector_mpq self, Vector_mpq other): return self + other*(-1) def copy(self): cdef int i cdef Vector_mpq w w = Vector_mpq(self.v.degree, self.v.num_nonzero) for i from 0 <= i < self.v.num_nonzero: mpq_set(w.v.entries[i], self.v.entries[i]) w.v.positions[i] = self.v.positions[i] return w def __mul__(x, y): if isinstance(x, Vector_mpq): self = x other = y elif isinstance(y, Vector_mpq): self = y other = x else: raise TypeError, "Invalid types." cdef object s, z cdef mpq_t t z = self.copy() mpq_init(t) s = str(other) mpq_set_str(t, s, 0) Vector_mpq_rescale(z, t) mpq_clear(t) return z def randomize(self, int sparcity, bound=3): """ randomize(self, int sparcity, exact=False): The sparcity is a bound on the number of nonzeros per row. """ cdef int i for i from 0 <= i < sparcity: self[random.randrange(self.v.degree)] = random.randrange(1,bound)
• ## sage/modules/vector_rational_sparse_c.pxi

 r3092 ############################################################# include "../ext/cdefs.pxi" cdef struct mpq_vector: mpq_t *entries      # array of nonzero entries int   *positions    # positions of those nonzero entries, starting at 0 int    degree       # the degree of this sparse vector int    num_nonzero  # the number of nonzero entries of this vector. cdef int allocate_mpq_vector(mpq_vector* v, int num_nonzero) except -1: include 'vector_rational_sparse_h.pxi' from sage.rings.rational cimport Rational cdef int allocate_mpq_vector(mpq_vector* v, Py_ssize_t num_nonzero) except -1: """ Allocate memory for a mpq_vector -- most user code won't call this. It does *not* clear the entries of v, if there are any. """ cdef int i cdef Py_ssize_t i v.entries = sage_malloc(num_nonzero*sizeof(mpq_t)) if v.entries == 0: if v.entries == NULL: raise MemoryError, "Error allocating memory" for i from 0 <= i < num_nonzero: mpq_init(v.entries[i]) v.positions = sage_malloc(num_nonzero*sizeof(int)) if v.positions == 0: v.positions = sage_malloc(num_nonzero*sizeof(Py_ssize_t)) if v.positions == NULL: for i from 0 <= i < num_nonzero: mpq_clear(v.entries[i]) sage_free(v.entries) v.entries = NULL raise MemoryError, "Error allocating memory" return 0 cdef int init_mpq_vector(mpq_vector* v, int degree, int num_nonzero) except -1: cdef int init_mpq_vector(mpq_vector* v, Py_ssize_t degree, Py_ssize_t num_nonzero) except -1: """ Initialize a mpq_vector -- most user code *will* call this. cdef void clear_mpq_vector(mpq_vector* v): cdef int i cdef Py_ssize_t i # Free all mpq objects allocated in creating v for i from 0 <= i < v.num_nonzero: sage_free(v.positions) cdef int mpq_binary_search0(mpq_t* v, int n, mpq_t x): cdef Py_ssize_t binary_search(Py_ssize_t* v, Py_ssize_t n, Py_ssize_t x, Py_ssize_t* ins): """ Find the position of the integer x in the array v, which has length n. Returns -1 if x is not in the array v, and in this case ins is set equal to the position where x should be inserted in order to obtain an ordered array. """ if n == 0: ins[0] = 0 return -1 cdef Py_ssize_t i, j, k i = 0 j = n-1 while i<=j: if i == j: if v[i] == x: ins[0] = i return i if v[i] < x: ins[0] = i + 1 else: ins[0] = i return -1 k = (i+j)/2 if v[k] > x: j = k-1 elif v[k] < x: i = k+1 else:   # only possibility is that v[k] == x ins[0] = k return k # end while ins[0] = j+1 return -1 cdef Py_ssize_t binary_search0(Py_ssize_t* v, Py_ssize_t n, Py_ssize_t x): """ Find the position of the int x in the array v, which has length n. Returns -1 if x is not in the array v. """ if n == 0: return -1 cdef Py_ssize_t i, j, k i = 0 j = n-1 while i<=j: if i == j: if v[i] == x: return i return -1 k = (i+j)/2 if v[k] > x: j = k-1 elif v[k] < x: i = k+1 else:   # only possibility is that v[k] == x return k return -1 cdef Py_ssize_t mpq_binary_search0(mpq_t* v, Py_ssize_t n, mpq_t x): """ Find the position of the rational x in the array v, which has length n. Returns -1 if x is not in the array v. """ cdef int i, j, k, c cdef Py_ssize_t i, j, k, c if n == 0: return -1 return -1 cdef int mpq_binary_search(mpq_t* v, int n, mpq_t x, int* ins): cdef Py_ssize_t mpq_binary_search(mpq_t* v, Py_ssize_t n, mpq_t x, Py_ssize_t* ins): """ Find the position of the integer x in the array v, which has length n. x -- mpq_t  (rational) OUTPUT: position of x (as an int) position of x (as an Py_ssize_t) ins -- (call be pointer), the insertion point if x is not found. """ cdef int i, j, k, c cdef Py_ssize_t i, j, k, c if n == 0: ins[0] = 0 return -1 cdef int mpq_vector_get_entry(mpq_t* ans, mpq_vector* v, int n) except -1: cdef int mpq_vector_get_entry(mpq_t* ans, mpq_vector* v, Py_ssize_t n) except -1: """ Returns the n-th entry of the sparse vector v.  This if n >= v.degree: raise IndexError, "Index must be between 0 and %s."%(v.degree - 1) cdef int m cdef Py_ssize_t m m = binary_search0(v.positions, v.num_nonzero, n) if m == -1: through the nonzero elements of x, in order. """ cdef X cdef sage.rings.rational.Rational a cdef int i cdef object X cdef Rational a cdef Py_ssize_t i X = [] for i from 0 <= i < v.num_nonzero: a = sage.rings.rational.Rational() a = Rational() a.set_from_mpq(v.entries[i]) X.append( (v.positions[i], a) ) cdef int mpq_vector_set_entry(mpq_vector* v, int n, mpq_t x) except -1: cdef int mpq_vector_set_entry(mpq_vector* v, Py_ssize_t n, mpq_t x) except -1: """ Set the n-th component of the sparse vector v equal to x. if n >= v.degree or n < 0: raise IndexError, "Index must be between 0 and the degree minus 1." return -1 cdef int i, m, ins cdef int m2, ins2 cdef int *pos cdef Py_ssize_t i, m, ins cdef Py_ssize_t m2, ins2 cdef Py_ssize_t *pos cdef mpq_t *e cdef mpq_t mpq_set_tmp mpq_init(mpq_set_tmp) cdef int mpq_vector_set_entry_str(mpq_vector* v, int n, char *x_str) except -1: cdef int mpq_vector_set_entry_str(mpq_vector* v, Py_ssize_t n, char *x_str) except -1: """ Set the n-th component of the sparse vector v equal to x. raise ArithmeticError, "The vectors must have the same degree." cdef int nz, i, j, k, do_multiply cdef Py_ssize_t nz, i, j, k, do_multiply cdef mpq_vector* z cdef mpq_t tmp init_mpq_vector(v, v.degree, 0) return 0 cdef int i cdef Py_ssize_t i for i from 0 <= i < v.num_nonzero: # v.entries[i] = scalar * v.entries[i] elif v.degree > w.degree: return 1 cdef int i, c cdef Py_ssize_t i cdef int c for i from 0 <= i < v.num_nonzero: c = mpq_cmp(v.entries[i], w.entries[i])
• ## setup.py

 r3104 libraries = ['gmp']) matrix_rational_sparse = Extension('sage.matrix.matrix_rational_sparse', ['sage/matrix/matrix_rational_sparse.pyx'], libraries = ['gmp']) matrix_integer_dense = Extension('sage.matrix.matrix_integer_dense', ['sage/matrix/matrix_integer_dense.pyx'], matrix_integer_dense, matrix_rational_dense, matrix_rational_sparse, matrix_integer_2x2, ##     matrix_integer_sparse,
Note: See TracChangeset for help on using the changeset viewer.