Don't materialize Q

src/continuations/continuations.jl

@@ -3,8 +3,8 @@ using LinearAlgebra
 using StaticArrays: StaticArray, SMatrix, SVector, SArray
 import DiffEqBase: solve, solve!, init, step!
 
-using ..ArrayUtils: _similar, _zeros, isalmostzero, zero_if_nan, _lq!,
-    _normalize!
+using ..ArrayUtils: _similar, _zeros, isalmostzero, zero_if_nan, _lq!, _det,
+    _normalize!, bottomrow
 
 include("interface.jl")
 include("euler_newton.jl")

src/continuations/euler_newton.jl

@@ -83,10 +83,21 @@ ContinuationCache(prob::AbstractContinuationProblem, args...) =
     verbose::Bool = false
 end
 
+rawtangent(Q) = vec(rawtangentmat(Q))
+function rawtangentmat(Q)
+    if Q isa StaticArray
+        return bottomrow(Q)
+    else
+        x = zeros(1, size(Q, 1))
+        x[1, end] = 1
+        rmul!(x, Q)
+        return x
+    end
+end
 
 function tangent(L, Q)
-    tJ = Q[end, :]
-    if det(Q) * det(@view L[1:end-1, 1:end-1]) < 0
+    tJ = rawtangent(Q)
+    if _det(Q) * det(@view L[1:end-1, 1:end-1]) < 0
         tJ *= -1
     end
     return tJ
@@ -102,7 +113,7 @@ function current_tangent(cache::ContinuationCache,
 
     H, J = residual_jacobian!(H, J, u, prob_cache)
     A = vcat(J, _zeros(J, 1, size(J, 2)))  # TODO: improve
-    L, Q = _lq!(Q, A)
+    L, Q = _lq!(A)
     return tangent(L, Q)
 end
 
@@ -118,15 +129,15 @@ function corrector_step!(H::HType,
                          }
     H, J = residual_jacobian!(H, J, v, prob_cache)
     A = vcat(J, _zeros(J, 1, size(J, 2)))  # TODO: improve
-    L, Q = _lq!(Q, A)
-    y = (@view L[1:end-1, 1:end-1]) \ H
-    dv = (@view Q[1:end-1, :])' * y
+    L, Q = _lq!(A)
+    y = vcat((@view L[1:end-1, 1:end-1]) \ H, false)
+    dv = Q' * y
     w = v - dv
     return (w :: vType,
             dv,
             H :: HType,
             L,
-            Q :: QType,
+            Q,
             J :: JType)
 end
 

src/continuations/zero_point.jl

@@ -1,7 +1,7 @@
 using Parameters: @with_kw
 
-calc_direction(_u, J, _L, Q) = det(vcat(J, (@view Q[end:end, :])))
-# Q[:, end] --- tangent w/o det-fixing
+calc_direction(_u, J, _L, Q) = det(vcat(J, rawtangentmat(Q)))
+# rawtangentmat(Q) == Q[end:end, 1] --- tangent w/o det-fixing
 
 find_simple_bifurcation!(cache, opts, sbint) =
     find_simple_bifurcation!(cache, opts, sbint.u0, sbint.u1,
@@ -57,12 +57,12 @@ function find_zero!(cache, opts, f, u0, u1, direction)
 
     H, J = residual_jacobian!(H, J, u1, prob_cache)
     A = vcat(J, _zeros(J, 1, size(J, 2)))  # TODO: improve
-    L, Q = _lq!(Q, A)
+    L, Q = _lq!(A)
     f1 = f(u1, J, L, Q)
 
     H, J = residual_jacobian!(H, J, u0, prob_cache)
     A = vcat(J, _zeros(J, 1, size(J, 2)))  # TODO: improve
-    L, Q = _lq!(Q, A)
+    L, Q = _lq!(A)
     tJ = tangent(L, Q)
     f0 = f(u0, J, L, Q)
 

src/utils/array_utils.jl

@@ -3,6 +3,7 @@ module ArrayUtils
 using Base: typename
 
 using LinearAlgebra
+using LinearAlgebra: QRPackedQ, LQPackedQ
 
 using ForwardDiff: Dual
 using StaticArrays: SVector, SMatrix, StaticArray, Size, similar_type
@@ -84,17 +85,32 @@ zero_if_nan(x) = isnan(x) ? zero(x) : x
     A
 end
 
-function _lq!(Q, A)
-    F = lq!(A)
-    lmul!(F.Q, eye!(Q))  # Q = Matrix(F[:Q])[...]; but lesser allocation
-    return (LowerTriangular(F.L), Q)
+function _lq!(A)
+    L, Q = lq!(A)
+    return (LowerTriangular(L), Q)
 end
 
-function _lq!(_, A::SMatrix)
+function _lq!(A::SMatrix)
     Q, R = qr(A')
     return (LowerTriangular(R'), Q')
 end
 
+# https://github.com/JuliaLang/julia/pull/32887
+_det(Q::Union{QRPackedQ{T}, LQPackedQ{T}}) where {T <: Real} =
+    isodd(count(!iszero, Q.τ)) ? -1 : 1
+_det(Q::StaticArray) = det(Q)
+
+@inline foldlargs(op, x) = x
+@inline foldlargs(op, x1, x2, xs...) = foldlargs(op, op(x1, x2), xs...)
+
+bottomrow(M::AbstractArray) = @view M[end:end, :]
+bottomrow(M::SMatrix{S1, S2}) where {S1, S2} =
+    SMatrix{1, S2}(
+        foldlargs((), ntuple(identity, S2)...) do xs, i
+            (xs..., @inbounds M[S1, i])
+        end
+    )
+
 function _normalize!(x)
     normalize!(x)
     return x

test/test_utils.jl

@@ -39,8 +39,7 @@ end
 
         L0, Q0 = lq(Array(A))
 
-        Q1 = copy(Q0)
-        L1, Q1 = _lq!(Q1, copy(A))
+        L1, Q1 = _lq!(copy(A))
 
         @test L0 ≈ L1
         @test Q0 ≈ Q1