Improve DynamicSparseNumberFoo performance

[roystgnr]

The good news: we now have a working class that gives O(N_nonzero) instead of O(N) performance. This could be a big deal for https://github.com/libantioch/antioch sensitivity applications where our eventual N/N_nonzero will be ~1000.

The bad news: @friedmud reports that the current implementation has a constant performance overhead of ~1000. (based on profiling idaholab/moose#5661, and according to the profiler the cost is mostly from new/delete). He fixed the obvious flaws (pass-by-value to the sparsity operations) to little effect.

So let's brainstorm ideas. I'll sort what I've got so far by increasing level of difficulty IMHO:

PBR: Pass-by-reference in the sparsity operations (even if it's not a big win it's still better than nothing)

RI: Use reverse iteration to do in-place sparsity operations rather than creating temporary merged_foo vectors.

RV: Add C++11 rvalue operations so we can steal doomed input arguments' allocations where available.

CA: Give a custom allocator to the underlying std::vector.

CC: Replace std::vector with a custom container that keeps the first O(N_nonzero) elements on the stack and only hits the heap for larger cases.

ET: Use expression templates to postpone evaluations until we can do more of them at once without creating so many intermediate temporaries.

[friedmud]

1c31b25 cuts residual computation time in half by simply trying to NOT compute derivatives during residual computation... but I'm sure it's not quite architected right (I'm still not great at template metaprogramming).

[friedmud]

If you completely turn off derivative calculation (by replacing the two if() statements in 1c31b25 with if(false)) the residual computation will only be double (7 vs 3.5) when using NumberArray.

It does also help when using DynamicSparseNumberArray... but unfortunately it just makes the time 49 seconds instead of 186!

[friedmud]

BTW: When using NumberArray the amount of storage set aside for derivatives has a HUGE impact on the runtime.

10: 5 seconds
100: 6 seconds
1000: 32 seconds
10000: Too long for me to wait

[friedmud]

That timing was with derivative computation off. Here it is with it on:

10: 5.5 seconds
100: 26
1000: 286

[friedmud]

Just trying to get to know the datastructures here... I implemented something called FriedmudArray that's backed by a std::map... that was an incredible failure (as you might imagine!). It came it at around 300 seconds.

(You can see this abomination here: c2d7dc8 )

Once again though: the problem is allocation / deallocation of the map. It's definitely just that we're making too many temporaries... but I don't know exactly what to do about it.

Like @roystgnr outlined... it's probably going to take either ET or RV. Neither of which is doable by me.

[roystgnr]

PBR and RI are implemented as of 1254b55 - the combination of which appears to do exactly jack for performance. Going to take a crack at RV next, but only after some perf runs on meta to see what I'm misunderstanding.

[friedmud]

What do you think about my attempts at removing derivative computation ( 1c31b25 )? I have a few more local commits that do even more... but to actually make it work well I think we need to change the architecture a bit.

I was thinking that we may need to store _derivatives as a pointer... and just not even initialize it at all unless the DualNumber is supposed to be doing derivative calculations. This would mean that we need to guard _derivatives a bit more (don't return a raw reference to it ever)... adding a bit more interface to DualNumber to allow you to get/set derivative entries (would would do nothing if that DualNumber is not supposed to be doing anything with derivatives).

This would keep us from ever creating (or destroying) the derivatives storage in the case where we're computing the residual....

[roystgnr]

Honestly? IMHO the way to avoid doing derivative computations with DualNumber is "don't pass in a DualNumber". If someone's code is sufficiently data type independent that you set T=DualNumber in the first place, it's usually generic enough that you can do one instantiation with T=DualNumber and another with T=double.

Otherwise we're tossing "if (_derivatives)" tests into every operation, which can't be good for performance. (although I admit it'd probably be slightly better than the implicit "if (!derivatives.empty())" equivalent that you'd get from dragging empty DynamicSparseNumberArray objects around) I'd much rather create a new "PotentiallyDualNumber" class than make such frequent branch instructions unavoidable.

[roystgnr]

We're still spending 52% of runtime in new/delete. More informatively, we're seeing 6 new/delete calls per computeQpResidual() call in Meta.

The dot product x*x+y*y+z*z would naively force construction of temporaries for xx, yy, zz, xxpyy, and xxpyypzz, then maybe a copy off the stack if the RVO can't be applied (although I'm not sure why it can't be, when it's just being used to initialize a temporary before temporary.value() is called...). This I guess accounts for all 6 calls.

With my move operations, on the other hand, the xxpyy temporary should swipe xx's memory, then be swiped in turn by xxpyypzz. If I understand C++11 correctly the return should be done via move constructor even if RVO doesn't apply...

But what I'm getting at here is that:

1. I clearly don't understand C++11 correctly or I'd be using only half as many copy constructions as valgrind is seeing. Fortunately this should be easy to hunt down.
2. Even if move constructiors were working they'd only be fixing half the problem here (thanks to the surfeit of non-rvalue inputs in this basic operation) and we really need to be using a custom allocator to fix the other half; that could make the difference between another 33% speedup and another 100%.

[friedmud]

Yes, yes.... all true computer scientists would say that the thing to do is template this up the wazoo to avoid the derivative computation....

I've gone over this in my head and I can't find a way to do that without completely infecting MOOSE and destroying the simple C++ we present to our users. This is the main reason why MOOSE hasn't had an AD capability from the beginning. Every time I looked at it, and what it did to the code, I ran screaming for the hills.

The biggest issue with templating the world is that it would necessitate two copies of everything. Two copies of every object (one instance with T=double one with T=DualNumber). That includes every material property... including stateful material properties (which are stored at every quadrature point). Objects that hold a lot of data (like Neutronics cross-sections) would have that data multiplied by two as well.

Yuck... just yuck. I would be ok with two copies of objects if it didn't mean two copies of data.

I mean, think about the Mesh. Like I mentioned to you, I would like to actually store the mesh coordinates as DualNumbers so that we could compute derivatives with respect to displacements properly... but if you go this templating route you would have to have two copies of your mesh! If you could just turn the derivatives "off" when you don't need them, that goes away.

Some of it is due to the way MOOSE is designed with inherited interfaces that serve up data references. It's not as simple as just templating computeQpResidual()... because all of the data feeding your residual calculation would also need to be templated... so the whole class would have to be templated and every class all the way down the inheritance tree would have to be templated as well. Yuck.

Interestingly this is actually one of those cases where a non-compiled code (like Python) can actually be more optimal. You could actually change what DualNumber does on the fly without needing to have an if() statement (you could just replace the routine with one that doesn't do anything with derivatives... at runtime). In a sense, this is similar to what the PETSc guys do... they have "dynamic objects" where they can replace function pointers on the fly to do things differently. It's one of the reasons why they don't like C++ and stick with C... because C++ is more rigid in this regard.

Doing this at runtime with C++ in MOOSE... I don't see any other option other than using an if(). But maybe I haven't thought about it long enough.

How difficult would a PotentiallyDualNumber be to implement?

I understand that it wouldn't be "optimal" in the computer science sense... but if we can get the runtime hit down to 2x or less over the "normal" MOOSE code then it will be viable.

On other topics... would it help if we use R-value reference returns for computeQpResidual()? That would only avoid one extra copy per call - but that is still a lot. We could definitely do that (and require C++11 for this new capability).

[friedmud]

Just to throw some more blasphemy in here... I'm actually thinking that a global (!) variable like MetaPhysicL::compute_derivatives would be useful. I could set it to true during Jacobian computation and false for residual computation.

Trying to set it per-variable is actually REALLY tricky.

I know that sounds insane... but I've been over this quite a bit in my head and it would actually be cleaner that way.

[roystgnr]

Updates:

I do understand C++11 moves correctly after all. I'm just seeing 6 allocations instead of 3 because I was forgetting that each of the three remaining DNSA copies does a separate allocation for _data and _indices. Trudging through in the debugger, the moves are all working fine.

There's no simple way to turn that 6 allocations into 3 via STL without using a vector-of-pairs, which would kill any chance at SIMD operations. I can't swear that'd be a bad idea, though, not when our memory allocation is currently twice as expensive as our floating point operations.

Using pointers and a custom static allocator instead of STL might be a better idea anyways. The extra complexity of having to replace "= default" methods with custom code doesn't look any worse than the extra complexity of the STL allocator interface.

The C++ way to change function pointers on the fly is with inheritance. A vtable itself isn't so bad, but the fact that you're using function pointers at all (whether in C or C++) prevents inlining, which for small function bodies makes baby inner loops cry. Great idea for PETSc vectors where each function body is O(N_dofs), not so great here where we were hoping for O(1).

Non-compiled codes can theoretically be better at this stuff with a smart enough JIT backend. Last I checked Python didn't qualify but if you sacrificed a goat you had a chance with Java.

We should test DualNumber<Real,empty-DSNA> and see how it performs before creating a PotentiallyDualNumber.

PotentiallyDualNumber would be little more than a cut-and-paste-and-tweak job; tedious but not hard.

We wouldn't want a MetaPhysicL::compute_derivatives global, but a static bool PotentiallyDualNumber::compute_derivatives is almost certainly a better idea than a non-static for your use case.