Adjustments to AdjointQ in LinearAlgebra.jl (#114)

Author: dkarrasch

Project.toml

@@ -1,7 +1,7 @@
 name = "ArrayLayouts"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.8.14"
+version = "0.8.15"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

src/ArrayLayouts.jl

@@ -23,6 +23,8 @@ using LinearAlgebra: AbstractTriangular, AbstractQ, QRCompactWYQ, QRPackedQ, che
 
 using LinearAlgebra.BLAS: BlasFloat, BlasReal, BlasComplex
 
+AdjointQtype{T} = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ{T} : Adjoint{T,<:AbstractQ}
+
 using FillArrays: AbstractFill, getindex_value, axes_print_matrix_row, _copy_oftype
 
 using Base: require_one_based_indexing

src/factorizations.jl

@@ -196,14 +196,14 @@ end
 
 
 ### QcB
-materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
-    (A = M.A.parent; LAPACK.ormqr!('L','T',A.factors,A.τ,M.B))
-materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasComplex =
-    (A = M.A.parent; LAPACK.ormqr!('L','C',A.factors,A.τ,M.B))
-materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
-    (A = M.A.parent; LAPACK.gemqrt!('L','T',A.factors,A.T,M.B))
-materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AbstractMatrix{T},<:AbstractVecOrMat{T}}) where T<:BlasComplex =
-    (A = M.A.parent; LAPACK.gemqrt!('L','C',A.factors,A.T,M.B))
+materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AdjointQtype{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
+    (A = parent(M.A); LAPACK.ormqr!('L','T',A.factors,A.τ,M.B))
+materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AdjointQtype{T},<:AbstractVecOrMat{T}}) where T<:BlasComplex =
+    (A = parent(M.A); LAPACK.ormqr!('L','C',A.factors,A.τ,M.B))
+materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AdjointQtype{T},<:AbstractVecOrMat{T}}) where T<:BlasFloat =
+    (A = parent(M.A); LAPACK.gemqrt!('L','T',A.factors,A.T,M.B))
+materialize!(M::Lmul{<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractStridedLayout,<:AdjointQtype{T},<:AbstractVecOrMat{T}}) where T<:BlasComplex =
+    (A = parent(M.A); LAPACK.gemqrt!('L','C',A.factors,A.T,M.B))
 function materialize!(M::Lmul{<:AdjQRPackedQLayout})
     adjA,B = M.A, M.B
     require_one_based_indexing(B)
@@ -264,14 +264,14 @@ function materialize!(M::Rmul{<:Any,<:QRPackedQLayout})
 end
 
 ### AQc
-materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasReal =
-    (B = M.B.parent; LAPACK.ormqr!('R','T',B.factors,B.τ,M.A))
-materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasComplex =
-    (B = M.B.parent; LAPACK.ormqr!('R','C',B.factors,B.τ,M.A))
-materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasReal =
-    (B = M.B.parent; LAPACK.gemqrt!('R','T',B.factors,B.T,M.A))
-materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AbstractMatrix{T}}) where T<:BlasComplex =
-    (B = M.B.parent; LAPACK.gemqrt!('R','C',B.factors,B.T,M.A))
+materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AdjointQtype{T}}) where T<:BlasReal =
+    (B = parent(M.B); LAPACK.ormqr!('R','T',B.factors,B.τ,M.A))
+materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRPackedQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AdjointQtype{T}}) where T<:BlasComplex =
+    (B = parent(M.B); LAPACK.ormqr!('R','C',B.factors,B.τ,M.A))
+materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AdjointQtype{T}}) where T<:BlasReal =
+    (B = parent(M.B); LAPACK.gemqrt!('R','T',B.factors,B.T,M.A))
+materialize!(M::Rmul{<:AbstractStridedLayout,<:AdjQRCompactWYQLayout{<:AbstractStridedLayout,<:AbstractStridedLayout},<:AbstractVecOrMat{T},<:AdjointQtype{T}}) where T<:BlasComplex =
+    (B = parent(M.B); LAPACK.gemqrt!('R','C',B.factors,B.T,M.A))
 function materialize!(M::Rmul{<:Any,<:AdjQRPackedQLayout})
     A,adjQ = M.A,M.B
     Q = parent(adjQ)
@@ -379,6 +379,8 @@ function _cholesky!(::SymmetricLayout{<:AbstractColumnMajor}, axes, A::AbstractM
     ::CRowMaximum; tol = 0.0, check::Bool = true)
     AA, piv, rank, info = LAPACK.pstrf!(A.uplo, A.data, tol)
     C = CholeskyPivoted{eltype(AA),typeof(AA)}(AA, A.uplo, piv, rank, tol, info)
+    # TODO: when lower bound on julia is above v1.8, replace the above line by the following 
+    # C = CholeskyPivoted{eltype(AA),typeof(AA),typeof(piv)}(AA, A.uplo, piv, rank, tol, info)
     check && chkfullrank(C)
     return C
 end

src/memorylayout.jl

@@ -283,6 +283,9 @@ MemoryLayout(::Type{Transpose{T,P}}) where {T,P} = transposelayout(MemoryLayout(
 MemoryLayout(::Type{Adjoint{T,P}}) where {T,P} = adjointlayout(T, MemoryLayout(P))
 MemoryLayout(::Type{AdjointAbsVec{T,P}}) where {T,P<:AbstractVector} = DualLayout{typeof(adjointlayout(T,MemoryLayout(P)))}()
 MemoryLayout(::Type{TransposeAbsVec{T,P}}) where {T,P<:AbstractVector} = DualLayout{typeof(transposelayout(MemoryLayout(P)))}()
+if isdefined(LinearAlgebra, :AdjointQ)
+    MemoryLayout(::Type{LinearAlgebra.AdjointQ{T,P}}) where {T,P} = adjointlayout(T, MemoryLayout(P))
+end
 
 transposelayout(_) = UnknownLayout()
 transposelayout(::StridedLayout) = StridedLayout()

src/mul.jl

@@ -118,6 +118,7 @@ mul!(dest::AbstractArray, A::AbstractArray, B::AbstractArray, α::Number, β::Nu
 
 broadcastable(M::Mul) = M
 
+const FlexibleLeftQs = Union{LinearAlgebra.HessenbergQ, LinearAlgebra.QRCompactWYQ, LinearAlgebra.QRPackedQ}
 
 macro veclayoutmul(Typ)
     ret = quote
@@ -134,6 +135,16 @@ macro veclayoutmul(Typ)
         (*)(A::LinearAlgebra.AbstractQ, B::$Typ) = ArrayLayouts.mul(A,B)
         (*)(A::$Typ, B::LinearAlgebra.LQPackedQ) = ArrayLayouts.mul(A,B)
     end
+    if isdefined(LinearAlgebra, :AdjointQ)
+        ret = quote
+            $ret
+
+            # disambiguation for flexible left-mul Qs
+            (*)(A::FlexibleLeftQs, B::$Typ) = ArrayLayouts.mul(A,B)
+            # flexible right-mul/adjoint left-mul Qs
+            (*)(A::LinearAlgebra.AdjointQ{<:Any,<:LinearAlgebra.LQPackedQ}, B::$Typ) = ArrayLayouts.mul(A,B)
+        end
+    end
     for Struc in (:AbstractTriangular, :Diagonal)
         ret = quote
             $ret
@@ -198,6 +209,18 @@ macro layoutmul(Typ)
         (*)(A::LinearAlgebra.AbstractQ, B::$Typ) = ArrayLayouts.mul(A,B)
         (*)(A::$Typ, B::LinearAlgebra.AbstractQ) = ArrayLayouts.mul(A,B)
     end
+    if isdefined(LinearAlgebra, :AdjointQ)
+        ret = quote
+            $ret
+
+            # disambiguation for flexible left-mul/adjoint right-mul Qs
+            (*)(A::FlexibleLeftQs, B::$Typ) = ArrayLayouts.mul(A,B)
+            (*)(A::$Typ, B::LinearAlgebra.AdjointQ{<:Any,<:FlexibleLeftQs}) = ArrayLayouts.mul(A,B)
+            # disambiguation for flexible right-mul/adjoint left-mul Qs
+            (*)(A::$Typ, B::LinearAlgebra.LQPackedQ) = ArrayLayouts.mul(A,B)
+            (*)(A::LinearAlgebra.AdjointQ{<:Any,<:LinearAlgebra.LQPackedQ}, B::$Typ) = ArrayLayouts.mul(A,B)
+        end
+    end
     for Struc in (:AbstractTriangular, :Diagonal)
         ret = quote
             $ret