We will present a scalable parallel finite difference algorithm for computing the equilibrium configurations, of the order-parameter tensor field for nematic liquid crystals, in rectangular and ellipsoidal regions, by mimization of the Landau - de Gennes free energy functional. In this formulation, we solve for a symmetric traceless 3x3 tensor at each point, thus giving rise to large non-linear algebraic systems. Due to the complexity of the expressions symbolic algebra techniques are used in the production of the code. The target architectures for our implementation are primarily SIMD machines, with 3 dimensional rectangular grid networks, such as the Wavetracer DTC as opposed to hypercube networks such as the Thinking Machines Corporation CM-2. We will describe certain problems encountered in implementing multigrid methods on such architectures and contrast the performance of multigrid with other iterative solvers. We remark that the more computationally intensive Landau - de Gennes formulation permits the modelling of phenomena such as line defects that do not have finite free energy in the more common vector-based Oseen-Zocher-Frank model. It also allows more general behaviors such as non-uniform order and biaxiality to be considered.