Permalink
Browse files

Change dot and the Faddeeva functions to use new complex return possi…

…bility and remove the Faddeeva wrapper.
  • Loading branch information...
1 parent ab3440a commit b962474fc0de35128069861ddbcc9d15a2c8310c @andreasnoack andreasnoack committed Jan 27, 2013
Showing with 38 additions and 92 deletions.
  1. +1 −1 Makefile
  2. +27 −4 base/blas.jl
  3. +2 −12 base/math.jl
  4. +5 −10 base/matmul.jl
  5. +0 −45 deps/Faddeeva_wrapper.c
  6. +3 −20 deps/Makefile
View
@@ -48,7 +48,7 @@ JL_LIBS = julia-release julia-debug
# private libraries, that are installed in $(PREFIX)/lib/julia
JL_PRIVATE_LIBS = amd arpack cholmod colamd fftw3 fftw3f fftw3_threads \
fftw3f_threads glpk glpk_wrapper gmp gmp_wrapper grisu \
- history Faddeeva_wrapper openlibm openlibm-extras pcre \
+ history openlibm openlibm-extras pcre \
random readline Rmath spqr suitesparse_wrapper \
tk_wrapper umfpack z openblas
View
@@ -72,19 +72,42 @@ for (fname, elty) in ((:dscal_,:Float64), (:sscal_,:Float32),
end
end
-# ccall is unable to return complex values (Issue #85)
-#@blas_dot :zdotc_ Complex128
-#@blas_dot :cdotc_ Complex64
-# DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
+# dot
for (fname, elty) in ((:ddot_,:Float64), (:sdot_,:Float32))
@eval begin
+# DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
+# * .. Scalar Arguments ..
+# INTEGER INCX,INCY,N
+# * ..
+# * .. Array Arguments ..
+# DOUBLE PRECISION DX(*),DY(*)
function dot(n::Integer, DX::Union(Ptr{$elty},Array{$elty}), incx::Integer, DY::Union(Ptr{$elty},Array{$elty}), incy::Integer)
ccall(($(string(fname)),libblas), $elty,
(Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
&n, DX, &incx, DY, &incy)
end
end
end
+for (fname, elty, relty) in ((:zdotc_,:Complex128,:Float64), (:cdotc_,:Complex64,:Float32))
+ @eval begin
+# DOUBLE COMPLEX FUNCTION ZDOTC(N,ZX,INCX,ZY,INCY)
+# * .. Scalar Arguments ..
+# INTEGER INCX,INCY,N
+# * ..
+# * .. Array Arguments ..
+# DOUBLE COMPLEX ZX(*),ZY(*)
+ function dot(n::Integer, DX::Union(Ptr{$elty},Array{$elty}), incx::Integer, DY::Union(Ptr{$elty},Array{$elty}), incy::Integer)
+ convert($elty, ccall(($(string(fname)),libblas), ComplexPair{$relty},
+ (Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
+ &n, DX, &incx, DY, &incy))
+ end
+ end
+end
+function dot{T<:BlasFloat}(DX::Array{T}, DY::Array{T})
+ n = length(DX)
+ if n != length(DY) throw(BlasDimMisMatch) end
+ return dot(n, DX, 1, DY, 1)
+end
# DOUBLE PRECISION FUNCTION DNRM2(N,X,INCX)
for (fname, elty, ret_type) in ((:dnrm2_,:Float64,:Float64),
View
@@ -667,21 +667,11 @@ function zeta(z::Number)
end
@vectorize_1arg Number zeta
-const Faddeeva_tmp = Array(Float64,2)
-
-# wrappers for complex Faddeeva functions; these will get a lot simpler,
-# and can call openopenlibm_extras directly, once ccall supports C99 complex types.
for f in (:erf, :erfc, :erfcx, :erfi, :Dawson)
fname = (f === :Dawson) ? :dawson : f
@eval begin
- function ($fname)(z::Complex128)
- ccall(($(string("wrapFaddeeva_",f)),:libFaddeeva_wrapper), Void, (Ptr{Complex128},Ptr{Complex128},Float64,), Faddeeva_tmp, &z, zero(Float64))
- return complex128(Faddeeva_tmp[1],Faddeeva_tmp[2])
- end
- function ($fname)(z::Complex64)
- ccall(($(string("wrapFaddeeva_",f)),:libFaddeeva_wrapper), Void, (Ptr{Complex128},Ptr{Complex128},Float64,), Faddeeva_tmp, &complex128(z), float64(eps(Float32)))
- return complex64(Faddeeva_tmp[1],Faddeeva_tmp[2])
- end
+ ($fname)(z::Complex128) = complex128(ccall(($(string("Faddeeva_",f)),openlibm_extras), ComplexPair{Float64}, (ComplexPair{Float64}, Float64), z, zero(Float64)))
+ ($fname)(z::Complex64) = complex64(ccall(($(string("Faddeeva_",f)),openlibm_extras), ComplexPair{Float64}, (ComplexPair{Float64}, Float64), complex128(z), float64(eps(Float32))))
($fname)(z::Complex) = ($fname)(complex128(z))
end
end
View
@@ -36,29 +36,24 @@ diagmm(b::Vector, A::Matrix) =
# Dot products
-function dot{T<:Union(Vector{Float64}, Vector{Float32})}(x::T, y::T)
- length(x) != length(y) ? error("Inputs should be of same length") : true
- BLAS.dot(length(x), x, 1, y, 1)
-end
-
-function dot{T<:Union(Float64, Float32), TI<:Integer}(x::Vector{T}, rx::Union(Range1{TI},Range{TI}), y::Vector{T}, ry::Union(Range1{TI},Range{TI}))
+dot{T<:BLAS.BlasFloat}(x::Vector{T}, y::Vector{T}) = BLAS.dot(x, y)
+function dot{T<:BLAS.BlasFloat, TI<:Integer}(x::Vector{T}, rx::Union(Range1{TI},Range{TI}), y::Vector{T}, ry::Union(Range1{TI},Range{TI}))
length(rx) != length(ry) ? error("Ranges should be of same length") : true
if min(rx) < 1 || max(rx) > length(x) || min(ry) < 1 || max(ry) > length(y)
throw(BoundsError())
end
BLAS.dot(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
end
-
-Ac_mul_B(x::Vector, y::Vector) = [dot(x, y)]
-At_mul_B{T<:Real}(x::Vector{T}, y::Vector{T}) = [dot(x, y)]
-
function dot(x::Vector, y::Vector)
s = zero(eltype(x))
for i=1:length(x)
s += conj(x[i])*y[i]
end
s
end
+dot(x::Number, y::Number) = conj(x) * y
+Ac_mul_B(x::Vector, y::Vector) = [dot(x, y)]
+At_mul_B{T<:Real}(x::Vector{T}, y::Vector{T}) = [dot(x, y)]
dot(x::Number, y::Number) = conj(x) * y
@@ -1,45 +0,0 @@
-/* Copyright (c) 2012 Massachusetts Institute of Technology
- *
- * Permission is hereby granted, free of charge, to any person obtaining
- * a copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject to
- * the following conditions:
- *
- * The above copyright notice and this permission notice shall be
- * included in all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
- * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
- * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
- * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
- * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
- * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
- * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
- */
-
-/* C89-compatible wrappers for complex Faddeeva functions, that
- pass and return complex values via a double* which points to
- the real and imaginary parts consecutively. (Note that this
- is binary-compatible with both C99 complex numbers and with
- C++ std::complex<double>.) */
-
-#include "Faddeeva.h"
-
-#define WRAP(func) \
-void wrap ## Faddeeva_ ## func(double *func, const double *z, double relerr) \
-{ \
- const double complex *zc = (const double complex *) z; \
- double complex c = Faddeeva_ ## func(*zc, relerr); \
- func[0] = creal(c); \
- func[1] = cimag(c); \
-}
-
-WRAP(w)
-WRAP(erfcx)
-WRAP(erf)
-WRAP(erfi)
-WRAP(erfc)
-WRAP(Dawson)
View
@@ -20,13 +20,13 @@ endif
#autoconf configure-driven scripts: llvm readline pcre arpack fftw unwind gmp glpk patchelf
#custom configure-driven script: zlib nginx
-#custom Makefile rules: openlibm Faddeeva-wrapper rmath double-conversion random gmp-wrapper suitesparse-wrapper suitesparse lapack openblas uv tk-wrapper
+#custom Makefile rules: openlibm rmath double-conversion random gmp-wrapper suitesparse-wrapper suitesparse lapack openblas uv tk-wrapper
# prevent installing libs into usr/lib64 on opensuse
unexport CONFIG_SITE
STAGE1_DEPS = uv openlibm-extras random rmath double-conversion
-STAGE2_DEPS = gmp-wrapper Faddeeva-wrapper
+STAGE2_DEPS = gmp-wrapper
STAGE3_DEPS = glpk-wrapper suitesparse-wrapper
ifeq ($(USE_SYSTEM_LIBUNWIND), 0)
@@ -551,24 +551,7 @@ get-openlibm-extras: openlibm/Makefile.extras
configure-openlibm-extras: get-openlibm-extras
compile-openlibm-extras: $(OPENLIBMEXT_OBJ_SOURCE)
check-openlibm-extras: compile-openlibm-extras
-install-openlibm-extras: $(OPENLIBMEXT_OBJ_TARGET) $(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT)
-
-# Wrapper for openlibm/Faddeeva since ccall doesn't support C99 complex:
-$(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT): Faddeeva_wrapper.c $(OPENLIBMEXT_OBJ_TARGET)
- $(CC) $(CPPFLAGS) $(CFLAGS) $(LDFLAGS) -O2 -shared $(fPIC) -I openlibm/Faddeeva Faddeeva_wrapper.c -o $@ -L $(BUILD)/lib -lopenlibm-extras $(RPATH_ORIGIN)
- $(INSTALL_NAME_CMD)libFaddeeva_wrapper.$(SHLIB_EXT) $@
- touch -c $@
-
-clean-Faddeeva-wrapper:
- -rm -f $(OPENLIBMEXT_OBJ_TARGET) $(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT)
-distclean-Faddeeva-wrapper: clean-Faddeeva-wrapper
-
-get-Faddeeva-wrapper: get-openlibm
-configure-Faddeeva-wrapper: get-openlibm
-compile-Faddeeva-wrapper: $(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT)
-check-Faddeeva-wrapper: compile-Faddeeva-wrapper
-install-Faddeeva-wrapper: $(BUILD)/lib/libFaddeeva_wrapper.$(SHLIB_EXT)
-
+install-openlibm-extras: $(OPENLIBMEXT_OBJ_TARGET)
## LIBRANDOM ##

0 comments on commit b962474

Please sign in to comment.