Rewrite forward_substitution and back_substitution

[Andlon]

I had some problems with these functions because they returned an error for empty matrices. In my opinion, solving a system Ax = b for empty A and b should merely return an empty x. I've had some major grievances in the past with this exact situation with Eigen, which doesn't handle this situation very well. My application generated linear systems from geometrical meshes, and sometimes these systems turned out to be "empty". Graceful handling is thus important so that users don't need to special-case these scenarios.

In any case, I took the liberty to completely rewrite the forward and backward substitution functions. I used some of the techniques I discovered when working on #150 to rewrite them in terms of the highly optimized utils::dot function, as well as leveraging our row iterators. They no longer contain any unsafe code, except for needing to call get_unchecked(i, i) from BaseMatrix. See #151.

The resulting code is significantly simpler, does not allocate a new vector like the previous implementation and is also significantly faster. Here are the benchmarks before making the changes:

test linalg::triangular::solve_l_triangular_10000x10000 ... bench:  51,194,669 ns/iter (+/- 323,709)
test linalg::triangular::solve_l_triangular_1000x1000   ... bench:     559,368 ns/iter (+/- 8,496)
test linalg::triangular::solve_l_triangular_100x100     ... bench:       5,394 ns/iter (+/- 550)
test linalg::triangular::solve_u_triangular_10000x10000 ... bench:  50,927,341 ns/iter (+/- 879,902)
test linalg::triangular::solve_u_triangular_1000x1000   ... bench:     550,613 ns/iter (+/- 12,704)
test linalg::triangular::solve_u_triangular_100x100     ... bench:       4,708 ns/iter (+/- 17)

Here are the benchmarks after making the changes:

test linalg::triangular::solve_l_triangular_10000x10000 ... bench:  29,930,075 ns/iter (+/- 263,948)
test linalg::triangular::solve_l_triangular_1000x1000   ... bench:     296,419 ns/iter (+/- 10,734)
test linalg::triangular::solve_l_triangular_100x100     ... bench:       2,574 ns/iter (+/- 3)
test linalg::triangular::solve_u_triangular_10000x10000 ... bench:  30,100,720 ns/iter (+/- 325,360)
test linalg::triangular::solve_u_triangular_1000x1000   ... bench:     295,012 ns/iter (+/- 3,668)
test linalg::triangular::solve_u_triangular_100x100     ... bench:       3,188 ns/iter (+/- 4)

[Andlon]

Update: I took the liberty of removing the Any trait bounds from these functions too. AFAIK, removing a trait bound is not a breaking change...?

[AtheMathmo]

This is awesome! I don't think that removing Any is a breaking change either. I'll try to check this over properly later today. (It's a busy week for me with a midterm)

[Andlon]

@AtheMathmo: please take as long as you want with the reviews. I'd rather you focus properly on your midterm :) I just found myself with some unexpected time to work on rulinalg, because my current temporary situation precludes me from doing my actual thesis work.

[AtheMathmo]

I massively appreciate all of the changes you've put in so far - they're big leaps in usability and performance. Thank you!

[Andlon]

Thanks for saying so! I'm learning a lot myself in the process, so I'm very happy to work on it :)

[AtheMathmo]

I merged in #142 which has given some conflicts here. I'll try to review regardless.

[AtheMathmo]

This is definitely fine - I think if we are within 1e-8 we can safely assume the implementation is correct. We might be missing some tricks to improve accuracy but these should be captured in a separate test where necessary.

[AtheMathmo]

This looks good to me. Once the conflicts are resolved I am happy to merge.

[Andlon]

Should be good to go now assuming tests complete successfully. Thanks for taking a look!