Coverage for local/lib/python2.7/site-packages/sage/matrix/strassen.pyx: 93%

sage: A = matrix(QQ, 7, [5, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 3, 1, 0, -1, 0, 0, -1, 0, 1, 2, -1, 1, 0, -1, 0, 1, 3, -1, 1, 0, 0, -2, 0, 2, 0, 1, 0, 0, -1, 0, 1, 0, 1])

sage: B = A.__copy__(); B._echelon_strassen(1); B

[ 1 0 0 0 0 0 0]

[ 0 1 0 -1 0 1 0]

[ 0 0 1 0 0 0 0]

[ 0 0 0 0 1 0 0]

[ 0 0 0 0 0 0 1]

[ 0 0 0 0 0 0 0]

sage: C = A.__copy__(); C._echelon_strassen(2); C == B

True

sage: C = A.__copy__(); C._echelon_strassen(4); C == B

True

sage: n = 32; A = matrix(Integers(389),n,range(n^2))

sage: B = A.__copy__(); B._echelon_in_place_classical()

sage: C = A.__copy__(); C._echelon_strassen(2)

sage: B == C

True

TESTS::

sage: A = matrix(Integers(7), 4, 4, [1,2,0,3,0,0,1,0,0,1,0,0,0,0,0,1])

sage: B = A.__copy__(); B._echelon_in_place_classical()

sage: C = A.__copy__(); C._echelon_strassen(2)

sage: B == C

True

sage: A = matrix(Integers(7), 4, 4, [1,0,5,0,2,0,3,6,5,1,2,6,4,6,1,1])

sage: B = A.__copy__(); B._echelon_in_place_classical()

sage: C = A.__copy__(); C._echelon_strassen(2) #indirect doctest

sage: B == C

True

AUTHORS:

- Robert Bradshaw

"""

if cutoff < 1:

raise ValueError("cutoff must be at least 1")

sig_on()

strassen_echelon_c(A, cutoff, A._matrix._strassen_default_cutoff(A._matrix))

sig_off()

cdef strassen_echelon_c(MatrixWindow A, Py_ssize_t cutoff, Py_ssize_t mul_cutoff):

# The following notation will be used in the comments below, which should be understood to give

# the general idea of what's going on, as if there were no inconvenient non-pivot columns.

# The original matrix is given by [ A B ]

# [ C D ]

# For compactness, let A' denote the inverse of A

# top_left, top_right, bottom_left, and bottom_right loosely correspond to A, B, C, and D respectively,

# however, the "cut" between the top and bottom rows need not be the same.

cdef Py_ssize_t nrows, ncols

nrows = A.nrows()

ncols = A.ncols()

if (nrows <= cutoff or ncols <= cutoff):

return list(A.echelon_in_place())

cdef Py_ssize_t top_h, bottom_cut, bottom_h, bottom_start, top_cut

cdef Py_ssize_t prev_pivot_count

cdef Py_ssize_t split

split = nrows / 2

cdef MatrixWindow top, bottom, top_left, top_right, bottom_left, bottom_right, clear

top = A.matrix_window(0, 0, split, ncols)

bottom = A.matrix_window(split, 0, nrows-split, ncols)

top_pivots = strassen_echelon_c(top, cutoff, mul_cutoff)

# effectively "multiplied" top row by A^{-1}

# [ I A'B ]

# [ C D ]

top_pivot_intervals = int_range(top_pivots)

top_h = len(top_pivots)

if top_h == 0:

# [ 0 0 ]

# the whole top is a zero matrix, [ C D ]. Run echelon on the bottom

bottom_pivots = strassen_echelon_c(bottom, cutoff, mul_cutoff)

# [ 0 0 ]

# we now have [ I C'D ], proceed to sorting

else:

bottom_cut = max(top_pivots) + 1

bottom_left = bottom.matrix_window(0, 0, nrows-split, bottom_cut)

if top_h == ncols:

bottom.set_to_zero()

# [ I ]

# [ 0 ]

# proceed to sorting

else:

if bottom_cut == top_h:

clear = bottom_left

else:

clear = bottom_left.to_matrix().matrix_from_columns(top_pivots).matrix_window() # TODO: read only, can I do this faster? Also below

# Subtract off C time top from the bottom_right

if bottom_cut < ncols:

bottom_right = bottom.matrix_window(0, bottom_cut, nrows-split, ncols-bottom_cut)

subtract_strassen_product(bottom_right, clear, top.matrix_window(0, bottom_cut, top_h, ncols-bottom_cut), mul_cutoff);

# [ I A'B ]

# [ * D - CA'B ]

# Now subtract off C times the top from the bottom_left (pivots -> 0)

if bottom_cut == top_h:

bottom_left.set_to_zero()

bottom_start = bottom_cut

else:

for cols in top_pivot_intervals:

bottom_left.matrix_window(0, cols[0], nrows-split, cols[1]).set_to_zero()

non_pivots = int_range(0, bottom_cut) - top_pivot_intervals

for cols in non_pivots:

if cols[0] == 0: continue

prev_pivot_count = len(top_pivot_intervals - int_range(cols[0]+cols[1], bottom_cut - cols[0]+cols[1]))

subtract_strassen_product(bottom_left.matrix_window(0, cols[0], nrows-split, cols[1]),

clear.matrix_window(0, 0, nrows-split, prev_pivot_count),

top.matrix_window(0, cols[0], prev_pivot_count, cols[1]),

mul_cutoff)

bottom_start = non_pivots._intervals[0][0]

# [ I A'B ]

# [ 0 D - CA'B ]

# Now recursively do echelon form on the bottom

bottom_pivots_rel = strassen_echelon_c(bottom.matrix_window(0, bottom_start, nrows-split, ncols-bottom_start), cutoff, mul_cutoff)

# [ I A'B ]

# [ 0 I F ]

bottom_pivots = []

for pivot in bottom_pivots_rel:

bottom_pivots.append(pivot + bottom_start)

bottom_h = len(bottom_pivots)

if bottom_h == 0:

pass

# [ I A'B ]

# [ 0 0 ]

# proceed to sorting

else:

# [ I A'B ] = [ I E G ]

# let [ 0 I F ] = [ 0 I F ]

top_cut = max(max(bottom_pivots) + 1, bottom_cut)

# Note: left with respect to leftmost non-zero column of bottom

top_left = top.matrix_window(0, bottom_start, top_h, top_cut - bottom_start)

if bottom_h + top_h < ncols:

if top_cut - bottom_start == bottom_h:

clear = top_left

else:

clear = top_left.to_matrix().matrix_from_columns(bottom_pivots_rel).matrix_window()

# subtract off E times bottom from top right

if top_cut < ncols:

top_right = top.matrix_window(0, top_cut, top_h, ncols - top_cut)

subtract_strassen_product(top_right, clear, bottom.matrix_window(0, top_cut, bottom_h, ncols - top_cut), mul_cutoff);

# [ I * G - EF ]

# [ 0 I F ]

# Now subtract of E times bottom from top left

if top_cut - bottom_start == bottom_h:

top_left.set_to_zero()

else:

bottom_pivot_intervals = int_range(bottom_pivots)

non_pivots = int_range(bottom_start, top_cut - bottom_start) - bottom_pivot_intervals - top_pivot_intervals

for cols in non_pivots:

if cols[0] == 0: continue

prev_pivot_count = len(bottom_pivot_intervals - int_range(cols[0]+cols[1], top_cut - cols[0]+cols[1]))

subtract_strassen_product(top.matrix_window(0, cols[0], top_h, cols[1]),

clear.matrix_window(0, 0, top_h, prev_pivot_count),

bottom.matrix_window(0, cols[0], prev_pivot_count, cols[1]),

mul_cutoff)

for cols in bottom_pivot_intervals:

top.matrix_window(0, cols[0], top_h, cols[1]).set_to_zero()

# [ I 0 G - EF ]

# [ 0 I F ]

# proceed to sorting

# subrows already sorted...maybe I could do this more efficiently in cases with few pivot columns (e.g. merge sort)

pivots = top_pivots

pivots.extend(bottom_pivots)

pivots.sort()

cdef Py_ssize_t i, cur_row

for cur_row from 0 <= cur_row < len(pivots):

pivot = pivots[cur_row]

for i from cur_row <= i < nrows:

if not A.element_is_zero(i, pivot):

break

if i > cur_row and i < nrows:

A.swap_rows(i, cur_row)

return pivots

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

# lots of room for optimization....

# eventually, should I just pass these around rather than lists of ints for pivots?

# would need new from_cols

class int_range:

r"""

Represent a list of integers as a list of integer intervals.

.. NOTE::

Repetitions are not considered.

Useful class for dealing with pivots in the strassen echelon, could

have much more general application

INPUT:

It can be one of the following:

- ``indices`` - integer, start of the unique interval

- ``range`` - integer, length of the unique interval

- ``indices`` - list of integers, the integers to wrap into intervals

- ``indices`` - None (default), shortcut for an empty list

OUTPUT:

An instance of ``int_range``, i.e. a list of pairs ``(start, length)``.

EXAMPLES:

From a pair of integers::

sage: from sage.matrix.strassen import int_range

sage: int_range(2, 4)

[(2, 4)]

Default::

sage: int_range()

[]

From a list of integers::

sage: int_range([1,2,3,4])

[(1, 4)]

sage: int_range([1,2,3,4,6,7,8])

[(1, 4), (6, 3)]

sage: int_range([1,2,3,4,100,101,102])

[(1, 4), (100, 3)]

sage: int_range([1,1000,2,101,3,4,100,102])

[(1, 4), (100, 3), (1000, 1)]

Repetitions are not considered::

sage: int_range([1,2,3])

[(1, 3)]

sage: int_range([1,1,1,1,2,2,2,3])

[(1, 3)]

AUTHORS:

- Robert Bradshaw

"""

def __init__(self, indices=None, range=None):

r"""

See ``sage.matrix.strassen.int_range`` for full documentation.

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: int_range(2, 4)

[(2, 4)]

"""

if indices is None:

self._intervals = []

return

elif not range is None:

self._intervals = [(int(indices), int(range))]

else:

self._intervals = []

if len(indices) == 0:

return

indices.sort()

start = None

last = None

for ix in indices:

if last is None:

start = ix

elif ix-last > 1:

self._intervals.append((start, last-start+1))

start = ix

last = ix

self._intervals.append((start, last-start+1))

def __repr__(self):

r"""

String representation.

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: int_range([4,5,6,20,21,22,23])

[(4, 3), (20, 4)]

sage: int_range([])

[]

"""

return str(self._intervals)

def intervals(self):

r"""

Return the list of intervals.

OUTPUT:

A list of pairs of integers.

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: I = int_range([4,5,6,20,21,22,23])

sage: I.intervals()

[(4, 3), (20, 4)]

sage: type(I.intervals())

<... 'list'>

"""

return self._intervals

def to_list(self):

r"""

Return the (sorted) list of integers represented by this object.

OUTPUT:

A list of integers.

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: I = int_range([6,20,21,4,5,22,23])

sage: I.to_list()

[4, 5, 6, 20, 21, 22, 23]

sage: I = int_range(34, 9)

sage: I.to_list()

[34, 35, 36, 37, 38, 39, 40, 41, 42]

Repetitions are not considered::

sage: I = int_range([1,1,1,1,2,2,2,3])

sage: I.to_list()

[1, 2, 3]

"""

all = []

for iv in self._intervals:

for i in range(iv[0], iv[0]+iv[1]):

all.append(i)

return all

def __iter__(self):

r"""

Return an iterator over the intervals.

OUTPUT:

iterator

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: I = int_range([6,20,21,4,5,22,23])

sage: it = iter(I)

sage: next(it)

(4, 3)

sage: next(it)

(20, 4)

sage: next(it)

Traceback (most recent call last):

...

StopIteration

"""

return iter(self._intervals)

def __len__(self):

r"""

Return the number of integers represented by this object.

OUTPUT:

Python integer

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: I = int_range([6,20,21,4,5,22,23])

sage: len(I)

sage: I = int_range([1,1,1,1,2,2,2,3])

sage: len(I)

"""

len = 0

for iv in self._intervals:

len = len + iv[1]

return int(len)

def __add__(self, right):

r"""

Return the union of ``self`` and ``right``.

INPUT:

- ``right`` - an instance of ``int_range``

OUTPUT:

An instance of ``int_range``

.. NOTE::

Yes, this two could be a lot faster...

Basically, this class is for abstracting away what I was trying

to do by hand in several places

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: I = int_range([1,1,1,1,2,2,2,3])

sage: J = int_range([6,20,21,4,5,22,23])

sage: I + J

[(1, 6), (20, 4)]

"""

all = self.to_list()

all.extend(right.to_list())

return int_range(all)

def __sub__(self, right):

r"""

Return the set difference of ``self`` and ``right``.

INPUT:

- ``right`` - an instance of ``int_range``.

OUTPUT:

An instance of ``int_range``.

.. NOTE::

Yes, this two could be a lot faster...

Basically, this class is for abstracting away what I was trying

to do by hand in several places

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: I = int_range([1,2,3,4,5])

sage: J = int_range([6,20,21,4,5,22,23])

sage: J - I

[(6, 1), (20, 4)]

"""

all = self.to_list()

for i in right.to_list():

if i in all:

all.remove(i)

return int_range(all)

def __mul__(self, right):

r"""

Return the intersection of ``self`` and ``right``.

INPUT:

- ``right`` - an instance of ``int_range``.

OUTPUT:

An instance of ``int_range``.

EXAMPLES::

sage: from sage.matrix.strassen import int_range

sage: I = int_range([1,2,3,4,5])

sage: J = int_range([6,20,21,4,5,22,23])

sage: J * I

[(4, 2)]

"""

intersection = []

all = self.to_list()

for i in right.to_list():

if i in all:

intersection.append(i)

return int_range(intersection)

# Useful test code:

def test(n, m, R, c=2):

r"""

INPUT:

- ``n`` - integer

- ``m`` - integer

- ``R`` - ring

- ``c`` - integer (optional, default:2)

EXAMPLES::

sage: from sage.matrix.strassen import test

sage: for n in range(5):

....: print("{} {}".format(n, test(2*n,n,Frac(QQ['x']),2)))

0 True

1 True

2 True

3 True

4 True

"""

from sage.matrix.all import matrix

A = matrix(R,n,m,range(n*m))

B = A.__copy__(); B._echelon_in_place_classical()

C = A.__copy__(); C._echelon_strassen(c)

return B == C

# This stuff gets tested extensively elsewhere, and the functions

# below aren't callable now without using Pyrex.

## todo: doc cutoff parameter as soon as I work out what it really means

## EXAMPLES:

## The following matrix dimensions are chosen especially to exercise the

## eight possible parity combinations that could occur while subdividing

## the matrix in the strassen recursion. The base case in both cases will

## be a (4x5) matrix times a (5x6) matrix.

## TODO -- the doctests below are currently not

## tested/enabled/working -- enable them when linear algebra

## restructing gets going.

## sage: dim1 = 64; dim2 = 83; dim3 = 101

## sage: R = MatrixSpace(QQ, dim1, dim2)

## sage: S = MatrixSpace(QQ, dim2, dim3)

## sage: T = MatrixSpace(QQ, dim1, dim3)

## sage: A = R.random_element(range(-30, 30))

## sage: B = S.random_element(range(-30, 30))

## sage: C = T(0)

## sage: D = T(0)

## sage: A_window = A.matrix_window(0, 0, dim1, dim2)

## sage: B_window = B.matrix_window(0, 0, dim2, dim3)

## sage: C_window = C.matrix_window(0, 0, dim1, dim3)

## sage: D_window = D.matrix_window(0, 0, dim1, dim3)

## sage: from sage.matrix.strassen import strassen_window_multiply

## sage: strassen_window_multiply(C_window, A_window, B_window, 2) # use strassen method

## sage: D_window.set_to_prod(A_window, B_window) # use naive method

## sage: C_window == D_window

## True

## sage: dim1 = 79; dim2 = 83; dim3 = 101

## sage: R = MatrixSpace(QQ, dim1, dim2)

## sage: S = MatrixSpace(QQ, dim2, dim3)

## sage: T = MatrixSpace(QQ, dim1, dim3)

## sage: A = R.random_element(range(30))

## sage: B = S.random_element(range(30))

## sage: C = T(0)

## sage: D = T(0)

## sage: A_window = A.matrix_window(0, 0, dim1, dim2)

## sage: B_window = B.matrix_window(0, 0, dim2, dim3)

## sage: C_window = C.matrix_window(0, 0, dim1, dim3)

## sage: strassen_window_multiply(C, A, B, 2) # use strassen method

## sage: D.set_to_prod(A, B) # use naive method

## sage: C == D

## True

Coverage for local/lib/python2.7/site-packages/sage/matrix/strassen.pyx : 93%

245 statements 228 run 17 missing 0 excluded