def __partial_trace_mat(mat, trace_systems, dimensions, reverse=True):
trace_systems = sorted(trace_systems, reverse=True)<EOL>for j in trace_systems:<EOL><INDENT>dimension_trace = int(dimensions[j])  <EOL>if reverse:<EOL><INDENT>left_dimensions = dimensions[j + <NUM_LIT:1>:]<EOL>right_dimensions = dimensions[:j]<EOL>dimensions = right_dimensions + left_dimensions<EOL><DEDENT>else:<EOL><INDENT>left_dimensions = dimensions[:j]<EOL>right_dimensions = dimensions[j + <NUM_LIT:1>:]<EOL>dimensions = left_dimensions + right_dimensions<EOL><DEDENT>dimension_left = int(np.prod(left_dimensions))<EOL>dimension_right = int(np.prod(right_dimensions))<EOL>mat = mat.reshape([dimension_left, dimension_trace, dimension_right,<EOL>dimension_left, dimension_trace, dimension_right])<EOL>mat = mat.trace(axis1=<NUM_LIT:1>, axis2=<NUM_LIT:4>).reshape(<EOL>dimension_left * dimension_right,<EOL>dimension_left * dimension_right)<EOL><DEDENT>return mat<EOL>
Partial trace over subsystems of multi-partite matrix.

Note that subsystems are ordered as rho012 = rho0(x)rho1(x)rho2.

Args:
    mat (matrix_like): a matrix NxN.
    trace_systems (list(int)): a list of subsystems (starting from 0) to
                              trace over.
    dimensions (list(int)): a list of the dimensions of the subsystems.
                            If this is not set it will assume all
                            subsystems are qubits.
    reverse (bool): ordering of systems in operator.
        If True system-0 is the right most system in tensor product.
        If False system-0 is the left most system in tensor product.

Returns:
    ndarray: A density matrix with the appropriate subsystems traced over.