diff --git a/numpy/lib/tests/test_twodim_base.py b/numpy/lib/tests/test_twodim_base.py
index 9e81cfe4b26a..9548cad30069 100644
--- a/numpy/lib/tests/test_twodim_base.py
+++ b/numpy/lib/tests/test_twodim_base.py
@@ -466,5 +466,30 @@ def test_dtypes(self):
         yield (assert_array_equal, v, expected)
 
 
+class TestElementary(TestCase):
+    def test_switch_rows(self):
+        G = np.elementary(4, 0, 3)
+        assert_equal(G, np.array([[0., 0., 0., 1.],
+                                        [0., 1., 0., 0.], 
+                                        [0., 0., 1., 0.],
+                                        [1., 0., 0., 0.]]))
+
+
+    def test_add_row_multiple(self):
+        H = np.elementary(5, 1, 3, 7, 'i')
+        assert_equal(H, np.array([[1, 0, 0, 0, 0],
+                                        [0, 1, 0, 0, 0],
+                                        [0, 0, 1, 0, 0],
+                                        [0, 7, 0, 1, 0],
+                                        [0, 0, 0, 0, 1]]))
+
+
+    def test_multiply_row(self):
+        I = np.elementary(3, 0, multiplier=6, dtype=int)
+        assert_equal(I, np.array([[6, 0, 0], 
+                                        [0, 1, 0], 
+                                        [0, 0, 1]]))
+
+
 if __name__ == "__main__":
     run_module_suite()
diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py
index 20b5cdd67a7b..3639c5f0caef 100644
--- a/numpy/lib/twodim_base.py
+++ b/numpy/lib/twodim_base.py
@@ -6,7 +6,7 @@
 __all__ = ['diag', 'diagflat', 'eye', 'fliplr', 'flipud', 'rot90', 'tri',
            'triu', 'tril', 'vander', 'histogram2d', 'mask_indices',
            'tril_indices', 'tril_indices_from', 'triu_indices',
-           'triu_indices_from',
+           'triu_indices_from', 'elementary',
            ]
 
 from numpy.core.numeric import (
@@ -998,3 +998,115 @@ def triu_indices_from(arr, k=0):
     if arr.ndim != 2:
         raise ValueError("input array must be 2-d")
     return triu_indices(arr.shape[-2], k=k, m=arr.shape[-1])
+
+
+def elementary(N, i, j=None, multiplier=None, dtype=float):
+    """
+    Return a 2-D array arranged as an elementary matrix.
+
+    Parameters
+    ----------
+    N : int
+        The size of the NxN array to be returned. Elementary matrices are 
+        square by definition.
+    i : int
+        The index of the first row on which operations are to be performed.
+    j : int
+        If set, the index of the second row on which operations are to be 
+        performed.
+    multiplier : scalar
+        If set, the factor by which a given row will be multiplied.
+
+    Returns
+    -------
+    m : ndarray of shape (NxN)
+        An identity array that has had a single row operation performed on it.
+
+    See also
+    --------
+    eye, identity
+
+    Examples
+    --------
+    To swap the the first and third rows of a 4x4 identity array:
+
+    >>> L = elementary(4, 0, 2)
+    >>> L
+    array([[ 0.,  0.,  1.,  0.],
+           [ 0.,  1.,  0.,  0.],
+           [ 1.,  0.,  0.,  0.],
+           [ 0.,  0.,  0.,  1.]])
+
+    When the elementary array is matrix-multiplied by a given array, the 
+    result is the given array with its first and third rows swapped.
+
+    >>> H = np.array([[ 2,  3,  5,  7],
+                      [11, 13, 17, 19],
+                      [23, 29, 31, 37],
+                      [41, 43, 47, 53]])
+    >>> L.dot(H)
+    array([[ 23.,  29.,  31.,  37.],
+           [ 11.,  13.,  17.,  19.],
+           [  2.,   3.,   5.,   7.],
+           [ 41.,  43.,  47.,  53.]])
+
+    If the array is matrix-multiplied by the elementary array (i.e., the
+    multiplication takes place in reverse order), the result is the given 
+    array with its first and third columns swapped.
+
+    >>> H.dot(L)
+    array([[  5.,   3.,   2.,   7.],
+           [ 17.,  13.,  11.,  19.],
+           [ 31.,  29.,  23.,  37.],
+           [ 47.,  43.,  41.,  53.]])
+
+    To add a multiple of an array's second row to the same array's fourth row:
+
+    >>> H
+    array([[ 2,  3,  5,  7],
+           [11, 13, 17, 19],
+           [23, 29, 31, 37],
+           [41, 43, 47, 53]])
+    >> M = elementary(4, 1, 3, 10000, int)
+    >> M
+    array([[    1,     0,     0,     0],
+           [    0,     1,     0,     0],
+           [    0,     0,     1,     0],
+           [    0, 10000,     0,     1]])
+    >> M.dot(H)
+    array([[     2,      3,      5,      7],
+           [    11,     13,     17,     19],
+           [    23,     29,     31,     37],
+           [110041, 130043, 170047, 190053]])
+
+    To multiply only the last column of an array by a scalar:
+
+    >>> H
+    array([[ 2,  3,  5,  7],
+           [11, 13, 17, 19],
+           [23, 29, 31, 37],
+           [41, 43, 47, 53]])
+    >>> G = elementary(len(H), len(H)-1, multiplier=2)
+    >>> H.dot(G)
+    array([[   2.,    3.,    5.,   14.],
+           [  11.,   13.,   17.,   38.],
+           [  23.,   29.,   31.,   74.],
+           [  41.,   43.,   47.,  106.]])
+
+    Any nonsingular array can be formed by taking the product of an array of 
+    elementary arrays. 
+
+    """
+    m = eye(N, dtype=dtype)
+    if j is None and multiplier is None:
+        raise ValueError("'j' and 'multiplier' are both set to None.")
+    elif multiplier is None:
+        m[[i, j]] = m[[j, i]]
+        return m
+    elif j is None:
+        m[i] *= multiplier
+        return m
+    else:
+        m[j] += (multiplier * m[i])
+        return m
+