$GL_3(\mathbb{F}_2)$ is a simple group

I'm trying to prove that $$G := GL_3(\mathbb{F}_2)$$, the group of $$3 \times 3$$ matrices with entries in $$\mathbb{F}_2$$ is a simple group. The steps outlined for me look like:

1) Construct a list of representatives for the conjugacy classes of $$G$$.

2) Compute the size of each of these conjugacy classes.

3) Show that $$G$$ is simple.

I was able to solve step (1) using the fact that every matrix in $$G$$ is conjugate to a unique block matrix, where each of the blocks are companion matrices of a list of invariant factors. That is, for each matrix $$A \in G$$ there exists unique (up to associates) $$\delta_1 \mid \cdots \mid \delta_n$$, $$\delta_i \in \mathbb{F}_2[x]$$ such that $$A \sim$$diag(Com($$\delta_1$$), $$\ldots$$, Com($$\delta_n$$)). Using this fact I was able to construct the following list of representatives for conjugacy classes of matrices: $$\begin{bmatrix} 0 & 0 & 1\\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix} , \begin{bmatrix} 0 & 0 & 1\\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix}, \begin{bmatrix} 0 & 0 & 1\\ 1 & 0 & 1 \\ 0 & 1 & 1 \end{bmatrix}, \begin{bmatrix} 0 & 0 & 1\\ 1 & 0 & 0 \\ 0 & 1 & 1 \end{bmatrix} , \begin{bmatrix} 1 & 0 & 0\\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix}, \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

I get a little stuck on part (2). I know the size of the conjugacy class of each of the above matrices is equal to the index of the centralizer. I could compute the centralizer of each of the above matrices directly and that would give me the answer but I'm a little hesitant to just multiply matrices for 10-15 minutes. This problem was on a practice qualifying exam so I suspect that there is a faster/more clever way to compute the sizes of these conjugacy classes. This is really what I want. One idea I have:

Two matrices are conjugate if and only if they have the same list of invariant factors. For many of the matrices the list of invariant factors is a single degree 3 polynomial. In this case I know both the minimal polynomial and characteristic polynomial of any matrix conjugate to my representative. These observations do not seem to make the computation much faster though.

I suspect that once I can do (2), (3) will follow relatively quickly.

Calling your six matrices $$a, b, c, d, e, f$$ you can find their centralizers by looking at the action of $$G$$ on the seven non-zero vectors in $$\mathbb{F}_2^3$$:
$$a$$ maps 100 to 010, 010 to 001 and 001 to 100 giving an orbit of length $$3$$. Another orbit of length 3 is $$\{011, 101, 110\}$$ and there is the fixed point $$111$$. As permutation $$a$$ is the product of two $$3$$-cycles. In the symmetric group $$S_7$$ its centralizer looks like $$Z_3^2\rtimes Z_2$$, as its elements have to respect the orbits of $$\langle a\rangle$$ by either fixing them (mapping its elements back into the orbit, no elementwise fixing) or exchanging orbits of the same size. Take a centralizing element $$x$$ that exchanges the orbits, i.e., $$xa=ax$$ and let's say without loss of generality that $$x\cdot100=110$$ (otherwise replace $$x$$ by $$ax$$ or $$a^2x$$). Then $$x\cdot010 = xa\cdot100=ax\cdot100=a\cdot110=011$$, and so by linearity of $$x$$ we get $$x\cdot110=x\cdot100 + x\cdot010=110 + 011=101$$ contradicting that $$x$$ exchanges both orbits of length $$3$$. A linear permutation fixing the orbit $$\{100,010,001\}$$ has to be a power of $$a$$, showing that $$\langle a\rangle$$ is self-centralizing in $$G$$.
Alternatively you can prove that an element $$x$$ of the centralizer of $$a$$ has to map each orbit of order $$3$$ to itself by observing that $$\{011,101,110\}$$ plus the zero vector is a subspace of dimension $$2$$ whose image under a linear $$x$$ has to intersect itself non-trivially. The restriction of $$x$$ to the other orbit (i.e., viewing its image in $$S_3$$) has to be a power of $$a$$, but as this orbit is a basis of the vector space, a linear $$x$$ has to be a power of $$a$$.
$$b$$ and $$d$$ have orbits of length $$7$$, and have as $$7$$-cycles already in $$S_7$$ the centralizers $$\langle b\rangle$$ rsp. $$\langle d\rangle$$.
$$c$$ has orbits $$\{100, 010, 001, 111\}$$, $$\{011, 110\}$$ and fixed point $$101$$. An element $$x\in S_7$$ centralizing $$c$$ has to map the orbit of length $$4$$ to itself, and restricted to this orbit be a power of $$c$$. As the orbit contains a basis of the vector space, a linear $$x$$ is therefor a power of $$c$$.
As you now know the sizes of the conjugacy classes of $$a,b,c,d,f$$ you skip this kind of argument for $$e$$, as its conjugacy classes consists of the remaining elements of $$G$$.