LowRankFactorizationNEP to NEPTypes

[jarlebring]

Moved the LowRankFactorizedNEP to NEPTypes. Made schrodinger_movebc use it + added a bit of functionality. Not in this PR: Some of the compute functions can be made more efficient by not constructing the full matrix, but only working with the factorizations. Some more Gallery examples could benefit from using this (as commented in the gallery).

I'm adding max as reviewer as courtesy. I know you're busy.

[jarlebring]

This should handle #157

[jarlebring]

I'm going to merge this now. A note for future: we can consider using the LowRankMatrix in the package LowRankApprox. They provide a type LowRankMatrix, which is an AbstractMatrix and can be handed into the SPMF, cf https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838.

[maxbennedich]

Sorry for the delayed response, I think these changes look good, I had just a couple of stylistic suggestions which shouldn't affect functionality or performance.

[maxbennedich]

I guess this could be removed now? Along with G, and the definition of F could optionally be done directly into Uv2.

[jarlebring]

Thanks. I've fixed it in master.

[maxbennedich]

We could consider simplifying these methods a bit, and even combining them to get rid of the duplicated code. I like to use higher order functions and fusing dots since it usually handles type stability and avoids unnecessary allocations. I think the above two methods could be written as a single method like this:

function LowRankFactorizedNEP(L::AbstractVector{S}, U::AbstractVector{S}, f, A = L .* adjoint.(U)) where {S<:AbstractMatrix}
    rank = mapreduce(u -> size(u, 2), +, U)
    return LowRankFactorizedNEP{S}(SPMF_NEP(A, f, align_sparsity_patterns=true), rank, L, U)
end

Alternatively, keeping the two methods if we don't want to compute A as part of the argument list.

[jarlebring]

Thx. I went for two functions.

[maxbennedich]

Alternatively, with higher order functions and dot notation:

function LowRankFactorizedNEP(Amf::AbstractVector{LowRankMatrixAndFunction{S}}) where {T<:Number, S<:AbstractMatrix{T}}
    # if A is not specified, create it from LU factors
    A = [isempty(M.A) ? M.L * M.U' : M.A for M in Amf]
    L = getfield.(Amf, :L)
    U = getfield.(Amf, :U)
    f = getfield.(Amf, :f)
    rank = mapreduce(M -> size(M.U, 2), +, Amf)
    return LowRankFactorizedNEP(SPMF_NEP(A, f, align_sparsity_patterns=true), rank, L, U)
end

[jarlebring]

Sweet. Fixed in master.

[maxbennedich]

Btw, does it make sense to use the conjugate transpose here? I believe that I did it this way since that's how it was implemented in the original NLEIGS code, but would this be preferred (and not confusing) over just A = L * U, without transpose?

[jarlebring]

I think it's just because linear algebra researcher tend to prefer tall skinny matrices over short and fat. It's the same in LowRankMatrix in LowRankApprox-package. https://github.com/JuliaMatrices/LowRankApprox.jl/blob/master/src/lowrankmatrix.jl . To me what is currently called low rank matrix (both in our code and LowRankApprox) has not so much to do with the rank, but just a representation of a product of two matrices. A better name is maybe just be MatrixProductNEP, where every term is product of two matrices. Such a change would be too big for right now.