  We propose a novel parallel numerical algorithm for calculating the smallest
eigenvalues of highly ill-conditioned matrices. It is based on the {\it LDLT}
decomposition and involves finding a $k \times k$ sub-matrix of the inverse of
the original $N \times N$ Hankel matrix $H_N^{-1}$ . The computation involves
extremely high precision arithmetic, message passing interface, and shared
memory parallelisation. We demonstrate that this approach achieves good
scalability on a high performance computing cluster (HPCC) which constitute a
major improvement of the earlier approaches. We use this method to study a
family of Hankel matrices generated by the weight $w(x)={\rm e}^{-x^{\beta}},$
supported on $[0,\infty)$ and $\beta>0.$ Such weight generates Hankel
determinant, a fundamental object in random matrix theory. In the situation
where $\beta>1/2,$ the smallest eigenvalue tend to 0, exponentially fast as $N$
gets large. If $\beta<1/2,$ the situation where the classical moment problem is
indeterminate, the smallest eigenvalue is bounded from below by a positive
number for all $N$, including infinity. If $\beta=1/2,$ it is conjectured that
the smallest eigenvalue tends to 0 algebraically, with a precise exponent. The
algorithm run on the HPCC producing fantastic match between the theoretical
value of $2/\pi$ and the numerical result.
