  Our problem is to accurately solve linear systems on a general purpose
graphics processing unit with double double and quad double arithmetic. The
linear systems originate from the application of Newton's method on polynomial
systems. Newton's method is applied as a corrector in a path following method,
so the linear systems are solved in sequence and not simultaneously. One
solution path may require the solution of thousands of linear systems. In
previous work we reported good speedups with our implementation to evaluate and
differentiate polynomial systems on the NVIDIA Tesla C2050. Although the cost
of evaluation and differentiation often dominates the cost of linear system
solving in Newton's method, because of the limited bandwidth of the
communication between CPU and GPU, we cannot afford to send the linear system
to the CPU for solving during path tracking.
  Because of large degrees, the Jacobian matrix may contain extreme values,
requiring extended precision, leading to a significant overhead. This overhead
of multiprecision arithmetic is our main motivation to develop a massively
parallel algorithm. To allow overdetermined linear systems we solve linear
systems in the least squares sense, computing the QR decomposition of the
matrix by the modified Gram-Schmidt algorithm. We describe our implementation
of the modified Gram-Schmidt orthogonalization method for the NVIDIA Tesla
C2050, using double double and quad double arithmetic. Our experimental results
show that the achieved speedups are sufficiently high to compensate for the
overhead of one extra level of precision.
