LDA is a supervised classifier, no question about it. The referred post states it clearly in its second part:
How does LDA work?
Linear discriminant analysis (LDA) is another linear transformation
technique that is used for dimensionality reduction. Unlike PCA,
however, LDA is a supervised learning method, which means it takes
class labels into account when finding directions of maximum variance.
It seems to me as if you're mixing some terms, so let's try to put some order in things:
A supervised learning algorithm means that we have data $\{x_i\}$ with proper corresponding tags $\{y_i\}$ and we're trying to find a function $y_i=f(x_i)$ such that it minimizes a loss function $L(f(x_i),y_i)$ in the parametric case (leave the nonparametric methods aside for this discussion).
An unsupervised learning algorithm means that we have data $\{x_i\}$ and no proper tags; Rather than finding a relation, we try to find a latent variable which explains an underlying property of the data, such as grouping to clusters, relations between covariates or simply division into classes.
A dimensionality reduction method aims to take a $p$-dimensional dataset and reduce it to $k$ dimensions, obviously $k<p$, usually aiming to minimize loss of variance. These are not learning algorithms per se.
However, there is some relation between learning and dimensionality reduction: Methods such as Factor Analysis use dimensionality reduction in order to find the underlying linear structure of the covariates (read more about SEM); The wavelet regression lets the user reduce a huge proportion of the data dimensionality while still working under the label of supervised learning (a common example is the MNIST data).
Back to the original question (and the post), assume we have two $p$-dimensional data classes $A,B$ with probabilities mixture $\pi_A,\pi_B$, and assume they have a normal distribution (if not, you can always normalize the data) such that $(x_i|y_i=A,\mu_A,\Sigma_A)\sim\mathcal{N}(\mu_A,\Sigma_A)$ and similarly $(x_i|y_i=B,\mu_B,\Sigma_B)\sim\mathcal{N}(\mu_B,\Sigma_B)$. That's the basis of Gaussian Discriminant Analysis, whose two most popular applications are Quadratic Discriminant Analysis and Linear Discriminant Analysis.
For LDA we further assume $\Sigma_A=\Sigma_B=\Lambda^{-1}$, so the discriminant line formula is
$$f(x)=(\mu_A-\mu_B)^T\Lambda x-\left( \log\frac{\pi_A}{\pi_B}-0.5\mu_A^T\Lambda\mu_A + 0.5\mu_B^T\Lambda\mu_B \right)$$
So in some essence - yes, we kinda reduce the $p$-dimensions into a linear formula - but this is not a dimensionality reduction method, in the same manner that a logistic regression on the same $p$-dimensional data (in which we have $\hat{\theta}_i=x_i^T\hat{\beta}$ and $P(\hat{y}_i=1)=sigmoid(\hat{\theta}_i)$ is not a dimensionality reduction method.
