Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

""" 

Matrix windows 

""" 

  

#***************************************************************************** 

# Distributed under the terms of the GNU General Public License (GPL) 

# as published by the Free Software Foundation; either version 2 of 

# the License, or (at your option) any later version. 

# http://www.gnu.org/licenses/ 

#***************************************************************************** 

  

from cpython.tuple cimport * 

  

######################################################################### 

# Generic matrix windows, which are used for block echelon and strassen # 

# algorithms. # 

######################################################################### 

cdef class MatrixWindow: 

############################ 

# creation and initialization 

############################ 

  

cpdef MatrixWindow new_matrix_window(MatrixWindow self, Matrix matrix, 

Py_ssize_t row, Py_ssize_t col, 

Py_ssize_t n_rows, Py_ssize_t n_cols): 

""" 

This method is here only to provide a fast cdef way of 

constructing new matrix windows. The only implicit assumption 

is that self._matrix and matrix are over the same base ring 

(so share the zero). 

""" 

cdef type t = type(self) 

cdef MatrixWindow M = <MatrixWindow>t.__new__(t) 

M._matrix = matrix 

M._row = row 

M._col = col 

M._nrows = n_rows 

M._ncols = n_cols 

M._cached_zero = self._cached_zero 

return M 

  

def __init__(MatrixWindow self, Matrix matrix, 

Py_ssize_t row, Py_ssize_t col, Py_ssize_t nrows, Py_ssize_t ncols): 

self._matrix = matrix 

self._row = row 

self._col = col 

self._nrows = nrows 

self._ncols = ncols 

  

cdef object _zero(self): 

if self._cached_zero is None: 

self._cached_zero = self._matrix.base_ring()(0) # expensive 

return self._cached_zero 

  

cpdef MatrixWindow matrix_window(MatrixWindow self, Py_ssize_t row, Py_ssize_t col, 

Py_ssize_t n_rows, Py_ssize_t n_cols): 

""" 

Returns a matrix window relative to this window of the 

underlying matrix. 

""" 

if row == 0 and col == 0 and n_rows == self._nrows and n_cols == self._ncols: 

return self 

return self.new_matrix_window(self._matrix, self._row + row, self._col + col, n_rows, n_cols) 

  

cpdef new_empty_window(MatrixWindow self, Py_ssize_t nrows, Py_ssize_t ncols): 

a = self._matrix.new_matrix(nrows, ncols) 

return self.new_matrix_window(a, 0, 0, nrows, ncols) 

  

def __repr__(self): 

return "Matrix window of size %s x %s at (%s,%s):\n%s"%( 

self._nrows, self._ncols, self._row, self._col, self._matrix) 

  

############################ 

# Getting and setting entries 

############################ 

  

def set(MatrixWindow self, MatrixWindow src): 

if self._matrix._parent.base_ring() is not src._matrix._parent.base_ring(): 

raise TypeError("Parents must be equal.") 

self.set_to(src) 

  

cpdef set_unsafe(self, Py_ssize_t i, Py_ssize_t j, x): 

self._matrix.set_unsafe(i + self._row, j + self._col, x) 

  

cpdef get_unsafe(self, Py_ssize_t i, Py_ssize_t j): 

return self._matrix.get_unsafe(i + self._row, j + self._col) 

  

def __setitem__(self, ij, x): 

cdef Py_ssize_t i, j 

if isinstance(ij, tuple): 

# ij is a tuple, so we get i and j efficiently, construct corresponding integer entry. 

if PyTuple_Size(ij) != 2: 

raise IndexError("index must be an integer or pair of integers") 

i = <object> PyTuple_GET_ITEM(ij, 0) 

j = <object> PyTuple_GET_ITEM(ij, 1) 

if i<0 or i >= self._nrows or j<0 or j >= self._ncols: 

raise IndexError("matrix index out of range") 

self.set_unsafe(i, j, x) 

else: 

# If ij is not a tuple, coerce to an integer and set the row. 

i = ij 

for j from 0 <= j < self._ncols: 

self.set_unsafe(i, j, x) 

  

def __getitem__(self, ij): 

cdef Py_ssize_t i, j 

cdef object x 

  

if isinstance(ij, tuple): 

# ij is a tuple, so we get i and j efficiently, construct corresponding integer entry. 

if PyTuple_Size(ij) != 2: 

raise IndexError("index must be an integer or pair of integers") 

i = <object> PyTuple_GET_ITEM(ij, 0) 

j = <object> PyTuple_GET_ITEM(ij, 1) 

if i<0 or i >= self._nrows or j<0 or j >= self._ncols: 

raise IndexError("matrix index out of range") 

return self.get_unsafe(i, j) 

else: 

# If ij is not a tuple, coerce to an integer and get the row. 

i = ij 

return self.row(i) 

  

cpdef matrix(MatrixWindow self): 

""" 

Returns the underlying matrix that this window is a view of. 

""" 

return self._matrix 

  

  

cpdef to_matrix(MatrixWindow self): 

""" 

Returns an actual matrix object representing this view. 

""" 

cdef MatrixWindow w 

a = self._matrix.new_matrix(self._nrows, self._ncols) # zero matrix 

w = self.new_matrix_window(a, 0, 0, self._nrows, self._ncols) 

w.set_to(self) 

return a 

  

def nrows(MatrixWindow self): 

return self._nrows 

  

def ncols(MatrixWindow self): 

return self._ncols 

  

cpdef set_to(MatrixWindow self, MatrixWindow A): 

""" 

Change self, making it equal A. 

""" 

cdef Py_ssize_t i, j 

if self._nrows != A._nrows or self._ncols != A._ncols: 

raise ArithmeticError("incompatible dimensions") 

for i from 0 <= i < self._nrows: 

for j from 0 <= j < self._ncols: 

self.set_unsafe(i, j, A.get_unsafe(i, j)) 

return 0 

  

cpdef set_to_zero(MatrixWindow self): 

cdef Py_ssize_t i, j 

z = self._zero() 

for i from 0 <= i < self._nrows: 

for j from 0 <= j < self._ncols: 

self.set_unsafe(i, j, z) 

  

cpdef add(MatrixWindow self, MatrixWindow A): 

cdef Py_ssize_t i, j 

if self._nrows != A._nrows or self._ncols != A._ncols: 

raise ArithmeticError("incompatible dimensions") 

for i from 0 <= i < self._nrows: 

for j from 0 <= j < self._ncols: 

self.set_unsafe(i, j, self.get_unsafe(i, j) + A.get_unsafe(i, j)) 

  

cpdef subtract(MatrixWindow self, MatrixWindow A): 

cdef Py_ssize_t i, j 

if self._nrows != A._nrows or self._ncols != A._ncols: 

raise ArithmeticError("incompatible dimensions") 

for i from 0 <= i < self._nrows: 

for j from 0 <= j < self._ncols: 

self.set_unsafe(i, j, self.get_unsafe(i, j) - A.get_unsafe(i, j)) 

  

cpdef set_to_sum(MatrixWindow self, MatrixWindow A, MatrixWindow B): 

cdef Py_ssize_t i, j 

if self._nrows != A._nrows or self._ncols != A._ncols: 

raise ArithmeticError("incompatible dimensions") 

if self._nrows != B._nrows or self._ncols != B._ncols: 

raise ArithmeticError("incompatible dimensions") 

for i from 0 <= i < self._nrows: 

for j from 0 <= j < self._ncols: 

self.set_unsafe(i, j, A.get_unsafe(i, j) + B.get_unsafe(i, j)) 

  

cpdef set_to_diff(MatrixWindow self, MatrixWindow A, MatrixWindow B): 

cdef Py_ssize_t i, j 

for i from 0 <= i < self._nrows: 

for j from 0 <= j < self._ncols: 

self.set_unsafe(i, j, A.get_unsafe(i, j) - B.get_unsafe(i, j)) 

  

cpdef set_to_prod(MatrixWindow self, MatrixWindow A, MatrixWindow B): 

cdef Py_ssize_t i, j, k 

if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols: 

raise ArithmeticError("incompatible dimensions") 

for i from 0 <= i < A._nrows: 

for j from 0 <= j < B._ncols: 

s = self._zero() 

for k from 0 <= k < A._ncols: 

s = s + A.get_unsafe(i, k) * B.get_unsafe(k, j) 

self.set_unsafe(i, j, s) 

  

cpdef add_prod(MatrixWindow self, MatrixWindow A, MatrixWindow B): 

cdef Py_ssize_t i, j, k 

if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols: 

raise ArithmeticError("incompatible dimensions") 

for i from 0 <= i < A._nrows: 

for j from 0 <= j < B._ncols: 

s = self.get_unsafe(i, j) 

for k from 0 <= k < A._ncols: 

s = s + A.get_unsafe(i, k) * B.get_unsafe(k, j) 

self.set_unsafe(i, j, s) 

  

cpdef subtract_prod(MatrixWindow self, MatrixWindow A, MatrixWindow B): 

cdef Py_ssize_t i, j, k 

if A._ncols != B._nrows or self._nrows != A._nrows or self._ncols != B._ncols: 

raise ArithmeticError("incompatible dimensions") 

for i from 0 <= i < A._nrows: 

for j from 0 <= j < B._ncols: 

s = self.get_unsafe(i, j) 

for k from 0 <= k < A._ncols: 

s = s - A.get_unsafe(i, k) * B.get_unsafe(k, j) 

self.set_unsafe(i, j, s) 

  

cpdef swap_rows(MatrixWindow self, Py_ssize_t a, Py_ssize_t b): 

self._matrix.swap_rows_c(self._row + a, self._row + b) 

  

def echelon_in_place(MatrixWindow self): 

""" 

Calculate the echelon form of this matrix, returning the list of pivot columns 

""" 

echelon = self.to_matrix() 

echelon.echelonize() # TODO: read only, only need to copy pointers 

self.set_to(echelon.matrix_window()) 

return echelon.pivots() 

  

cpdef bint element_is_zero(MatrixWindow self, Py_ssize_t i, Py_ssize_t j): 

return self._matrix.get_unsafe(i+self._row, j+self._col) == self._zero()