WIP: TFQMR classes and reference kernels

[hartwiganzt]

this is the initial tfqmr pull request. Core functionality + reference kernels

[pratikvn]

Maybe this if clause can go into an helper function after combining the two steps 1 and 5 so that you just have to pass in your different inputs to the helper function.

[gflegar]

I didn't go through the whole thing, but I did find some serious problems with it - it seem that something is not right here. See inline comments. (That's why I haven't look at the rest in more detail.)

Also, I think we have a general problem with reviewing Krylov solvers, see #29.

Spent the whole day on this, but at least something useful came out of it: found a bug here, found several possible improvements, started a discussion on how we review Krylov solvers, reminded myself of the details of CG, and created a draft for better CG documentation.

[gflegar]

step_2 description says sigma = (theta^2 / alpha) * eta

[gflegar]

description says c_mp1 = 1 / (1 + theta)

[gflegar]

Trying to mask two clearly different iterations into one has serious drawbacks:

1. It is virtually impossible to understand (I just spent several hours trying to figure out what's going on here!). Krylov methods already have extremely complicated loop invariants, adding control statements to them just to decrease code size makes reasoning about them difficult. You need to mentally unroll the loop anyway, so why not just write it like that in the first place?
2. You can merge more kernels than you have currently merged.
3. You have quite a few additional operations that are completely unnecessary.

After figuring out what this does (all without understanding the Krylov method, just by looking at the sequence of operations), I found that there must be at least 1 bug in the algorithm here, that you actually need only 5 merged kernels for 2 iterations (instead of 7 for 1) and that you can avoid copying vectors 3 times and scalars 1 time over the course of 2 iterations. This is all by taking only a superficial look, and it might be possible to get even more savings. Of course, it didn't make sense to look into it in more detail, since there's a bug in the algorithm anyway.

See the following file for details: tfqmr.txt

[hartwiganzt]

I tried to keep the code as close as possible to the MATLAB code. The idea was that this could be used by users to match the functionality. I you outlined in #29 there may be another, more efficient way to go.

You can merge more kernels than you have currently merged.

Yes, I think we had this discussion on the phone. I would like to see timings to decide whether it is a good idea. But sure, I can do that.

You have quite a few additional operations that are completely unnecessary.

Not sure about that. For example, I think

[u_mp1] = u_m - [alpha] * v ! ERROR: the result is never used

is needed in the update step

pu_m = M^-1 [u_mp1] <- apply

for the even iteration counts.

you actually need only 5 merged kernels for 2 iterations (instead of 7 for 1)

Yes, we have more merging in MAGMA-sparse. And I wanted to see the performance difference because less merging is easier to understand.

[gflegar]

I tried to keep the code as close as possible to the MATLAB code. The idea was that this could be used by users to match the functionality.

Functionality isn't dictated by implementation details, we shouldn't blindly stick to something if it doesn't work well.

Yes, I think we had this discussion on the phone. I would like to see timings to decide whether it is a good idea. But sure, I can do that.

I am aware of our discussion:

1. Even if you keep scalar updates separate from vector updates, you can still merge more.
2. Either way, you should merge Ginkgo kernels, and have some of them call two CUDA kernels. This performance issue is an implementation detail, which can differ for different devices, and shouldn't leak to the top-level solver implementation.

Not sure about that. For example, I think

[u_mp1] = u_m - [alpha] * v ! ERROR: the result is never used

is needed in the update step

pu_m = M^-1 [u_mp1] <- apply

for the even iteration counts.

u_mp1 gets updated before it's used here: [u_mp1] = w + [beta] * u_m. Just search for u_mp1 in any of the codes (the original file, my transcript of the code, or the unrolled one), and you'll see that there are two consecutive assigns to u_mp1.
This is something that's hard to notice in the original version of the code due to the convoluted implementation. If that's based on MATLAB's code, then their implementation is quite horrible.
Please take some more time to look through my comments (at least the error, and opt ones), and point out exactly what's wrong if something is wrong with them.

less merging is easier to understand.

I disagree, we have comments below, so you don't look at the merged call anyway. Where's the difference for the user if he sees calls to step_1 ..., step_5 as oposed to step_1, ...., step_7?
Also this code has far greater issues when it comes to "understanding".

[gflegar]

Closing this due to inactivity and interface changes done to the main Ginkgo branch that make this implementation obsolete (e.g. updated EnableLinOp mixin implementation, stopping criteria, loggers).

We should probably keep the branch for reference, until a proper TFQMR implementation is made.