1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187
|
# This file is a part of Julia. License is MIT: https://julialang.org/license
module TestLQ
using Test, LinearAlgebra, Random
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, rmul!, lmul!
n = 10
# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1
Random.seed!(1234321)
areal = randn(n,n)/2
aimg = randn(n,n)/2
a2real = randn(n,n)/2
a2img = randn(n,n)/2
breal = randn(n,2)/2
bimg = randn(n,2)/2
# helper functions to unambiguously recover explicit forms of an LQPackedQ
squareQ(Q::LinearAlgebra.LQPackedQ) = (n = size(Q.factors, 2); lmul!(Q, Matrix{eltype(Q)}(I, n, n)))
rectangularQ(Q::LinearAlgebra.LQPackedQ) = convert(Array, Q)
@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64)
a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
asym = a' + a # symmetric indefinite
apd = a' * a # symmetric positive-definite
ε = εa = eps(abs(float(one(eltya))))
@testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int)
b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
εb = eps(abs(float(one(eltyb))))
ε = max(εa,εb)
α = rand(eltya)
aα = fill(α,1,1)
@test lq(α).L*lq(α).Q ≈ lq(aα).L*lq(aα).Q
@test abs(lq(α).Q[1,1]) ≈ one(eltya)
tab = promote_type(eltya,eltyb)
for i = 1:2
let a = i == 1 ? a : view(a, 1:n - 1, 1:n - 1), b = i == 1 ? b : view(b, 1:n - 1), n = i == 1 ? n : n - 1
lqa = lq(a)
l,q = lqa.L, lqa.Q
qra = qr(a)
@testset "Basic ops" begin
@test size(lqa,1) == size(a,1)
@test size(lqa,3) == 1
@test size(lqa.Q,3) == 1
@test_throws ErrorException lqa.Z
@test Array(copy(adjoint(lqa))) ≈ a'
@test q*squareQ(q)' ≈ Matrix(I, n, n)
@test l*q ≈ a
@test Array(lqa) ≈ a
@test Array(copy(lqa)) ≈ a
lstring = sprint(show,l)
qstring = sprint(show,q)
@test sprint(show,lqa) == "$(typeof(lqa)) with factors L and Q:\n$lstring\n$qstring"
end
@testset "Binary ops" begin
@test a*(lqa\b) ≈ b atol=3000ε
@test lqa*b ≈ qra.Q*qra.R*b atol=3000ε
@test (sq = size(q.factors, 2); *(Matrix{eltyb}(I, sq, sq), adjoint(q))*squareQ(q)) ≈ Matrix(I, n, n) atol=5000ε
if eltya != Int
@test Matrix{eltyb}(I, n, n)*q ≈ convert(AbstractMatrix{tab},q)
end
@test q*b ≈ squareQ(q)*b atol=100ε
@test transpose(q)*b ≈ transpose(squareQ(q))*b atol=100ε
@test q'*b ≈ squareQ(q)'*b atol=100ε
@test a*q ≈ a*squareQ(q) atol=100ε
@test a*transpose(q) ≈ a*transpose(squareQ(q)) atol=100ε
@test a*q' ≈ a*squareQ(q)' atol=100ε
@test a'*q ≈ a'*squareQ(q) atol=100ε
@test a'*q' ≈ a'*squareQ(q)' atol=100ε
@test_throws DimensionMismatch q*b[1:n1 + 1]
@test_throws DimensionMismatch adjoint(q) * Matrix{eltya}(undef,n+2,n+2)
@test_throws DimensionMismatch Matrix{eltyb}(undef,n+2,n+2)*q
if isa(a, DenseArray) && isa(b, DenseArray)
# use this to test 2nd branch in mult code
pad_a = vcat(I, a)
pad_b = hcat(I, b)
@test pad_a*q ≈ pad_a*squareQ(q) atol=100ε
@test transpose(q)*pad_b ≈ transpose(squareQ(q))*pad_b atol=100ε
@test q'*pad_b ≈ squareQ(q)'*pad_b atol=100ε
end
end
end
end
@testset "Matmul with LQ factorizations" begin
lqa = lq(a[:,1:n1])
l,q = lqa.L, lqa.Q
@test rectangularQ(q)*rectangularQ(q)' ≈ Matrix(I, n1, n1)
@test squareQ(q)'*squareQ(q) ≈ Matrix(I, n1, n1)
@test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1),q)
@test lmul!(adjoint(q), rectangularQ(q)) ≈ Matrix(I, n1, n1)
@test_throws DimensionMismatch rmul!(Matrix{eltya}(I, n+1, n+1), adjoint(q))
@test_throws BoundsError size(q,-1)
end
end
end
@testset "getindex on LQPackedQ (#23733)" begin
local m, n
function getqs(F::LinearAlgebra.LQ)
implicitQ = F.Q
sq = size(implicitQ.factors, 2)
explicitQ = lmul!(implicitQ, Matrix{eltype(implicitQ)}(I, sq, sq))
return implicitQ, explicitQ
end
m, n = 3, 3 # reduced Q 3-by-3, full Q 3-by-3
implicitQ, explicitQ = getqs(lq(randn(m, n)))
@test implicitQ[1, 1] == explicitQ[1, 1]
@test implicitQ[m, 1] == explicitQ[m, 1]
@test implicitQ[1, n] == explicitQ[1, n]
@test implicitQ[m, n] == explicitQ[m, n]
m, n = 3, 4 # reduced Q 3-by-4, full Q 4-by-4
implicitQ, explicitQ = getqs(lq(randn(m, n)))
@test implicitQ[1, 1] == explicitQ[1, 1]
@test implicitQ[m, 1] == explicitQ[m, 1]
@test implicitQ[1, n] == explicitQ[1, n]
@test implicitQ[m, n] == explicitQ[m, n]
@test implicitQ[m+1, 1] == explicitQ[m+1, 1]
@test implicitQ[m+1, n] == explicitQ[m+1, n]
m, n = 4, 3 # reduced Q 3-by-3, full Q 3-by-3
implicitQ, explicitQ = getqs(lq(randn(m, n)))
@test implicitQ[1, 1] == explicitQ[1, 1]
@test implicitQ[n, 1] == explicitQ[n, 1]
@test implicitQ[1, n] == explicitQ[1, n]
@test implicitQ[n, n] == explicitQ[n, n]
end
@testset "size on LQPackedQ (#23780)" begin
# size(Q::LQPackedQ) yields the shape of Q's full/square form
for ((mA, nA), nQ) in (
((3, 3), 3), # A 3-by-3 => full/square Q 3-by-3
((3, 4), 4), # A 3-by-4 => full/square Q 4-by-4
((4, 3), 3) )# A 4-by-3 => full/square Q 3-by-3
@test size(lq(randn(mA, nA)).Q) == (nQ, nQ)
end
end
@testset "postmultiplication with / right-application of LQPackedQ (#23779)" begin
function getqs(F::LinearAlgebra.LQ)
implicitQ = F.Q
explicitQ = lmul!(implicitQ, Matrix{eltype(implicitQ)}(I, size(implicitQ)...))
return implicitQ, explicitQ
end
# for any shape m-by-n of LQ-factored matrix, where Q is an LQPackedQ
# A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) operations should work for
# *-by-n (n-by-*) C, which we test below via n-by-n C
for (mA, nA) in ((3, 3), (3, 4), (4, 3))
implicitQ, explicitQ = getqs(lq(randn(mA, nA)))
C = randn(nA, nA)
@test *(C, implicitQ) ≈ *(C, explicitQ)
@test *(C, adjoint(implicitQ)) ≈ *(C, adjoint(explicitQ))
@test *(adjoint(C), implicitQ) ≈ *(adjoint(C), explicitQ)
@test *(adjoint(C), adjoint(implicitQ)) ≈ *(adjoint(C), adjoint(explicitQ))
end
# where the LQ-factored matrix has at least as many rows m as columns n,
# Q's full/square and reduced/rectangular forms have the same shape (n-by-n). hence we expect
# _only_ *-by-n (n-by-*) C to work in A_mul_B*(C, Q) (Ac_mul_B*(C, Q)) ops.
# and hence the n-by-n C tests above suffice.
#
# where the LQ-factored matrix has more columns n than rows m,
# Q's full/square form is n-by-n whereas its reduced/rectangular form is m-by-n.
# hence we need also test *-by-m C with
# A*_mul_B(C, Q) ops, as below via m-by-m C.
mA, nA = 3, 4
implicitQ, explicitQ = getqs(lq(randn(mA, nA)))
C = randn(mA, mA)
zeroextCright = hcat(C, zeros(eltype(C), mA))
zeroextCdown = vcat(C, zeros(eltype(C), (1, mA)))
@test *(C, implicitQ) ≈ *(zeroextCright, explicitQ)
@test *(adjoint(C), implicitQ) ≈ *(adjoint(zeroextCdown), explicitQ)
@test_throws DimensionMismatch C * adjoint(implicitQ)
@test_throws DimensionMismatch adjoint(C) * adjoint(implicitQ)
end
end # module TestLQ
|