Modify inverters to return residuals

[jpfoley]

At the moment, as I understand it, the QUDA inverter routines do not return residuals. So currently, if you want to pass the inverter residual back to the client program, you call the inverter routine and then use the host source and solution vectors to work out the residual. This seems wasteful. Why not return the residual calculated on the GPU from the inverter routine?
Also, the residual calculated at the moment can be a poor measure of convergence for heavy quark masses.
Since the heavy-quark propagator drops off quickly, you can have a small relative residual and still have large relative error in a heavy-quark propagator. The Fermilab folk use an alternate error estimate. They compute the per-site residual divided by the norm of the solution vector at that site, and average this ratio over the whole lattice. I don't propose that we use the Fermilab residual as a stopping criterion in the inverters, but it would be nice to return this quantity with the standard relative residual from the inverters.

[maddyscientist]

I imagine we can this fairly easily.

1. To return the residual we could make an optional extra argument which stores a pointer to the residual vector. This would only be computed and returned if non-null.
2. The error measure you propose will require a new kernel. Most convenient to locate in blas I guess. We could make the choice of solver termination an interface option.

On Mar 21, 2012, at 16:49, jpfoleyreply@reply.github.com wrote:

At the moment, as I understand it, the QUDA inverter routines do not return residuals. So currently, if you want to pass the inverter residual back to the client program, you call the inverter routine and then use the host source and solution vectors to work out the residual. This seems wasteful. Why not return the residual calculated on the GPU from the inverter routine?
Also, the residual calculated at the moment can be a poor measure of convergence for heavy quark masses.
Since the heavy-quark propagator drops off quickly, you can have a small relative residual and still have large relative error in a heavy-quark propagator. The Fermilab folk use an alternate error estimate. They compute the per-site residual divided by the norm of the solution vector at that site, and average this ratio over the whole lattice. I don't propose that we use the Fermilab residual as a stopping criterion in the inverters, but it would be nice to return this quantity with the standard relative residual from the inverters.

---

Reply to this email directly or view it on GitHub:
#55

[jpfoley]

I am in the process of changing the inverter code to optionally return the residuals as well as solution vectors. Will push this afternoon.

[jpfoley]

Slight oversight on my part. I had assumed that I could overload functions like invertQuda to return the residual on request, but, since the interface is C, I cannot. Therefore, I would like to replace the interface inverter functions with new functions invertQuda(void h_x, void *h_b, QudaInvertParam *param) --> invertQuda(double final_residual, void *h_x, void *h_b, QudaInvertParam *param), and the residual will only be returned if final_residual is not NULL. Initially, I had thought to copy the residual to param, but I think I prefer the transparency of having an additional parameter. How does this sit with people?

[jpfoley]

Correction to the last post. I meant of course invertQuda(double* final_residual...)

[maddyscientist]

I'm slightly confused. Do you mean to return the L2 norm of the residual, or the residual vector itself? I presume the latter, but I was confused by "the residual".

I'm not in favour of modifying the interface this way, I would think it makes much more sense to return the scalar norm in the struct.

[jpfoley]

Apologies for the confusing post. I was referring to the residual norm. Well, really the rescaled norm that is used to
to test for convergence. So, the norm of the residual divided by the norm of the source, for example.
As for returning the residual norm (ratio, etc, etc) in the param struct, that's fine with me.
Personally, I like to separate input parameters and output parameters in a function call, so one would use
references to const to pass input parameters, and constant pointers for output variables. But that is a question of taste.

[maddyscientist]

Since the interface has to be ripped up anyway, I'm loath to change the interface in such a small way at moment. I guess we should (finally) start to prototype what the new interface will be exactly…..

This email message is for the sole use of the intended recipient(s) and may contain
confidential information. Any unauthorized review, use, disclosure or distribution
is prohibited. If you are not the intended recipient, please contact the sender by

reply email and destroy all copies of the original message.

[jpfoley]

Agreed. There are some other changes, however, that Carleton and the MILC folk want now, which will impact on the interface. For example, the ability to use different stopping criteria in the inverters. It's not the most efficient way to proceed, I know...