Mixed precision for xGELS

[jgpallero]

Hello,

for xGESV and xPOSV exist (DS/ZC)GESV and (DS/ZC)DSPOSV, but I've not found the mixed precision routines for xGELS. Exist any technical or numerical reason because mixed precision for xGELS functions are not implemented?

Thanks

[Julien Langou]

The algorithm for iterative refinement and least square is a little more tricky, and we never got enough request to quite get to write the code. See: Ake Bjorck and Gene Golub. "Iterative refinement of linear least squares solution by Householder transformation (Algol Programming)". BIT, 7:322–337, 1967. A good recent study with a good bibliography and survey of previous algorithms was done by Demmel, Hida, Riedy, Li ("Extra-Precise Iterative Refinement for Overdetermined Least Squares Problems", ACM TOMS, Volume 35 Issue 4, February 2009, Article No. 28). This is in the context of Extra-Precise Iterative Refinement, so, for example, you work with the low precision and use iterative refinement with the higher precision to get good forward error in with the lower precision. This is not what the "DSGESV" are about. In DSGESV, you work in the high precision and factor in low precision to get speed and then you try to recover the accuracy of the high precision with iterative refinement. So one is about forward error, the other is about speed. This changes stopping criteria. Cheers, Julien.

[jgpallero]

Thank you for your answer.

So, I understand that about the xGELS algorithm the Bjork's method is about speed and the Dmmels et al. one is about precision. Am I right? For example, in the PLASMA software there is a function called DSUNGESV that uses the Bjork's results to perform a single precision QR decomposition and then make an iterative refinement. In MAGMA there is the DSGEQRSV that apparently does also the GELS in mixed (single/double) precision.

[Julien Langou]

I understand that about the xGELS algorithm the Bjork's method is about speed and the Dmmels et al. one is about precision. Am I right?

Well, yep kind of. This is kind of the same method everywhere, the difference is subtle and is simply in the stopping criteria. Assume you do single-double iterative refinement. When do you stop? Do you stop at the lower accuracy (single precision)? In which case I would say you are shooting for accuracy. That is your natural working precision is the lower precision (single) and you are using the higher precision (double) to have a good accuracy on the lower precision (single). Or do you stop at the higher precision (double)? In which case, I would say you are shooting for speed. That is your natural working precision is the higher precision (double) and you are using the lower precision (single) to factor the matrix faster (than the higher).

For example, in the PLASMA software there is a function called DSUNGESV that uses the Bjork's results to perform a single precision QR decomposition and then make an iterative refinement. In MAGMA there is the DSGEQRSV that apparently does also the GELS in mixed (single/double) precision.

Yup. And all these are shooting for speed.