source: sage/matrix/matrix_dense.pyx @ 1637:abb74e940771

Revision 1637:abb74e940771, 11.4 KB checked in by William Stein <wstein@…>, 7 years ago (diff)

(1) Pyrex introspection (this also required modifying pyrex!); (2) matrix reorganization and implementation.
Line 
1#############################################
2## Generic *DENSE* matrices over any commutative ring
3#############################################
4
5def _convert_dense_entries_to_list(entries):
6    # Create matrix from a list of vectors
7    e = []
8    for v in entries:
9        e = e+ v.list()
10    copy = False
11    return e
12
13include "../ext/interrupt.pxi"
14
15cimport matrix_generic
16import matrix_generic
17
18cdef class Matrix_dense(matrix_generic.Matrix):
19    """
20    The \\class{Matrix_dense} class derives from
21    \\class{Matrix}, and defines functionality for dense matrices over
22    any base ring.  Matrices are represented by a list of elements in
23    the base ring, and element access operations are implemented in
24    this class.
25    """
26    def __init__(self, parent,
27                 entries = None,
28                 copy = True,
29                 coerce = True):
30        matrix_generic.Matrix.__init__(self, parent)
31        self._nrows = parent.nrows()
32        self._ncols = parent.ncols()
33        cdef int i, n
34        if entries:
35            if not isinstance(entries, list):
36                raise TypeError, 'entries must be a list'
37            if not (coerce or copy):
38                self.__entries = entries
39            else:
40                self.__entries = [None]*(self._nrows*self._ncols)                       
41                n = len(entries)
42                if coerce:
43                    R = parent.base_ring()
44                    for i from 0 <= i < n:
45                        self.__entries[i] = R(entries[i])
46                else:
47                    for i from 0 <= i < n:
48                        self.__entries[i] = entries[i]
49        else:
50            self.__entries = [None]*(self._nrows*self._ncols)
51
52        self._init_row_indices()
53
54    def _init_row_indices(self):
55        self._row_indices = <int*> PyMem_Malloc(sizeof(int*) * self._nrows)
56        n = 0
57        for i from 0 <= i < self._nrows:
58            self._row_indices[i] = n
59            n = n + self._ncols
60       
61    def  __dealloc__(self):
62        PyMem_Free(self._row_indices)
63       
64    def __getitem__(self, ij):
65        """
66        INPUT:
67            A[i, j] -- the i,j of A, and
68            A[i]    -- the i-th row of A.
69        """
70        #self._require_mutable()   # todo
71        cdef int i, j
72        i, j = ij
73        if i < 0 or i >= self._nrows:
74            raise IndexError
75        return self.__entries[self._row_indices[i] + j]
76   
77    def __setitem__(self, ij, value):
78        """
79        INPUT:
80            A[i, j] = value -- set the (i,j) entry of A
81            A[i] = value    -- set the ith row of A
82        """
83        #self._require_mutable()   # todo
84        cdef int i, j
85        i, j = ij
86        if i < 0 or i >= self._nrows:
87            raise IndexError
88        self.__entries[self._row_indices[i] + j] = value
89
90    def matrix_window(self, int row=0, int col=0, int nrows=-1, int ncols=-1):
91        if nrows == -1:
92            nrows = self._nrows
93            ncols = self._ncols
94        return MatrixWindow(self, row, col, nrows, ncols)
95
96    cdef classical_multiply_cdef(self, matrix_generic.Matrix right):
97        """
98        Multiply the matrices self and right using the classical $O(n^3)$
99        algorithm.
100
101        EXAMPLES
102
103            sage: include the 0 rows and 0 columns cases  -- do dense examples
104        """
105        cdef int i, j, k, m, n, r, nr, nc, snc
106        cdef object v
107
108        if self._ncols != right._nrows:
109            raise IndexError, "Number of columns of self must equal number of rows of other."
110
111        cdef Matrix_dense A
112        nr = self._nrows
113        nc = right._ncols
114        snc = self._ncols
115
116        R = self.base_ring()
117        P = self.matrix_space(nr, nc)
118        A = Matrix_dense(P)
119        v = A.__entries
120        zero = R(0)
121        r = 0
122        for i from 0 <= i < nr:
123            m = self._row_indices[i]
124            for j from 0 <= j < nc:
125                z = zero
126                for k from 0 <= k < snc:
127                    z = z + self.__entries[m + k] * right.__entries[right._row_indices[k]+j]
128                v[r] = z
129                r = r + 1
130        return A
131
132    def _entries(self):
133        return self.__entries
134
135    def list(self):
136        return self.__entries
137   
138    def antitranspose(self):
139        f = []
140        e = self.list()
141        (nc, nr) = (self.ncols(), self.nrows())
142        for j in reversed(xrange(nc)):
143            for i in reversed(xrange(nr)):
144                f.append(e[i*nc + j])
145        return self.new_matrix(nrows = nc, ncols = nr,
146                               entries = f, copy=False, coerce=False)
147
148    def transpose(self):
149        """
150        Returns the transpose of self, without changing self.
151       
152        EXAMPLES:
153        We create a matrix, compute its transpose, and note that the
154        original matrix is not changed.
155            sage: M = MatrixSpace(QQ,  2)
156            sage: A = M([1,2,3,4])
157            sage: B = A.transpose()
158            sage: print B
159            [1 3]
160            [2 4]
161            sage: print A
162            [1 2]
163            [3 4]
164        """
165        f = []
166        e = self.list()
167        (nc, nr) = (self.ncols(), self.nrows())
168        for j in xrange(nc):       
169            for i in xrange(nr):
170                f.append(e[i*nc + j])
171        return self.new_matrix(nrows = nc, ncols = nr,
172                               entries = f, copy=False,
173                               coerce=False)
174
175
176
177cdef class MatrixWindow:
178
179    def __init__(MatrixWindow self, matrix, int row, int col, int nrows, int ncols):
180        self._matrix = matrix
181        self._row = row
182        self._col = col
183        self._nrows = nrows
184        self._ncols = ncols
185
186    def __repr__(self):
187        return "Matrix window of size %s x %s at (%s,%s):\n%s"%(
188            self._nrows, self._ncols, self._row, self._col, self._matrix)
189   
190    def matrix(MatrixWindow self):
191        """
192        Returns the underlying matrix that this window is a view of.
193
194        EXAMPLES:
195       
196        """
197        return self._matrix
198       
199       
200    def to_matrix(MatrixWindow self):
201        """
202        Returns an actual matrix object representing this view. (Copy)
203        """
204        entries = self.list()
205        return self._matrix.new_matrix(self._nrows, self._ncols, entries=entries,
206                                       coerce=False, copy=False)
207
208    def list(MatrixWindow self):
209        v = self._matrix.__entries
210        w = [None]*(self._nrows*self._ncols)
211        cdef int i, j, k, l
212        k = 0
213        for i from 0 <= i < self._nrows:
214            l = self._matrix._row_indices[i]
215            for j from 0 <= j < self._ncols:
216                w[k] = v[l + j]
217                k = k + 1
218        return w
219
220    def matrix_window(MatrixWindow self, int row, int col, int n_rows, int n_cols):
221        """
222        Returns a matrix window relative to this window of the underlying matrix.
223        """
224        return self._matrix.matrix_window(self._matrix, _row + row, _col + col, n_rows, n_cols)
225
226    def nrows(MatrixWindow self):
227        return self._nrows
228   
229    def ncols(MatrixWindow self):
230        return self._ncols
231   
232   
233    def set_to(MatrixWindow self, MatrixWindow A):
234        cdef int i, j
235        cdef int start, self_ix
236        cdef int A_ix
237        for i from 0 <= i < self._nrows:
238            A_ix = self._matrix._row_indices[i+A._row] + a._col
239            for self_ix from self._row_indices[i+self._row] + self._col <= self_ix < self_start + self._ncols:
240                self._matrix.__entries[self_ix] = A._matrix.__entries[A_ix]
241                A_ix = A_ix + 1
242               
243               
244    def set_to_zero(MatrixWindow self):
245        cdef int i, j
246        cdef int start, self_ix
247        zero = self._matrix.base_ring(0)
248        for i from 0 <= i < self._nrows:
249            for self_ix from self._row_indices[i+self._row] + self._col <= self_ix < self_start + self._ncols:
250                self._matrix.__entries[self_ix] = zero
251
252
253    def add(MatrixWindow self, MatrixWindow A):
254        cdef int i, j
255        cdef int start, self_ix
256        cdef int A_ix
257        for i from 0 <= i < self._nrows:
258            A_ix = A._matrix._row_indices[i+A._row] + A._col
259            for self_ix from self._row_indices[i+self._row] + self._col <= self_ix < self_start + self._ncols:
260                self._matrix.__entries[self_ix] = slef._matrix.__entries[self_ix] + A._matrix.__entries[A_ix]
261                A_ix = A_ix + 1
262   
263   
264    def subtract(MatrixWindow self, MatrixWindow A):
265        cdef int i, j
266        cdef int start, self_ix
267        cdef int A_ix
268        for i from 0 <= i < self._nrows:
269            A_ix = A._matrix._row_indices[i+A._row] + a._col
270            for self_ix from self._row_indices[i+self._row] + self._col <= self_ix < self_start + self._ncols:
271                self._matrix.__entries[self_ix] = self._matrix.__entries[self_ix] - A._matrix.__entries[A_ix]
272                A_ix = A_ix + 1
273               
274
275    def set_to_sum(MatrixWindow self, MatrixWindow A, MatrixWindow B):
276        cdef int i, j
277        cdef int start, self_ix
278        cdef int A_ix
279        for i from 0 <= i < self._nrows:
280            A_ix = A._matrix._row_indices[i+A._row] + A._col
281            B_ix = B._matrix._row_indices[i+B._row] + B._col
282            for self_ix from self._row_indices[i+self._row] + self._col <= self_ix < self_start + self._ncols:
283                self._matrix.__entries[self_ix] = A._matrix.__entries[A_ix] + B._matrix.__entries[B_ix]
284                A_ix = A_ix + 1
285                B_ix = B_ix + 1
286
287
288    def set_to_diff(MatrixWindow self, MatrixWindow A, MatrixWindow B):
289        cdef int i, j
290        cdef int start, self_ix
291        cdef int A_ix
292        for i from 0 <= i < self._nrows:
293            A_ix = A._matrix._row_indices[i+A._row] + A._col
294            B_ix = B._matrix._row_indices[i+B._row] + B._col
295            for self_ix from self._row_indices[i+self._row] + self._col <= self_ix < self_start + self._ncols:
296                self._matrix.__entries[self_ix] = A._matrix.__entries[A_ix] - B._matrix.__entries[B_ix]
297                A_ix = A_ix + 1
298                B_ix = B_ix + 1
299   
300   
301    def set_to_prod(MatrixWindow self, MatrixWindow A, MatrixWindow B):
302        cdef int i, j, k
303        cdef int start, self_ix
304        for i from 0 <= i < A._nrows:
305            for j from 0 <= j < B._ncols:
306                sum = A._matrix.__entries[ A._row_indices[A._row+i]+A._col ] *B._matrix.__entries[ B._row_indices[B._row]+B._col+j ]
307                for k from 1 <= k < A._ncols:
308                    sum = sum + A._matrix.__entries[ A._row_indices[A._row+i]+A._col + k ] * B._matrix.__entries[ B._row_indices[B._row+k]+B._col+j ]
309                self._matrix.__entries[ self._row_indices[self_.row+i]+self._col+j ] = sum
310               
311               
312    def add_prod(MatrixWindow self, MatrixWindow A, MatrixWindow B):
313        cdef int i, j, k
314        cdef int start, self_ix
315        for i from 0 <= i < A._nrows:
316            for j from 0 <= j < B._ncols:
317                sum = A._matrix.__entries[ A._row_indices[A._row+i]+A._col ] *B._matrix.__entries[ B._row_indices[B._row]+B._col+j ]
318                for k from 1 <= k < A._ncols:
319                    sum = sum + A._matrix.__entries[ A._row_indices[A._row+i]+A._col + k ] * B._matrix.__entries[ B._row_indices[B._row+k]+B._col+j ]
320                self._matrix.__entries[ self._row_indices[self_.row+i]+self._col+j ] = self._matrix.__entries[ self._row_indices[self_.row+i]+self._col+j ] + sum
321
322
323    def new_empty_window(MatrixWindow self, int nrows, int ncols, zero=True):
324        return self._matrix.new_matrix(nrows,ncols, zero=zero).matrix_window()
Note: See TracBrowser for help on using the repository browser.