Add alloc-free in-place Tridiagonal solves and lu!

[dkarrasch]

This implements @stevengj's suggestions, lu!(F::LU{<:Tridiagonal}, ::Tridiagonal), which essentially copies the data from the second argument to the first and applies the in-place lu!. This avoids the allocation of ipiv (and potentially du2) and is therefore completely allocation-free. Second, this adds an ldiv!(::Tridiagonal, VecOrMat) which combines the in-place lu and the in-place solve. Naturally, this needs to overwrite the tridiagonal matrix, but this seems to be fine for some use cases, see CliMA/ClimaAtmos.jl#715. It seems that this ldiv! is even slightly faster than the Thomas algorithm(must have been a measurement artifact), but comes with the advantage of stability due to partial pivoting.

Closes #46099.

[dkarrasch]

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.11.0-DEV.64 (2023-07-12)
 _/ |\__'_|_|_|\__'_|  |  dk/lutridiag/90fbd6300c* (fork: 1 commits, 1 day)
|__/                   |

julia> using LinearAlgebra, BenchmarkTools

julia> T = Tridiagonal(randn(1000), randn(1001), randn(1000));

julia> B = randn(1001,100);

julia> L = lu(T);

julia> T \ B ≈ ldiv!(copy(T), copy(B))
true

julia> T \ B ≈ LAPACK.gtsv!(copy(T.dl), copy(T.d), copy(T.du), copy(B))
true

julia> @btime LAPACK.gtsv!($(copy(T.dl)), $(copy(T.d)), $(copy(T.du)), $(copy(B)));
  911.000 μs (0 allocations: 0 bytes)

julia> @btime $T \ $B;
  1.159 ms (9 allocations: 822.02 KiB)

julia> @btime lu($T);
  11.318 μs (7 allocations: 39.91 KiB)

julia> @btime ldiv!($L, $(copy(B)));
  1.060 ms (0 allocations: 0 bytes)

julia> @btime ldiv!($(copy(T)), $(copy(B)));
  841.458 μs (0 allocations: 0 bytes)

Comparing to ClimaAtmos.jl's Thomas algorithm implementation (defined only for vector rhs's):

function thomas_algorithm!(A, b)
    nrows = size(A, 1)
    # first row
    @inbounds A[1, 2] /= A[1, 1]
    @inbounds b[1] /= A[1, 1]
    # interior rows
    for row in 2:(nrows - 1)
        @inbounds fac = A[row, row] - (A[row, row - 1] * A[row - 1, row])
        @inbounds A[row, row + 1] /= fac
        @inbounds b[row] = (b[row] - A[row, row - 1] * b[row - 1]) / fac
    end
    # last row
    @inbounds fac = A[nrows, nrows] - A[nrows - 1, nrows] * A[nrows, nrows - 1]
    @inbounds b[nrows] = (b[nrows] - A[nrows, nrows - 1] * b[nrows - 1]) / fac
    # back substitution
    for row in (nrows - 1):-1:1
        @inbounds b[row] -= b[row + 1] * A[row, row + 1]
    end
    return b
end

julia> b = randn(1001);

julia> @btime thomas_algorithm!($(copy(T)), $(copy(b)));
  8.957 μs (0 allocations: 0 bytes)

julia> @btime ldiv!($(copy(T)), $(copy(b)));
  14.488 μs (0 allocations: 0 bytes)

[dkarrasch]

I'd merge this soon if there are no objections.