From dadd6ca842d345f89e015a83a1a091294e9c65ad Mon Sep 17 00:00:00 2001 From: Andreas Noack Date: Sat, 11 Oct 2014 20:05:06 -0400 Subject: [PATCH] Add function for machine precision number to avoid expensive calculation in the Givens rotation algorithms. Remove subtype relation to AbstractMatrix and therefore also the now unnecessary ambiguity warning avoiding method definitions. Use @simd and @inbounds in methods for applying Givens rotations. --- NEWS.md | 2 ++ base/linalg/givens.jl | 46 +++++++++++++++++++++---------------------- 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/NEWS.md b/NEWS.md index 5d658a27bec7b..565f01e0d5aff 100644 --- a/NEWS.md +++ b/NEWS.md @@ -75,6 +75,8 @@ Library improvements * `deepcopy` recurses through immutable types and makes copies of their mutable fields ([#8560]). + * Givens type doesn't have a size anymore and is no longer a subtype of AbstractMatrix ([#8660]) + Deprecated or removed --------------------- diff --git a/base/linalg/givens.jl b/base/linalg/givens.jl index 4da127377d26c..4ad93cd3a13d5 100644 --- a/base/linalg/givens.jl +++ b/base/linalg/givens.jl @@ -1,5 +1,4 @@ -immutable Givens{T} <: AbstractMatrix{T} - size::Int +immutable Givens{T} i1::Int i2::Int c::T @@ -11,6 +10,10 @@ type Rotation{T} end typealias AbstractRotation Union(Givens, Rotation) +realmin2(::Type{Float32}) = reinterpret(Float32, 0x26000000) +realmin2(::Type{Float64}) = reinterpret(Float64, 0x21a0000000000000) +realmin2{T}(::Type{T}) = (twopar = 2one(T); twopar^itrunc(log(realmin(T)/eps(T))/log(twopar)/twopar)) + function givensAlgorithm{T<:FloatingPoint}(f::T, g::T) zeropar = zero(T) onepar = one(T) @@ -18,7 +21,7 @@ function givensAlgorithm{T<:FloatingPoint}(f::T, g::T) safmin = realmin(T) epspar = eps(T) - safmn2 = twopar^itrunc(log(safmin/epspar)/log(twopar)/twopar) + safmn2 = realmin2(T) safmx2 = onepar/safmn2 if g == 0 @@ -84,7 +87,7 @@ function givensAlgorithm{T<:FloatingPoint}(f::Complex{T}, g::Complex{T}) abs1(ff) = max(abs(real(ff)), abs(imag(ff))) safmin = realmin(T) epspar = eps(T) - safmn2 = twopar^itrunc(log(safmin/epspar)/log(twopar)/twopar) + safmn2 = realmin2(T) safmx2 = onepar/safmn2 scalepar = max(abs1(f), abs1(g)) fs = f @@ -181,33 +184,25 @@ function givensAlgorithm{T<:FloatingPoint}(f::Complex{T}, g::Complex{T}) return cs, sn, r end -function givens{T}(f::T, g::T, i1::Integer, i2::Integer, size::Integer) - i2 <= size || error("indices cannot be larger than size Givens rotation matrix") +function givens{T}(f::T, g::T, i1::Integer, i2::Integer) i1 < i2 || error("second index must be larger than the first") - h = givensAlgorithm(f, g) - Givens(size, i1, i2, convert(T, h[1]), convert(T, h[2]), convert(T, h[3])) + c, s, r = givensAlgorithm(f, g) + Givens(i1, i2, convert(T, c), convert(T, s), convert(T, r)) end function givens{T}(A::AbstractMatrix{T}, i1::Integer, i2::Integer, col::Integer) i1 < i2 || error("second index must be larger than the first") - h = givensAlgorithm(A[i1,col], A[i2,col]) - Givens(size(A, 1), i1, i2, convert(T, h[1]), convert(T, h[2]), convert(T, h[3])) + c, s, r = givensAlgorithm(A[i1,col], A[i2,col]) + Givens(i1, i2, convert(T, c), convert(T, s), convert(T, r)) end -*{T}(G1::Givens{T}, G2::Givens{T}) = Rotation(push!(push!(Givens{T}[], G2), G1)) -*(G::Givens, B::BitArray{2}) = error("method not defined") -*{TBf,TBi}(G::Givens, B::SparseMatrixCSC{TBf,TBi}) = error("method not defined") -*(R::AbstractRotation, A::AbstractMatrix) = A_mul_B!(R, copy(A)) - -A_mul_Bc(A::AbstractMatrix, R::Rotation) = A_mul_Bc!(copy(A), R) - getindex(G::Givens, i::Integer, j::Integer) = i == j ? (i == G.i1 || i == G.i2 ? G.c : one(G.c)) : (i == G.i1 && j == G.i2 ? G.s : (i == G.i2 && j == G.i1 ? -G.s : zero(G.s))) -size(G::Givens) = (G.size, G.size) -size(G::Givens, i::Integer) = 1 <= i <= 2 ? G.size : 1 + A_mul_B!(G1::Givens, G2::Givens) = error("Operation not supported. Consider *") function A_mul_B!(G::Givens, A::AbstractMatrix) m, n = size(A) - for i = 1:n + G.i2 <= m || throw(DimensionMismatch("column indices for rotation are outside the matrix")) + @inbounds @simd for i = 1:n tmp = G.c*A[G.i1,i] + G.s*A[G.i2,i] A[G.i2,i] = G.c*A[G.i2,i] - conj(G.s)*A[G.i1,i] A[G.i1,i] = tmp @@ -216,7 +211,8 @@ function A_mul_B!(G::Givens, A::AbstractMatrix) end function A_mul_Bc!(A::AbstractMatrix, G::Givens) m, n = size(A) - for i = 1:m + G.i2 <= n || throw(DimensionMismatch("column indices for rotation are outside the matrix")) + @inbounds @simd for i = 1:m tmp = G.c*A[i,G.i1] + conj(G.s)*A[i,G.i2] A[i,G.i2] = G.c*A[i,G.i2] - G.s*A[i,G.i1] A[i,G.i1] = tmp @@ -228,14 +224,18 @@ function A_mul_B!(G::Givens, R::Rotation) return R end function A_mul_B!(R::Rotation, A::AbstractMatrix) - for i = 1:length(R.rotations) + @inbounds for i = 1:length(R.rotations) A_mul_B!(R.rotations[i], A) end return A end function A_mul_Bc!(A::AbstractMatrix, R::Rotation) - for i = 1:length(R.rotations) + @inbounds for i = 1:length(R.rotations) A_mul_Bc!(A, R.rotations[i]) end return A end +*{T}(G1::Givens{T}, G2::Givens{T}) = Rotation(push!(push!(Givens{T}[], G2), G1)) +*(R::AbstractRotation, A::AbstractMatrix) = A_mul_B!(R, copy(A)) + +A_mul_Bc(A::AbstractMatrix, R::AbstractRotation) = A_mul_Bc!(copy(A), R)