Ticket #10876: trac_10876-elementary-matrices-v2.patch

File trac_10876-elementary-matrices-v2.patch, 13.6 KB (added by rbeezer, 10 years ago)
  • sage/matrix/all.py

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1299281755 28800
    # Node ID 50434db8e46c5c37045765bf3f2af6a92f5fa8fd
    # Parent  6ea10f8b8ce2e45bf4587ceafa59eddc844a2389
    10876: constructor for elementary matrices
    
    diff -r 6ea10f8b8ce2 -r 50434db8e46c sage/matrix/all.py
    a b  
    11from matrix_space import MatrixSpace, is_MatrixSpace
    2 from constructor import matrix, Matrix, column_matrix, random_matrix, diagonal_matrix, identity_matrix, block_matrix, block_diagonal_matrix, jordan_block, zero_matrix, ones_matrix
     2from constructor import matrix, Matrix, column_matrix, random_matrix, diagonal_matrix, identity_matrix, block_matrix, block_diagonal_matrix, jordan_block, zero_matrix, ones_matrix, elementary_matrix
    33from matrix import is_Matrix
    44from berlekamp_massey import berlekamp_massey
    55
  • sage/matrix/constructor.py

    diff -r 6ea10f8b8ce2 -r 50434db8e46c sage/matrix/constructor.py
    a b  
    16411641    one = ring(1)
    16421642    return matrix_space.MatrixSpace(ring, nrows, ncols, sparse).matrix([one]*nents)
    16431643
     1644def elementary_matrix(arg0, arg1=None, row1=None, row2=None, col1=None, col2=None, scale=None, sparse=False):
     1645    r"""
     1646    Creates a square matrix that corresponds to a row operation or a column operation.
     1647
     1648    FORMATS:
     1649
     1650    In each case, ``R`` is the base ring, and is optional (the default
     1651    is the integers). ``n`` is the size of the square matrix created.
     1652    We describe the matrices by their associated row operation, see
     1653    the output description for more.
     1654
     1655    - ``elementary_matrix(R, n, row1=i, row2=j)``
     1656
     1657      The matrix which swaps rows ``i`` and ``j``.
     1658
     1659    - ``elementary_matrix(R, n, row1=i, scale=s)``
     1660
     1661      The matrix which multiplies row ``i`` by ``s``.
     1662
     1663    - ``elementary_matrix(R, n, row1=i, row2=j, scale=s)``
     1664
     1665      The matrix which multiplies row ``j`` by ``s``
     1666      and adds it to row ``i``.
     1667
     1668    Elementary matrices representing column operations are created
     1669    in an entirely analogous way, replacing ``row1`` by ``col1`` and
     1670    replacing ``row2`` by ``col2``.
     1671
     1672    Keywords for the rows or columns and the scale must be used,
     1673    since positional arguments are not supported (and will almost surely
     1674    lead to unexpected results).
     1675
     1676    OUTPUT:
     1677
     1678    An elementary matrix is a square matrix that is very close to being
     1679    an identity matrix.  If ``E`` is an elementary matrix and ``A`` is any
     1680    matrix with the same number of rows, then ``E*A`` is the result of
     1681    applying a row operation to ``A``.  This is how the three types
     1682    created by this function are described.  Similarly, an elementary matrix
     1683    can be associated with a column operation, so if ``E`` has the same number
     1684    of columns as ``A`` then ``A*E`` is the result of performing a column
     1685    operation on ``A``.
     1686
     1687    An elementary matrix representing a row operation is created if ``row1``
     1688    is specified, while an elementary matrix representing a column operation
     1689    is created if ``col1`` is specified.
     1690
     1691    EXAMPLES:
     1692
     1693    Over the integers, creating row operations. Recall that row
     1694    and column numbering begins at zero.  ::
     1695
     1696        sage: A = matrix(ZZ, 4, 10, range(40)); A
     1697        [ 0  1  2  3  4  5  6  7  8  9]
     1698        [10 11 12 13 14 15 16 17 18 19]
     1699        [20 21 22 23 24 25 26 27 28 29]
     1700        [30 31 32 33 34 35 36 37 38 39]
     1701
     1702        sage: E = elementary_matrix(4, row1=1, row2=3); E
     1703        [1 0 0 0]
     1704        [0 0 0 1]
     1705        [0 0 1 0]
     1706        [0 1 0 0]
     1707        sage: E*A
     1708        [ 0  1  2  3  4  5  6  7  8  9]
     1709        [30 31 32 33 34 35 36 37 38 39]
     1710        [20 21 22 23 24 25 26 27 28 29]
     1711        [10 11 12 13 14 15 16 17 18 19]
     1712
     1713        sage: E = elementary_matrix(4, row1=2, scale=10); E
     1714        [ 1  0  0  0]
     1715        [ 0  1  0  0]
     1716        [ 0  0 10  0]
     1717        [ 0  0  0  1]
     1718        sage: E*A
     1719        [  0   1   2   3   4   5   6   7   8   9]
     1720        [ 10  11  12  13  14  15  16  17  18  19]
     1721        [200 210 220 230 240 250 260 270 280 290]
     1722        [ 30  31  32  33  34  35  36  37  38  39]
     1723
     1724        sage: E = elementary_matrix(4, row1=2, row2=1, scale=10); E
     1725        [ 1  0  0  0]
     1726        [ 0  1  0  0]
     1727        [ 0 10  1  0]
     1728        [ 0  0  0  1]
     1729        sage: E*A
     1730        [  0   1   2   3   4   5   6   7   8   9]
     1731        [ 10  11  12  13  14  15  16  17  18  19]
     1732        [120 131 142 153 164 175 186 197 208 219]
     1733        [ 30  31  32  33  34  35  36  37  38  39]
     1734
     1735    Over the rationals, now as column operations. Recall that row
     1736    and column numbering begins at zero.  Checks now have the
     1737    elementary matrix on the right.  ::
     1738
     1739        sage: A = matrix(QQ, 5, 4, range(20)); A
     1740        [ 0  1  2  3]
     1741        [ 4  5  6  7]
     1742        [ 8  9 10 11]
     1743        [12 13 14 15]
     1744        [16 17 18 19]
     1745
     1746        sage: E = elementary_matrix(QQ, 4, col1=1, col2=3); E
     1747        [1 0 0 0]
     1748        [0 0 0 1]
     1749        [0 0 1 0]
     1750        [0 1 0 0]
     1751        sage: A*E
     1752        [ 0  3  2  1]
     1753        [ 4  7  6  5]
     1754        [ 8 11 10  9]
     1755        [12 15 14 13]
     1756        [16 19 18 17]
     1757
     1758        sage: E = elementary_matrix(QQ, 4, col1=2, scale=10); E
     1759        [ 1  0  0  0]
     1760        [ 0  1  0  0]
     1761        [ 0  0 10  0]
     1762        [ 0  0  0  1]
     1763        sage: A*E
     1764        [  0   1  20   3]
     1765        [  4   5  60   7]
     1766        [  8   9 100  11]
     1767        [ 12  13 140  15]
     1768        [ 16  17 180  19]
     1769
     1770        sage: E = elementary_matrix(QQ, 4, col1=2, col2=1, scale=10); E
     1771        [ 1  0  0  0]
     1772        [ 0  1 10  0]
     1773        [ 0  0  1  0]
     1774        [ 0  0  0  1]
     1775        sage: A*E
     1776        [  0   1  12   3]
     1777        [  4   5  56   7]
     1778        [  8   9 100  11]
     1779        [ 12  13 144  15]
     1780        [ 16  17 188  19]
     1781
     1782    An elementary matrix is always nonsingular.  Then repeated row
     1783    operations can be represented by products of elementary matrices,
     1784    and this product is again nonsingular.  If row operations are to
     1785    preserve fundamental properties of a matrix (like rank), we do not
     1786    allow scaling a row by zero.  Similarly, the corresponding elementary
     1787    matrix is not constructed.  Also, we do not allow adding a multiple
     1788    of a row to itself, since this could also lead to a new zero row.  ::
     1789
     1790        sage: A = matrix(QQ, 4, 10, range(40)); A
     1791        [ 0  1  2  3  4  5  6  7  8  9]
     1792        [10 11 12 13 14 15 16 17 18 19]
     1793        [20 21 22 23 24 25 26 27 28 29]
     1794        [30 31 32 33 34 35 36 37 38 39]
     1795
     1796        sage: E1 = elementary_matrix(QQ, 4, row1=0, row2=1)
     1797        sage: E2 = elementary_matrix(QQ, 4, row1=3, row2=0, scale=100)
     1798        sage: E = E2*E1
     1799        sage: E.is_singular()
     1800        False
     1801        sage: E*A
     1802        [  10   11   12   13   14   15   16   17   18   19]
     1803        [   0    1    2    3    4    5    6    7    8    9]
     1804        [  20   21   22   23   24   25   26   27   28   29]
     1805        [1030 1131 1232 1333 1434 1535 1636 1737 1838 1939]
     1806
     1807        sage: E3 = elementary_matrix(QQ, 4, row1=3, scale=0)
     1808        Traceback (most recent call last):
     1809        ...
     1810        ValueError: scale parameter of row of elementary matrix must be non-zero
     1811
     1812        sage: E4 = elementary_matrix(QQ, 4, row1=3, row2=3, scale=12)
     1813        Traceback (most recent call last):
     1814        ...
     1815        ValueError: cannot add a multiple of a row to itself
     1816
     1817    Returned matrices have a dense implementation by default,
     1818    but a sparse implementation may be requested.  ::
     1819
     1820        sage: E = elementary_matrix(4, row1=0, row2=1)
     1821        sage: E.is_dense()
     1822        True
     1823
     1824        sage: E = elementary_matrix(4, row1=0, row2=1, sparse=True)
     1825        sage: E.is_sparse()
     1826        True
     1827
     1828    And the ridiculously small cases.  The zero-row matrix fails
     1829    since there are no rows to manipulate. ::
     1830
     1831        sage: elementary_matrix(QQ, 1, 0, 0)
     1832        [1]
     1833        sage: elementary_matrix(QQ, 0, 0, 0)
     1834        Traceback (most recent call last):
     1835        ...
     1836        ValueError: row of elementary matrix must be positive and smaller than 0, not 0
     1837
     1838    TESTS::
     1839
     1840        sage: E = elementary_matrix(QQ, row1=3, scale=12)
     1841        Traceback (most recent call last):
     1842        ...
     1843        ValueError: size of elementary matrix must be given
     1844
     1845        sage: E = elementary_matrix(ZZ, 4/3, row1=3, row2=1, scale=12)
     1846        Traceback (most recent call last):
     1847        ...
     1848        TypeError: size of elementary matrix must be an integer, not 4/3
     1849
     1850        sage: E = elementary_matrix(ZZ, -3, row1=3, row2=1, scale=12)
     1851        Traceback (most recent call last):
     1852        ...
     1853        ValueError: size of elementary matrix must be positive, not -3
     1854
     1855        sage: E = elementary_matrix(ZZ, 5, row2=1, scale=12)
     1856        Traceback (most recent call last):
     1857        ...
     1858        ValueError: row1 or col1 must be specified
     1859
     1860        sage: E = elementary_matrix(ZZ, 5, row1=3, col1=3, scale=12)
     1861        Traceback (most recent call last):
     1862        ...
     1863        ValueError: cannot specify both row1 and col1
     1864
     1865        sage: E = elementary_matrix(ZZ, 5, row1=4/3, row2=1, scale=12)
     1866        Traceback (most recent call last):
     1867        ...
     1868        TypeError: row of elementary matrix must be an integer, not 4/3
     1869
     1870        sage: E = elementary_matrix(ZZ, 5, col1=5, col2=1, scale=12)
     1871        Traceback (most recent call last):
     1872        ...
     1873        ValueError: column of elementary matrix must be positive and smaller than 5, not 5
     1874
     1875        sage: E = elementary_matrix(ZZ, 5, col1=3, col2=4/3, scale=12)
     1876        Traceback (most recent call last):
     1877        ...
     1878        TypeError: column of elementary matrix must be an integer, not 4/3
     1879
     1880        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=-1, scale=12)
     1881        Traceback (most recent call last):
     1882        ...
     1883        ValueError: row of elementary matrix must be positive and smaller than 5, not -1
     1884
     1885        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=1, scale=4/3)
     1886        Traceback (most recent call last):
     1887        ...
     1888        TypeError: scale parameter of elementary matrix must an element of Integer Ring, not 4/3
     1889
     1890        sage: E = elementary_matrix(ZZ, 5, row1=3)
     1891        Traceback (most recent call last):
     1892        ...
     1893        ValueError: insufficient parameters provided to construct elementary matrix
     1894
     1895        sage: E = elementary_matrix(ZZ, 5, row1=3, row2=3, scale=12)
     1896        Traceback (most recent call last):
     1897        ...
     1898        ValueError: cannot add a multiple of a row to itself
     1899
     1900        sage: E = elementary_matrix(ZZ, 5, col1=3, scale=0)
     1901        Traceback (most recent call last):
     1902        ...
     1903        ValueError: scale parameter of column of elementary matrix must be non-zero
     1904
     1905    AUTHOR:
     1906
     1907    - Rob Beezer (2011-03-04)
     1908    """
     1909    # determine ring and matrix size
     1910    if rings.is_Ring(arg0):
     1911        R = arg0
     1912        arg0 = arg1
     1913    else:
     1914        R = rings.ZZ
     1915    if arg0 is None:
     1916        raise ValueError('size of elementary matrix must be given')
     1917    try:
     1918        n = rings.Integer(arg0)
     1919    except:
     1920        raise TypeError('size of elementary matrix must be an integer, not {0}'.format(arg1))
     1921    if n < 0:
     1922        raise ValueError('size of elementary matrix must be positive, not {0}'.format(n))
     1923
     1924    # row operations or column operations?
     1925    # column operation matrix will be transpose of a row operation matrix
     1926    if row1 is None and col1 is None:
     1927        raise ValueError('row1 or col1 must be specified')
     1928    if not row1 is None and not col1 is None:
     1929        raise ValueError('cannot specify both row1 and col1')
     1930    rowop = not row1 is None
     1931    if rowop:
     1932        opstring = "row"
     1933    else:
     1934        opstring = "column"
     1935        row1 = col1
     1936        row2 = col2
     1937
     1938    # analyze parameters
     1939    try:
     1940        row1 = rings.Integer(row1)
     1941    except:
     1942        raise TypeError('{0} of elementary matrix must be an integer, not {1}'.format(opstring, row1))
     1943    if row1 < 0 or row1 >= n :
     1944        raise ValueError('{0} of elementary matrix must be positive and smaller than {1}, not {2}'.format(opstring, n, row1))
     1945    if not row2 is None:
     1946        try:
     1947            row2 = rings.Integer(row2)
     1948        except:
     1949            raise TypeError('{0} of elementary matrix must be an integer, not {1}'.format(opstring, row2))
     1950        if row2 < 0 or row2 >= n :
     1951            raise ValueError('{0} of elementary matrix must be positive and smaller than {1}, not {2}'.format(opstring, n, row2))
     1952    if not scale is None:
     1953        try:
     1954            scale = R(scale)
     1955        except:
     1956            raise TypeError('scale parameter of elementary matrix must an element of {0}, not {1}'.format(R, scale))
     1957
     1958    # determine type of matrix and adjust an identity matrix
     1959    # return row operation matrix or the transpose as a column operation matrix
     1960    elem = identity_matrix(R, n, sparse=sparse)
     1961    if row2 is None and scale is None:
     1962        raise ValueError('insufficient parameters provided to construct elementary matrix')
     1963    elif not row2 is None and not scale is None:
     1964        if row1 == row2:
     1965            raise ValueError('cannot add a multiple of a {0} to itself'.format(opstring))
     1966        elem[row1, row2] = scale
     1967    elif not row2 is None and scale is None:
     1968        elem[row1, row1] = 0
     1969        elem[row2, row2] = 0
     1970        elem[row1, row2] = 1
     1971        elem[row2, row1] = 1
     1972    elif row2 is None and not scale is None:
     1973        if scale == 0:
     1974            raise ValueError('scale parameter of {0} of elementary matrix must be non-zero'.format(opstring))
     1975        elem[row1, row1] = scale
     1976    if rowop:
     1977        return elem
     1978    else:
     1979        return elem.transpose()
     1980
    16441981def _determine_block_matrix_grid(sub_matrices):
    16451982    r"""
    16461983    For internal use. This tries to determine the dimensions