source: sage/matrix/matrix_window_modn_dense.pyx @ 2847:5077101063da

Revision 2847:5077101063da, 10.8 KB checked in by Robert Bradshaw <robertwb@…>, 6 years ago (diff)

Multi-modular matrix window optimization, gather
Line 
1include "../ext/cdefs.pxi"
2include "../ext/stdsage.pxi"
3
4import matrix_window
5cimport matrix_window
6
7from sage.matrix.matrix_modn_dense import Matrix_modn_dense
8from sage.matrix.matrix_modn_dense cimport Matrix_modn_dense
9
10# TODO: what if we used the actual pointers for our loop variables?
11
12cdef class MatrixWindow_modn_dense(matrix_window.MatrixWindow):
13     
14    cdef new_empty_window(self, Py_ssize_t nrows, Py_ssize_t ncols):
15# the current code is all python, goes through inits, and possibly creates a parent...
16# can we get away with something faster?
17#       a = MatrixWindow_modn_dense.__new__(self._parent, )?
18       a = self._matrix.new_matrix(nrows, ncols)
19       return self.new_matrix_window(a, 0, 0, nrows, ncols)
20       
21    cdef set_to(MatrixWindow_modn_dense self, MatrixWindow A):
22        """
23        Change self, making it equal A.
24        """
25        cdef Py_ssize_t i, j
26        cdef mod_int** self_rows
27        cdef mod_int** A_rows
28        if self._nrows != A._nrows or self._ncols != A._ncols:
29            raise ArithmeticError, "incompatible dimensions"
30        self_rows = ( <Matrix_modn_dense> self._matrix ).matrix
31        A_rows = ( <Matrix_modn_dense> A._matrix ).matrix
32        for i from 0 <= i < self._nrows:
33            memcpy(self_rows[i+self._row] + self._col, A_rows[i+A._row] + A._col, self._ncols * sizeof(mod_int*))
34               
35    cdef set_to_zero(MatrixWindow_modn_dense self):
36        cdef Py_ssize_t i, j
37        cdef mod_int** rows
38        rows = ( <Matrix_modn_dense> self._matrix ).matrix
39        for i from self._row <= i < self._row + self._nrows:
40            memset(rows[i] + self._col, 0, self._ncols * sizeof(mod_int*))
41
42    cdef add(self, MatrixWindow A):
43        cdef Py_ssize_t i, j
44        cdef mod_int p
45        cdef mod_int* self_row
46        cdef mod_int* A_row
47        if self._nrows != A._nrows or self._ncols != A._ncols:
48            raise ArithmeticError, "incompatible dimensions"
49        p = ( <Matrix_modn_dense> self._matrix ).p
50        for i from 0 <= i < self._nrows:
51            self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col
52            A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col
53            for j from 0 <= j < self._ncols:
54                self_row[j] += A_row[j]
55                if self_row[j] >= p:
56                    self_row[j] -= p
57   
58    cdef subtract(self, MatrixWindow A):
59        cdef Py_ssize_t i, j
60        cdef mod_int p
61        cdef mod_int* self_row
62        cdef mod_int* A_row
63        if self._nrows != A._nrows or self._ncols != A._ncols:
64            raise ArithmeticError, "incompatible dimensions"
65        p = ( <Matrix_modn_dense> self._matrix ).p
66        for i from 0 <= i < self._nrows:
67            self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col
68            A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col
69            for j from 0 <= j < self._ncols:
70                if self_row[j] >= A_row[j]:
71                    self_row[j] -= A_row[j]
72                else:
73                    self_row[j] += p - A_row[j]
74       
75    cdef set_to_sum(self, MatrixWindow A, MatrixWindow B):
76        cdef Py_ssize_t i, j
77        cdef mod_int p
78        cdef mod_int* self_row
79        cdef mod_int* A_row
80        cdef mod_int* B_row
81        if self._nrows != A._nrows or self._ncols != A._ncols:
82            raise ArithmeticError, "incompatible dimensions"
83        if self._nrows != B._nrows or self._ncols != B._ncols:
84            raise ArithmeticError, "incompatible dimensions"
85        p = ( <Matrix_modn_dense> self._matrix ).p
86        for i from 0 <= i < self._nrows:
87            self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col
88            A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col
89            B_row    = ( <Matrix_modn_dense>    B._matrix ).matrix[i +    B._row] + B._col
90            for j from 0 <= j < self._ncols:
91                self_row[j] = A_row[j] + B_row[j]
92                if self_row[j] >= p:
93                    self_row[j] -= p
94
95    cdef set_to_diff(self, MatrixWindow A, MatrixWindow B):
96        cdef Py_ssize_t i, j
97        cdef mod_int p
98        cdef mod_int* self_row
99        cdef mod_int* A_row
100        cdef mod_int* B_row
101        if self._nrows != A._nrows or self._ncols != A._ncols:
102            raise ArithmeticError, "incompatible dimensions"
103        if self._nrows != B._nrows or self._ncols != B._ncols:
104            raise ArithmeticError, "incompatible dimensions"
105        p = ( <Matrix_modn_dense> self._matrix ).p
106        for i from 0 <= i < self._nrows:
107            self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col
108            A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col
109            B_row    = ( <Matrix_modn_dense>    B._matrix ).matrix[i +    B._row] + B._col
110            for j from 0 <= j < self._ncols:
111                if A_row[j] > B_row[j]:
112                    self_row[j] = A_row[j] - B_row[j]
113                else:
114                    self_row[j] = p + A_row[j] - B_row[j]
115               
116    cdef set_to_prod(self, MatrixWindow A, MatrixWindow B):
117        cdef Py_ssize_t i, j, k, gather, top, A_ncols
118        cdef mod_int p, s
119        cdef mod_int* self_row
120        cdef mod_int* A_row
121        cdef mod_int** B_matrix_off
122        if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols:
123            raise ArithmeticError, "incompatible dimensions"
124        B_matrix_off = ( <Matrix_modn_dense> B._matrix ).matrix + B._row
125        p = ( <Matrix_modn_dense> self._matrix ).p
126        gather = ( <Matrix_modn_dense> self._matrix ).gather
127        A_ncols = A._ncols
128       
129        if gather <= 1:
130            for i from 0 <= i < A._nrows:
131                self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col
132                A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col
133                for j from 0 <= j < B._ncols:
134                    s = 0
135                    for k from 0 <= k < A._ncols:
136                        s = (s + A_row[k] * B_matrix_off[k][j+B._col]) % p
137                    self_row[j] = s
138        else:
139            for i from 0 <= i < A._nrows:
140                self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i + self._row] + self._col
141                A_row    = ( <Matrix_modn_dense>    A._matrix ).matrix[i +    A._row] + A._col
142                for j from 0 <= j < B._ncols:
143                    s = 0
144                    k = 0
145                    while k < A_ncols:
146                        top = k + gather
147                        if top > A_ncols:
148                            top = A_ncols
149                        for k from k <= k < top: # = min(k+gather, A._ncols)
150                            s += A_row[k] * B_matrix_off[k][j+B._col]
151                        s %= p
152                    self_row[j] = s
153       
154    cdef add_prod(self, MatrixWindow A, MatrixWindow B):
155        # TODO: gather
156        cdef Py_ssize_t i, j, k, A_ncols
157        cdef mod_int p, s
158        cdef mod_int* self_row
159        cdef mod_int* A_row
160        cdef mod_int* B_row
161        if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols:
162            raise ArithmeticError, "incompatible dimensions"
163        p = ( <Matrix_modn_dense> self._matrix ).p
164        gather = ( <Matrix_modn_dense> self._matrix ).gather
165        A_ncols = A._ncols
166       
167        if gather <= 1:
168            for i from 0 <= i < A._nrows:
169                self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i+self._row] + self._col
170                A_row = ( <Matrix_modn_dense> A._matrix ).matrix[i+B._row] + B._col
171                B_row = ( <Matrix_modn_dense> B._matrix ).matrix[i+A._row] + A._col
172                for j from 0 <= j < B._ncols:
173                    s = self_row[j]
174                    for k from 0 <= k < A._ncols:
175                        s = ( s + A_row[k] * B_row[j] ) % p
176                    self_row[j] = s
177                   
178        else:
179            for i from 0 <= i < A._nrows:
180                self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i+self._row] + self._col
181                A_row = ( <Matrix_modn_dense> A._matrix ).matrix[i+B._row] + B._col
182                B_row = ( <Matrix_modn_dense> B._matrix ).matrix[i+A._row] + A._col
183                for j from 0 <= j < B._ncols:
184                    s = self_row[j]
185                    k = 0
186                    while k < A_ncols:
187                        top = k + gather
188                        if top > A_ncols:
189                            top = A_ncols
190                        for k from k <= k < top: # = min(k+gather, A._ncols)
191                            s += A_row[k] * B_matrix_off[k][j+B._col]
192                        s %= p
193                    self_row[j] = s
194                   
195               
196    cdef subtract_prod(self, MatrixWindow A, MatrixWindow B):
197        # TODO: gather
198        cdef Py_ssize_t i, j, k, A_ncols
199        cdef mod_int p, s
200        cdef mod_int* self_row
201        cdef mod_int* A_row
202        cdef mod_int* B_row
203        if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols:
204            raise ArithmeticError, "incompatible dimensions"
205        p = ( <Matrix_modn_dense> self._matrix ).p
206        A_ncols = A._ncols
207       
208        if gather <= 1:
209            for i from 0 <= i < A._nrows:
210                self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i+self._row] + self._col
211                A_row = ( <Matrix_modn_dense> A._matrix ).matrix[i+B._row] + B._col
212                B_row = ( <Matrix_modn_dense> B._matrix ).matrix[i+A._row] + A._col
213                for j from 0 <= j < B._ncols:
214                    s = self_row[j]
215                    for k from 0 <= k < A._ncols:
216                        s = s + p - ( A_row[k] * B_row[j] ) % p
217                    self_row[j] = s
218                   
219        else:
220            for i from 0 <= i < A._nrows:
221                self_row = ( <Matrix_modn_dense> self._matrix ).matrix[i+self._row] + self._col
222                A_row = ( <Matrix_modn_dense> A._matrix ).matrix[i+B._row] + B._col
223                B_row = ( <Matrix_modn_dense> B._matrix ).matrix[i+A._row] + A._col
224                for j from 0 <= j < B._ncols:
225                    s = self_row[j] + p * gather # because they're unsigned
226                    k = 0
227                    while k < A_ncols:
228                        top = k + gather
229                        if top > A_ncols:
230                            top = A_ncols
231                        for k from k <= k < top: # = min(k+gather, A._ncols)
232                            s -= A_row[k] * B_matrix_off[k][j+B._col]
233                        s = s % p  +  p * gather # because they're unsigne
234                    self_row[j] = s % p
235
236    cdef int element_is_zero(self, Py_ssize_t i, Py_ssize_t j):
237        return (<Matrix_modn_dense>self._matrix).matrix[i+self._row][j+self._col] == 0
238   
239
Note: See TracBrowser for help on using the repository browser.