From bcc6bcc4d9b44aa6e2422abf445ebf3f82ac7cc5 Mon Sep 17 00:00:00 2001 From: Andreas Noack Jensen Date: Mon, 22 Oct 2012 22:36:31 +0200 Subject: [PATCH] Added some eps methods for complex floats. Fixed a bug in dgesld for complex types and added input argument for singularity value. Updated \ accordingly. Added some methods to handle right and left hand side of different type. Restricted methods in operators apply for Numbers only. --- base/complex.jl | 4 +++- base/lapack.jl | 53 +++++++++++++++++++++++++++++++------------- base/linalg_dense.jl | 6 +++-- base/operators.jl | 22 +++++++++--------- 4 files changed, 55 insertions(+), 30 deletions(-) diff --git a/base/complex.jl b/base/complex.jl index 1e3ca874661f5..82d3b9ef8daf6 100644 --- a/base/complex.jl +++ b/base/complex.jl @@ -195,7 +195,9 @@ isequal(x::Real, z::Complex) = real_valued(z) && isequal(real(z),x) hash(z::Complex) = (r = hash(real(z)); real_valued(z) ? r : bitmix(r,hash(imag(z)))) -eps(z::Complex) = eps(abs(z)) +eps(z::Complex) = eps(real(z)) +eps(::Type{Complex64}) = eps(Float32) +eps(::Type{Complex128}) = eps(Float64) conj(z::Complex) = complex(real(z),-imag(z)) abs(z::Complex) = hypot(real(z), imag(z)) diff --git a/base/lapack.jl b/base/lapack.jl index 00927d3dc067e..e0e4de7eb025f 100644 --- a/base/lapack.jl +++ b/base/lapack.jl @@ -383,30 +383,51 @@ for (gels, gelsd, gesv, getrs, getri, elty) in # * .. Array Arguments .. # INTEGER IWORK( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) - function gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::FloatingPoint) Lapack.chkstride1(A, B) - m, n = size(A) + m, n = size(A) if size(B,1) != m; throw(Lapack.LapackDimMisMatch("gelsd!")); end - s = Array($elty, min(m, n)) - rcond = eps(Float32) - rnk = Array(Int32, 1) - info = Array(Int32, 1) - work = Array($elty, 1) - lwork = int32(-1) - iwork = Array(Int32, 1) + Rtyp = typeof(real(A[1])) + s = Array(Rtyp, min(m, n)) + cmplx = iscomplex(A) + if cmplx + rwork = Array(Rtyp, 1) + end + rnk = Array(Int32, 1) + info = Array(Int32, 1) + work = Array($elty, 1) + lwork = int32(-1) + iwork = Array(Int32, 1) for i in 1:2 - ccall(dlsym(Base.liblapack, $(string(gelsd))), Void, - (Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, - Ptr{$elty}, Ptr{Int32}, Ptr{$elty}, Ptr{Int32}, - Ptr{$elty}, Ptr{$elty}, Ptr{Int32}, Ptr{$elty}, - Ptr{Int32}, Ptr{Int32}, Ptr{Int32}), - &m, &n, &size(B,2), A, &max(1,stride(A,2)), - B, &max(1,stride(B,2)), s, &rcond, rnk, work, &lwork, iwork, info) + if cmplx + ccall(dlsym(Base.liblapack, $(string(gelsd))), Void, + (Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, + Ptr{Int32}, Ptr{$elty}, Ptr{Int32}, Ptr{Rtyp}, + Ptr{Rtyp}, Ptr{Int32}, Ptr{$elty}, Ptr{Int32}, + Ptr{Rtyp}, Ptr{Int32}, Ptr{Int32}), + &m, &n, &size(B,2), A, + &max(1,stride(A,2)), B, &max(1,stride(B,2)), s, + &rcond, rnk, work, &lwork, + rwork, iwork, info) + else + ccall(dlsym(Base.liblapack, $(string(gelsd))), Void, + (Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{$elty}, + Ptr{Int32}, Ptr{$elty}, Ptr{Int32}, Ptr{$elty}, + Ptr{$elty}, Ptr{Int32}, Ptr{$elty}, Ptr{Int32}, + Ptr{Int32}, Ptr{Int32}), + &m, &n, &size(B,2), A, + &max(1,stride(A,2)), B, &max(1,stride(B,2)), s, + &rcond, rnk, work, &lwork, + iwork, info) + end if info[1] != 0 throw(LapackException(info[1])) end if lwork < 0 lwork = int32(real(work[1])) work = Array($elty, lwork) iwork = Array(Int32, iwork[1]) + if cmplx + rwork = Array(Rtyp, int(rwork[1])) + end end end isa(B, Vector) ? B[1:n] : B[1:n,:], rnk[1] diff --git a/base/linalg_dense.jl b/base/linalg_dense.jl index 90d8295bc7336..cc4d22cdfe2b7 100644 --- a/base/linalg_dense.jl +++ b/base/linalg_dense.jl @@ -617,12 +617,14 @@ function (\){T<:LapackType}(A::StridedMatrix{T}, B::StridedVecOrMat{T}) if ishermitian(A) return Lapack.sysv!('U', Acopy, X)[1] end return Lapack.gesv!(Acopy, X)[3] end - Lapack.gelsd!(Acopy, X)[1] + Lapack.gelsd!(Acopy, X, -1.0)[1] end +(\){T1<:LapackType, T2<:Real}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(A, convert(Array{T1}, B)) +(\){T1<:Real, T2<:LapackType}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(convert(Array{T2}, A), B) (\){T1<:Real, T2<:Real}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(float64(A), float64(B)) +(\){T1<:Complex, T2<:Complex}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(convert(Array{Complex128}, A), convert(Array{Complex128}, B)) -# TODO: use *gels transpose argument (/)(A::StridedVecOrMat, B::StridedVecOrMat) = (B' \ A')' ##TODO: Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly diff --git a/base/operators.jl b/base/operators.jl index 0bddce2c06ef5..cadb524985fa1 100644 --- a/base/operators.jl +++ b/base/operators.jl @@ -14,6 +14,8 @@ isequal(x,y) = is(x,y) > (x,y) = y < x <=(x,y) = !(y < x) >=(x,y) = (y <= x) +.> (x,y) = y.=(x,y) = y.<=x # these definitions allow Number types to implement # == and < instead of isequal and isless, which is more idiomatic: @@ -52,22 +54,20 @@ for op = (:+, :*, :&, :|, :$, :min, :max) end end -\(x,y) = y/x +\(x::Number,y::Number) = y/x # . defaults to -./(x,y) = x/y -.\(x,y) = y./x -.*(x,y) = x*y -.^(x,y) = x^y +./(x::Number,y::Number) = x/y +.\(x::Number,y::Number) = y./x +.*(x::Number,y::Number) = x*y +.^(x::Number,y::Number) = x^y .+(x,y) = x+y .-(x,y) = x-y -.==(x,y) = x==y -.!=(x,y) = x!=y -.< (x,y) = x (x,y) = y.=(x,y) = y.<=x +.==(x::Number,y::Number) = x==y +.!=(x::Number,y::Number) = x!=y +.< (x::Real,y::Real) = x> and >>> takes Int32 as second arg <<(x,y::Integer) = x << convert(Int32,y)