[proofplan]
We use the matrix-entry coordinates on $M(n,\mathbb F)$, viewing the complex case as real coordinates given by real and imaginary parts. Since $GL(n,\mathbb F)$ is an open submanifold of $M(n,\mathbb F)$ by [citetheorem:8772], it is enough to prove that the coordinate components of multiplication and inversion are smooth functions on these open coordinate domains. Multiplication is polynomial in the entries, while inversion is given by the adjugate formula divided by the nonvanishing determinant.
[/proofplan]
[step:Use matrix-entry coordinates on the open general linear group]
For $\mathbb F=\mathbb R$, identify $M(n,\mathbb R)$ with $\mathbb R^{n^2}$ by the coordinate functions
\begin{align*}
x_{ij}:M(n,\mathbb R)\to\mathbb R,\qquad A\mapsto A_{ij},
\end{align*}
where $1\le i,j\le n$.
For $\mathbb F=\mathbb C$, identify $M(n,\mathbb C)$ with $\mathbb R^{2n^2}$ by writing
\begin{align*}
A_{ij}=x_{ij}(A)+i\,y_{ij}(A),
\end{align*}
where
\begin{align*}
x_{ij}:M(n,\mathbb C)\to\mathbb R,\qquad A\mapsto \operatorname{Re}(A_{ij})
\end{align*}
and
\begin{align*}
y_{ij}:M(n,\mathbb C)\to\mathbb R,\qquad A\mapsto \operatorname{Im}(A_{ij}).
\end{align*}
By [citetheorem:8772], $GL(n,\mathbb F)$ is open in $M(n,\mathbb F)$, so these ambient matrix-entry coordinates restrict to smooth coordinates on $GL(n,\mathbb F)$. Therefore, by the coordinate criterion for smoothness of maps between open submanifolds of finite-dimensional real vector spaces, it suffices to prove smoothness of the coordinate component functions of $m$ and $i$.
[/step]
[step:Write multiplication as polynomial functions of the entries]
Let $A,B\in GL(n,\mathbb F)$. For each pair $1\le i,j\le n$, the $(i,j)$-entry of $m(A,B)=AB$ is
\begin{align*}
(AB)_{ij}=\sum_{k=1}^{n} A_{ik}B_{kj}.
\end{align*}
If $\mathbb F=\mathbb R$, this is a polynomial function of the real coordinate functions $A_{ik}$ and $B_{kj}$.
If $\mathbb F=\mathbb C$, write $A_{ik}=a_{ik}+i\alpha_{ik}$ and $B_{kj}=b_{kj}+i\beta_{kj}$, where $a_{ik},\alpha_{ik},b_{kj},\beta_{kj}\in\mathbb R$. Then
\begin{align*}
\operatorname{Re}(AB)_{ij}=\sum_{k=1}^{n}(a_{ik}b_{kj}-\alpha_{ik}\beta_{kj})
\end{align*}
and
\begin{align*}
\operatorname{Im}(AB)_{ij}=\sum_{k=1}^{n}(a_{ik}\beta_{kj}+\alpha_{ik}b_{kj}).
\end{align*}
Thus every real coordinate component of $m$ is a polynomial function on the real coordinate space. Polynomial functions on finite-dimensional real vector spaces are smooth, hence $m$ is smooth.
[/step]
[step:Express inversion by adjugates and nonvanishing determinants]
For $A\in GL(n,\mathbb F)$, define the adjugate matrix $\operatorname{adj}(A)\in M(n,\mathbb F)$ by
\begin{align*}
\operatorname{adj}(A)_{ij}=C_{ji}(A),
\end{align*}
where $C_{ji}(A)$ denotes the $(j,i)$-cofactor of $A$. Each cofactor is a signed determinant of an $(n-1)\times(n-1)$ submatrix of $A$, so each function
\begin{align*}
A\mapsto \operatorname{adj}(A)_{ij}
\end{align*}
is polynomial in the entries of $A$. The determinant function
\begin{align*}
\det:M(n,\mathbb F)\to\mathbb F,\qquad A\mapsto \det A
\end{align*}
is also polynomial in the entries of $A$. Since $A\in GL(n,\mathbb F)$, $\det A\ne 0$.
By the adjugate inverse formula,
\begin{align*}
A^{-1}=\frac{1}{\det A}\operatorname{adj}(A).
\end{align*}
For $\mathbb F=\mathbb R$, each coordinate function of $A^{-1}$ is therefore a quotient of two smooth real-valued functions with nonzero denominator on $GL(n,\mathbb R)$, hence is smooth.
For $\mathbb F=\mathbb C$, write
\begin{align*}
\operatorname{adj}(A)_{ij}=p_{ij}(A)+i q_{ij}(A)
\end{align*}
and
\begin{align*}
\det A=u(A)+i v(A),
\end{align*}
where $p_{ij},q_{ij},u,v:GL(n,\mathbb C)\to\mathbb R$ are real polynomial functions in the real and imaginary matrix-entry coordinates. Since $\det A\ne 0$, one has $u(A)^2+v(A)^2>0$. Hence
\begin{align*}
\operatorname{Re}(A^{-1})_{ij}=\frac{p_{ij}(A)u(A)+q_{ij}(A)v(A)}{u(A)^2+v(A)^2}
\end{align*}
and
\begin{align*}
\operatorname{Im}(A^{-1})_{ij}=\frac{q_{ij}(A)u(A)-p_{ij}(A)v(A)}{u(A)^2+v(A)^2}.
\end{align*}
Both are quotients of smooth real-valued functions with positive denominator, so every real coordinate component of $i$ is smooth.
[guided]
The point of this step is to turn inversion into ordinary coordinate functions. Let
\begin{align*}
\operatorname{adj}:M(n,\mathbb F)\to M(n,\mathbb F)
\end{align*}
denote the adjugate map, whose $(i,j)$-entry is the $(j,i)$-cofactor of the matrix. A cofactor is a signed determinant of an $(n-1)\times(n-1)$ submatrix, and determinants are finite sums of products of entries. Therefore each coordinate function
\begin{align*}
A\mapsto \operatorname{adj}(A)_{ij}
\end{align*}
is polynomial in the entries of $A$. The determinant function
\begin{align*}
\det:M(n,\mathbb F)\to\mathbb F,\qquad A\mapsto \det A
\end{align*}
is polynomial for the same reason.
Now restrict to $GL(n,\mathbb F)$. The defining property of $GL(n,\mathbb F)$ is invertibility, equivalently $\det A\ne 0$. Thus the adjugate formula applies:
\begin{align*}
A^{-1}=\frac{1}{\det A}\operatorname{adj}(A).
\end{align*}
When $\mathbb F=\mathbb R$, the $(i,j)$-entry of $A^{-1}$ is a polynomial function divided by the nonzero smooth function $\det A$. Quotients by nonvanishing smooth functions are smooth, so every coordinate function of inversion is smooth.
For $\mathbb F=\mathbb C$, smoothness means real smoothness after identifying $M(n,\mathbb C)$ with $\mathbb R^{2n^2}$. Write the relevant complex numerator and denominator as
\begin{align*}
\operatorname{adj}(A)_{ij}=p_{ij}(A)+i q_{ij}(A)
\end{align*}
and
\begin{align*}
\det A=u(A)+i v(A),
\end{align*}
where $p_{ij},q_{ij},u,v$ are real polynomial functions of the real and imaginary parts of the entries of $A$. Because $A\in GL(n,\mathbb C)$, the complex number $u(A)+i v(A)$ is nonzero, so
\begin{align*}
u(A)^2+v(A)^2>0.
\end{align*}
Dividing complex numbers by multiplying numerator and denominator by $u(A)-i v(A)$ gives
\begin{align*}
\frac{p_{ij}(A)+iq_{ij}(A)}{u(A)+iv(A)}
=
\frac{(p_{ij}(A)+iq_{ij}(A))(u(A)-iv(A))}{u(A)^2+v(A)^2}.
\end{align*}
Therefore the real and imaginary parts are
\begin{align*}
\operatorname{Re}(A^{-1})_{ij}=\frac{p_{ij}(A)u(A)+q_{ij}(A)v(A)}{u(A)^2+v(A)^2}
\end{align*}
and
\begin{align*}
\operatorname{Im}(A^{-1})_{ij}=\frac{q_{ij}(A)u(A)-p_{ij}(A)v(A)}{u(A)^2+v(A)^2}.
\end{align*}
Each numerator and denominator is a smooth real-valued function in the real coordinates, and the denominator is positive on $GL(n,\mathbb C)$. Hence each real coordinate component of inversion is smooth.
[/guided]
[/step]
[step:Conclude smoothness of both group operations]
The multiplication map $m$ has smooth coordinate components in the matrix-entry coordinates, and the inversion map $i$ also has smooth coordinate components in those same coordinates. Since $GL(n,\mathbb F)$ is an open smooth submanifold of $M(n,\mathbb F)$ and the complex case is interpreted through the real coordinate identification $M(n,\mathbb C)\cong\mathbb R^{2n^2}$, the coordinate criterion for smoothness gives that
\begin{align*}
m:GL(n,\mathbb F)\times GL(n,\mathbb F)\to GL(n,\mathbb F)
\end{align*}
and
\begin{align*}
i:GL(n,\mathbb F)\to GL(n,\mathbb F)
\end{align*}
are smooth maps of real smooth manifolds.
[/step]
