def get_P_matrix(bravais_lattice):
    """
    Return a tuple of length 2 with the P matrix and its inverse::

      (P, invP)

    with :math:`invP = P^{-1}`.

    These :math:`P` matrices are obtained from Table 3 of the HPKOT 
    paper.

    The P matrix is a :math:`3\\times 3` matrix is the matrix that converts
    the lattice vectors from crystallographic conventional
    :math:`(a,b,c)` to crystallographic primitive :math:`(a_P, b_P, c_P)`
    as follows: :math:`(a_P, b_P, c_P) = (a,b,c) P`

    The change of (real space) coordinate triples follows instead: 
    :math:`(x_P, y_P, z_P)^T = (P^{-1}) (x,y,z)^T`

    .. note:: the :math:`invP = P^{-1}` matrix is always integer (with values
        only :math:`-1, 0, 1`) while :math:`P` is rational (non-integer values can be
        :math:`\\pm \\frac 1 2` and :math:`\\pm \\frac 1 3`).
    """
    import numpy as np
    if bravais_lattice in ['cP', 'tP', 'hP', 'oP', 'mP']:
        P = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
        invP = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
    elif bravais_lattice in ['cF', 'oF']:
        P = 1.0 / 2.0 * np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
        invP = np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]])
    elif bravais_lattice in ['cI', 'tI', 'oI']:
        P = 1.0 / 2.0 * np.array([[-1, 1, 1], [1, -1, 1], [1, 1, -1]])
        invP = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
    elif bravais_lattice == 'hR':
        P = 1.0 / 3.0 * np.array([[2, -1, -1], [1, 1, -2], [1, 1, 1]])
        invP = np.array([[1, 0, 1], [-1, 1, 1], [0, -1, 1]])
    elif bravais_lattice == 'oC':
        P = 1.0 / 2.0 * np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 2]])
        invP = np.array([[1, -1, 0], [1, 1, 0], [0, 0, 1]])
    elif bravais_lattice == 'oA':
        P = 1.0 / 2.0 * np.array([[0, 0, 2], [1, 1, 0], [-1, 1, 0]])
        invP = np.array([[0, 1, -1], [0, 1, 1], [1, 0, 0]])
    elif bravais_lattice == 'mC':
        P = 1.0 / 2.0 * np.array([[1, -1, 0], [1, 1, 0], [0, 0, 2]])
        invP = np.array([[1, 1, 0], [-1, 1, 0], [0, 0, 1]])
    elif bravais_lattice == 'aP':
        P = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
        invP = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
    else:
        raise ValueError('Invalid bravais_lattice {}'.format(bravais_lattice))
    return P, invP
