How does one calculate the position eigenvalues of the matrix corresponding to the position operator? The matrix representation corresponding to the position operator is:
$$x = \sqrt{\frac{\hbar}{2 m \omega}} \left[
\begin{array}{ccccc}
0 & \sqrt{1} & 0 & 0 & \cdots \\
\sqrt{1} & 0 & \sqrt{2} & 0 & \cdots \\
0 & \sqrt{2} & 0 & \sqrt{3} & \cdots \\
0 & 0 & \sqrt{3} & 0 & \cdots \\
\vdots & \vdots & \vdots & \vdots & \ddots
\end{array}\right] \, .$$
I found this matrix is this article by Eric Weisstein (equation (3)).  Now every symmetric, finite dimensional is diagonalizable (see 5.20 in this article). But in this case, the matrix corresponding to $x$ is symmetric but infinite dimensional. So this matrix can't be normalized in the "normal" way by solving the characteristic polynomial.
I know the position eigenfunctions are the delta Dirac functions (actually, distibutions). But how to get the eigenvalues (which have a continuous spectrum) from this matrix. I.e. how to diagonalize this matrix in which (the diagonalized matrix) the diagonal elements resemble the collection of real numbers?
 A: Indirectly. 
To start with, set the pesky and useless constants $\hbar/m\omega \mapsto 1 $, to work in natural units. You then recognize the matrix, as Eric emphasizes by identifying to his (5), as written in the number basis of oscillator discrete energy eigenstates $|n\rangle$. Now for a celebrated basis change. Infinity should not bother you. 
Consider the state 
$$
|q\rangle\equiv \sum_n \psi^*_n(x) |n\rangle ,
$$
where $\psi_n(x)= \langle x|n\rangle $, are (real) Hermite functions, eigenfunctions of the oscillator Schroedinger operator in the x representation-- dot with $\langle m|$: get it?
Dotting with $\langle y|$ and exploiting the completeness of said Hermite functions, 
$$\langle y|q\rangle =\sum_n \psi^*_n(x)  \psi_n(y)=\delta(x-y),
 $$
so you identify this state with $|q\rangle=|x\rangle$.
The pros normally use simpler operators instead of complete sums; that is, they sum the sum to 
$$
|x\rangle =\frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2}|0\rangle ,   
$$
so that , as you should check,
$$
\hat{x} |x\rangle= \frac{(a+a^\dagger)}{\sqrt 2} ~ \frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle= x~\left (  \frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle \right ) =x|x\rangle ~.
$$ 
So the answer to your question is the eigenvalues are all continuous values x, for  eigenvectors $e^{x^2/2}  e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle  /\pi^{1/4} $.
To check your algebra, try the evident unnormalized eigenvector corresponding to x=0, $(1,0,-1/\sqrt{2},0,\sqrt{3/8},0,-\sqrt{5}/4,0,...)$, associated to Hermite numbers.
For a quick crib-sheet on the number basis ops cf Messiah, Quantum Mechanics vI Ch XII §5.
A: *

*That "matrix" form of the position operator, infinite-dimensional or not, is meaningless since neither you nor the page you link say what basis it is supposed to be with respect to.

*Since the $\delta$-distributions which one might take as the eigenvectors of the position operator are not part of the Hilbert space of square-integrable functions the position operator acts on, it cannot be diagonalized on that space. For an operator on an infinite-dimensional separable Hilbert space to be diagonalizable in a meaningful sense, its spectrum has to be discrete, since the usual meaning of a "basis" of such a space is that of a Hilbert basis, which is countable in the separable case.
A: If we approach this problem as you state it, without assuming any context, we have the following: we are given a matrix in an unspecified basis, which Cosmas identified as the energy eigenstates of a harmonic oscillator, and which we are told corresponds to the position operator.
In this (Hilbert) basis, the delta functions are the standard basis elements, which are not position eigenfunctions, since they are not eigenfunctions of the given matrix.
Let's call our basis elements $|0\rangle,|1\rangle,|2\rangle,\ldots,|k\rangle,\ldots$. Their linear span is dense in the state space, and a general element can be written as a limit w.r.t. the obvious inner product that makes this into an orthonormal basis, which is a square integrable sequence $(a_k)_{k\ge0} = (a_0,a_2,a_3,\ldots)$. 
Your question is how to find the eigenvalues if you cannot write down a characteristic polynomial. When the matrix has a nice form like the present one, in which the matrix entries are an expression in number $n$ of the row they're in, we get a recurrence that allows you to find the coefficients one by one. 
In our example, for $(a_k)$ to form an eigenvector with eigenvalue $\lambda$ we must have 
$$\sqrt1a_1 = \lambda a_0$$
and
$$\sqrt{n}a_{n-1} + \sqrt{n+1}a_{n+1} = \lambda a_n.$$
This defines a recurrence that can be simplified by writing $b_k = \sqrt{k!}a_k$, for which we get
$$b_1 = \lambda b_0$$
and
$$nb_{n-1} + b_{n+1} = \lambda b_n.$$
The solution with $b_0 = 0$ gives the $0$ vector, so we can assume that $b_0 = 1$. It is clear that for every $\lambda$ we get a sequence $b_k$ that is a solution ($\lambda = 0$ gives Cosmas' solution, as it should). This gives an infinite series in the original basis, whose limit, if it existed, would be the position eigenstate for the position $\lambda$. Note that going back to the harmonic oscillator, this is an explicit sum in energy eigenstates whose explicit form is known.
What is not so clear is if the associated sequence $(a_k)$ is actually square integrable, i.e. if $\langle a,a\rangle = \sum a_k^2 = \sum \frac{b_k^2}{k!} < \infty$, for which we need to bound the growth of the $b_k$. Even better would be to find a closed form expression for the $b_k$.
We can start by setting $f(x) = \sum_{k = 0}^\infty b_kx^k$ be the generating function of the $b_k$. Then the recurrence relation implies that $f - \lambda xf + x^2(xf)' = f - \lambda xf + x^2f + x^3f'$ is a power series whose coefficient of $x^{n+1}$ is $nb_{n-1} + b_{n+1} - \lambda b_n$ for $n\ge 0$, so this is 0, while the constant term is $b_0 = 1$, so finally we can obtain 
$$(1 - \lambda x + x^2)f(x) + x^3f'(x) = 1.$$
This is general as well. In our case, this doesn't look very easy at all, but sympy actually finds the solution
$$f{\left (x \right )} = \frac{1}{x} \left(C_{1} + \int \frac{e^{\frac{\lambda}{x}}}{x^{2}} e^{- \frac{1}{2 x^{2}}}\, dx\right) e^{- \frac{1}{2 x^{2}} \left(2 \lambda x - 1\right)}$$
For nice $f$, you might find a closed form expression for the $b_n$, hence the $a_n$ from its Taylor coefficients. In this case, if anyone knows how to find a closed form expression for the coefficients from this (or in any other way) or how to bound the coefficients $b_n$, please comment. 
As a final remark, the sum of the series may exist in some formal sense even if it lies outside the state space. If the sum of the squares diverges, that means that it isn't normalizable, so that the standard probabilistic interpretation of the state vector will present some difficulties. As ACuriousMind remarked, the fact that there is an uncountable number of position eigenvalues should imply that none of these sums actually is square integrable.
