Ticket #7723: trac_7723-double_sparse.patch

File trac_7723-double_sparse.patch, 41.9 KB (added by dagss, 11 years ago)

  • doc/en/reference/matrices.rst 

    # HG changeset patch
# User Dag Sverre Seljebotn <dagss@student.matnat.uio.no>
# Date 1261052817 -3600
# Node ID 25f413c99ea3becff7701c12da9feda75323f831
# Parent  977f52f47b09cdf5cfbca81ff74d5fd9f5c257fe
[mq]: double_sparse

diff -r 977f52f47b09 -r 25f413c99ea3 doc/en/reference/matrices.rst
    a b  
    8080   sage/matrix/matrix_real_double_dense
    8181
    8282   sage/matrix/matrix_complex_double_dense
     83
     84   sage/matrix/matrix_double_sparse

  • module_list.py 

    diff -r 977f52f47b09 -r 25f413c99ea3 module_list.py
    a b  
    687687              libraries=[BLAS, BLAS2],
    688688              include_dirs = numpy_include_dirs),
    689689
     690    Extension('sage.matrix.matrix_double_sparse',
     691              sources = ['sage/matrix/matrix_double_sparse.pyx'],
     692              include_dirs = numpy_include_dirs),
     693
    690694    Extension('sage.matrix.matrix_generic_dense',
    691695              sources = ['sage/matrix/matrix_generic_dense.pyx']),
    692696   

  • sage/matrix/docs.py 

    diff -r 977f52f47b09 -r 25f413c99ea3 sage/matrix/docs.py
    a b  
    179179       * cdef _list -- list of underlying elements (need not be a copy)
    180180       * cdef _dict -- sparse dictionary of underlying elements
    181181       * cdef _add_ -- add two matrices with identical parents
    182        * _matrix_times_matrix_c_impl -- multiply two matrices with compatible dimensions and
    183                                         identical base rings (both sparse or both dense)
     182       * _matrix_times_matrix_ -- multiply two matrices with compatible dimensions and
     183                                  identical base rings (both sparse or both dense)
     184       * _generic_matrix_times_matrix_ -- multiply two matrices with compatible dimensions but
     185                                          different base rings (for optimized code paths
     186                                          that works between different base rings or
     187                                          sparse-dense-multiplication)
    184188       * cdef _cmp_c_impl -- compare two matrices with identical parents
    185        * cdef _lmul_c_impl -- multiply this matrix on the right by a scalar, i.e., self * scalar
    186        * cdef _rmul_c_impl -- multiply this matrix on the left by a scalar, i.e., scalar * self
     189       * cdef _lmul_ -- multiply this matrix on the right by a scalar, i.e., self * scalar
     190       * cdef _rmul_ -- multiply this matrix on the left by a scalar, i.e., scalar * self
    187191       * __copy__
    188192       * __neg__
    189193   

  • sage/matrix/matrix_double_dense.pxd 

    diff -r 977f52f47b09 -r 25f413c99ea3 sage/matrix/matrix_double_dense.pxd
    a b  
    1313    cdef Matrix_double_dense _new(self, int nrows=*, int ncols=*)
    1414    cdef _c_compute_LU(self)
    1515    cdef cnumpy.ndarray _matrix_numpy
     16
     17cdef Matrix_double_dense create_uninitialized_Matrix_double_dense(parent, cls=*)

  • sage/matrix/matrix_double_dense.pyx 

    diff -r 977f52f47b09 -r 25f413c99ea3 sage/matrix/matrix_double_dense.pyx
    a b  
    4040import sage.rings.complex_double
    4141from matrix cimport Matrix
    4242from sage.structure.element cimport ModuleElement,Vector
     43from sage.structure.parent cimport Parent
    4344from constructor import matrix
    4445from sage.modules.free_module_element import vector
    4546cimport sage.structure.element
    4647from matrix_space import MatrixSpace
     48from matrix_double_sparse cimport Matrix_double_sparse
    4749
    4850cimport numpy as cnumpy
    4951
    … …  
    5557
    5658
    5759from matrix_complex_double_dense import Matrix_complex_double_dense
     60from matrix_real_double_dense import Matrix_real_double_dense
     61
     62cdef Matrix_double_dense create_uninitialized_Matrix_double_dense(parent, cls=None):
     63    """
     64    Creates an uninitialized Matrix_double_dense with the given parent.
     65    The _matrix_numpy field of the result must be filled in by the
     66    caller.
     67    """
     68    if cls is None:
     69        if parent._base == sage.rings.real_double.RDF:
     70            cls = Matrix_real_double_dense
     71        elif parent._base == sage.rings.complex_double.CDF:
     72            cls = Matrix_complex_double_dense
     73        else:
     74            raise ValueError("Invalid cls")
     75    return cls.__new__(cls,parent,None,None,None)
    5876
    5977cdef class Matrix_double_dense(matrix_dense.Matrix_dense):
    6078    """
    … …  
    291309        if ncols == -1:
    292310            ncols = self._ncols
    293311        parent = self.matrix_space(nrows, ncols)
    294         m = self.__class__.__new__(self.__class__,parent,None,None,None)
     312        m = create_uninitialized_Matrix_double_dense(parent, self.__class__)
    295313        return m
    296314
    297315
    … …  
    374392    # def _pickle(self):                        #unsure how to implement
    375393    # def _unpickle(self, data, int version):   # use version >= 0 #unsure how to implement
    376394    ######################################################################
    377     cdef sage.structure.element.Matrix _matrix_times_matrix_(self, sage.structure.element.Matrix right):
     395    cdef sage.structure.element.Matrix _generic_matrix_times_matrix_(self,
     396            sage.structure.element.Matrix right, Parent codomain):
    378397        """
    379398        Multiply self*right as matrices.
    380399
    381         EXAMPLES
     400        Seperate code path taken for dense time sparse.
     401
     402        TESTS::
    382403            sage: A = matrix(RDF,3,range(1,10))
    383404            sage: B = matrix(RDF,3,range(1,13))
    384405            sage: A*B
    385406            [ 38.0  44.0  50.0  56.0]
    386407            [ 83.0  98.0 113.0 128.0]
    387408            [128.0 152.0 176.0 200.0]
    388         """
    389         if self._ncols!=right._nrows:
    390             raise IndexError, "Number of columns of self must equal number of rows of right"
    391        
     409
     410            sage: C=A*B.sparse_matrix(); C
     411            [ 38.0  44.0  50.0  56.0]
     412            [ 83.0  98.0 113.0 128.0]
     413            [128.0 152.0 176.0 200.0]
     414            sage: parent(C)
     415            Full MatrixSpace of 3 by 4 dense matrices over Real Double Field
     416
     417            sage: C = A.change_ring(CDF)*B.sparse_matrix(); C
     418            [ 38.0  44.0  50.0  56.0]
     419            [ 83.0  98.0 113.0 128.0]
     420            [128.0 152.0 176.0 200.0]
     421            sage: parent(C)
     422            Full MatrixSpace of 3 by 4 dense matrices over Complex Double Field
     423
     424            sage: C = A.change_ring(CDF)*B; C
     425            [ 38.0  44.0  50.0  56.0]
     426            [ 83.0  98.0 113.0 128.0]
     427            [128.0 152.0 176.0 200.0]
     428            sage: parent(C)
     429            Full MatrixSpace of 3 by 4 dense matrices over Complex Double Field
     430
     431            sage: A*B.change_ring(ZZ)
     432            [ 38.0  44.0  50.0  56.0]
     433            [ 83.0  98.0 113.0 128.0]
     434            [128.0 152.0 176.0 200.0]
     435
     436            sage: A.change_ring(ZZ)*B
     437            [ 38.0  44.0  50.0  56.0]
     438            [ 83.0  98.0 113.0 128.0]
     439            [128.0 152.0 176.0 200.0]
     440       """
    392441        if self._nrows == 0 or self._ncols == 0 or right._nrows == 0 or right._ncols == 0:
    393             return self.matrix_space(self._nrows, right._ncols).zero_matrix()
     442            return codomain.zero_matrix()
    394443
    395         cdef Matrix_double_dense M,_right,_left
    396         M = self._new(self._nrows, right._ncols)
    397         _right = right
    398         _left = self
    399         global numpy
     444        cdef Matrix_double_dense left, M
     445
     446        # There's no real advantage to have a mixed complex/real code path here,
     447        # change both operands to same base ring
     448        left = self
     449        base = codomain._base
     450        if left._parent._base is not base:
     451            left = left.change_ring(base)
     452        if right._parent._base is not base:
     453            right = right.change_ring(base)
     454
    400455        if numpy is None:
    401456            import numpy
    402 
    403         M._matrix_numpy = numpy.dot(_left._matrix_numpy, _right._matrix_numpy)
     457        M = create_uninitialized_Matrix_double_dense(codomain, left.__class__)
     458        if right.is_dense_c():
     459            M._matrix_numpy = numpy.dot(left._matrix_numpy,
     460                                        (<Matrix_double_dense?>right)._matrix_numpy)
     461        else:
     462            # Do a transposed multiplication -- scipy.sparse only defines right-mul.
     463            # Avoid using Sage matrices to avoid some data copying -- transposing
     464            # a scipy.sparse CSC matrix results in a CSR matrix of the same data,
     465            # and then specific CSR multiplication code is used.
     466            right_T = (<Matrix_double_sparse?>right)._entries_csc.transpose(copy=False)
     467            left_T = left._matrix_numpy.T
     468           
     469            M._matrix_numpy = numpy.ascontiguousarray((right_T * left_T).T)
    404470        return M
    405471
    406472    # cdef int _cmp_c_impl(self, Matrix right) except -2:

  • new file sage/matrix/matrix_double_sparse.pxd 

    diff -r 977f52f47b09 -r 25f413c99ea3 sage/matrix/matrix_double_sparse.pxd
    - +  
     1cimport matrix_sparse
     2cimport numpy as cnumpy
     3
     4cdef class Matrix_double_sparse(matrix_sparse.Matrix_sparse):
     5    cdef object _entries_csc
     6    cdef object _numpy_dtype
     7    cdef object _python_dtype
     8    cdef object _sage_dtype
     9    cdef object _sage_vector_dtype
     10    cdef cnumpy.ndarray _L_M, _U_M, _P_M
     11    cdef int _LU_valid
     12    cdef Matrix_double_sparse _new(self, int nrows=*, int ncols=*)
     13
     14cdef Matrix_double_sparse create_uninitialized_Matrix_double_sparse(parent)

  • new file sage/matrix/matrix_double_sparse.pyx 

    diff -r 977f52f47b09 -r 25f413c99ea3 sage/matrix/matrix_double_sparse.pyx
    - +  
     1"""
     2Sparse matrices over Real Double and Complex Double fields
     3
     4EXAMPLES::
     5
     6    sage: matrix(RDF, {(0, 0): 1, (1, 1): 2})
     7    [1.0 0.0]
     8    [0.0 2.0]
     9    sage: matrix(RDF, 3, 3, 1, sparse=True)
     10    [1.0 0.0 0.0]
     11    [0.0 1.0 0.0]
     12    [0.0 0.0 1.0]
     13
     14::
     15
     16    sage: b = MatrixSpace(RDF,2,3,sparse=True).basis()
     17    sage: b[0]
     18    [1.0 0.0 0.0]
     19    [0.0 0.0 0.0]
     20
     21
     22We deal with the case of zero rows or zero columns::
     23
     24    sage: m = MatrixSpace(RDF,0,3,sparse=True)
     25    sage: m.zero_matrix()
     26    []
     27
     28AUTHORS:
     29    - Dag Sverre Seljebotn
     30"""
     31
     32##############################################################################
     33#    Copyright (C) 2009 Dag Sverre Seljebotn <dagss@student.matnat.uio.no>
     34#  Distributed under the terms of the GNU General Public License (GPL)
     35#  The full text of the GPL is available at:
     36#                  http://www.gnu.org/licenses/
     37##############################################################################
     38
     39# TODO list:
     40# - Efficient cmp?
     41# - Efficient get/set_unsafe when not changing sparsity
     42# - Cholesky
     43# - Matrix inversion
     44
     45import math
     46
     47import sage.rings.real_double
     48from sage.rings.real_double import RDF
     49import sage.rings.complex_double
     50from sage.rings.complex_double import CDF
     51from matrix cimport Matrix
     52from sage.structure.element cimport ModuleElement, Vector, RingElement, Element
     53from constructor import matrix
     54from sage.modules.free_module_element import vector
     55cimport sage.structure.element
     56from matrix_space import MatrixSpace
     57cimport matrix_sparse
     58import warnings
     59from sage.structure.parent cimport Parent
     60from matrix_double_dense cimport Matrix_double_dense, create_uninitialized_Matrix_double_dense
     61from matrix_generic_sparse cimport _convert_sparse_entries_to_dict
     62from sage.modules.vector_double_dense cimport Vector_double_dense
     63
     64numpy = None
     65scipy_sparse = None
     66
     67cdef bint VERBOSE = False
     68
     69times = 0
     70
     71cdef class _verbose:
     72    """
     73    The _verbose context manager allows unit testing of different
     74    code paths being taken. Used by in-module doctests *only*.
     75    """
     76    def __enter__(self):
     77        global VERBOSE
     78        VERBOSE = True
     79        return None
     80    def __exit__(self, a, b, c):
     81        global VERBOSE
     82        VERBOSE = False
     83        return False # do not suppress exception
     84
     85cdef Matrix_double_sparse create_uninitialized_Matrix_double_sparse(parent):
     86    """
     87    Creates an uninitialized Matrix_double_dense with the given parent.
     88    The _entries_csc field of the result must be filled in by the
     89    caller.
     90    """
     91    return Matrix_double_sparse.__new__(Matrix_double_sparse, parent,
     92                                        None, None, None)
     93
     94cdef class Matrix_double_sparse(matrix_sparse.Matrix_sparse):
     95    """
     96    Sparse, compressed matrix over double floating point fields (RDF
     97    or CDF).
     98
     99    The compressed format supports fast matrix operations and
     100    algorithms. However, changing the sparsity structure (assigning to
     101    an element which was previously zero, or setting a non-zero
     102    element to zero) is slow. Changes which do not affect sparsity are
     103    faster but still slower than with a hash-based or dense matrix.
     104    It is recommended to construct the matrix with all non-zero
     105    positions filled in and keep the sparsity pattern. Remember to
     106    call :meth:`set_immutable` after construction if possible.
     107   
     108    Operations that are not directly supported will still work, but
     109    may use algorithms which are slow, use a lot of memory and are
     110    numerically unstable. Currently supported efficient operations:
     111
     112     - Addition, subtraction, negation, transpose, pickling
     113
     114     - Multiplications (scalar, vector, dense and sparse
     115       matrices). Multiplying with a dense matrix or vector does not
     116       require more memory than the size of the (dense) result.
     117
     118    The storage requirement is the same as the in-memory requirement,
     119    which is about ``4*nrows + 4*nnz + element_size*nnz`` bytes where
     120    ``nnz`` is the number of non-zero elements, and ``element_size``
     121    is 8 for RDF and 16 for CDF. In addition comes any cached
     122    information (which can be cleared with a call to
     123    meth:`_clear_cache`).
     124
     125    The intention is to eventually support more efficient sparse
     126    matrix algorithms, such as LU factorization and Cholesky
     127    decomposition.  The exact format is currently CSC, this might
     128    change in the future (and may dynamically change in respose to
     129    requested algorithms), but it will remain a compressed format with
     130    inefficient single-element access.
     131
     132    EXAMPLES::
     133   
     134        sage: m = matrix(RDF, [[1,2],[3,4]], sparse=True)
     135        sage: m**2
     136        [ 7.0 10.0]
     137        [15.0 22.0]
     138        sage: n= m^(-1); n
     139        [-2.0  1.0]
     140        [ 1.5 -0.5]   
     141
     142        sage: m = matrix(CDF, [[1,2*I],[3+I,4]], sparse=True)
     143        sage: m**2
     144        [-1.0 + 6.0*I       10.0*I]
     145        [15.0 + 5.0*I 14.0 + 6.0*I]
     146        sage: n= m^(-1); n
     147        [  0.333333333333 + 0.333333333333*I   0.166666666667 - 0.166666666667*I]
     148        [ -0.166666666667 - 0.333333333333*I 0.0833333333333 + 0.0833333333333*I]
     149    """
     150   
     151    def __cinit__(self, Parent parent, entries, copy, coerce):
     152        """
     153        Set up a new matrix
     154
     155        TESTS:
     156       
     157            sage: matrix(RDF, 2**31, 2, sparse=True)
     158            Traceback (most recent call last):
     159                ...
     160            NotImplementedError: Number of rows and columns must currently be smaller than 2^31, even on 64-bit systems. This only concerns RDF/CDF sparse matrices.
     161       
     162            sage: matrix(RDF, 2**31-1, 2, sparse=True)
     163            2147483647 x 2 sparse matrix over Real Double Field
     164        """
     165       
     166        global scipy_sparse, numpy
     167        if scipy_sparse is None:
     168            from scipy import sparse as scipy_sparse
     169        if numpy is None:
     170            import numpy
     171        matrix_sparse.Matrix_sparse.__init__(self, parent)
     172        if self._nrows >= 2**31 or self._ncols >= 2**31:
     173            raise NotImplementedError(
     174                'Number of rows and columns must currently be smaller than 2^31, even on 64-bit systems. This only concerns RDF/CDF sparse matrices.')
     175            # This is a limitation in csc_matrix, even on 64-bit systems
     176        self._entries_csc = None # this will be set either by __init__ or _unpickle
     177
     178        self._sage_dtype = parent._base
     179        if parent._base == RDF:
     180            self._numpy_dtype = numpy.dtype('float64')
     181            self._python_dtype = float
     182        elif parent._base == CDF:
     183            self._numpy_dtype = numpy.dtype('complex128')
     184            self._python_dtype = complex
     185        else:
     186            raise ValueError("Parent has invalid base ring")
     187
     188    def __dealloc__(self):
     189        """ Deallocate any memory that was initialized."""
     190        return
     191
     192    def __richcmp__(Matrix self, right, int op):  # always need for mysterious reasons.
     193        return self._richcmp(right, op)
     194
     195    def __init__(self, parent, entries, copy, coerce):
     196        """
     197        Fill the matrix with entries. 
     198
     199        EXAMPLES::
     200       
     201            sage: matrix(RDF,3,sparse=True)
     202            [0.0 0.0 0.0]
     203            [0.0 0.0 0.0]
     204            [0.0 0.0 0.0]
     205            sage: matrix(RDF,3,range(9),sparse=True)
     206            [0.0 1.0 2.0]
     207            [3.0 4.0 5.0]
     208            [6.0 7.0 8.0]
     209            sage: matrix(CDF,3,3,2,sparse=True)
     210            [2.0   0   0]
     211            [  0 2.0   0]
     212            [  0   0 2.0]
     213
     214        TESTS::
     215       
     216            sage: matrix(RDF,3,0,sparse=True)
     217            []
     218            sage: matrix(RDF,3,3,0,sparse=True)
     219            [0.0 0.0 0.0]
     220            [0.0 0.0 0.0]
     221            [0.0 0.0 0.0]
     222            sage: matrix(RDF,3,3,1,sparse=True)
     223            [1.0 0.0 0.0]
     224            [0.0 1.0 0.0]
     225            [0.0 0.0 1.0]
     226            sage: matrix(RDF,3,3,2,sparse=True)
     227            [2.0 0.0 0.0]
     228            [0.0 2.0 0.0]
     229            [0.0 0.0 2.0]
     230            sage: matrix(CDF,3,0,sparse=True)
     231            []
     232            sage: matrix(CDF,3,3,0,sparse=True)
     233            [0 0 0]
     234            [0 0 0]
     235            [0 0 0]
     236            sage: matrix(CDF,3,3,1,sparse=True)
     237            [1.0   0   0]
     238            [  0 1.0   0]
     239            [  0   0 1.0]
     240            sage: matrix(CDF,3,3,range(9),sparse=True)
     241            [  0 1.0 2.0]
     242            [3.0 4.0 5.0]
     243            [6.0 7.0 8.0]
     244            sage: matrix(CDF,2,2,[CDF(1+I)*j for j in range(4)],sparse=True)
     245            [          0 1.0 + 1.0*I]
     246            [2.0 + 2.0*I 3.0 + 3.0*I]
     247
     248        The following invokes the constructor with list arguments::
     249       
     250            sage: M=MatrixSpace(RDF, 2, 2, sparse=True)
     251            sage: M(matrix(RDF, 2, 2, 1, sparse=False))
     252            [1.0 0.0]
     253            [0.0 1.0]
     254         
     255        """
     256        cdef Py_ssize_t i, j
     257        global scipy_sparse
     258       
     259        if self._nrows == 0 or self._ncols == 0:
     260            return # leave _entries_csc as None
     261        shape = (self._nrows, self._ncols)
     262
     263        # Below is taken verbatim from matrix_generic_sparse.pyx.
     264        # Something should be done to avoid copying code like this
     265        # across all the matrix classes and a couple of places in
     266        # matrix_space.py, but I'm too lazy and make another copy...
     267        if isinstance(entries, list) and len(entries) > 0 and \
     268                hasattr(entries[0], "is_vector"):
     269            entries = _convert_sparse_entries_to_dict(entries)
     270
     271        if isinstance(entries, list):
     272            if len(entries) != self.nrows() * self.ncols():
     273                raise TypeError, "entries has the wrong length"
     274            x = entries
     275            entries = {}
     276            k = 0
     277            for i from 0 <= i < self._nrows:
     278                for j from 0 <= j < self._ncols:
     279                    if x[k] != 0:
     280                        entries[(int(i),int(j))] = x[k]
     281                    k = k + 1
     282            copy = False
     283        # End copy
     284       
     285        if isinstance(entries, dict):
     286            M = scipy_sparse.dok_matrix(shape, dtype=self._numpy_dtype)
     287            if coerce:
     288                for key, value in entries.iteritems():
     289                    M[key] = self._python_dtype(value)
     290                    z = self._python_dtype(value)
     291                    if z != 0:
     292                        M[key] = z
     293            else:
     294                # TODO: Is the coerce flag present for speed or to catch nonsense?
     295                # Could drop check below if the purpose is speed
     296                for key, value in entries.iteritems():
     297                    if value.parent() != self._sage_dtype:
     298                        raise TypeError, ("all entries must be in %s" % self._sage_dtype)
     299                    z = self._python_dtype(value)
     300                    if z != 0:
     301                        M[key] = z
     302            self._entries_csc = M.tocsc()
     303        else:
     304            if entries is None:
     305                z = 0
     306            else:
     307                try:
     308                    z = self._python_dtype(entries)
     309                except TypeError:
     310                    raise TypeError, ("entries must be coercible to dict or %s" %
     311                                      self._python_dtype.__name__)
     312            if z == 0:
     313                self._entries_csc = scipy_sparse.csc_matrix(shape, dtype=self._numpy_dtype)
     314            else:
     315                # Diagonal matrix of constant.
     316                # TODO: Faster direct construction.
     317                M = scipy_sparse.eye(self._nrows, self._ncols, dtype=self._numpy_dtype)
     318                M *= z
     319                self._entries_csc = M.tocsc()
     320       
     321    cdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, object value):
     322        """
     323        Set the (i,j) entry to value without any bounds checking,
     324        mutability checking, etc.
     325        """
     326        # TODO: More efficient setter if not changing sparsity pattern
     327       
     328        # We call the self._python_dtype function on the value since
     329        # numpy does not know how to deal with complex numbers other
     330        # than the built-in complex number type.
     331        with warnings.catch_warnings():
     332            warnings.simplefilter('ignore', scipy_sparse.SparseEfficiencyWarning)
     333            self._entries_csc[i, j] = self._python_dtype(value)
     334
     335    cdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j):
     336        """
     337        Get the (i,j) entry without any bounds checking, etc.
     338        """
     339        # TODO: More efficient getter
     340       
     341        return self._sage_dtype(self._entries_csc[i, j])
     342
     343    def _pickle(self):
     344        """
     345        See _unpickle.
     346        """
     347        version = 0
     348        data = self._entries_csc.data
     349        indices = self._entries_csc.indices
     350        indptr = self._entries_csc.indptr
     351        return (data, indices, indptr), version
     352       
     353    def _unpickle(self, data, int version):
     354        """
     355        Format version 0: Returns a tuple (data, indices, indptr),
     356        three NumPy arrays containing the matrix in CSC format.
     357       
     358        TESTS::
     359       
     360            sage: a = matrix(RDF,2,range(4), sparse=True)
     361            sage: loads(dumps(a)) == a
     362            True
     363            sage: a = matrix(CDF,2,range(4), sparse=True)
     364            sage: loads(dumps(a)) == a
     365            True
     366        """       
     367        if version == 0:
     368            self._entries_csc = scipy_sparse.csc_matrix(data, (self._nrows, self._ncols))
     369        else:
     370            raise RuntimeError, "unknown matrix version (=%s)"%version
     371       
     372    def _dict(self):
     373        """
     374        Unsafe version of the dict method, mainly for internal use.
     375        See ancestor docstring.
     376
     377        TESTS::
     378       
     379            sage: matrix(RDF, 3, 0, sparse=True)._dict()
     380            {}
     381       
     382            sage: x=matrix(RDF, 3, 3, 2, sparse=True)._dict().items(); x.sort(); x
     383            [((0, 0), 2.0), ((1, 1), 2.0), ((2, 2), 2.0)]
     384        """
     385        # TODO: Lots of potential for optimization
     386        if self._ncols == 0 or self._nrows == 0:
     387            return {}       
     388        d = self.fetch('dict')
     389        if d is not None:
     390            return d
     391        d = {}
     392        num_d = self._entries_csc.todok()
     393        for key, value in num_d.iteritems():
     394            d[key] = self._sage_dtype(value)
     395        self.cache('dict', d)
     396        return d
     397
     398    def _nonzero_positions_by_row(self, copy=True):
     399        """
     400        See ancestor docstring.
     401       
     402        TESTS::
     403       
     404            sage: x=matrix(RDF, 3, 3, {(1, 2) : 1, (2, 1) : 1, (0, 1): 1}, sparse=True)
     405            sage: x._nonzero_positions_by_row()
     406            [(0, 1), (1, 2), (2, 1)]
     407        """
     408        # TODO: Lots of potential for optimization
     409        x = self.fetch('nonzero_positions')
     410        if not x is None:
     411            if copy:
     412                return list(x)
     413            return x
     414        x = self._entries_csc.todok().keys()
     415        x.sort()
     416        self.cache('nonzero_positions', x)
     417        return x
     418       
     419    def _nonzero_positions_by_column(self, copy=True):
     420        """
     421        See ancestor docstring.
     422       
     423        TESTS::
     424       
     425            sage: x=matrix(RDF, 3, 3, {(1, 2) : 1, (2, 1) : 1, (0, 1): 1}, sparse=True)
     426            sage: x._nonzero_positions_by_column()
     427            [(0, 1), (2, 1), (1, 2)]
     428        """
     429        # TODO: Lots of potential for optimization   
     430        x = self.fetch('nonzero_positions_by_column')
     431        if not x is None:
     432            if copy:
     433                return list(x)
     434            return x
     435        x = self._dict().keys()
     436        x.sort(key=_transpose_index)
     437        self.cache('nonzero_positions_by_column', x)
     438        return x
     439
     440    cdef Matrix_double_sparse _new(self, int nrows=-1, int ncols=-1):
     441        """
     442        Return a new uninitialized matrix with same parent as self.
     443        Note that _entries_csc will be None in the returned matrix
     444        and *must* be set by the caller (unless nrows or ncols == 0).
     445
     446        INPUT
     447            nrows -- (default self._nrows) number of rows in returned matrix
     448            ncols -- (default self._ncols) number of columns in returned matrix
     449        """
     450        cdef Matrix_double_sparse m
     451        if nrows == -1:
     452            nrows = self._nrows
     453        if ncols == -1:
     454            ncols = self._ncols
     455        parent = self.matrix_space(nrows, ncols)
     456        m = self.__class__.__new__(self.__class__,parent,None,None,None)
     457        return m
     458
     459    def __copy__(self):
     460        r"""
     461        Returns a new copy of this matrix.
     462
     463        EXAMPLES::
     464       
     465            sage: a = matrix(RDF,1,3, [1,2,-3],sparse=True)
     466            sage: a
     467            [ 1.0  2.0 -3.0]
     468            sage: b = a.__copy__()
     469            sage: b
     470            [ 1.0  2.0 -3.0]
     471            sage: b is a
     472            False
     473            sage: b == a
     474            True
     475            sage: b[0,0] = 3
     476            sage: a[0,0] # note that a hasn't changed
     477            1.0
     478        """
     479        if self._nrows == 0 or self._ncols == 0:
     480            return self.new_matrix(self._nrows, self._ncols)
     481
     482        cdef Matrix_double_sparse A
     483        A = self._new(self._nrows, self._ncols)
     484        A._entries_csc = self._entries_csc.copy()
     485        if self.subdivisions is not None:
     486            A.subdivide(*self.get_subdivisions())
     487        return A
     488
     489    def transpose(self):
     490        """
     491        Returns the transpose of self, without changing self.
     492       
     493        EXAMPLES:
     494
     495        We create a matrix and compute its transpose. Subdivisions
     496        are preserved::
     497       
     498            sage: A = matrix(RDF, 2, 3, range(6), sparse=True)
     499            sage: A.subdivide(None, 1)
     500            sage: A
     501            [0.0|1.0 2.0]
     502            [3.0|4.0 5.0]
     503            sage: B = A.transpose(); B
     504            [0.0 3.0]
     505            [-------]
     506            [1.0 4.0]
     507            [2.0 5.0]
     508
     509        Note that changing the transpose does not affect the
     510        original matrix::
     511
     512            sage: B[0,0] = 10
     513            sage: A
     514            [0.0|1.0 2.0]
     515            [3.0|4.0 5.0]
     516       
     517        """
     518        cdef Matrix_double_sparse trans = self._new(self._ncols, self._nrows)
     519        trans._entries_csc = self._entries_csc.transpose(copy=True).tocsc()
     520        if self.subdivisions is not None:
     521            # TODO: Consistently move this up in the hierarchy everywhere
     522            # and introduce _transpose_, so that e.g. other sparse matrices
     523            # still have subdivisions applied
     524            row_divs, col_divs = self.get_subdivisions()
     525            trans.subdivide(col_divs, row_divs)
     526        return trans
     527       
     528
     529
     530    cdef int _cmp_c_impl(self, Element right) except -2:
     531        """
     532        Compare two matrices. See parent docstring for details.
     533
     534        TESTS::
     535
     536            sage: matrix(RDF,2,[1,2,3,4],sparse=True) < matrix(RDF,2,[3,2,3,4],sparse=True)
     537            True
     538            sage: matrix(RDF,2,[1,2,3,4],sparse=True) > matrix(RDF,2,[3,2,3,4],sparse=True)
     539            False
     540            sage: matrix(RDF,2,[0,2,3,4],sparse=True) < matrix(RDF,2,[0,3,3,4],sparse=True)
     541            True
     542            sage: matrix(RDF,2,{(0,0): 1, (1,1): 2},sparse=True) == matrix(RDF,2,[1,0,0,2])
     543            True
     544        """
     545        # TODO: I intended to do this, realized that it requires more than 30 minutes,
     546        # but have already written the tests...so cop out for now.
     547        # Also the tests uncovered a serious bug in __init__ and should be left in.
     548        return cmp(self._dict(), right._dict())
     549
     550    ########################################################################
     551    # Arithmetic operations. Just use scipy.sparse. The result is converted
     552    # using tocsc() in all cases, but that just returns "self" in the
     553    # current implementation (but the API technically appears free to
     554    # to return other forms)
     555    #
     556    # - _lmul_
     557    # - _add_
     558    # - __neg__
     559    # - _generic_matrix_times_matrix_
     560    ########################################################################
     561   
     562    cpdef ModuleElement _lmul_(self, RingElement right):
     563        """
     564        EXAMPLES::
     565       
     566            sage: A = matrix(RDF, 2, 3, range(6), sparse=True)
     567            sage: A * RDF(0.5)
     568            [0.0 0.5 1.0]
     569            [1.5 2.0 2.5]
     570            sage: 1/2 * A
     571            [0.0 0.5 1.0]
     572            [1.5 2.0 2.5]
     573        """
     574        cdef Matrix_double_sparse M = self._new()
     575        M._entries_csc = self._entries_csc * self._python_dtype(right)
     576        return M
     577
     578    cpdef ModuleElement _add_(self, ModuleElement right):
     579        """
     580        Add two matrices together.
     581
     582        EXAMPLES:
     583            sage: A = matrix(RDF,3,range(1,10),sparse=True)
     584            sage: A+A
     585            [ 2.0  4.0  6.0]
     586            [ 8.0 10.0 12.0]
     587            [14.0 16.0 18.0]
     588        """
     589        if self._nrows == 0 or self._ncols == 0:
     590            return self.__copy__()
     591
     592        cdef Matrix_double_sparse M, _left, _right
     593        _left = self
     594        _right = right
     595
     596        M = self._new()
     597        M._entries_csc = (_left._entries_csc + _right._entries_csc).tocsc()
     598        return M
     599
     600    cpdef ModuleElement _sub_(self, ModuleElement right):
     601        """
     602        Return self - right
     603
     604        EXAMPLES:
     605            sage: A = matrix(RDF,3,range(1,10),sparse=True)
     606            sage: (A-A).is_zero()
     607            True
     608        """
     609        if self._nrows == 0 or self._ncols == 0:
     610            return self.__copy__()
     611
     612        cdef Matrix_double_sparse M,_right,_left
     613        _right = right
     614        _left = self
     615
     616        M = self._new()
     617        M._entries_csc = (_left._entries_csc - _right._entries_csc).tocsc()
     618        return M
     619
     620    def __neg__(self):
     621        """
     622        Negate this matrix
     623
     624        EXAMPLES:
     625            sage: A = matrix(RDF,3,range(1,10),sparse=True)
     626            sage: -A
     627            [-1.0 -2.0 -3.0]
     628            [-4.0 -5.0 -6.0]
     629            [-7.0 -8.0 -9.0]
     630            sage: B = -A ; (A+B).is_zero()
     631            True
     632        """
     633        if self._nrows == 0 or self._ncols == 0:
     634            return self.__copy__()
     635
     636        cdef Matrix_double_sparse M
     637        M = self._new()
     638        M._entries_csc = (-self._entries_csc).tocsc()
     639        return M
     640
     641    cdef sage.structure.element.Matrix _generic_matrix_times_matrix_(self,
     642            sage.structure.element.Matrix right, Parent codomain):
     643        """
     644        Multiply ``self`` with ``right``.
     645
     646        Overrides generic version to provide:
     647       
     648         - faster code paths for sparse times dense
     649         - complex times real uses less memory (indices not copied)
     650
     651        TESTS:
     652
     653            sage: from sage.matrix.matrix_double_sparse import _verbose
     654            sage: A = matrix(RDF,3,range(1,10),sparse=True)
     655            sage: B = matrix(RDF,3,range(1,13),sparse=True)
     656            sage: AC = A.change_ring(CDF); BC = B.change_ring(CDF)
     657           
     658            sage: with _verbose(): A*B
     659            Sparse times sparse
     660            [ 38.0  44.0  50.0  56.0]
     661            [ 83.0  98.0 113.0 128.0]
     662            [128.0 152.0 176.0 200.0]
     663
     664        Different sparseness takes different code path::
     665       
     666            sage: with _verbose(): A*B.dense_matrix()
     667            Sparse times dense
     668            [ 38.0  44.0  50.0  56.0]
     669            [ 83.0  98.0 113.0 128.0]
     670            [128.0 152.0 176.0 200.0]
     671
     672        Multiplying with matrix with non-double base ring works::
     673       
     674            sage: A*B.change_ring(ZZ)
     675            [ 38.0  44.0  50.0  56.0]
     676            [ 83.0  98.0 113.0 128.0]
     677            [128.0 152.0 176.0 200.0]
     678
     679        Different base rings takes different code path::
     680       
     681            sage: with _verbose(): A.change_ring(CDF) * B
     682            Sparse times sparse
     683            [ 38.0  44.0  50.0  56.0]
     684            [ 83.0  98.0 113.0 128.0]
     685            [128.0 152.0 176.0 200.0]
     686            sage: with _verbose(): A.change_ring(CDF) * B.dense_matrix()
     687            Sparse times dense
     688            [ 38.0  44.0  50.0  56.0]
     689            [ 83.0  98.0 113.0 128.0]
     690            [128.0 152.0 176.0 200.0]
     691           
     692        A.dense_matrix()*B is handled in ``Matrix_double_dense``.
     693        """
     694        cdef sage.structure.element.Matrix left
     695        cdef Matrix_double_sparse result_sparse, right_sparse
     696        cdef Matrix_double_dense result_dense, right_dense
     697
     698        if self._nrows == 0 or self._ncols == 0 or right._nrows == 0 or right._ncols == 0:
     699            return codomain.zero_matrix()
     700
     701        base = codomain._base
     702        if base is not RDF and base is not CDF:
     703            # Don't know how to handle base ring. Convert left matrix and
     704            # forward call.
     705            left = self.change_ring(codomain._base)
     706            return left._generic_matrix_times_matrix(right, codomain)
     707
     708        # Note: Codomain base ring may be either of RDF or CDF.
     709        # Codomain is dense if and only if right.is_dense().
     710        if right.is_dense_c():
     711            if VERBOSE: print 'Sparse times dense'
     712            result_dense = create_uninitialized_Matrix_double_dense(codomain)
     713            result_dense._matrix_numpy = (self._entries_csc *
     714                                          (<Matrix_double_dense?>right)._matrix_numpy)
     715            assert result_dense._matrix_numpy.dtype == result_dense._numpy_dtype
     716            return result_dense
     717        else:
     718            if VERBOSE: print 'Sparse times sparse'
     719            if right._parent._base is not base:
     720                right = right.change_ring(codomain._base)
     721            result_sparse = create_uninitialized_Matrix_double_sparse(codomain)
     722            result_sparse._entries_csc = (self._entries_csc *
     723                                          (<Matrix_double_sparse?>right)._entries_csc).tocsc()
     724            assert result_sparse._entries_csc.dtype == result_sparse._numpy_dtype
     725            return result_sparse
     726
     727    cdef Vector _matrix_times_vector_(self, Vector v):
     728        """
     729        TESTS::
     730
     731            sage: A = matrix(RDF, [[1,0,1],[0,1,1]], sparse=True)
     732            sage: v = vector(RDF, [1, 2, 3],  sparse=True)
     733            sage: z = A * v; z; parent(z)
     734            (4.0, 5.0)
     735            Sparse vector space of dimension 2 over Real Double Field
     736
     737            sage: z = A * v.dense_vector(); z; parent(z)
     738            (4.0, 5.0)
     739            Sparse vector space of dimension 2 over Real Double Field
     740
     741            sage: z = vector(CDF, [1+1j, 1-1j], sparse=True) * A; z; parent(z)
     742            (1.0 + 1.0*I, 1.0 - 1.0*I, 2.0)
     743            Sparse vector space of dimension 3 over Complex Double Field
     744
     745        """
     746        # This keeps existing Sage semantics, but is wrong from a
     747        # numerical POV -- should discuss whether it is OK to return a
     748        # dense vector when multiplying with a dense vector (that must
     749        # be fixed in the coercions though, not just here.)
     750        from sage.modules.free_module_element import vector
     751        cdef Vector_double_dense v_dense = v.dense_vector()
     752        result = self._entries_csc * v_dense._vector_numpy
     753        result = vector(self._sage_dtype, result)
     754        return result.sparse_vector()
     755
     756
     757    ########################################################################
     758    # Solvers and factorizations.
     759    #
     760    # Re-implement right_solve using LU factorizations -- the default
     761    # implementation is useless in a numerical setting (begins by finding
     762    # matrix rank...), and since solving is so different for numerical
     763    # matrices, we redefine docstring rather than introduce a
     764    # _right_solve_
     765    ########################################################################
     766
     767    # SuperLU is a great library with lots of cool options (like doing
     768    # iterative refinement to get a solution within a specified error tolerance,
     769    # and caching many parts of the computation). I need to think more about how
     770    # I want to expose it properly to user-space in a friendly manner which doesn't
     771    # make it awkward to maintain backwards compatability later.
     772    # The below docstring is a start though, so I keep it in the source.
     773    # -- dagss
     774   
     775   
     776##     def solve_right(self, B, check=True, algorithm="lu", **kwds):
     777##         r"""
     778##         If self is a matrix `A`, then this function returns a
     779##         vector or matrix `X` such that `A X = B`. If
     780##         `B` is a vector then `X` is a vector and if
     781##         `B` is a matrix, then `X` is a matrix.
     782
     783##         Multiple algorithms are available and work well for different
     784##         matrices.  Temporary results such as computed decompositions
     785##         are cached and available to subsequent calls to
     786##         :meth:`solve_right` or :meth:`solve_left`. Consider making the
     787##         matrix immutable (:meth:`set_immutable`) prior to solving.
     788
     789##         Currently supported algorithms:
     790           
     791##           - "lu" - (the default) (P)LU factorization. Usually needs no tuning
     792##             arguments, see :meth:`LU` for the arguments which are
     793##             accepted.
     794           
     795##           - "naive" - Row reduction to echelon form. This is almost always a bad
     796##             idea. Takes no tuning arguments.
     797
     798##         .. note::
     799
     800##            In Sage one can also write ``A \backslash  B`` for
     801##            ``A.solve_right(B)``, i.e., Sage implements the "the
     802##            MATLAB/Octave backslash operator".
     803       
     804##         INPUT:
     805               
     806##         -  ``B`` - a matrix or vector
     807##         -  ``check`` - bool (default: True) - if False and self
     808##            is nonsquare, may not raise an error message even if there is no
     809##            solution. This is faster but more dangerous.
     810##         -  ``algorithm`` - string (default: "lu"). Which algorithm to use.
     811##            See above list.
     812##         - ``**kwds`` - Tuning arguments can optionally supplied. The arguments
     813##            varies depending on selected algorithm; see above.
     814       
     815##         OUTPUT: a matrix or vector
     816       
     817##         .. seealso::
     818
     819##            :meth:`solve_left`
     820##            :meth:`solve_left`
     821
     822           
     823##         """
     824##         # Ideas:
     825##         # - Cholesky solver (CHOLMOD/SuiteSparse)
     826##         # - QR (through SuiteSparse)
     827##         # - For LU: UMFPACK in addition to SuperLU
     828##         # - Perhaps make available particularily versatile algorithms from
     829##         #   the endless list of iterative methods...and perhaps not
     830       
     831   
     832       
     833
     834def _transpose_index(idx):
     835    return (idx[1], idx[0])

  • sage/matrix/matrix_generic_sparse.pxd 

    diff -r 977f52f47b09 -r 25f413c99ea3 sage/matrix/matrix_generic_sparse.pxd
    a b  
    44    cdef object _entries
    55    cdef object _zero
    66
     7cdef _convert_sparse_entries_to_dict(entries)

  • sage/matrix/matrix_space.py 

    diff -r 977f52f47b09 -r 25f413c99ea3 sage/matrix/matrix_space.py
    a b  
    821821            <type 'sage.matrix.matrix_integer_sparse.Matrix_integer_sparse'>
    822822            sage: type(matrix(SR, 2, 2, 0))
    823823            <type 'sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense'>
     824            sage: type(matrix(RDF, 2, 2, 0))
     825            <type 'sage.matrix.matrix_real_double_dense.Matrix_real_double_dense'>
     826            sage: type(matrix(CDF, 2, 2, 0))
     827            <type 'sage.matrix.matrix_complex_double_dense.Matrix_complex_double_dense'>
     828            sage: type(matrix(RDF, 2, 2, 0, sparse=True))
     829            <type 'sage.matrix.matrix_double_sparse.Matrix_double_sparse'>
     830            sage: type(matrix(CDF, 2, 2, 0, sparse=True))
     831            <type 'sage.matrix.matrix_double_sparse.Matrix_double_sparse'>
    824832        """
    825833        R = self.base_ring()
    826834        if self.is_dense():
    … …  
    860868                return matrix_rational_sparse.Matrix_rational_sparse
    861869            elif sage.rings.integer_ring.is_IntegerRing(R):
    862870                return matrix_integer_sparse.Matrix_integer_sparse
     871            elif R == sage.rings.real_double.RDF or R == sage.rings.complex_double.CDF:
     872                import matrix_double_sparse
     873                return matrix_double_sparse.Matrix_double_sparse
    863874            # the default
    864875            return matrix_generic_sparse.Matrix_generic_sparse
    865876           
    … …  
    10831094                    j = k // self.__ncols
    10841095                    new_x.append( x[ i*self.__nrows + j ] )
    10851096                x = new_x
    1086            
    10871097        return self.__matrix_class(self, entries=x, copy=copy, coerce=coerce)
    10881098     
    10891099    def matrix_space(self, nrows=None, ncols=None, sparse=False):