  This paper describes the application of the code generated by the CAMPARY
software to accelerate the solving of linear systems in the least squares sense
on Graphics Processing Units (GPUs), in double double, quad double, and octo
double precision. The goal is to use accelerators to offset the cost overhead
caused by multiple double precision arithmetic. For the blocked Householder QR
and the back substitution, of interest are those dimensions at which teraflop
performance is attained. The other interesting question is the cost overhead
factor that appears each time the precision is doubled.
  Experimental results are reported on five different NVIDIA GPUs, with a
particular focus on the P100 and the V100, both capable of teraflop
performance. Thanks to the high Compute to Global Memory Access (CGMA) ratios
of multiple double arithmetic, teraflop performance is already attained running
the double double QR on 1,024-by-1,024 matrices, both on the P100 and the V100.
For the back substitution, the dimension of the upper triangular system must be
as high as 17,920 to reach one teraflops on the V100, in quad double precision,
and then taking only the times spent by the kernels into account. The lower
performance of the back substitution in small dimensions does not prevent
teraflop performance of the solver at dimension 1,024, as the time for the QR
decomposition dominates.
  In doubling the precision from double double to quad double and from quad
double to octo double, the observed cost overhead factors are lower than the
factors predicted by the arithmetical operation counts. This observation
correlates with the increased performance for increased precision, which can
again be explained by the high CGMA ratios.
