@@ -35,6 +35,7 @@ convert(::Type{AbstractQ{T}}, adjQ::AdjointQ{T}) where {T} = adjQ
3535convert (:: Type{AbstractQ{T}} , adjQ:: AdjointQ ) where {T} = convert (AbstractQ{T}, adjQ. Q)'
3636
3737# ... to matrix
38+ collect (Q:: AbstractQ ) = copyto! (Matrix {eltype(Q)} (undef, size (Q)), Q)
3839Matrix {T} (Q:: AbstractQ ) where {T} = convert (Matrix{T}, Q* I) # generic fallback, yields square matrix
3940Matrix {T} (adjQ:: AdjointQ{S} ) where {T,S} = convert (Matrix{T}, lmul! (adjQ, Matrix {S} (I, size (adjQ))))
4041Matrix (Q:: AbstractQ{T} ) where {T} = Matrix {T} (Q)
@@ -56,6 +57,15 @@ function size(Q::AbstractQ, dim::Integer)
5657end
5758size (adjQ:: AdjointQ ) = reverse (size (adjQ. Q))
5859
60+ # comparison
61+ (== )(Q:: AbstractQ , A:: AbstractMatrix ) = lmul! (Q, Matrix {eltype(Q)} (I, size (A))) == A
62+ (== )(A:: AbstractMatrix , Q:: AbstractQ ) = Q == A
63+ (== )(Q:: AbstractQ , P:: AbstractQ ) = Matrix (Q) == Matrix (P)
64+ isapprox (Q:: AbstractQ , A:: AbstractMatrix ; kwargs... ) =
65+ isapprox (lmul! (Q, Matrix {eltype(Q)} (I, size (A))), A, kwargs... )
66+ isapprox (A:: AbstractMatrix , Q:: AbstractQ ; kwargs... ) = isapprox (Q, A, kwargs... )
67+ isapprox (Q:: AbstractQ , P:: AbstractQ ; kwargs... ) = isapprox (Matrix (Q), Matrix (P), kwargs... )
68+
5969# pseudo-array behaviour, required for indexing with `begin` or `end`
6070axes (Q:: AbstractQ ) = map (Base. oneto, size (Q))
6171axes (Q:: AbstractQ , d:: Integer ) = d in (1 , 2 ) ? axes (Q)[d] : Base. OneTo (1 )
@@ -125,36 +135,60 @@ function show(io::IO, ::MIME{Symbol("text/plain")}, Q::AbstractQ)
125135end
126136
127137# multiplication
138+ # generically, treat AbstractQ like a matrix with its definite size
139+ qsize_check (Q:: AbstractQ , B:: AbstractVecOrMat ) =
140+ size (Q, 2 ) == size (B, 1 ) ||
141+ throw (DimensionMismatch (" second dimension of Q, $(size (Q,2 )) , must coincide with first dimension of B, $(size (B,1 )) " ))
142+ qsize_check (A:: AbstractVecOrMat , Q:: AbstractQ ) =
143+ size (A, 2 ) == size (Q, 1 ) ||
144+ throw (DimensionMismatch (" second dimension of A, $(size (A,2 )) , must coincide with first dimension of Q, $(size (Q,1 )) " ))
145+ qsize_check (Q:: AbstractQ , P:: AbstractQ ) =
146+ size (Q, 2 ) == size (P, 1 ) ||
147+ throw (DimensionMismatch (" second dimension of A, $(size (Q,2 )) , must coincide with first dimension of B, $(size (P,1 )) " ))
148+
128149(* )(Q:: AbstractQ , J:: UniformScaling ) = Q* J. λ
129150function (* )(Q:: AbstractQ , b:: Number )
130151 T = promote_type (eltype (Q), typeof (b))
131152 lmul! (convert (AbstractQ{T}, Q), Matrix {T} (b* I, size (Q)))
132153end
133- function (* )(A:: AbstractQ , B:: AbstractVecOrMat )
134- T = promote_type (eltype (A), eltype (B))
135- lmul! (convert (AbstractQ{T}, A), copy_similar (B, T))
154+ function (* )(Q:: AbstractQ , B:: AbstractVector )
155+ T = promote_type (eltype (Q), eltype (B))
156+ qsize_check (Q, B)
157+ mul! (similar (B, T, size (Q, 1 )), convert (AbstractQ{T}, Q), B)
158+ end
159+ function (* )(Q:: AbstractQ , B:: AbstractMatrix )
160+ T = promote_type (eltype (Q), eltype (B))
161+ qsize_check (Q, B)
162+ mul! (similar (B, T, (size (Q, 1 ), size (B, 2 ))), convert (AbstractQ{T}, Q), B)
136163end
137164
138165(* )(J:: UniformScaling , Q:: AbstractQ ) = J. λ* Q
139166function (* )(a:: Number , Q:: AbstractQ )
140167 T = promote_type (typeof (a), eltype (Q))
141168 rmul! (Matrix {T} (a* I, size (Q)), convert (AbstractQ{T}, Q))
142169end
143- * (a:: AbstractVector , Q:: AbstractQ ) = reshape (a, length (a), 1 ) * Q
170+ function (* )(A:: AbstractVector , Q:: AbstractQ )
171+ T = promote_type (eltype (A), eltype (Q))
172+ qsize_check (A, Q)
173+ return mul! (similar (A, T, length (A)), A, convert (AbstractQ{T}, Q))
174+ end
144175function (* )(A:: AbstractMatrix , Q:: AbstractQ )
145176 T = promote_type (eltype (A), eltype (Q))
146- return rmul! (copy_similar (A, T), convert (AbstractQ{T}, Q))
177+ qsize_check (A, Q)
178+ return mul! (similar (A, T, (size (A, 1 ), size (Q, 2 ))), A, convert (AbstractQ{T}, Q))
147179end
148180(* )(u:: AdjointAbsVec , Q:: AbstractQ ) = (Q' u' )'
149181
150182# ## Q*Q (including adjoints)
151- * (Q:: AbstractQ , P:: AbstractQ ) = Q * (P* I)
183+ ( * ) (Q:: AbstractQ , P:: AbstractQ ) = Q * (P* I)
152184
153185# ## mul!
154- function mul! (C:: AbstractVecOrMat{T} , Q:: AbstractQ{T} , B:: Union{AbstractVecOrMat{T} ,AbstractQ{T} } ) where {T}
186+ function mul! (C:: AbstractVecOrMat{T} , Q:: AbstractQ{T} , B:: Union{AbstractVecOrMat,AbstractQ} ) where {T}
155187 require_one_based_indexing (C, B)
156- mB = size (B, 1 )
157- mC = size (C, 1 )
188+ mB, nB = size (B, 1 ), size (B, 2 )
189+ mC, nC = size (C, 1 ), size (C, 2 )
190+ qsize_check (Q, B)
191+ nB != nC && throw (DimensionMismatch ())
158192 if mB < mC
159193 inds = CartesianIndices (axes (B))
160194 copyto! (view (C, inds), B)
@@ -164,9 +198,21 @@ function mul!(C::AbstractVecOrMat{T}, Q::AbstractQ{T}, B::Union{AbstractVecOrMat
164198 return lmul! (Q, copyto! (C, B))
165199 end
166200end
167- mul! (C:: AbstractVecOrMat{T} , A:: AbstractVecOrMat{T} , Q:: AbstractQ{T} ) where {T} = rmul! (copyto! (C, A), Q)
168- mul! (C:: AbstractVecOrMat{T} , adjQ:: AdjointQ{T} , B:: AbstractVecOrMat{T} ) where {T} = lmul! (adjQ, copyto! (C, B))
169- mul! (C:: AbstractVecOrMat{T} , A:: AbstractVecOrMat{T} , adjQ:: AdjointQ{T} ) where {T} = rmul! (copyto! (C, A), adjQ)
201+ function mul! (C:: AbstractVecOrMat{T} , A:: AbstractVecOrMat , Q:: AbstractQ{T} ) where {T}
202+ require_one_based_indexing (C, A)
203+ mA, nA = size (A, 1 ), size (A, 2 )
204+ mC, nC = size (C, 1 ), size (C, 2 )
205+ mA != mC && throw (DimensionMismatch ())
206+ qsize_check (A, Q)
207+ if nA < nC
208+ inds = CartesianIndices (axes (A))
209+ copyto! (view (C, inds), A)
210+ C[CartesianIndices ((axes (C, 1 ), nA+ 1 : nC))] .= zero (T)
211+ return rmul! (C, Q)
212+ else
213+ return rmul! (copyto! (C, A), Q)
214+ end
215+ end
170216
171217# ## division
172218\ (Q:: AbstractQ , A:: AbstractVecOrMat ) = Q' * A
@@ -319,7 +365,7 @@ rmul!(A::StridedVecOrMat{T}, B::QRCompactWYQ{T,<:StridedMatrix}) where {T<:BlasF
319365 LAPACK. gemqrt! (' R' , ' N' , B. factors, B. T, A)
320366rmul! (A:: StridedVecOrMat{T} , B:: QRPackedQ{T,<:StridedMatrix} ) where {T<: BlasFloat } =
321367 LAPACK. ormqr! (' R' , ' N' , B. factors, B. τ, A)
322- function rmul! (A:: AbstractMatrix , Q:: QRPackedQ )
368+ function rmul! (A:: AbstractVecOrMat , Q:: QRPackedQ )
323369 require_one_based_indexing (A)
324370 mQ, nQ = size (Q. factors)
325371 mA, nA = size (A,1 ), size (A,2 )
@@ -354,7 +400,7 @@ rmul!(A::StridedVecOrMat{T}, adjQ::AdjointQ{<:Any,<:QRPackedQ{T}}) where {T<:Bla
354400 (Q = adjQ. Q; LAPACK. ormqr! (' R' , ' T' , Q. factors, Q. τ, A))
355401rmul! (A:: StridedVecOrMat{T} , adjQ:: AdjointQ{<:Any,<:QRPackedQ{T}} ) where {T<: BlasComplex } =
356402 (Q = adjQ. Q; LAPACK. ormqr! (' R' , ' C' , Q. factors, Q. τ, A))
357- function rmul! (A:: AbstractMatrix , adjQ:: AdjointQ{<:Any,<:QRPackedQ} )
403+ function rmul! (A:: AbstractVecOrMat , adjQ:: AdjointQ{<:Any,<:QRPackedQ} )
358404 require_one_based_indexing (A)
359405 Q = adjQ. Q
360406 mQ, nQ = size (Q. factors)
@@ -459,42 +505,12 @@ lmul!(adjQ::AdjointQ{<:Any,<:HessenbergQ{T}}, X::Adjoint{T,<:StridedVecOrMat{T}}
459505rmul! (X:: Adjoint{T,<:StridedVecOrMat{T}} , adjQ:: AdjointQ{<:Any,<:HessenbergQ{T}} ) where {T} = lmul! (adjQ' , X' )'
460506
461507# flexible left-multiplication (and adjoint right-multiplication)
462- function (* )(Q:: Union{QRPackedQ,QRCompactWYQ,HessenbergQ} , b:: AbstractVector )
463- T = promote_type (eltype (Q), eltype (b))
464- if size (Q. factors, 1 ) == length (b)
465- bnew = copy_similar (b, T)
466- elseif size (Q. factors, 2 ) == length (b)
467- bnew = [b; zeros (T, size (Q. factors, 1 ) - length (b))]
468- else
469- throw (DimensionMismatch (" vector must have length either $(size (Q. factors, 1 )) or $(size (Q. factors, 2 )) " ))
470- end
471- lmul! (convert (AbstractQ{T}, Q), bnew)
472- end
473- function (* )(Q:: Union{QRPackedQ,QRCompactWYQ,HessenbergQ} , B:: AbstractMatrix )
474- T = promote_type (eltype (Q), eltype (B))
475- if size (Q. factors, 1 ) == size (B, 1 )
476- Bnew = copy_similar (B, T)
477- elseif size (Q. factors, 2 ) == size (B, 1 )
478- Bnew = [B; zeros (T, size (Q. factors, 1 ) - size (B,1 ), size (B, 2 ))]
479- else
480- throw (DimensionMismatch (" first dimension of matrix must have size either $(size (Q. factors, 1 )) or $(size (Q. factors, 2 )) " ))
481- end
482- lmul! (convert (AbstractQ{T}, Q), Bnew)
483- end
484- function (* )(A:: AbstractMatrix , adjQ:: AdjointQ{<:Any,<:Union{QRPackedQ,QRCompactWYQ,HessenbergQ}} )
485- Q = adjQ. Q
486- T = promote_type (eltype (A), eltype (adjQ))
487- adjQQ = convert (AbstractQ{T}, adjQ)
488- if size (A, 2 ) == size (Q. factors, 1 )
489- AA = copy_similar (A, T)
490- return rmul! (AA, adjQQ)
491- elseif size (A, 2 ) == size (Q. factors, 2 )
492- return rmul! ([A zeros (T, size (A, 1 ), size (Q. factors, 1 ) - size (Q. factors, 2 ))], adjQQ)
493- else
494- throw (DimensionMismatch (" matrix A has dimensions $(size (A)) but Q-matrix B has dimensions $(size (adjQ)) " ))
495- end
496- end
497- (* )(u:: AdjointAbsVec , Q:: AdjointQ{<:Any,<:Union{QRPackedQ,QRCompactWYQ,HessenbergQ}} ) = (Q' u' )'
508+ qsize_check (Q:: Union{QRPackedQ,QRCompactWYQ,HessenbergQ} , B:: AbstractVecOrMat ) =
509+ size (B, 1 ) in size (Q. factors) ||
510+ throw (DimensionMismatch (" first dimension of B, $(size (B,1 )) , must equal one of the dimensions of Q, $(size (Q. factors)) " ))
511+ qsize_check (A:: AbstractVecOrMat , adjQ:: AdjointQ{<:Any,<:Union{QRPackedQ,QRCompactWYQ,HessenbergQ}} ) =
512+ (Q = adjQ. Q; size (A, 2 ) in size (Q. factors) ||
513+ throw (DimensionMismatch (" second dimension of A, $(size (A,2 )) , must equal one of the dimensions of Q, $(size (Q. factors)) " )))
498514
499515det (Q:: HessenbergQ ) = _det_tau (Q. τ)
500516
@@ -518,104 +534,41 @@ convert(::Type{AbstractQ{T}}, Q::LQPackedQ) where {T} = LQPackedQ{T}(Q)
518534size (Q:: LQPackedQ ) = (n = size (Q. factors, 2 ); return n, n)
519535
520536# # Multiplication
521- # ## QB / QcB
522- lmul! (A:: LQPackedQ{T} , B:: StridedVecOrMat{T} ) where {T<: BlasFloat } = LAPACK. ormlq! (' L' ,' N' ,A. factors,A. τ,B)
523- lmul! (adjA:: AdjointQ{<:Any,<:LQPackedQ{T}} , B:: StridedVecOrMat{T} ) where {T<: BlasReal } =
524- (A = adjA. Q; LAPACK. ormlq! (' L' , ' T' , A. factors, A. τ, B))
525- lmul! (adjA:: AdjointQ{<:Any,<:LQPackedQ{T}} , B:: StridedVecOrMat{T} ) where {T<: BlasComplex } =
526- (A = adjA. Q; LAPACK. ormlq! (' L' , ' C' , A. factors, A. τ, B))
537+ # out-of-place right application of LQPackedQs
538+ #
539+ # these methods: (1) check whether the applied-to matrix's (A's) appropriate dimension
540+ # (columns for A_*, rows for Ac_*) matches the number of columns (nQ) of the LQPackedQ (Q),
541+ # and if so effectively apply Q's square form to A without additional shenanigans; and
542+ # (2) if the preceding dimensions do not match, check whether the appropriate dimension of
543+ # A instead matches the number of rows of the matrix of which Q is a factor (i.e.
544+ # size(Q.factors, 1)), and if so implicitly apply Q's truncated form to A by zero extending
545+ # A as necessary for check (1) to pass (if possible) and then applying Q's square form
527546
528- function (* )(adjA:: AdjointQ{<:Any,<:LQPackedQ} , B:: AbstractVector )
529- A = adjA. Q
530- T = promote_type (eltype (A), eltype (B))
531- if length (B) == size (A. factors, 2 )
532- C = copy_similar (B, T)
533- elseif length (B) == size (A. factors, 1 )
534- C = [B; zeros (T, size (A. factors, 2 ) - size (A. factors, 1 ), size (B, 2 ))]
535- else
536- throw (DimensionMismatch (" length of B, $(length (B)) , must equal one of the dimensions of A, $(size (A)) " ))
537- end
538- lmul! (convert (AbstractQ{T}, adjA), C)
539- end
540- function (* )(adjA:: AdjointQ{<:Any,<:LQPackedQ} , B:: AbstractMatrix )
541- A = adjA. Q
542- T = promote_type (eltype (A), eltype (B))
543- if size (B,1 ) == size (A. factors,2 )
544- C = copy_similar (B, T)
545- elseif size (B,1 ) == size (A. factors,1 )
546- C = [B; zeros (T, size (A. factors, 2 ) - size (A. factors, 1 ), size (B, 2 ))]
547- else
548- throw (DimensionMismatch (" first dimension of B, $(size (B,1 )) , must equal one of the dimensions of A, $(size (A)) " ))
549- end
550- lmul! (convert (AbstractQ{T}, adjA), C)
551- end
547+ qsize_check (adjQ:: AdjointQ{<:Any,<:LQPackedQ} , B:: AbstractVecOrMat ) =
548+ size (B, 1 ) in size (adjQ. Q. factors) ||
549+ throw (DimensionMismatch (" first dimension of B, $(size (B,1 )) , must equal one of the dimensions of Q, $(size (adjQ. Q. factors)) " ))
550+ qsize_check (A:: AbstractVecOrMat , Q:: LQPackedQ ) =
551+ size (A, 2 ) in size (Q. factors) ||
552+ throw (DimensionMismatch (" second dimension of A, $(size (A,2 )) , must equal one of the dimensions of Q, $(size (Q. factors)) " ))
552553
553554# in-place right-application of LQPackedQs
554555# these methods require that the applied-to matrix's (A's) number of columns
555556# match the number of columns (nQ) of the LQPackedQ (Q) (necessary for in-place
556557# operation, and the underlying LAPACK routine (ormlq) treats the implicit Q
557558# as its (nQ-by-nQ) square form)
558- rmul! (A:: StridedMatrix {T} , B:: LQPackedQ{T} ) where {T<: BlasFloat } =
559+ rmul! (A:: StridedVecOrMat {T} , B:: LQPackedQ{T} ) where {T<: BlasFloat } =
559560 LAPACK. ormlq! (' R' , ' N' , B. factors, B. τ, A)
560- rmul! (A:: StridedMatrix {T} , adjB:: AdjointQ{<:Any,<:LQPackedQ{T}} ) where {T<: BlasReal } =
561+ rmul! (A:: StridedVecOrMat {T} , adjB:: AdjointQ{<:Any,<:LQPackedQ{T}} ) where {T<: BlasReal } =
561562 (B = adjB. Q; LAPACK. ormlq! (' R' , ' T' , B. factors, B. τ, A))
562- rmul! (A:: StridedMatrix {T} , adjB:: AdjointQ{<:Any,<:LQPackedQ{T}} ) where {T<: BlasComplex } =
563+ rmul! (A:: StridedVecOrMat {T} , adjB:: AdjointQ{<:Any,<:LQPackedQ{T}} ) where {T<: BlasComplex } =
563564 (B = adjB. Q; LAPACK. ormlq! (' R' , ' C' , B. factors, B. τ, A))
564565
565- # out-of-place right application of LQPackedQs
566- #
567- # these methods: (1) check whether the applied-to matrix's (A's) appropriate dimension
568- # (columns for A_*, rows for Ac_*) matches the number of columns (nQ) of the LQPackedQ (Q),
569- # and if so effectively apply Q's square form to A without additional shenanigans; and
570- # (2) if the preceding dimensions do not match, check whether the appropriate dimension of
571- # A instead matches the number of rows of the matrix of which Q is a factor (i.e.
572- # size(Q.factors, 1)), and if so implicitly apply Q's truncated form to A by zero extending
573- # A as necessary for check (1) to pass (if possible) and then applying Q's square form
574- #
575- function (* )(A:: AbstractVector , Q:: LQPackedQ )
576- T = promote_type (eltype (A), eltype (Q))
577- if 1 == size (Q. factors, 2 )
578- C = copy_similar (A, T)
579- elseif 1 == size (Q. factors, 1 )
580- C = zeros (T, length (A), size (Q. factors, 2 ))
581- copyto! (C, 1 , A, 1 , length (A))
582- else
583- _rightappdimmismatch (" columns" )
584- end
585- return rmul! (C, convert (AbstractQ{T}, Q))
586- end
587- function (* )(A:: AbstractMatrix , Q:: LQPackedQ )
588- T = promote_type (eltype (A), eltype (Q))
589- if size (A, 2 ) == size (Q. factors, 2 )
590- C = copy_similar (A, T)
591- elseif size (A, 2 ) == size (Q. factors, 1 )
592- C = zeros (T, size (A, 1 ), size (Q. factors, 2 ))
593- copyto! (C, 1 , A, 1 , length (A))
594- else
595- _rightappdimmismatch (" columns" )
596- end
597- return rmul! (C, convert (AbstractQ{T}, Q))
598- end
599- function (* )(adjA:: AdjointAbsMat , Q:: LQPackedQ )
600- A = adjA. parent
601- T = promote_type (eltype (A), eltype (Q))
602- if size (A, 1 ) == size (Q. factors, 2 )
603- C = copy_similar (adjA, T)
604- elseif size (A, 1 ) == size (Q. factors, 1 )
605- C = zeros (T, size (A, 2 ), size (Q. factors, 2 ))
606- adjoint! (view (C, :, 1 : size (A, 1 )), A)
607- else
608- _rightappdimmismatch (" rows" )
609- end
610- return rmul! (C, convert (AbstractQ{T}, Q))
611- end
612- (* )(u:: AdjointAbsVec , Q:: LQPackedQ ) = (Q' u' )'
613-
614- _rightappdimmismatch (rowsorcols) =
615- throw (DimensionMismatch (string (" the number of $(rowsorcols) of the matrix on the left " ,
616- " must match either (1) the number of columns of the (LQPackedQ) matrix on the right " ,
617- " or (2) the number of rows of that (LQPackedQ) matrix's internal representation " ,
618- " (the factorization's originating matrix's number of rows)" )))
566+ # ## QB / QcB
567+ lmul! (A:: LQPackedQ{T} , B:: StridedVecOrMat{T} ) where {T<: BlasFloat } = LAPACK. ormlq! (' L' ,' N' ,A. factors,A. τ,B)
568+ lmul! (adjA:: AdjointQ{<:Any,<:LQPackedQ{T}} , B:: StridedVecOrMat{T} ) where {T<: BlasReal } =
569+ (A = adjA. Q; LAPACK. ormlq! (' L' , ' T' , A. factors, A. τ, B))
570+ lmul! (adjA:: AdjointQ{<:Any,<:LQPackedQ{T}} , B:: StridedVecOrMat{T} ) where {T<: BlasComplex } =
571+ (A = adjA. Q; LAPACK. ormlq! (' L' , ' C' , A. factors, A. τ, B))
619572
620573# In LQ factorization, `Q` is expressed as the product of the adjoint of the
621574# reflectors. Thus, `det` has to be conjugated.
0 commit comments