Fixing Julia for Thomas algorithm

[leios]

Yeah. That last version of Julia was complete and utter garbage. This fixes that and is modelled after the C code.

[Butt4cak3]

This new implementation that returns the solution vector still modifies the c vector. Should we just make it so that the algorithm doesn't modify the input vectors at all? The C implementation is an example of how that could happen. It creates an intermediate double y[size].

[leios]

to be fair, this implementation seems to modify both c and d and the C solution modifies what I called d in this case, but C called x.

This is a bit of a mess all around

[jiegillet]

+1 on not modifying any input. Functional programming FTW :)

[Butt4cak3]

It's not that difficult to change the algorithm to not modify c and d. Create two more arrays, one is the result and the other is there to be modified instead of c and you're good to go.

[leios]

To be fair, we are heavily using the mathematical notation of the chapter. While we are revising the implementation, maybe we should also change the variable names to diagonal, lower_diagonal, and upper_diagonal?

Then the modified temp variable can be solution and temp_coefficient (I don't know about that last one). All the texts usually use a, b, c, and d, which I am fine with, and then call the new variables c' and d', but we could do c_prime and d_prime if we stick with that notation.

[jiegillet]

We are not questioning the notations a, b and c I think, just the fact that the inputs should not be modified by the algorithm.

[leios]

I was bringing up the topic because a, b, c, and d are really unclear and it's something I thought we could clean up while we are at it. Otherwise, we can do precisely what you guys are saying by adding the c_prime and d_prime variables (which is common mathematical notation).

[jiegillet]

I like c_prime and d_prime, but the chapter also mentions x with is kind of d_prime_prime, so maybe just using c_prime and x would be good.

[leios]

x is the set of solutions, which we can use d_prime for. I'll rework this with c_prime and x.

[berquist]

I usually suggest things like a = a / b get converted to a /= b, but what do you think? No opinion here.

[berquist]

These need to be

x = copy(d)
c_prime = copy(d)

[berquist]

(otherwise you're still overwriting the input.)

[berquist]

I was not here for the modifying input vs. returning conversation, but my 2 cents:

I don't think the average numerical algorithm is more understandable by being made pure, and most people who want to understand the algorithm so they can implement it should rather implement the version that mutates a preallocated input argument. Maybe it should go as far as using other input arguments as scratch, as this and other implementations used to. This is how many calls to BLAS and LAPACK work. Passing matrices around with copy-by-value semantics is a death sentence for performance.

Perhaps the solution is to say "there should be no surprises". This is well-documented behavior in real code, but not in the AAA, so being pure is good. Sorry for the inconclusive language, but I never read this rebuttal against "make everything more functional" even though writing pure code prevents necessary memory optimizations in numerical code.

/rant over

TL;DR looks good.

[leios]

I agree with @berquist , though. I am not sure what is the best way to do this. I know we are returning values for readability, but in some languages, this does kill performance and code shouldn't be written in general like that.

[jiegillet]

We were never concerned with performance though, so I'm not worried about that.
I'm just thinking a function like this one should be pure in case you need to use your matrix for other purposes later :)

[berquist]

That is possible, yes. It is also possible that the people who would be affected already know about this sort of thing. I am fine with making everything pure, but only pointing out something I have never seen written down that I think should be revisited in some form later (chapter?).