function mysum(A)
    acc = zero(eltype(A))
    for a in A
       acc += a
    end
    acc
end

julia> @code_lowered mysum(ones(1))
CodeInfo(:(begin
        nothing
        acc = (Main.zero)((Main.eltype)(A))
        #= line 3 =#
        SSAValue(0) = A
        #temp# = (Base.start)(SSAValue(0))
        6: 
        unless !((Base.done)(SSAValue(0), #temp#)) goto 15
        SSAValue(1) = (Base.next)(SSAValue(0), #temp#)
        a = (Core.getfield)(SSAValue(1), 1)
        #temp# = (Core.getfield)(SSAValue(1), 2)
        #= line 4 =#
        acc = acc + a
        13: 
        goto 6
        15: 
        #= line 6 =#
        return acc
    end))

julia> @code_typed mysum(ones(1))
CodeInfo(:(begin
        acc = (Base.sitofp)(Float64, 0)::Float64
        #= line 3 =#
        #temp# = 1
        4: 
        unless (Base.not_int)(
            (#temp# === (Base.add_int)((Base.arraylen)(A)::Int64, 1)::Int64)::Bool
        )::Bool goto 14
        SSAValue(2) = (Base.arrayref)(true, A, #temp#)::Float64
        SSAValue(3) = (Base.add_int)(#temp#, 1)::Int64
        a = SSAValue(2)
        #temp# = SSAValue(3)
        #= line 4 =#
        acc = (Base.add_float)(acc, a)::Float64
        12: 
        goto 4
        14: 
        #= line 6 =#
        return acc
    end)) => Float64

julia> @code_warntype mysum_unstable(ones(1))
Variables:
  A::Array{Float64,1}
  acc::Union{Float64, Int64}
  #temp#@_4::Int64
  a::Float64
  #temp#@_6::Float64

Body:
  begin
      acc::Union{Float64, Int64} = 0
      #= line 3 =#
      #temp#@_4::Int64 = 1
      4: 
      unless (Base.not_int)((#temp#@_4::Int64 === (Base.add_int)((Base.arraylen)(A::Array{Float64,1})::Int64, 1)::Int64)::Bool)::Bool goto 26
      SSAValue(3) = (Base.arrayref)(true, A::Array{Float64,1}, #temp#@_4::Int64)::Float64
      SSAValue(4) = (Base.add_int)(#temp#@_4::Int64, 1)::Int64
      a::Float64 = SSAValue(3)
      #temp#@_4::Int64 = SSAValue(4)
      #= line 4 =#
      unless (acc::Union{Float64, Int64} isa Float64)::Bool goto 14
      #temp#@_6::Float64 = $(Expr(:invoke, MethodInstance for +(::Float64, ::Float64), :(Main.+), :(acc), :(a)))::Float64
      goto 22
      14: 
      unless (acc::Union{Float64, Int64} isa Int64)::Bool goto 18
      #temp#@_6::Float64 = $(Expr(:invoke, MethodInstance for +(::Int64, ::Float64), :(Main.+), :(acc), :(a)))::Float64
      goto 22
      18: 
      goto 20
      20: 
      (Base.error)("fatal error in type inference (type bound)")::Any
      22: 
      acc::Union{Float64, Int64} = #temp#@_6::Float64
      24: 
      goto 4
      26: 
      #= line 6 =#
      return acc::Union{Float64, Int64}
  end::Union{Float64, Int64}

julia> @code_llvm mysum(ones(1))

; Function mysum
; Location: REPL[52]:2
define double @julia_mysum_64348(%jl_value_t addrspace(10)* nonnull dereferenceable(40)) {
top:
; Location: REPL[52]:3
  %1 = addrspacecast %jl_value_t addrspace(10)* %0 to %jl_value_t addrspace(11)*
  %2 = bitcast %jl_value_t addrspace(11)* %1 to %jl_array_t addrspace(11)*
  %3 = getelementptr inbounds %jl_array_t, %jl_array_t addrspace(11)* %2, i64 0, i32 1
  %4 = load i64, i64 addrspace(11)* %3, align 8
  %5 = icmp eq i64 %4, 0
  br i1 %5, label %L14, label %if.lr.ph

if.lr.ph:                                         ; preds = %top
  %6 = bitcast %jl_value_t addrspace(11)* %1 to double* addrspace(11)*
  br label %if

if:                                               ; preds = %if.lr.ph, %idxend
  %"#temp#.04" = phi i64 [ 1, %if.lr.ph ], [ %14, %idxend ]
  %acc.03 = phi double [ 0.000000e+00, %if.lr.ph ], [ %15, %idxend ]
  %7 = add i64 %"#temp#.04", -1
  %8 = icmp ult i64 %7, %4
  br i1 %8, label %idxend, label %oob

L14.loopexit:                                     ; preds = %idxend
; Location: REPL[52]:6
  br label %L14

L14:                                              ; preds = %L14.loopexit, %top
  %acc.0.lcssa = phi double [ 0.000000e+00, %top ], [ %15, %L14.loopexit ]
  ret double %acc.0.lcssa

oob:                                              ; preds = %if
; Location: REPL[52]:3
  %9 = alloca i64, align 8
  store i64 %"#temp#.04", i64* %9, align 8
  %10 = addrspacecast %jl_value_t addrspace(10)* %0 to %jl_value_t addrspace(12)*
  call void @jl_bounds_error_ints(%jl_value_t addrspace(12)* %10, i64* nonnull %9, i64 1)
  unreachable

idxend:                                           ; preds = %if
  %11 = load double*, double* addrspace(11)* %6, align 8
  %12 = getelementptr double, double* %11, i64 %7
  %13 = load double, double* %12, align 8
  %14 = add i64 %"#temp#.04", 1
; Location: REPL[52]:4
  %15 = fadd double %acc.03, %13
; Location: REPL[52]:3
  %16 = icmp eq i64 %"#temp#.04", %4
  br i1 %16, label %L14.loopexit, label %if
}

julia> @code_native mysum(ones(1))
	.text
; Function mysum {
; Location: REPL[52]:2
	movq	8(%rdi), %rcx
	vxorpd	%xmm0, %xmm0, %xmm0
	testq	%rcx, %rcx
	je	L54
	vxorpd	%xmm0, %xmm0, %xmm0
	xorl	%eax, %eax
	nopw	%cs:(%rax,%rax)
L32:
	cmpq	%rcx, %rax
	jae	L55
	movq	(%rdi), %rdx
; Location: REPL[52]:4
	vaddsd	(%rdx,%rax,8), %xmm0, %xmm0
; Location: REPL[52]:3
	addq	$1, %rax
	cmpq	%rax, %rcx
	jne	L32
; Location: REPL[52]:6
L54:
	retq
L55:
	pushq	%rbp
	movq	%rsp, %rbp
; Location: REPL[52]:3
	movq	%rsp, %rcx
	leaq	-16(%rcx), %rsi
	movq	%rsi, %rsp
	addq	$1, %rax
	movq	%rax, -16(%rcx)
	movabsq	$jl_bounds_error_ints, %rax
	movl	$1, %edx
	callq	*%rax
	nop
;}

@profile fun() # Profile a specific function
Profile.clear() # Clear the recorded profile
Profile.print() # Print the recorded profile
Profile.print(C=true) # Print the profile including calls into C

# The textual output of the profiler can be hard to understand
# ProfileView.jl gives you a graphical representation
using ProfileView
ProfileView.view()

# LLVM IR uses <2 x f64> to denote vector operations

Pkg.clone("https://github.com/vchuravy/IACA.jl")

# Similar to code_llvm
IACA.code_llvm(func, (Type,))

# What has my compiler ever done for me
IACA.code_llvm(func, (Type,), optimize=false)

# What would my code do on a different platform
IACA.code_llvm(func, (Type,), cpu="knl")

# Switching optimization level
IACA.optlevel[] = 3
IACA.code_llvm(func, (Type,))

struct RGBA
    r :: Float64
    g :: Float64
    b :: Float64
    a :: Float64
end

Base.rand(::Type{RGBA}) = RGBA(rand(), rand(), rand(), rand())

function myscale!(colors::AbstractVector{RGBA}, α::Float64)
    @inbounds for i in 1:length(colors)
        c = colors[i]
        r = α * c.r
        g = α * c.g
        b = α * c.b
        a = α * c.a
        colors[i] = RGBA(r, g, b, a)
    end
end

IACA.code_llvm(myscale!, (typeof([rand(RGBA)]),Float64))

IACA.optlevel[] = 3

IACA.code_llvm(myscale!, (typeof([rand(RGBA)]),Float64))

> lscpu
...
Thread(s) per core:  2
Core(s) per socket:  2
...
Model name:          Intel(R) Core(TM) i7-7660U CPU @ 2.50GHz
...
L1d cache:           32K
L1i cache:           32K
L2 cache:            256K
L3 cache:            4096K
...
Flags: ... avx2 ...

> lscpu
...
Thread(s) per core:  2
Core(s) per socket:  12
...
Model name:          AMD Ryzen Threadripper 1920X 12-Core Processor
...
L1d cache:           32K
L1i cache:           64K
L2 cache:            512K
L3 cache:            8192K
...
Flags: ... avx2 ...

using BenchmarkTools
using Base.Threads
println("Number of threads: ", nthreads())

function myfun(rng::MersenneTwister)
    s = 0.0
    N = 10000
    for i=1:N
        s+=det(randn(rng, 3,3))
    end
    s/N
end

rgi   = [MersenneTwister(rand(UInt)) for _ in 1:nthreads()]
function bench(rgi)
    a  = zeros(1000)
    @threads for i=1:length(a)
        a[i] = myfun(rgi[threadid()])
    end
end

result = @benchmark bench($rgi)
display(result)

@edit det(zeros(3, 3)) -> det(lufact(A))

det!(A) = det(lufact!(A))
det!(A) = det(LinAlg.generic_lufact!(A))

using BenchmarkTools
using StaticArrays
using Base.Threads
println("Number of threads: ", nthreads())

function myfun(rng::MersenneTwister)
    s = 0.0
    N = 10000
    for i=1:N
        s += det(randn(rng, SMatrix{3, 3}))
    end
    s/N
end

rgi   = [MersenneTwister(abs(rand(Int))) for s in 1:nthreads()]

function bench(rgi)
    a  = zeros(1000)
    @threads for i=1:length(a)
        @inbounds a[i] = myfun(rgi[threadid()])
    end
end

result = @benchmark bench($rgi)
display(result)

BenchmarkTools.Trial: 
  memory estimate:  7.98 KiB
  allocs estimate:  2
  --------------
  minimum time:     588.229 ms (0.00% GC)
  median time:      654.259 ms (0.00% GC)
  mean time:        641.654 ms (0.00% GC)
  maximum time:     669.961 ms (0.00% GC)
  --------------
  samples:          8
  evals/sample:     1

BenchmarkTools.Trial: 
  memory estimate:  3.52 GiB
  allocs estimate:  34212074
  --------------
  minimum time:     3.346 s (0.00% GC)
  median time:      3.960 s (10.85% GC)
  mean time:        3.960 s (10.85% GC)
  maximum time:     4.574 s (18.78% GC)
  --------------
  samples:          2
  evals/sample:     1

numactl --hardware
available: 1 nodes (0)
node 0 cpus: 0 1 2 3
node 0 size: 15934 MB
node 0 free: 2410 MB
node distances:
node   0 
  0:  10 

vchuravy@helios ~ $ numactl --hardware
available: 2 nodes (0-1)
node 0 cpus: 0 1 2 3 4 5 12 13 14 15 16 17
node 0 size: 16000 MB
node 0 free: 15394 MB
node 1 cpus: 6 7 8 9 10 11 18 19 20 21 22 23
node 1 size: 16102 MB
node 1 free: 15403 MB
node distances:
node   0   1 
  0:  10  16 
  1:  16  10 

Machine (31GB total) + Package L#0
  NUMANode L#0 (P#0 16GB)
    L3 L#0 (8192KB)
      L2 L#0 (512KB) + L1d L#0 (32KB) + L1i L#0 (64KB) + Core L#0
        PU L#0 (P#0)
        PU L#1 (P#12)
      L2 L#1 (512KB) + L1d L#1 (32KB) + L1i L#1 (64KB) + Core L#1
        PU L#2 (P#1)
        PU L#3 (P#13)
      L2 L#2 (512KB) + L1d L#2 (32KB) + L1i L#2 (64KB) + Core L#2
        PU L#4 (P#2)
        PU L#5 (P#14)
    L3 L#1 (8192KB)
      L2 L#3 (512KB) + L1d L#3 (32KB) + L1i L#3 (64KB) + Core L#3
        PU L#6 (P#3)
        PU L#7 (P#15)
      L2 L#4 (512KB) + L1d L#4 (32KB) + L1i L#4 (64KB) + Core L#4
        PU L#8 (P#4)
        PU L#9 (P#16) 
   L2 L#5 (512KB) + L1d L#5 (32KB) + L1i L#5 (64KB) + Core L#5
      PU L#10 (P#5)
      PU L#11 (P#17)

numactl --cpunodebind=1 --membind=1 ./julia

numactl --physcpubind=0-5 --membind=0

julia -p 4 # Start Julia with 4 local processes

julia> addprocs(4) # Add 4 local processes

julia> nprocs()
5

julia> nworkers()
4

julia> myid()
1

julia> workers()
[2, 3, 4, 5]

# remotecall will not block
julia> f = remotecall(myid, 2)
Future(2, 1, 6, Nullable{Any}())

julia> wait(f)
# Note: `fetch` will `wait` on the value being available
julia> fetch(f)
2

# `wait` and `fetch` require multiple messages
# to be send over the network. For efficiency
# two variants are provided.
# * remotecall_wait = wait(remotecall())
# * remotecall_fetch = fetch(remotecall())

julia> remotecall_wait(x->2x, 4, 2)

# Closures get serialized

julia> g(x) = 2x
julia> remotecall_fetch(g, 3, 2)
ERROR: On worker 3:
UndefVarError: #g not defined
...

# Methods are not.
julia> @everywhere g(x) = 2x
julia> remotecall_fetch(g, 3, 2)
4

@sync begin
    @async (sleep(3); println("Today is reverse day!"))
    @async (sleep(2); println(" class!"))
    @async print("Hello")
end

@sync begin
    for wid in workers()
        @async remotecall_wait(println, wid, "Hello from worker ", string(wid))
    end
end

# Use `@verywhere` to execute a top-level block on each process

@everywhere begin
    using Test
    include("src.jl")
end

# Define on all processes the variable bar
@everywhere bar = 1

julia> foo =  10
julia> @everywhere baz = $foo

help?> @parallel
  A parallel for loop of the form :

  @parallel [reducer] for var = range
      body
  end

  The specified range is partitioned and locally executed
  across all workers. In case an optional reducer function
  is specified, @parallel performs local reductions on each
  worker with a final reduction on the calling process.

  Note that without a reducer function, @parallel 
  executes asynchronously, i.e. it spawns independent tasks
  on all available workers and returns immediately without waiting
  for completion. To wait for completion, prefix the call with
  @sync, like :

  @sync @parallel for var = range
      body
  end

julia> pmap(f, data) # chunk up data and distribute work across workers
# Better use DistributedArrays.jl

 addprocs([("anubis.juliacomputing.io", 4)], 
          exename="/data/vchuravy/julia/julia",
          dir="/data/vchuravy")

@everywhere using DistributedArrays

# Distribute local array across processes
distribute(a::Array)

# Convert a distributedArray to a local one
convert(Array, a::DArray)

# access the process local part of the DArray
localpart(dA::DArray)
dA[:L]

# Which indexes are local?
localindexes(a::DArray)

# Construct a distributed array
julia> dA = dfill(0.0, 500, 500)
julia> remotecall_fetch(()->size(localpart(dA)), 2)
(250, 250)
julia> remotecall_fetch(localindexes, 3, dA)
(251:500, 1:250)

# Darray(init, dims, [procs, dist])
julia> dA = dfill(1.0, (128, 128), workers(), [1, nworkers()]);
julia> remotecall_fetch(()->size(localpart(dA)), 2)
(128, 32)

function Base.map!{F}(f::F, dest::DArray, src::DArray)
    @sync for p in procs(dest)
        @async remotecall_wait(p) do
            map!(f, localpart(dest), src[localindexes(dest)...])
            nothing
        end
    end
    return dest
end

import MPI

MPI.Init()

comm = MPI.COMM_WORLD
MPI.Barrier(comm)

root = 0
r = MPI.Comm_rank(comm)

sr = MPI.Reduce(r, MPI.SUM, root, comm)

if(MPI.Comm_rank(comm) == root)
   @printf("sum of ranks: %s\n", sr)
end

MPI.Finalize()

mpirun -np 4 ./julia example.jl

function cpu_vadd(a, b, c)
  for i = 1:length(a)
    c[i] = a[i] + b[i]
  end
end

len = 16
a = rand(len)
b = rand(len)
c = similar(a)




cpu_vadd(a, b, c)

function gpu_vadd(a, b, c)
  i = get_my_index()
  c[i] = a[i] + b[i]
end


len = 16
a = rand(len)
b = rand(len)

gpu_a = CuArray(a)
gpu_b = CuArray(b)
gpu_c = CuArray(Float64, len)

@cuda (...) gpu_vadd(gpu_a, gpu_b, gpu_c)

c = Array(d_c)

function gpu_vadd(a, b, c)
  i, j = threadIdx().x, threadIdx().y
  c[i,j] = a[i,j] + b[i,j]
end

threads = (size(a,1), size(a,2))
@cuda (1,threads) gpu_vadd(...)

get_my_index() ?
threadIdx() & blockDim()
@cuda (1, threads) kernel(...)

threadIdx(), blockDim()
blockIdx(), gridDim()
@cuda (blocks, threads) kernel(...)

using CUDAdrv

@show CUDAdrv.devices

ctx0 = CuContext(CuDevice(0))
ctx1 = CuContext(CuDevice(1))

activate(ctx0)
# Async schedule kernel calls
activate(ctx1)
# Async schedule kernel calls