Merge pull request #2146 from andreasnoackjensen/master

Change dot and the Faddeeva functions to use new complex return possibility and remove the Faddeeva wrapper.

Author: ViralBShah

Makefile

@@ -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
 

base/blas.jl

@@ -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),

base/math.jl

@@ -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

base/matmul.jl

@@ -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
 

deps/Faddeeva_wrapper.c

@@ -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 ##