  We extend a systematic renormalization procedure for quantum field theory to
include particle masses and present several applications. We use a Hamiltonian
formulation and light-front quantization because this may produce a convergent
Fock-space expansion. The QCD Hamiltonian is systematically renormalized to
second order in the strong coupling and the Fock-space expansion is truncated
to lowest order to produce a finite-dimensional Hamiltonian matrix. The
renormalized Hamiltonian is used to calculate the spectra of the $b\bar{b}$ and
$c\bar{c}$ mesons as a lowest-order test of our procedure for full QCD.
  The analytic determination of the renormalized Hamiltonian matrix generates
expressions that must be numerically integrated to generate quantitative
results. The efficiency of the numerical calculation depends on how well the
basis functions can approximate the real state. We examine the effectiveness of
using Basis-Splines (B-Splines) to represent QCD states. After briefly
describing these functions, we test them using the one- and two-dimensional
harmonic oscillator problems. We test their ability to represent realistic
wavefunctions by using them to find the glueball mass spectrum.
  An efficient algorithm for numerically calculating the matrix elements in the
glueball and meson problems is necessary because the calculation is numerically
intensive. We describe our algorithm and discuss its parallel-cpu
implementation.
