Merge pull request JuliaLang#1324 from JuliaLang/vs/linalg

WIP: Linear Algebra reorganization

base/array.jl

@@ -1334,6 +1334,30 @@ function nonzeros{T}(A::StridedArray{T})
     return V
 end
 
+function findmax(a::Array)
+    m = typemin(eltype(a))
+    mi = 0
+    for i=1:length(a)
+        if a[i] > m
+            m = a[i]
+            mi = i
+        end
+    end
+    return (m, mi)
+end
+
+function findmin(a::Array)
+    m = typemax(eltype(a))
+    mi = 0
+    for i=1:length(a)
+        if a[i] < m
+            m = a[i]
+            mi = i
+        end
+    end
+    return (m, mi)
+end
+
 ## Reductions ##
 
 areduce{T}(f::Function, A::StridedArray{T}, region::Dimspec, v0) =

base/blas.jl

@@ -1,4 +1,4 @@
-typealias LapackScalar Union(Float64,Float32,Complex128,Complex64)
+typealias LapackType Union(Float64,Float32,Complex128,Complex64)
 
 module Blas
 import Base.*
@@ -54,17 +54,38 @@ end
 for (fname, elty) in ((:daxpy_,:Float64), (:saxpy_,:Float32),
                       (:zaxpy_,:Complex128), (:caxpy_,:Complex64))
     @eval begin
-        function axpy!(n::Integer, a::($elty),
-                       DX::Union(Ptr{$elty},Array{$elty}), incx::Integer,
-                       DY::Union(Ptr{$elty},Array{$elty}), incy::Integer)
+        function axpy!(n::Integer, alpha::($elty),
+                       dx::Union(Ptr{$elty},Array{$elty}), incx::Integer,
+                       dy::Union(Ptr{$elty},Array{$elty}), incy::Integer)
             ccall(dlsym(Base.libblas, $(string(fname))), Void,
                   (Ptr{Int32}, Ptr{$elty}, Ptr{$elty}, Ptr{Int32}, Ptr{$elty}, Ptr{Int32}),
-                  &n, &a, DX, &incx, DY, &incy)
-            DY
+                  &n, &alpha, dx, &incx, dy, &incy)
+            return dy
         end
     end
 end
 
+function axpy!{T,Ta<:Number}(alpha::Ta, x::Array{T}, y::Array{T})
+    if length(x) != length(y)
+        error("Inputs should be of the same length")
+    end
+    return axpy!(length(x), convert(T, alpha), x, 1, y, 1)
+end
+
+function axpy!{T,Ta<:Number,Ti<:Integer}(alpha::Ta, x::Array{T}, rx::Union(Range1{Ti},Range{Ti}), y::Array{T}, ry::Union(Range1{Ti},Range{Ti}))
+
+    if length(rx) != length(ry)
+        error("Ranges should be of the same length")
+    end
+
+    if min(rx) < 1 || max(rx) > length(x) || min(ry) < 1 || max(ry) > length(y)
+        throw(BoundsError())
+    end
+
+    return axpy!(length(rx), convert(T, alpha), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
+end
+
+
 # SUBROUTINE DSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
 # *     .. Scalar Arguments ..
 #       REAL ALPHA,BETA
@@ -286,3 +307,34 @@ for (fname, elty) in ((:dgemv_,:Float64), (:sgemv_,:Float32),
 end
 
 end # module
+
+# Use BLAS copy for small arrays where it is faster than memcpy, and for strided copying
+
+function copy_to{T<:LapackType}(dest::Ptr{T}, src::Ptr{T}, n::Integer)
+    if n < 200
+        Blas.copy!(n, src, 1, dest, 1)
+    else
+        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Uint), dest, src, n*sizeof(T))
+    end
+    return dest
+end
+
+function copy_to{T<:LapackType}(dest::Array{T}, src::Array{T})
+    n = numel(src)
+    if n > numel(dest); throw(BoundsError()); end
+    copy_to(pointer(dest), pointer(src), n)
+    return dest
+end
+
+function copy_to{T<:LapackType,Ti<:Integer}(dest::Array{T}, rdest::Union(Range1{Ti},Range{Ti}), 
+                                            src::Array{T}, rsrc::Union(Range1{Ti},Range{Ti}))
+    if min(rdest) < 1 || max(rdest) > length(dest) || min(rsrc) < 1 || max(rsrc) > length(src)
+        throw(BoundsError())
+    end
+    if length(rdest) != length(rsrc)
+        error("Ranges must be of the same length")
+    end
+    Blas.copy!(length(rsrc), pointer(src)+(first(rsrc)-1)*sizeof(T), step(rsrc),
+              pointer(dest)+(first(rdest)-1)*sizeof(T), step(rdest))
+    return dest
+end

base/client.jl

@@ -198,11 +198,11 @@ function _start()
     global const stderr_stream = make_stderr_stream()
     global OUTPUT_STREAM = stdout_stream
 
+    librandom_init()
+
     # set CPU core count
     global const CPU_CORES = ccall(:jl_cpu_cores, Int32, ())
 
-    _jl_librandom_init()
-
     atexit(()->flush(stdout_stream))
     try
         ccall(:jl_register_toplevel_eh, Void, ())

base/export.jl

@@ -8,6 +8,8 @@ export
     DSP,
     Lapack,
     Blas,
+    LibRandom,
+    RNG,
 
 # Types
     AbstractMatrix,
@@ -868,6 +870,9 @@ export
     quantile,
     quartile,
     quintile,
+    librandom_init,
+    Rng_MT,
+    Rng,
     rand,
     rand!,
     randbeta,