text:
@testset "Fornberg FD weights tests:" begin

    # all stencils provide derivatives at point labelled with 0

    # centered 1st order derivatives i.e. ...,-1,0,1,...
    @test Jecco.calculate_weights(1,2,1) == [-1//2, 0//1, 1//2]
    @test Jecco.calculate_weights(1,4,2) == [1//12, -2//3, 0//1, 2//3, -1//12]
    @test Jecco.calculate_weights(1,6,3) == [-1//60, 3//20, -3//4, 0//1, 3//4, -3//20, 1//60]
    @test Jecco.calculate_weights(1,8,4) == [1//280, -4//105, 1//5, -4//5, 0//1,
                                       4//5, -1//5, 4//105, -1//280 ]

    # forward 1st order derivatives i.e. 0,1,...
    @test Jecco.calculate_weights(1,1,0) == [-1//1, 1//1]
    @test Jecco.calculate_weights(1,2,0) == [-3//2, 2//1, -1//2]
    @test Jecco.calculate_weights(1,3,0) == [-11//6, 3//1, -3//2, 1//3]
    @test Jecco.calculate_weights(1,4,0) == [-25//12, 4//1, -3//1, 4//3, -1//4]
    @test Jecco.calculate_weights(1,5,0) == [-137//60, 5//1, -5//1, 10//3, -5//4, 1//5]
    @test Jecco.calculate_weights(1,6,0) == [-49//20, 6//1, -15//2, 20//3, -15//4, 6//5, -1//6]

    # forward 1st order derivatives i.e. -1,0,1,2,...
    @test Jecco.calculate_weights(1,3,1) == [-1//3, -1//2, 1//1, -1//6]
    @test Jecco.calculate_weights(1,4,1) == [-1//4, -5//6, 3//2, -1//2, 1//12]
    @test Jecco.calculate_weights(1,5,1) == [-1//5, -13//12, 2//1, -1//1, 1//3, -1//20]
    @test Jecco.calculate_weights(1,6,1) == [-1//6, -77//60, 5//2, -5//3, 5//6, -1//4, 1//30]
    @test Jecco.calculate_weights(1,7,1) == [-1//7, -29//20, 3//1, -5//2, 5//3, -3//4, 1//5, -1//42]

    # forward 1st order derivatives i.e. -2,-1,0,1,2,3,...
    @test Jecco.calculate_weights(1,5,2) == [1//20, -1//2, -1//3, 1//1, -1//4, 1//30]
    @test Jecco.calculate_weights(1,6,2) == [1//30, -2//5, -7//12, 4//3, -1//2, 2//15, -1//60]

    # backward 1st order derivatives i.e. ...,-1,0
    @test Jecco.calculate_weights(1,1,1) == [-1//1, 1//1]
    @test Jecco.calculate_weights(1,2,2) == [1//2, -2//1, 3//2]
    @test Jecco.calculate_weights(1,3,3) == [-1//3, 3//2, -3//1, 11//6]
    @test Jecco.calculate_weights(1,4,4) == [1//4, -4//3, 3//1, -4//1, 25//12]

    # backward 1st order derivatives i.e. ...,-2,-1,0,1
    @test Jecco.calculate_weights(1,3,2) == [1//6, -1//1, 1//2, 1//3]
    @test Jecco.calculate_weights(1,4,3) == [-1//12, 1//2, -3//2, 5//6, 1//4]
    @test Jecco.calculate_weights(1,5,4) == [1//20, -1//3, 1//1, -2//1, 13//12, 1//5]

    # backward 1st order derivatives i.e. ...,-3,-2,-1,0,1,2
    @test Jecco.calculate_weights(1,5,3) == [-1//30, 1//4, -1//1, 1//3, 1//2, -1//20]
    @test Jecco.calculate_weights(1,6,4) == [1//60, -2//15, 1//2, -4//3, 7//12, 2//5, -1//30]

    # centered 2nd order derivatives i.e. ...,-1,0,1,...
    @test Jecco.calculate_weights(2,2,1) == [1//1, -2//1, 1//1]
    @test Jecco.calculate_weights(2,4,2) == [-1//12, 4//3, -5//2, 4//3, -1//12]
    @test Jecco.calculate_weights(2,6,3) == [1//90, -3//20, 3//2, -49//18,
                                       3//2, -3//20, 1//90]
    @test Jecco.calculate_weights(2,8,4) == [-1//560, 8//315, -1//5, 8//5, -205//72,
                                       8//5, -1//5, 8//315, -1//560]

    # forward 2nd order derivatives i.e. 0,1,...
    @test Jecco.calculate_weights(2,2,0) == [1//1, -2//1, 1//1]
    @test Jecco.calculate_weights(2,3,0) == [2//1, -5//1, 4//1, -1//1]
    @test Jecco.calculate_weights(2,4,0) == [35//12, -26//3, 19//2, -14//3, 11//12]
    @test Jecco.calculate_weights(2,5,0) == [15//4, -77//6, 107//6, -13//1, 61//12, -5//6]
    @test Jecco.calculate_weights(2,6,0) == [203//45, -87//5, 117//4, -254//9, 33//2, -27//5, 137//180]
    @test Jecco.calculate_weights(2,7,0) == [469//90, -223//10, 879//20, -949//18, 41//1, -201//10,
                                       1019//180, -7//10]

    # forward 2nd order derivatives i.e. -1,0,1,2,...
    @test Jecco.calculate_weights(2,3,1) == [1//1, -2//1, 1//1, 0//1]
    @test Jecco.calculate_weights(2,4,1) == [11//12, -5//3, 1//2, 1//3, -1//12]
    @test Jecco.calculate_weights(2,5,1) == [5//6, -5//4, -1//3, 7//6, -1//2, 1//12]
    @test Jecco.calculate_weights(2,6,1) == [137//180, -49//60, -17//12, 47//18, -19//12,
                                       31//60, -13//180]

    # forward 2nd order derivatives i.e. -2,-1,0,1,2,3,...
    @test Jecco.calculate_weights(2,5,2) == [-1//12, 4//3, -5//2, 4//3, -1//12, 0//1]
    @test Jecco.calculate_weights(2,6,2) == [-13//180, 19//15, -7//3, 10//9,
                                       1//12, -1//15, 1//90]

    # backward 2nd order derivatives i.e. ...,-1,0
    @test Jecco.calculate_weights(2,2,2) == [1//1, -2//1, 1//1]
    @test Jecco.calculate_weights(2,3,3) == [-1//1, 4//1, -5//1, 2//1]
    @test Jecco.calculate_weights(2,4,4) == [11//12, -14//3, 19//2, -26//3, 35//12]

    # backward 2nd order derivatives i.e. ...,-2,-1,0,1
    @test Jecco.calculate_weights(2,3,2) == [0//1, 1//1, -2//1, 1//1]
    @test Jecco.calculate_weights(2,4,3) == [-1//12, 1//3, 1//2, -5//3, 11//12]
    @test Jecco.calculate_weights(2,5,4) == [1//12, -1//2, 7//6, -1//3, -5//4, 5//6]

    # backward 2nd order derivatives i.e. ...,-3,-2,-1,0,1,2
    @test Jecco.calculate_weights(2,5,3) == [0//1, -1//12, 4//3, -5//2, 4//3, -1//12]
    @test Jecco.calculate_weights(2,6,4) == [1//90, -1//15, 1//12, 10//9, -7//3, 19//15, -13//180]

    # centered 3rd order derivatives i.e ...,-3,-2,-1,0,1,2,3
    @test Jecco.calculate_weights(3,6,3) == [1//8, -1//1, 13//8, 0//1, -13//8, 1//1, -1//8]

    # forward 3rd order derivatives i.e 0,1,2,3,...
    @test Jecco.calculate_weights(3,3,0) == [-1//1, 3//1, -3//1, 1//1]
    @test Jecco.calculate_weights(3,4,0) == [-5//2, 9//1, -12//1, 7//1, -3//2]

    # forward 3rd order derivatives i.e -1,0,1,2,...
    @test Jecco.calculate_weights(3,3,1) == [-1//1, 3//1, -3//1, 1//1]
    @test Jecco.calculate_weights(3,4,1) == [-3//2, 5//1, -6//1, 3//1, -1//2]

    # backward 3rd order derivatives i.e -3,-2,-1,0,...
    @test Jecco.calculate_weights(3,3,3) == [-1//1, 3//1, -3//1, 1//1]
    @test Jecco.calculate_weights(3,4,4) == [3//2, -7//1, 12//1, -9//1, 5//2]

    # backward 3rd order derivatives i.e -2,-1,0,1...
    @test Jecco.calculate_weights(3,3,2) == [-1//1, 3//1, -3//1, 1//1]
    @test Jecco.calculate_weights(3,4,3) == [1//2, -3//1, 6//1, -5//1, 3//2]

end

@testset "EqualSizeStencilFD tests:" begin

    Dx = EqualSizeStencilFD{1}(1, 4, 1.0, 10)
    @test 12*Dx[1,:]  ≈ [-25.0, 48.0, -36.0, 16.0, -3, 0.0, 0.0, 0.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[2,:]  ≈ [-3.0, -10.0, 18.0, -6.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[3,:]  ≈ [1.0, -8.0, 0.0, 8.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[4,:]  ≈ [0.0, 1.0, -8.0, 0.0, 8.0, -1.0, 0.0, 0.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[5,:]  ≈ [0.0, 0.0, 1.0, -8.0, 0.0, 8.0, -1.0, 0.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[6,:]  ≈ [0.0, 0.0, 0.0, 1.0, -8.0, 0.0, 8.0, -1.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[7,:]  ≈ [0.0, 0.0, 0.0, 0.0, 1.0, -8.0, 0.0, 8.0, -1.0, 0.0] atol=1e-15
    @test 12*Dx[8,:]  ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -8.0, 0.0, 8.0, -1.0] atol=1e-15
    @test 12*Dx[9,:]  ≈ [0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 6.0, -18.0, 10.0, 3.0] atol=1e-15
    @test 12*Dx[10,:] ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 3.0, -16.0, 36.0, -48.0, 25.0] atol=1e-15

end

@testset "Periodic Derivative tests:" begin

    Dx = CenteredDiff{1}(1, 4, 1.0, 10)
    @test 12*Dx[1,:]  ≈ [0.0, 8.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -8.0] atol=1e-15
    @test 12*Dx[2,:]  ≈ [-8.0, 0.0, 8.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0] atol=1e-15
    @test 12*Dx[3,:]  ≈ [1.0, -8.0, 0.0, 8.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[4,:]  ≈ [0.0, 1.0, -8.0, 0.0, 8.0, -1.0, 0.0, 0.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[5,:]  ≈ [0.0, 0.0, 1.0, -8.0, 0.0, 8.0, -1.0, 0.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[6,:]  ≈ [0.0, 0.0, 0.0, 1.0, -8.0, 0.0, 8.0, -1.0, 0.0, 0.0] atol=1e-15
    @test 12*Dx[7,:]  ≈ [0.0, 0.0, 0.0, 0.0, 1.0, -8.0, 0.0, 8.0, -1.0, 0.0] atol=1e-15
    @test 12*Dx[8,:]  ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -8.0, 0.0, 8.0, -1.0] atol=1e-15
    @test 12*Dx[9,:]  ≈ [-1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -8.0, 0.0, 8.0] atol=1e-15
    @test 12*Dx[10,:] ≈ [8.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, -8.0, 0.0] atol=1e-15

    # 1D case

    xmin   = -2.0*pi
    xmax   =  2.0*pi
    xnodes =  600
    ord    =  4

    hx     = (xmax - xmin) / xnodes

    x  = collect(xmin:hx:xmax-hx)
    f  = sin.(x)

    D1 = CenteredDiff(1, ord, hx, length(x))
    D2 = CenteredDiff(2, ord, hx, length(x))

    df  = D1 * f
    d2f = D2 * f

    @test df  ≈ cos.(x) rtol=1e-8
    @test d2f ≈ -f rtol=1e-8

    # now for the callable, point-wise, methods
    df0 = similar(df)
    for i in eachindex(f)
        df0[i] = D1(f,i)
    end
    @test df0 == df

    # 3D case

    ymin   = -1.0*pi
    ymax   =  1.0*pi
    ynodes =  20
    zmin   = -1.0*pi
    zmax   =  1.0*pi
    znodes =  300

    hy     = (ymax - ymin) / ynodes
    hz     = (zmax - zmin) / znodes

    y      = collect(ymin:hy:ymax-hy)
    z      = collect(zmin:hz:zmax-hz)

    f      = [sin.(x1) .* sin.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dxf0   = [cos.(x1) .* sin.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dyf0   = [sin.(x1) .* cos.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dzf0   = [sin.(x1) .* sin.(x2) .* cos.(x3) for x1 in x, x2 in y, x3 in z]
    dxzf0  = [cos.(x1) .* sin.(x2) .* cos.(x3) for x1 in x, x2 in y, x3 in z]

    Dx     = CenteredDiff{1}(1, ord, hx, length(x))
    Dz     = CenteredDiff{3}(1, ord, hz, length(z))

    Dxx    = CenteredDiff{1}(2, ord, hx, length(x))
    Dzz    = CenteredDiff{3}(2, ord, hz, length(z))

    dxf    = Dx * f
    dzf    = Dz * f
    dxzf   = Dx * (Dz * f)

    @test dxf  ≈ dxf0
    @test dzf  ≈ dzf0
    @test dxzf ≈ dxzf0

    d2xf   = Dxx * f
    d2zf   = Dzz * f

    @test d2xf ≈ -f
    @test d2zf ≈ -f

    # now for callable, point-wise, methods
    @test Dx(f,2,10,120)  == dxf[2,10,120]
    @test Dx(f,42,20,300) == dxf[42,20,300]
    @test Dz(f,2,10,120)  == dzf[2,10,120]
    @test Dz(f,42,20,300) == dzf[42,20,300]
end

@testset "EqualSizeStencilFD Derivative tests:" begin

    # 1D case

    xmin   = -2.0*pi
    xmax   =  2.0*pi
    xnodes =  600
    ord    =  6

    hx     = (xmax - xmin) / xnodes

    x  = collect(xmin:hx:xmax-hx)
    f  = sin.(x)

    D1 = EqualSizeStencilFD(1, ord, hx, length(x))
    D2 = EqualSizeStencilFD(2, ord, hx, length(x))

    df  = D1 * f
    d2f = D2 * f

    @test df  ≈ cos.(x)
    @test d2f ≈ -f

    # now for the callable, point-wise, methods
    df0 = similar(df)
    for i in eachindex(f)
        df0[i] = D1(f,i)
    end
    @test df0 == df

    # 3D case

    ymin   = -1.0*pi
    ymax   =  1.0*pi
    ynodes =  20
    zmin   = -1.0*pi
    zmax   =  1.0*pi
    znodes =  300

    hy     = (ymax - ymin) / ynodes
    hz     = (zmax - zmin) / znodes

    y      = collect(ymin:hy:ymax-hy)
    z      = collect(zmin:hz:zmax-hz)

    f      = [sin.(x1) .* sin.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dxf0   = [cos.(x1) .* sin.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dyf0   = [sin.(x1) .* cos.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dzf0   = [sin.(x1) .* sin.(x2) .* cos.(x3) for x1 in x, x2 in y, x3 in z]
    dxzf0  = [cos.(x1) .* sin.(x2) .* cos.(x3) for x1 in x, x2 in y, x3 in z]

    Dx     = EqualSizeStencilFD{1}(1, ord, hx, length(x))
    Dz     = EqualSizeStencilFD{3}(1, ord, hz, length(z))

    Dxx    = EqualSizeStencilFD{1}(2, ord, hx, length(x))
    Dzz    = EqualSizeStencilFD{3}(2, ord, hz, length(z))

    dxf    = Dx * f
    dzf    = Dz * f
    dxzf   = Dx * (Dz * f)

    @test dxf  ≈ dxf0
    @test dzf  ≈ dzf0
    @test dxzf ≈ dxzf0

    d2xf   = Dxx * f
    d2zf   = Dzz * f

    @test d2xf ≈ -f
    @test d2zf ≈ -f

    # now for callable, point-wise, methods
    @test Dx(f,2,10,120)  == dxf[2,10,120]
    @test Dx(f,42,20,300) == dxf[42,20,300]
    @test Dz(f,2,10,120)  == dzf[2,10,120]
    @test Dz(f,42,20,300) == dzf[42,20,300]
end

@testset "Spectral Derivative tests:" begin

    # 1D case
    xmin   = -2.0
    xmax   =  2.0
    xnodes =  32

    x, = Jecco.cheb(xmin, xmax, xnodes)
    f = 0.5 * x.^2

    Dx  = ChebDeriv(1, xmin, xmax, xnodes)
    Dxx = ChebDeriv(2, xmin, xmax, xnodes)

    dxf  = Dx * f
    dxxf = Dxx * f
    @test dxf  ≈ x
    @test dxxf ≈ fill(1.0, size(dxxf))

    # now for the callable, point-wise, methods
    dxf0 = similar(dxf)
    for i in eachindex(f)
        dxf0[i] = Dx(f,i)
    end
    @test dxf0 ≈ dxf

    # 3D case

    ymin   = -1.0
    ymax   =  1.0
    ynodes =  8

    zmin   = -1.0
    zmax   =  1.0
    znodes =  16

    y, = Jecco.cheb(ymin, ymax, ynodes)
    z, = Jecco.cheb(zmin, zmax, znodes)

    Dy  = ChebDeriv{2}(1, ymin, ymax, ynodes)
    Dyy = ChebDeriv{2}(2, ymin, ymax, ynodes)
    Dz  = ChebDeriv{3}(1, zmin, zmax, znodes)
    Dzz = ChebDeriv{3}(2, zmin, zmax, znodes)

    f     = [0.5 * x1.^2 .* cos.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dxf0  = [x1 .* cos.(x2) .* sin.(x3)          for x1 in x, x2 in y, x3 in z]
    dzf0  = [0.5 * x1.^2 .* cos.(x2) .* cos.(x3) for x1 in x, x2 in y, x3 in z]
    dxxf0 = [cos.(x2) .* sin.(x3)                for x1 in x, x2 in y, x3 in z]
    dzzf0 = -copy(f)

    dxf = Dx * f
    dzf = Dz * f
    d2xf = Dxx * f
    d2zf = Dzz * f

    @test dxf ≈ dxf0
    @test dzf ≈ dzf0

    @test d2xf ≈ dxxf0
    @test d2zf ≈ dzzf0

    # now for callable, point-wise, methods
    @test Dx(f,2,4,16)  ≈ dxf[2,4,16]
    @test Dx(f,1,6,12)  ≈ dxf[1,6,12]
    @test Dz(f,2,3,8)   ≈ dzf[2,3,8]
    @test Dz(f,16,8,1)  ≈ dzf[16,8,1]
end

@testset "Cross FD derivative tests:" begin
    # 2D FD case
    xmin   = -2.0*pi
    xmax   =  2.0*pi
    xnodes =  600
    ymin   = -1.0*pi
    ymax   =  1.0*pi
    ynodes =  300
    ord    =  4

    hx     = (xmax - xmin) / xnodes
    hy     = (ymax - ymin) / ynodes

    x      = collect(xmin:hx:xmax-hx)
    y      = collect(ymin:hy:ymax-hy)

    f      = [sin.(x1) .* sin.(x2) for x1 in x, x2 in y]

    Dx     = CenteredDiff{1}(1, ord, hx, length(x))
    Dy     = CenteredDiff{2}(1, ord, hy, length(y))

    dxyf   = Dx * (Dy * f)

    dxyf0  = similar(dxyf)
    for j in 1:ynodes
        for i in 1:xnodes
            dxyf0[i,j] = Dx(Dy, f,i, j)
        end
    end
    @test dxyf0 ≈ dxyf rtol=1e-12 atol=1e-12
end

@testset "PeriodicFD cross derivative for general arrays tests:" begin
    xmin   = -2.0*pi
    xmax   =  2.0*pi
    xnodes =  600
    ord    =  4

    ymin   = -1.0
    ymax   =  1.0
    ynodes =  16

    zmin   = -1.0*pi
    zmax   =  1.0*pi
    znodes =  300

    hx     = (xmax - xmin) / xnodes
    hz     = (zmax - zmin) / znodes

    x      = collect(xmin:hx:xmax-hx)
    y,     = Jecco.cheb(ymin, ymax, ynodes)
    z      = collect(zmin:hz:zmax-hz)

    f   = [0.5 * sin.(x1) .* x2.^2 .* sin.(x3)  for x1 in x, x2 in y, x3 in z]

    Dx  = CenteredDiff{1}(1, ord, hx, length(x))
    Dz  = CenteredDiff{3}(1, ord, hz, length(z))

    dxzf  = Dx * (Dz * f)

    dxzf0  = similar(dxzf)
    for k in 1:znodes
        for j in 1:ynodes
            for i in 1:xnodes
                dxzf0[i,j,k] = Dx(Dz, f, i, j, k)
            end
        end
    end
    @test dxzf0 ≈ dxzf rtol=1e-12 atol=1e-12

    g   = [sin.(x3) .* 0.5 * sin.(x2) .* x1.^2 for x3 in z, x1 in y, x2 in x]

    Dz  = CenteredDiff{1}(1, ord, hz, length(z))
    Dx  = CenteredDiff{3}(1, ord, hx, length(x))

    dxzg  = Dx * (Dz * g)

    @test Dz(Dx, g, 100,2,1)     ≈ dxzg[100,2,1]
    @test Dz(Dx, g, 1,1,1)       ≈ dxzg[1,1,1]
    @test Dz(Dx, g, 300,16,600)  ≈ dxzg[300,16,600]
    @test Dz(Dx, g, 10,8,150)    ≈ dxzg[10,8,150]
end

@testset "Fourier Derivative tests:" begin
    # 1D case
    xmin   = -2.0*pi
    xmax   =  2.0*pi
    xnodes =  600
    ord    =  4

    hx     = (xmax - xmin) / xnodes

    x  = collect(xmin:hx:xmax-hx)
    f  = sin.(x)

    D1 = FourierDeriv(1, xmax-xmin, length(x))
    D2 = FourierDeriv(2, xmax-xmin, length(x))

    df  = D1 * f
    d2f = D2 * f

    @test df  ≈ cos.(x)
    @test d2f ≈ -f

    # now for the callable, point-wise, methods
    df0  = similar(df)
    d2f0 = similar(d2f)
    for i in eachindex(f)
        df0[i]  = D1(f,i)
        d2f0[i] = D2(f,i)
    end
    @test df0  ≈ df  atol=1e-12
    @test d2f0 ≈ d2f atol=1e-10

    # 3D case
    ymin   = -1.0*pi
    ymax   =  1.0*pi
    ynodes =  20
    zmin   = -1.0*pi
    zmax   =  1.0*pi
    znodes =  300

    hy     = (ymax - ymin) / ynodes
    hz     = (zmax - zmin) / znodes

    y      = collect(ymin:hy:ymax-hy)
    z      = collect(zmin:hz:zmax-hz)

    f      = [sin.(x1) .* sin.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dxf0   = [cos.(x1) .* sin.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dyf0   = [sin.(x1) .* cos.(x2) .* sin.(x3) for x1 in x, x2 in y, x3 in z]
    dzf0   = [sin.(x1) .* sin.(x2) .* cos.(x3) for x1 in x, x2 in y, x3 in z]
    dxzf0  = [cos.(x1) .* sin.(x2) .* cos.(x3) for x1 in x, x2 in y, x3 in z]

    Dx     = FourierDeriv{1}(1, xmax-xmin, length(x))
    Dz     = FourierDeriv{3}(1, zmax-zmin, length(z))

    Dxx    = FourierDeriv{1}(2, xmax-xmin, length(x))
    Dzz    = FourierDeriv{3}(2, zmax-zmin, length(z))

    dxf    = Dx * f
    dzf    = Dz * f
    dxzf   = Dx * (Dz * f)

    @test dxf  ≈ dxf0
    @test dzf  ≈ dzf0
    @test dxzf ≈ dxzf0

    d2xf   = Dxx * f
    d2zf   = Dzz * f

    @test d2xf ≈ -f
    @test d2zf ≈ -f

    # now for callable, point-wise, methods
    @test Dx(f,2,10,120)  ≈ dxf[2,10,120]
    @test Dx(f,42,20,300) ≈ dxf[42,20,300]
    @test Dz(f,2,10,120)  ≈ dzf[2,10,120]
    @test Dz(f,42,20,300) ≈ dzf[42,20,300]
end

@testset "EqualSizeStencilFD cross derivative for general arrays tests:" begin
    xmin   = -2.0*pi
    xmax   =  2.0*pi
    xnodes =  600
    ord    =  6

    ymin   = -1.0
    ymax   =  1.0
    ynodes =  16

    zmin   = -1.0*pi
    zmax   =  1.0*pi
    znodes =  300

    hx     = (xmax - xmin) / xnodes
    hz     = (zmax - zmin) / znodes

    x      = collect(xmin:hx:xmax-hx)
    y,     = Jecco.cheb(ymin, ymax, ynodes)
    z      = collect(zmin:hz:zmax-hz)

    f   = [0.5 * sin.(x1) .* x2.^2 .* sin.(x3)  for x1 in x, x2 in y, x3 in z]

    Dx  = EqualSizeStencilFD{1}(1, ord, hx, length(x))
    Dz  = EqualSizeStencilFD{3}(1, ord, hz, length(z))

    dxzf  = Dx * (Dz * f)

    dxzf0  = similar(dxzf)
    for k in 1:znodes
        for j in 1:ynodes
            for i in 1:xnodes
                dxzf0[i,j,k] = Dx(Dz, f, i, j, k)
            end
        end
    end
    @test dxzf0 ≈ dxzf rtol=1e-12 atol=1e-12

    g   = [sin.(x3) .* 0.5 * sin.(x2) .* x1.^2 for x3 in z, x1 in y, x2 in x]

    Dz  = EqualSizeStencilFD{1}(1, ord, hz, length(z))
    Dx  = EqualSizeStencilFD{3}(1, ord, hx, length(x))

    dxzg  = Dx * (Dz * g)

    @test Dz(Dx, g, 100,2,1)     ≈ dxzg[100,2,1]
    @test Dz(Dx, g, 1,1,1)       ≈ dxzg[1,1,1]
    @test Dz(Dx, g, 300,16,600)  ≈ dxzg[300,16,600]
    @test Dz(Dx, g, 10,8,150)    ≈ dxzg[10,8,150]
end

@testset "PeriodicFD and EqualSizeStencilFD cross derivative for general arrays tests:" begin
    xmin   = -2.0*pi
    xmax   =  2.0*pi
    xnodes =  600
    ord    =  6

    ymin   = -1.0
    ymax   =  1.0
    ynodes =  16

    zmin   = -1.0*pi
    zmax   =  1.0*pi
    znodes =  300

    hx     = (xmax - xmin) / xnodes
    hz     = (zmax - zmin) / znodes

    x      = collect(xmin:hx:xmax-hx)
    y,     = Jecco.cheb(ymin, ymax, ynodes)
    z      = collect(zmin:hz:zmax-hz)

    f   = [0.5 * sin.(x1) .* x2.^2 .* sin.(x3)  for x1 in x, x2 in y, x3 in z]

    Dx  = EqualSizeStencilFD{1}(1, ord, hx, length(x))
    Dz  = CenteredDiff{3}(1, ord, hz, length(z))

    dxzf  = Dx * (Dz * f)

    dxzf0  = similar(dxzf)
    for k in 1:znodes
        for j in 1:ynodes
            for i in 1:xnodes
                dxzf0[i,j,k] = Dx(Dz, f, i, j, k)
            end
        end
    end
    @test dxzf0 ≈ dxzf rtol=1e-12 atol=1e-12

    g   = [sin.(x3) .* 0.5 * sin.(x2) .* x1.^2 for x3 in z, x1 in y, x2 in x]

    Dz  = CenteredDiff{1}(1, ord, hz, length(z))
    Dx  = EqualSizeStencilFD{3}(1, ord, hx, length(x))

    dxzg  = Dx * (Dz * g)

    @test Dz(Dx, g, 100,2,1)     ≈ dxzg[100,2,1]
    @test Dz(Dx, g, 1,1,1)       ≈ dxzg[1,1,1]
    @test Dz(Dx, g, 300,16,600)  ≈ dxzg[300,16,600]
    @test Dz(Dx, g, 10,8,150)    ≈ dxzg[10,8,150]
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module BLAS
@doc """
Interface to BLAS subroutines.
""" BLAS

import ..axpy!, ..axpby!
import Base: copyto!

export
# Level 1
    asum,
    axpy!,
    axpby!,
    blascopy!,
    dot,
    dotc,
    dotu,
    scal!,
    scal,
    nrm2,
    iamax,
# Level 2
    gbmv!,
    gbmv,
    gemv!,
    gemv,
    hemv!,
    hemv,
    sbmv!,
    sbmv,
    symv!,
    symv,
    trsv!,
    trsv,
    trmv!,
    trmv,
    ger!,
    syr!,
    her!,
# Level 3
    herk!,
    herk,
    her2k!,
    her2k,
    gemm!,
    gemm,
    symm!,
    symm,
    hemm!,
    hemm,
    syrk!,
    syrk,
    syr2k!,
    syr2k,
    trmm!,
    trmm,
    trsm!,
    trsm

const libblas = Base.libblas_name
const liblapack = Base.liblapack_name

import LinearAlgebra
import LinearAlgebra: BlasReal, BlasComplex, BlasFloat, BlasInt, DimensionMismatch, checksquare, stride1, chkstride1, axpy!

import Libdl

# utility routines
function vendor()
    lib = Libdl.dlopen_e(Base.libblas_name)
    vend = :unknown
    if lib != C_NULL
        if Libdl.dlsym_e(lib, :openblas_set_num_threads) != C_NULL
            vend = :openblas
        elseif Libdl.dlsym_e(lib, :openblas_set_num_threads64_) != C_NULL
            vend = :openblas64
        elseif Libdl.dlsym_e(lib, :MKL_Set_Num_Threads) != C_NULL
            vend = :mkl
        end
        Libdl.dlclose(lib)
    end
    return vend
end

if vendor() == :openblas64
    macro blasfunc(x)
        return Expr(:quote, Symbol(x, "64_"))
    end
    openblas_get_config() = strip(unsafe_string(ccall((:openblas_get_config64_, Base.libblas_name), Ptr{UInt8}, () )))
else
    macro blasfunc(x)
        return Expr(:quote, x)
    end
    openblas_get_config() = strip(unsafe_string(ccall((:openblas_get_config, Base.libblas_name), Ptr{UInt8}, () )))
end

"""
    set_num_threads(n)

Set the number of threads the BLAS library should use.
"""
function set_num_threads(n::Integer)
    blas = vendor()
    if blas == :openblas
        return ccall((:openblas_set_num_threads, Base.libblas_name), Cvoid, (Int32,), n)
    elseif blas == :openblas64
        return ccall((:openblas_set_num_threads64_, Base.libblas_name), Cvoid, (Int32,), n)
    elseif blas == :mkl
        # MKL may let us set the number of threads in several ways
        return ccall((:MKL_Set_Num_Threads, Base.libblas_name), Cvoid, (Cint,), n)
    end

    # OSX BLAS looks at an environment variable
    @static if Sys.isapple()
        ENV["VECLIB_MAXIMUM_THREADS"] = n
    end

    return nothing
end

function check()
    blas = vendor()
    if blas == :openblas || blas == :openblas64
        openblas_config = openblas_get_config()
        openblas64 = contains(openblas_config, r".*USE64BITINT.*")
        if Base.USE_BLAS64 != openblas64
            if !openblas64
                @error """
                    OpenBLAS was not built with 64bit integer support.
                    You're seeing this error because Julia was built with USE_BLAS64=1.
                    Please rebuild Julia with USE_BLAS64=0"""
            else
                @error """
                    Julia was not built with support for OpenBLAS with 64bit integer support.
                    You're seeing this error because Julia was built with USE_BLAS64=0.
                    Please rebuild Julia with USE_BLAS64=1"""
            end
            println("Quitting.")
            exit()
        end
    elseif blas == :mkl
        if Base.USE_BLAS64
            ENV["MKL_INTERFACE_LAYER"] = "ILP64"
        end
    end

    #
    # Check if BlasInt is the expected bitsize, by triggering an error
    #
    (_, info) = LinearAlgebra.LAPACK.potrf!('U', [1.0 0.0; 0.0 -1.0])
    if info != 2 # mangled info code
        if info == 2^33
            error("BLAS and LAPACK are compiled with 32-bit integer support, but Julia expects 64-bit integers. Please build Julia with USE_BLAS64=0.")
        elseif info == 0
            error("BLAS and LAPACK are compiled with 64-bit integer support but Julia expects 32-bit integers. Please build Julia with USE_BLAS64=1.")
        else
            error("The LAPACK library produced an undefined error code. Please verify the installation of BLAS and LAPACK.")
        end
    end

end

# Level 1
## copy

"""
    blascopy!(n, X, incx, Y, incy)

Copy `n` elements of array `X` with stride `incx` to array `Y` with stride `incy`. Returns `Y`.
"""
function blascopy! end

for (fname, elty) in ((:dcopy_,:Float64),
                      (:scopy_,:Float32),
                      (:zcopy_,:ComplexF64),
                      (:ccopy_,:ComplexF32))
    @eval begin
        # SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
        function blascopy!(n::Integer, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},AbstractArray{$elty}}, incy::Integer)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 n, DX, incx, DY, incy)
            DY
        end
    end
end

## scal

"""
    scal!(n, a, X, incx)

Overwrite `X` with `a*X` for the first `n` elements of array `X` with stride `incx`. Returns `X`.
"""
function scal! end

"""
    scal(n, a, X, incx)

Return `X` scaled by `a` for the first `n` elements of array `X` with stride `incx`.
"""
function scal end

for (fname, elty) in ((:dscal_,:Float64),
                      (:sscal_,:Float32),
                      (:zscal_,:ComplexF64),
                      (:cscal_,:ComplexF32))
    @eval begin
        # SUBROUTINE DSCAL(N,DA,DX,INCX)
        function scal!(n::Integer, DA::$elty, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer)
            ccall((@blasfunc($fname), libblas), Cvoid,
                  (Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                  n, DA, DX, incx)
            DX
        end
    end
end
scal(n, DA, DX, incx) = scal!(n, DA, copy(DX), incx)

## dot

"""
    dot(n, X, incx, Y, incy)

Dot product of two vectors consisting of `n` elements of array `X` with stride `incx` and
`n` elements of array `Y` with stride `incy`.

# Examples
```jldoctest
julia> dot(10, fill(1.0, 10), 1, fill(1.0, 20), 2)
10.0
```
"""
function dot end

"""
    dotc(n, X, incx, U, incy)

Dot function for two complex vectors, consisting of `n` elements of array `X`
with stride `incx` and `n` elements of array `U` with stride `incy`,
conjugating the first vector.

# Examples
```jldoctest
julia> BLAS.dotc(10, fill(1.0im, 10), 1, fill(1.0+im, 20), 2)
10.0 - 10.0im
```
"""
function dotc end

"""
    dotu(n, X, incx, Y, incy)

Dot function for two complex vectors consisting of `n` elements of array `X`
with stride `incx` and `n` elements of array `Y` with stride `incy`.

# Examples
```jldoctest
julia> BLAS.dotu(10, fill(1.0im, 10), 1, fill(1.0+im, 20), 2)
-10.0 + 10.0im
```
"""
function dotu end

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},AbstractArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},AbstractArray{$elty}}, incy::Integer)
            ccall((@blasfunc($fname), libblas), $elty,
                (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 n, DX, incx, DY, incy)
        end
    end
end
for (fname, elty) in ((:cblas_zdotc_sub,:ComplexF64),
                      (:cblas_cdotc_sub,:ComplexF32))
    @eval begin
                #       DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
                # *     .. Scalar Arguments ..
                #       INTEGER INCX,INCY,N
                # *     ..
                # *     .. Array Arguments ..
                #       DOUBLE PRECISION DX(*),DY(*)
        function dotc(n::Integer, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},AbstractArray{$elty}}, incy::Integer)
            result = Ref{$elty}()
            ccall((@blasfunc($fname), libblas), Cvoid,
                (BlasInt, Ptr{$elty}, BlasInt, Ptr{$elty}, BlasInt, Ptr{$elty}),
                 n, DX, incx, DY, incy, result)
            result[]
        end
    end
end
for (fname, elty) in ((:cblas_zdotu_sub,:ComplexF64),
                      (:cblas_cdotu_sub,:ComplexF32))
    @eval begin
                #       DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
                # *     .. Scalar Arguments ..
                #       INTEGER INCX,INCY,N
                # *     ..
                # *     .. Array Arguments ..
                #       DOUBLE PRECISION DX(*),DY(*)
        function dotu(n::Integer, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},AbstractArray{$elty}}, incy::Integer)
            result = Ref{$elty}()
            ccall((@blasfunc($fname), libblas), Cvoid,
                (BlasInt, Ptr{$elty}, BlasInt, Ptr{$elty}, BlasInt, Ptr{$elty}),
                 n, DX, incx, DY, incy, result)
            result[]
        end
    end
end

function dot(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasReal
    n = length(DX)
    if n != length(DY)
        throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    end
    GC.@preserve DX DY dot(n, pointer(DX), stride(DX, 1), pointer(DY), stride(DY, 1))
end
function dotc(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasComplex
    n = length(DX)
    if n != length(DY)
        throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    end
    GC.@preserve DX DY dotc(n, pointer(DX), stride(DX, 1), pointer(DY), stride(DY, 1))
end
function dotu(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasComplex
    n = length(DX)
    if n != length(DY)
        throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    end
    GC.@preserve DX DY dotu(n, pointer(DX), stride(DX, 1), pointer(DY), stride(DY, 1))
end

## nrm2

"""
    nrm2(n, X, incx)

2-norm of a vector consisting of `n` elements of array `X` with stride `incx`.

# Examples
```jldoctest
julia> BLAS.nrm2(4, fill(1.0, 8), 2)
2.0

julia> BLAS.nrm2(1, fill(1.0, 8), 2)
1.0
```
"""
function nrm2 end

for (fname, elty, ret_type) in ((:dnrm2_,:Float64,:Float64),
                                (:snrm2_,:Float32,:Float32),
                                (:dznrm2_,:ComplexF64,:Float64),
                                (:scnrm2_,:ComplexF32,:Float32))
    @eval begin
        # SUBROUTINE DNRM2(N,X,INCX)
        function nrm2(n::Integer, X::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer)
            ccall((@blasfunc($fname), libblas), $ret_type,
                (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 n, X, incx)
        end
    end
end
nrm2(x::Union{AbstractVector,DenseArray}) = GC.@preserve x nrm2(length(x), pointer(x), stride1(x))

## asum

"""
    asum(n, X, incx)

Sum of the absolute values of the first `n` elements of array `X` with stride `incx`.

# Examples
```jldoctest
julia> BLAS.asum(5, fill(1.0im, 10), 2)
5.0

julia> BLAS.asum(2, fill(1.0im, 10), 5)
2.0
```
"""
function asum end

for (fname, elty, ret_type) in ((:dasum_,:Float64,:Float64),
                                (:sasum_,:Float32,:Float32),
                                (:dzasum_,:ComplexF64,:Float64),
                                (:scasum_,:ComplexF32,:Float32))
    @eval begin
        # SUBROUTINE ASUM(N, X, INCX)
        function asum(n::Integer, X::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer)
            ccall((@blasfunc($fname), libblas), $ret_type,
                (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 n, X, incx)
        end
    end
end
asum(x::Union{AbstractVector,DenseArray}) = GC.@preserve x asum(length(x), pointer(x), stride1(x))

## axpy

"""
    axpy!(a, X, Y)

Overwrite `Y` with `a*X + Y`, where `a` is a scalar. Return `Y`.

# Examples
```jldoctest
julia> x = [1; 2; 3];

julia> y = [4; 5; 6];

julia> BLAS.axpy!(2, x, y)
3-element Array{Int64,1}:
  6
  9
 12
```
"""
function axpy! end

for (fname, elty) in ((:daxpy_,:Float64),
                      (:saxpy_,:Float32),
                      (:zaxpy_,:ComplexF64),
                      (:caxpy_,:ComplexF32))
    @eval begin
                # SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
                # DY <- DA*DX + DY
                #*     .. Scalar Arguments ..
                #      DOUBLE PRECISION DA
                #      INTEGER INCX,INCY,N
                #*     .. Array Arguments ..
                #      DOUBLE PRECISION DX(*),DY(*)
        function axpy!(n::Integer, alpha::($elty), dx::Union{Ptr{$elty}, AbstractArray{$elty}}, incx::Integer, dy::Union{Ptr{$elty}, AbstractArray{$elty}}, incy::Integer)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 n, alpha, dx, incx, dy, incy)
            dy
        end
    end
end
function axpy!(alpha::Number, x::Union{DenseArray{T},StridedVector{T}}, y::Union{DenseArray{T},StridedVector{T}}) where T<:BlasFloat
    if length(x) != length(y)
        throw(DimensionMismatch("x has length $(length(x)), but y has length $(length(y))"))
    end
    GC.@preserve x y axpy!(length(x), convert(T,alpha), pointer(x), stride(x, 1), pointer(y), stride(y, 1))
    y
end

function axpy!(alpha::Number, x::Array{T}, rx::Union{UnitRange{Ti},AbstractRange{Ti}},
               y::Array{T}, ry::Union{UnitRange{Ti},AbstractRange{Ti}}) where {T<:BlasFloat,Ti<:Integer}
    if length(rx) != length(ry)
        throw(DimensionMismatch("ranges of differing lengths"))
    end
    if minimum(rx) < 1 || maximum(rx) > length(x)
        throw(ArgumentError("range out of bounds for x, of length $(length(x))"))
    end
    if minimum(ry) < 1 || maximum(ry) > length(y)
        throw(ArgumentError("range out of bounds for y, of length $(length(y))"))
    end
    GC.@preserve x y axpy!(length(rx), convert(T, alpha), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry))
    y
end

"""
    axpby!(a, X, b, Y)

Overwrite `Y` with `X*a + Y*b`, where `a` and `b` are scalars. Return `Y`.

# Examples
```jldoctest
julia> x = [1., 2, 3];

julia> y = [4., 5, 6];

julia> BLAS.axpby!(2., x, 3., y)
3-element Array{Float64,1}:
14.0
19.0
24.0
```
"""
function axpby! end

for (fname, elty) in ((:daxpby_,:Float64), (:saxpby_,:Float32),
                      (:zaxpby_,:ComplexF64), (:caxpby_,:ComplexF32))
    @eval begin
        # SUBROUTINE DAXPBY(N,DA,DX,INCX,DB,DY,INCY)
        # DY <- DA*DX + DB*DY
        #*     .. Scalar Arguments ..
        #      DOUBLE PRECISION DA,DB
        #      INTEGER INCX,INCY,N
        #*     .. Array Arguments ..
        #      DOUBLE PRECISION DX(*),DY(*)
        function axpby!(n::Integer, alpha::($elty), dx::Union{Ptr{$elty},
                        AbstractArray{$elty}}, incx::Integer, beta::($elty),
                        dy::Union{Ptr{$elty}, AbstractArray{$elty}}, incy::Integer)
            ccall((@blasfunc($fname), libblas), Cvoid, (Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                n, alpha, dx, incx, beta, dy, incy)
            dy
        end
    end
end

function axpby!(alpha::Number, x::Union{DenseArray{T},AbstractVector{T}}, beta::Number, y::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasFloat
    if length(x) != length(y)
        throw(DimensionMismatch("x has length $(length(x)), but y has length $(length(y))"))
    end
    GC.@preserve x y axpby!(length(x), convert(T,alpha), pointer(x), stride(x, 1), convert(T,beta), pointer(y), stride(y, 1))
    y
end

## iamax
for (fname, elty) in ((:idamax_,:Float64),
                      (:isamax_,:Float32),
                      (:izamax_,:ComplexF64),
                      (:icamax_,:ComplexF32))
    @eval begin
        function iamax(n::Integer, dx::Union{Ptr{$elty}, AbstractArray{$elty}}, incx::Integer)
            ccall((@blasfunc($fname), libblas),BlasInt,
                (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                n, dx, incx)
        end
    end
end
iamax(dx::Union{AbstractVector,DenseArray}) = GC.@preserve dx iamax(length(dx), pointer(dx), stride1(dx))

# Level 2
## mv
### gemv
for (fname, elty) in ((:dgemv_,:Float64),
                      (:sgemv_,:Float32),
                      (:zgemv_,:ComplexF64),
                      (:cgemv_,:ComplexF32))
    @eval begin
             #SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
             #*     .. Scalar Arguments ..
             #      DOUBLE PRECISION ALPHA,BETA
             #      INTEGER INCX,INCY,LDA,M,N
             #      CHARACTER TRANS
             #*     .. Array Arguments ..
             #      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
        function gemv!(trans::Char, alpha::($elty), A::AbstractVecOrMat{$elty}, X::AbstractVector{$elty}, beta::($elty), Y::AbstractVector{$elty})
            m,n = size(A,1),size(A,2)
            if trans == 'N' && (length(X) != n || length(Y) != m)
                throw(DimensionMismatch("A has dimensions $(size(A)), X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'C' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("the adjoint of A has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            elseif trans == 'T' && (length(X) != m || length(Y) != n)
                throw(DimensionMismatch("the transpose of A has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))"))
            end
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                 trans, size(A,1), size(A,2), alpha,
                 A, max(1,stride(A,2)), X, stride(X,1),
                 beta, Y, stride(Y,1))
            Y
        end
        function gemv(trans::Char, alpha::($elty), A::AbstractMatrix{$elty}, X::AbstractVector{$elty})
            gemv!(trans, alpha, A, X, zero($elty), similar(X, $elty, size(A, (trans == 'N' ? 1 : 2))))
        end
        function gemv(trans::Char, A::AbstractMatrix{$elty}, X::AbstractVector{$elty})
            gemv!(trans, one($elty), A, X, zero($elty), similar(X, $elty, size(A, (trans == 'N' ? 1 : 2))))
        end
    end
end

"""
    gemv!(tA, alpha, A, x, beta, y)

Update the vector `y` as `alpha*A*x + beta*y` or `alpha*A'x + beta*y`
according to [`tA`](@ref stdlib-blas-trans).
`alpha` and `beta` are scalars. Return the updated `y`.
"""
gemv!

"""
    gemv(tA, alpha, A, x)

Return `alpha*A*x` or `alpha*A'x` according to [`tA`](@ref stdlib-blas-trans).
`alpha` is a scalar.
"""
gemv(tA, alpha, A, x)

"""
    gemv(tA, A, x)

Return `A*x` or `A'x` according to [`tA`](@ref stdlib-blas-trans).
"""
gemv(tA, A, x)

### (GB) general banded matrix-vector multiplication

"""
    gbmv!(trans, m, kl, ku, alpha, A, x, beta, y)

Update vector `y` as `alpha*A*x + beta*y` or `alpha*A'*x + beta*y` according to [`trans`](@ref stdlib-blas-trans).
The matrix `A` is a general band matrix of dimension `m` by `size(A,2)` with `kl`
sub-diagonals and `ku` super-diagonals. `alpha` and `beta` are scalars. Return the updated `y`.
"""
function gbmv! end

"""
    gbmv(trans, m, kl, ku, alpha, A, x)

Return `alpha*A*x` or `alpha*A'*x` according to [`trans`](@ref stdlib-blas-trans).
The matrix `A` is a general band matrix of dimension `m` by `size(A,2)` with `kl` sub-diagonals and `ku`
super-diagonals, and `alpha` is a scalar.
"""
function gbmv end

for (fname, elty) in ((:dgbmv_,:Float64),
                      (:sgbmv_,:Float32),
                      (:zgbmv_,:ComplexF64),
                      (:cgbmv_,:ComplexF32))
    @eval begin
             # SUBROUTINE DGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
             # *     .. Scalar Arguments ..
             #       DOUBLE PRECISION ALPHA,BETA
             #       INTEGER INCX,INCY,KL,KU,LDA,M,N
             #       CHARACTER TRANS
             # *     .. Array Arguments ..
             #       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
        function gbmv!(trans::Char, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty})
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}),
                 trans, m, size(A,2), kl,
                 ku, alpha, A, max(1,stride(A,2)),
                 x, stride(x,1), beta, y, stride(y,1))
            y
        end
        function gbmv(trans::Char, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            n = size(A,2)
            leny = trans == 'N' ? m : n
            gbmv!(trans, m, kl, ku, alpha, A, x, zero($elty), similar(x, $elty, leny))
        end
        function gbmv(trans::Char, m::Integer, kl::Integer, ku::Integer, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            gbmv(trans, m, kl, ku, one($elty), A, x)
        end
    end
end

### symv

"""
    symv!(ul, alpha, A, x, beta, y)

Update the vector `y` as `alpha*A*x + beta*y`. `A` is assumed to be symmetric.
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
`alpha` and `beta` are scalars. Return the updated `y`.
"""
function symv! end

for (fname, elty, lib) in ((:dsymv_,:Float64,libblas),
                           (:ssymv_,:Float32,libblas),
                           (:zsymv_,:ComplexF64,liblapack),
                           (:csymv_,:ComplexF32,liblapack))
    # Note that the complex symv are not BLAS but auiliary functions in LAPACK
    @eval begin
             #      SUBROUTINE DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
             #     .. Scalar Arguments ..
             #      DOUBLE PRECISION ALPHA,BETA
             #      INTEGER INCX,INCY,LDA,N
             #      CHARACTER UPLO
             #     .. Array Arguments ..
             #      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
        function symv!(uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty})
            m, n = size(A)
            if m != n
                throw(DimensionMismatch("matrix A is $m by $n but must be square"))
            end
            if n != length(x)
                throw(DimensionMismatch("A has size $(size(A)), and x has length $(length(x))"))
            end
            if m != length(y)
                throw(DimensionMismatch("A has size $(size(A)), and y has length $(length(y))"))
            end
            chkstride1(A)
            ccall((@blasfunc($fname), $lib), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}),
                 uplo, n, alpha, A,
                 max(1,stride(A,2)), x, stride(x,1), beta,
                 y, stride(y,1))
            y
        end
        function symv(uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
                symv!(uplo, alpha, A, x, zero($elty), similar(x))
        end
        function symv(uplo::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            symv(uplo, one($elty), A, x)
        end
    end
end

"""
    symv(ul, alpha, A, x)

Return `alpha*A*x`. `A` is assumed to be symmetric.
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
`alpha` is a scalar.
"""
symv(ul, alpha, A, x)

"""
    symv(ul, A, x)

Return `A*x`. `A` is assumed to be symmetric.
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
"""
symv(ul, A, x)

### hemv
for (fname, elty) in ((:zhemv_,:ComplexF64),
                      (:chemv_,:ComplexF32))
    @eval begin
        function hemv!(uplo::Char, α::$elty, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, β::$elty, y::AbstractVector{$elty})
            m, n = size(A)
            if m != n
                throw(DimensionMismatch("matrix A is $m by $n but must be square"))
            end
            if n != length(x)
                throw(DimensionMismatch("A has size $(size(A)), and x has length $(length(x))"))
            end
            if m != length(y)
                throw(DimensionMismatch("A has size $(size(A)), and y has length $(length(y))"))
            end
            chkstride1(A)
            lda = max(1, stride(A, 2))
            incx = stride(x, 1)
            incy = stride(y, 1)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}),
                uplo, n, α, A,
                lda, x, incx, β,
                y, incy)
            y
        end
        function hemv(uplo::Char, α::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            hemv!(uplo, α, A, x, zero($elty), similar(x))
        end
        function hemv(uplo::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            hemv(uplo, one($elty), A, x)
        end
    end
end

### sbmv, (SB) symmetric banded matrix-vector multiplication
for (fname, elty) in ((:dsbmv_,:Float64),
                      (:ssbmv_,:Float32))
    @eval begin
             #       SUBROUTINE DSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
             # *     .. Scalar Arguments ..
             #       DOUBLE PRECISION ALPHA,BETA
             #       INTEGER INCX,INCY,K,LDA,N
             #       CHARACTER UPLO
             # *     .. Array Arguments ..
             #       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
        function sbmv!(uplo::Char, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty})
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                 uplo, size(A,2), k, alpha,
                 A, max(1,stride(A,2)), x, stride(x,1),
                 beta, y, stride(y,1))
            y
        end
        function sbmv(uplo::Char, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            n = size(A,2)
            sbmv!(uplo, k, alpha, A, x, zero($elty), similar(x, $elty, n))
        end
        function sbmv(uplo::Char, k::Integer, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            sbmv(uplo, k, one($elty), A, x)
        end
    end
end

"""
    sbmv(uplo, k, alpha, A, x)

Return `alpha*A*x` where `A` is a symmetric band matrix of order `size(A,2)` with `k`
super-diagonals stored in the argument `A`.
Only the [`uplo`](@ref stdlib-blas-uplo) triangle of `A` is used.
"""
sbmv(uplo, k, alpha, A, x)

"""
    sbmv(uplo, k, A, x)

Return `A*x` where `A` is a symmetric band matrix of order `size(A,2)` with `k`
super-diagonals stored in the argument `A`.
Only the [`uplo`](@ref stdlib-blas-uplo) triangle of `A` is used.
"""
sbmv(uplo, k, A, x)

"""
    sbmv!(uplo, k, alpha, A, x, beta, y)

Update vector `y` as `alpha*A*x + beta*y` where `A` is a a symmetric band matrix of order
`size(A,2)` with `k` super-diagonals stored in the argument `A`. The storage layout for `A`
is described the reference BLAS module, level-2 BLAS at
<http://www.netlib.org/lapack/explore-html/>.
Only the [`uplo`](@ref stdlib-blas-uplo) triangle of `A` is used.

Return the updated `y`.
"""
sbmv!

### hbmv, (HB) Hermitian banded matrix-vector multiplication
for (fname, elty) in ((:zhbmv_,:ComplexF64),
                      (:chbmv_,:ComplexF32))
    @eval begin
             #       SUBROUTINE ZHBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
             # *     .. Scalar Arguments ..
             #       DOUBLE PRECISION ALPHA,BETA
             #       INTEGER INCX,INCY,K,LDA,N
             #       CHARACTER UPLO
             # *     .. Array Arguments ..
             #       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
        function hbmv!(uplo::Char, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty})
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                 uplo, size(A,2), k, alpha,
                 A, max(1,stride(A,2)), x, stride(x,1),
                 beta, y, stride(y,1))
            y
        end
        function hbmv(uplo::Char, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            n = size(A,2)
            hbmv!(uplo, k, alpha, A, x, zero($elty), similar(x, $elty, n))
        end
        function hbmv(uplo::Char, k::Integer, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            hbmv(uplo, k, one($elty), A, x)
        end
    end
end

### trmv, Triangular matrix-vector multiplication

"""
    trmv(ul, tA, dA, A, b)

Return `op(A)*b`, where `op` is determined by [`tA`](@ref stdlib-blas-trans).
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
[`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or
are assumed to be all ones.
"""
function trmv end

"""
    trmv!(ul, tA, dA, A, b)

Return `op(A)*b`, where `op` is determined by [`tA`](@ref stdlib-blas-trans).
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
[`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or
are assumed to be all ones.
The multiplication occurs in-place on `b`.
"""
function trmv! end

for (fname, elty) in ((:dtrmv_,:Float64),
                        (:strmv_,:Float32),
                        (:ztrmv_,:ComplexF64),
                        (:ctrmv_,:ComplexF32))
    @eval begin
                #       SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
                # *     .. Scalar Arguments ..
                #       INTEGER INCX,LDA,N
                #       CHARACTER DIAG,TRANS,UPLO
                # *     .. Array Arguments ..
                #       DOUBLE PRECISION A(LDA,*),X(*)
        function trmv!(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            n = checksquare(A)
            if n != length(x)
                throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))"))
            end
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 uplo, trans, diag, n,
                 A, max(1,stride(A,2)), x, max(1,stride(x, 1)))
            x
        end
        function trmv(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            trmv!(uplo, trans, diag, A, copy(x))
        end
    end
end

### trsv, Triangular matrix-vector solve

"""
    trsv!(ul, tA, dA, A, b)

Overwrite `b` with the solution to `A*x = b` or one of the other two variants determined by
[`tA`](@ref stdlib-blas-trans) and [`ul`](@ref stdlib-blas-uplo).
[`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or
are assumed to be all ones.
Return the updated `b`.
"""
function trsv! end

"""
    trsv(ul, tA, dA, A, b)

Return the solution to `A*x = b` or one of the other two variants determined by
[`tA`](@ref stdlib-blas-trans) and [`ul`](@ref stdlib-blas-uplo).
[`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or
are assumed to be all ones.
"""
function trsv end

for (fname, elty) in ((:dtrsv_,:Float64),
                        (:strsv_,:Float32),
                        (:ztrsv_,:ComplexF64),
                        (:ctrsv_,:ComplexF32))
    @eval begin
                #       SUBROUTINE DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
                #       .. Scalar Arguments ..
                #       INTEGER INCX,LDA,N
                #       CHARACTER DIAG,TRANS,UPLO
                #       .. Array Arguments ..
                #       DOUBLE PRECISION A(LDA,*),X(*)
        function trsv!(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            n = checksquare(A)
            if n != length(x)
                throw(DimensionMismatch("size of A is $n != length(x) = $(length(x))"))
            end
            chkstride1(A)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 uplo, trans, diag, n,
                 A, max(1,stride(A,2)), x, stride(x, 1))
            x
        end
        function trsv(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty})
            trsv!(uplo, trans, diag, A, copy(x))
        end
    end
end

### ger

"""
    ger!(alpha, x, y, A)

Rank-1 update of the matrix `A` with vectors `x` and `y` as `alpha*x*y' + A`.
"""
function ger! end

for (fname, elty) in ((:dger_,:Float64),
                      (:sger_,:Float32),
                      (:zgerc_,:ComplexF64),
                      (:cgerc_,:ComplexF32))
    @eval begin
        function ger!(α::$elty, x::AbstractVector{$elty}, y::AbstractVector{$elty}, A::AbstractMatrix{$elty})
            m, n = size(A)
            if m != length(x) || n != length(y)
                throw(DimensionMismatch("A has size ($m,$n), x has length $(length(x)), y has length $(length(y))"))
            end
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                 Ref{BlasInt}),
                 m, n, α, x,
                 stride(x, 1), y, stride(y, 1), A,
                 max(1,stride(A,2)))
            A
        end
    end
end

### syr

"""
    syr!(uplo, alpha, x, A)

Rank-1 update of the symmetric matrix `A` with vector `x` as `alpha*x*transpose(x) + A`.
[`uplo`](@ref stdlib-blas-uplo) controls which triangle of `A` is updated. Returns `A`.
"""
function syr! end

for (fname, elty, lib) in ((:dsyr_,:Float64,libblas),
                           (:ssyr_,:Float32,libblas),
                           (:zsyr_,:ComplexF64,liblapack),
                           (:csyr_,:ComplexF32,liblapack))
    @eval begin
        function syr!(uplo::Char, α::$elty, x::AbstractVector{$elty}, A::AbstractMatrix{$elty})
            n = checksquare(A)
            if length(x) != n
                throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))"))
            end
            ccall((@blasfunc($fname), $lib), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 uplo, n, α, x,
                 stride(x, 1), A, max(1,stride(A, 2)))
            A
        end
    end
end

### her

"""
    her!(uplo, alpha, x, A)

Methods for complex arrays only. Rank-1 update of the Hermitian matrix `A` with vector `x`
as `alpha*x*x' + A`.
[`uplo`](@ref stdlib-blas-uplo) controls which triangle of `A` is updated. Returns `A`.
"""
function her! end

for (fname, elty, relty) in ((:zher_,:ComplexF64, :Float64),
                             (:cher_,:ComplexF32, :Float32))
    @eval begin
        function her!(uplo::Char, α::$relty, x::AbstractVector{$elty}, A::AbstractMatrix{$elty})
            n = checksquare(A)
            if length(x) != n
                throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))"))
            end
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{$relty}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 uplo, n, α, x,
                 stride(x, 1), A, max(1,stride(A,2)))
            A
        end
    end
end

# Level 3
## (GE) general matrix-matrix multiplication

"""
    gemm!(tA, tB, alpha, A, B, beta, C)

Update `C` as `alpha*A*B + beta*C` or the other three variants according to
[`tA`](@ref stdlib-blas-trans) and `tB`. Return the updated `C`.
"""
function gemm! end

for (gemm, elty) in
        ((:dgemm_,:Float64),
         (:sgemm_,:Float32),
         (:zgemm_,:ComplexF64),
         (:cgemm_,:ComplexF32))
    @eval begin
             # SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
             # *     .. Scalar Arguments ..
             #       DOUBLE PRECISION ALPHA,BETA
             #       INTEGER K,LDA,LDB,LDC,M,N
             #       CHARACTER TRANSA,TRANSB
             # *     .. Array Arguments ..
             #       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
        function gemm!(transA::Char, transB::Char, alpha::($elty), A::AbstractVecOrMat{$elty}, B::AbstractVecOrMat{$elty}, beta::($elty), C::AbstractVecOrMat{$elty})
#           if any([stride(A,1), stride(B,1), stride(C,1)] .!= 1)
#               error("gemm!: BLAS module requires contiguous matrix columns")
#           end  # should this be checked on every call?
            m = size(A, transA == 'N' ? 1 : 2)
            ka = size(A, transA == 'N' ? 2 : 1)
            kb = size(B, transB == 'N' ? 1 : 2)
            n = size(B, transB == 'N' ? 2 : 1)
            if ka != kb || m != size(C,1) || n != size(C,2)
                throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
            end
            chkstride1(A)
            chkstride1(B)
            chkstride1(C)
            ccall((@blasfunc($gemm), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}),
                 transA, transB, m, n,
                 ka, alpha, A, max(1,stride(A,2)),
                 B, max(1,stride(B,2)), beta, C,
                 max(1,stride(C,2)))
            C
        end
        function gemm(transA::Char, transB::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            gemm!(transA, transB, alpha, A, B, zero($elty), similar(B, $elty, (size(A, transA == 'N' ? 1 : 2), size(B, transB == 'N' ? 2 : 1))))
        end
        function gemm(transA::Char, transB::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            gemm(transA, transB, one($elty), A, B)
        end
    end
end

"""
    gemm(tA, tB, alpha, A, B)

Return `alpha*A*B` or the other three variants according to [`tA`](@ref stdlib-blas-trans) and `tB`.
"""
gemm(tA, tB, alpha, A, B)

"""
    gemm(tA, tB, A, B)

Return `A*B` or the other three variants according to [`tA`](@ref stdlib-blas-trans) and `tB`.
"""
gemm(tA, tB, A, B)

## (SY) symmetric matrix-matrix and matrix-vector multiplication
for (mfname, elty) in ((:dsymm_,:Float64),
                       (:ssymm_,:Float32),
                       (:zsymm_,:ComplexF64),
                       (:csymm_,:ComplexF32))
    @eval begin
             #     SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
             #     .. Scalar Arguments ..
             #     DOUBLE PRECISION ALPHA,BETA
             #     INTEGER LDA,LDB,LDC,M,N
             #     CHARACTER SIDE,UPLO
             #     .. Array Arguments ..
             #     DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
        function symm!(side::Char, uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, beta::($elty), C::AbstractMatrix{$elty})
            m, n = size(C)
            j = checksquare(A)
            if j != (side == 'L' ? m : n)
                throw(DimensionMismatch("A has size $(size(A)), C has size ($m,$n)"))
            end
            if size(B,2) != n
                throw(DimensionMismatch("B has second dimension $(size(B,2)) but needs to match second dimension of C, $n"))
            end
            chkstride1(A)
            chkstride1(B)
            chkstride1(C)
            ccall((@blasfunc($mfname), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                 Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                 side, uplo, m, n,
                 alpha, A, max(1,stride(A,2)), B,
                 max(1,stride(B,2)), beta, C, max(1,stride(C,2)))
            C
        end
        function symm(side::Char, uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            symm!(side, uplo, alpha, A, B, zero($elty), similar(B))
        end
        function symm(side::Char, uplo::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            symm(side, uplo, one($elty), A, B)
        end
    end
end

"""
    symm(side, ul, alpha, A, B)

Return `alpha*A*B` or `alpha*B*A` according to [`side`](@ref stdlib-blas-side).
`A` is assumed to be symmetric. Only
the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
"""
symm(side, ul, alpha, A, B)

"""
    symm(side, ul, A, B)

Return `A*B` or `B*A` according to [`side`](@ref stdlib-blas-side).
`A` is assumed to be symmetric. Only the [`ul`](@ref stdlib-blas-uplo)
triangle of `A` is used.
"""
symm(side, ul, A, B)

"""
    symm!(side, ul, alpha, A, B, beta, C)

Update `C` as `alpha*A*B + beta*C` or `alpha*B*A + beta*C` according to [`side`](@ref stdlib-blas-side).
`A` is assumed to be symmetric. Only the [`ul`](@ref stdlib-blas-uplo) triangle of
`A` is used. Return the updated `C`.
"""
symm!

## (HE) Hermitian matrix-matrix and matrix-vector multiplication
for (mfname, elty) in ((:zhemm_,:ComplexF64),
                       (:chemm_,:ComplexF32))
    @eval begin
             #     SUBROUTINE DHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
             #     .. Scalar Arguments ..
             #     DOUBLE PRECISION ALPHA,BETA
             #     INTEGER LDA,LDB,LDC,M,N
             #     CHARACTER SIDE,UPLO
             #     .. Array Arguments ..
             #     DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
        function hemm!(side::Char, uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, beta::($elty), C::AbstractMatrix{$elty})
            m, n = size(C)
            j = checksquare(A)
            if j != (side == 'L' ? m : n)
                throw(DimensionMismatch("A has size $(size(A)), C has size ($m,$n)"))
            end
            if size(B,2) != n
                throw(DimensionMismatch("B has second dimension $(size(B,2)) but needs to match second dimension of C, $n"))
            end
            chkstride1(A)
            chkstride1(B)
            chkstride1(C)
            ccall((@blasfunc($mfname), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},
                 Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}),
                 side, uplo, m, n,
                 alpha, A, max(1,stride(A,2)), B,
                 max(1,stride(B,2)), beta, C, max(1,stride(C,2)))
            C
        end
        function hemm(side::Char, uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            hemm!(side, uplo, alpha, A, B, zero($elty), similar(B))
        end
        function hemm(side::Char, uplo::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            hemm(side, uplo, one($elty), A, B)
        end
    end
end

## syrk

"""
    syrk!(uplo, trans, alpha, A, beta, C)

Rank-k update of the symmetric matrix `C` as `alpha*A*transpose(A) + beta*C` or `alpha*transpose(A)*A +
beta*C` according to [`trans`](@ref stdlib-blas-trans).
Only the [`uplo`](@ref stdlib-blas-uplo) triangle of `C` is used. Returns `C`.
"""
function syrk! end

"""
    syrk(uplo, trans, alpha, A)

Returns either the upper triangle or the lower triangle of `A`,
according to [`uplo`](@ref stdlib-blas-uplo),
of `alpha*A*transpose(A)` or `alpha*transpose(A)*A`,
according to [`trans`](@ref stdlib-blas-trans).
"""
function syrk end

for (fname, elty) in ((:dsyrk_,:Float64),
                      (:ssyrk_,:Float32),
                      (:zsyrk_,:ComplexF64),
                      (:csyrk_,:ComplexF32))
   @eval begin
       # SUBROUTINE DSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
       # *     .. Scalar Arguments ..
       #       REAL ALPHA,BETA
       #       INTEGER K,LDA,LDC,N
       #       CHARACTER TRANS,UPLO
       # *     .. Array Arguments ..
       #       REAL A(LDA,*),C(LDC,*)
       function syrk!(uplo::Char, trans::Char,
                      alpha::($elty), A::AbstractVecOrMat{$elty},
                      beta::($elty), C::AbstractMatrix{$elty})
           n = checksquare(C)
           nn = size(A, trans == 'N' ? 1 : 2)
           if nn != n throw(DimensionMismatch("C has size ($n,$n), corresponding dimension of A is $nn")) end
           k  = size(A, trans == 'N' ? 2 : 1)
           chkstride1(A)
           chkstride1(C)
           ccall((@blasfunc($fname), libblas), Cvoid,
                 (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                  Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty},
                  Ptr{$elty}, Ref{BlasInt}),
                 uplo, trans, n, k,
                 alpha, A, max(1,stride(A,2)), beta,
                 C, max(1,stride(C,2)))
            C
        end
    end
end
function syrk(uplo::Char, trans::Char, alpha::Number, A::AbstractVecOrMat)
    T = eltype(A)
    n = size(A, trans == 'N' ? 1 : 2)
    syrk!(uplo, trans, convert(T,alpha), A, zero(T), similar(A, T, (n, n)))
end
syrk(uplo::Char, trans::Char, A::AbstractVecOrMat) = syrk(uplo, trans, one(eltype(A)), A)

"""
    herk!(uplo, trans, alpha, A, beta, C)

Methods for complex arrays only. Rank-k update of the Hermitian matrix `C` as `alpha*A*A' +
beta*C` or `alpha*A'*A + beta*C` according to [`trans`](@ref stdlib-blas-trans).
Only the [`uplo`](@ref stdlib-blas-uplo) triangle of `C` is updated.
Returns `C`.
"""
function herk! end

"""
    herk(uplo, trans, alpha, A)

Methods for complex arrays only.
Returns the [`uplo`](@ref stdlib-blas-uplo) triangle of `alpha*A*A'` or `alpha*A'*A`,
according to [`trans`](@ref stdlib-blas-trans).
"""
function herk end

for (fname, elty, relty) in ((:zherk_, :ComplexF64, :Float64),
                             (:cherk_, :ComplexF32, :Float32))
   @eval begin
       # SUBROUTINE CHERK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
       # *     .. Scalar Arguments ..
       #       REAL ALPHA,BETA
       #       INTEGER K,LDA,LDC,N
       #       CHARACTER TRANS,UPLO
       # *     ..
       # *     .. Array Arguments ..
       #       COMPLEX A(LDA,*),C(LDC,*)
       function herk!(uplo::Char, trans::Char, α::$relty, A::AbstractVecOrMat{$elty},
                      β::$relty, C::AbstractMatrix{$elty})
           n = checksquare(C)
           nn = size(A, trans == 'N' ? 1 : 2)
           if nn != n
               throw(DimensionMismatch("the matrix to update has dimension $n but the implied dimension of the update is $(size(A, trans == 'N' ? 1 : 2))"))
           end
           chkstride1(A)
           chkstride1(C)
           k  = size(A, trans == 'N' ? 2 : 1)
           ccall((@blasfunc($fname), libblas), Cvoid,
                 (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                  Ref{$relty}, Ptr{$elty}, Ref{BlasInt}, Ref{$relty},
                  Ptr{$elty}, Ref{BlasInt}),
                 uplo, trans, n, k,
                 α, A, max(1,stride(A,2)), β,
                 C, max(1,stride(C,2)))
           C
       end
       function herk(uplo::Char, trans::Char, α::$relty, A::AbstractVecOrMat{$elty})
           n = size(A, trans == 'N' ? 1 : 2)
           herk!(uplo, trans, α, A, zero($relty), similar(A, (n,n)))
       end
       herk(uplo::Char, trans::Char, A::AbstractVecOrMat{$elty}) = herk(uplo, trans, one($relty), A)
   end
end

## syr2k
for (fname, elty) in ((:dsyr2k_,:Float64),
                      (:ssyr2k_,:Float32),
                      (:zsyr2k_,:ComplexF64),
                      (:csyr2k_,:ComplexF32))
    @eval begin
             #       SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
             #
             #       .. Scalar Arguments ..
             #       REAL PRECISION ALPHA,BETA
             #       INTEGER K,LDA,LDB,LDC,N
             #       CHARACTER TRANS,UPLO
             #       ..
             #       .. Array Arguments ..
             #       REAL PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
        function syr2k!(uplo::Char, trans::Char,
                        alpha::($elty), A::AbstractVecOrMat{$elty}, B::AbstractVecOrMat{$elty},
                        beta::($elty), C::AbstractMatrix{$elty})
            n = checksquare(C)
            nn = size(A, trans == 'N' ? 1 : 2)
            if nn != n throw(DimensionMismatch("C has size ($n,$n), corresponding dimension of A is $nn")) end
            k  = size(A, trans == 'N' ? 2 : 1)
            chkstride1(A)
            chkstride1(B)
            chkstride1(C)
            ccall((@blasfunc($fname), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ref{BlasInt}),
                 uplo, trans, n, k,
                 alpha, A, max(1,stride(A,2)), B, max(1,stride(B,2)), beta,
                 C, max(1,stride(C,2)))
            C
        end
    end
end
function syr2k(uplo::Char, trans::Char, alpha::Number, A::AbstractVecOrMat, B::AbstractVecOrMat)
    T = eltype(A)
    n = size(A, trans == 'N' ? 1 : 2)
    syr2k!(uplo, trans, convert(T,alpha), A, B, zero(T), similar(A, T, (n, n)))
end
syr2k(uplo::Char, trans::Char, A::AbstractVecOrMat, B::AbstractVecOrMat) = syr2k(uplo, trans, one(eltype(A)), A, B)

for (fname, elty1, elty2) in ((:zher2k_,:ComplexF64,:Float64), (:cher2k_,:ComplexF32,:Float32))
   @eval begin
       # SUBROUTINE CHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
       #
       #       .. Scalar Arguments ..
       #       COMPLEX ALPHA
       #       REAL BETA
       #       INTEGER K,LDA,LDB,LDC,N
       #       CHARACTER TRANS,UPLO
       #       ..
       #       .. Array Arguments ..
       #       COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
       function her2k!(uplo::Char, trans::Char, alpha::($elty1),
                       A::AbstractVecOrMat{$elty1}, B::AbstractVecOrMat{$elty1},
                       beta::($elty2), C::AbstractMatrix{$elty1})
           n = checksquare(C)
           nn = size(A, trans == 'N' ? 1 : 2)
           if nn != n throw(DimensionMismatch("C has size ($n,$n), corresponding dimension of A is $nn")) end
           chkstride1(A)
           chkstride1(B)
           chkstride1(C)
           k  = size(A, trans == 'N' ? 2 : 1)
           ccall((@blasfunc($fname), libblas), Cvoid,
                 (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                  Ref{$elty1}, Ptr{$elty1}, Ref{BlasInt}, Ptr{$elty1}, Ref{BlasInt},
                  Ref{$elty2},  Ptr{$elty1}, Ref{BlasInt}),
                 uplo, trans, n, k,
                 alpha, A, max(1,stride(A,2)), B, max(1,stride(B,2)),
                 beta, C, max(1,stride(C,2)))
           C
       end
       function her2k(uplo::Char, trans::Char, alpha::($elty1), A::AbstractVecOrMat{$elty1}, B::AbstractVecOrMat{$elty1})
           n = size(A, trans == 'N' ? 1 : 2)
           her2k!(uplo, trans, alpha, A, B, zero($elty2), similar(A, $elty1, (n,n)))
       end
       her2k(uplo::Char, trans::Char, A::AbstractVecOrMat{$elty1}, B::AbstractVecOrMat{$elty1}) = her2k(uplo, trans, one($elty1), A, B)
   end
end

## (TR) Triangular matrix and vector multiplication and solution

"""
    trmm!(side, ul, tA, dA, alpha, A, B)

Update `B` as `alpha*A*B` or one of the other three variants determined by
[`side`](@ref stdlib-blas-side) and [`tA`](@ref stdlib-blas-trans).
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
[`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or
are assumed to be all ones.
Returns the updated `B`.
"""
function trmm! end

"""
    trmm(side, ul, tA, dA, alpha, A, B)

Returns `alpha*A*B` or one of the other three variants determined by
[`side`](@ref stdlib-blas-side) and [`tA`](@ref stdlib-blas-trans).
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
[`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or
are assumed to be all ones.
"""
function trmm end

"""
    trsm!(side, ul, tA, dA, alpha, A, B)

Overwrite `B` with the solution to `A*X = alpha*B` or one of the other three variants
determined by [`side`](@ref stdlib-blas-side) and [`tA`](@ref stdlib-blas-trans).
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
[`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or
are assumed to be all ones.
Returns the updated `B`.
"""
function trsm! end

"""
    trsm(side, ul, tA, dA, alpha, A, B)

Return the solution to `A*X = alpha*B` or one of the other three variants determined by
determined by [`side`](@ref stdlib-blas-side) and [`tA`](@ref stdlib-blas-trans).
Only the [`ul`](@ref stdlib-blas-uplo) triangle of `A` is used.
[`dA`](@ref stdlib-blas-diag) determines if the diagonal values are read or
are assumed to be all ones.
"""
function trsm end

for (mmname, smname, elty) in
        ((:dtrmm_,:dtrsm_,:Float64),
         (:strmm_,:strsm_,:Float32),
         (:ztrmm_,:ztrsm_,:ComplexF64),
         (:ctrmm_,:ctrsm_,:ComplexF32))
    @eval begin
        #       SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
        # *     .. Scalar Arguments ..
        #       DOUBLE PRECISION ALPHA
        #       INTEGER LDA,LDB,M,N
        #       CHARACTER DIAG,SIDE,TRANSA,UPLO
        # *     .. Array Arguments ..
        #       DOUBLE PRECISION A(LDA,*),B(LDB,*)
        function trmm!(side::Char, uplo::Char, transa::Char, diag::Char, alpha::Number,
                       A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            m, n = size(B)
            nA = checksquare(A)
            if nA != (side == 'L' ? m : n)
                throw(DimensionMismatch("size of A, $(size(A)), doesn't match $side size of B with dims, $(size(B))"))
            end
            chkstride1(A)
            chkstride1(B)
            ccall((@blasfunc($mmname), libblas), Cvoid,
                  (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                   Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                  side, uplo, transa, diag, m, n,
                  alpha, A, max(1,stride(A,2)), B, max(1,stride(B,2)))
            B
        end
        function trmm(side::Char, uplo::Char, transa::Char, diag::Char,
                      alpha::$elty, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            trmm!(side, uplo, transa, diag, alpha, A, copy(B))
        end
        #       SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
        # *     .. Scalar Arguments ..
        #       DOUBLE PRECISION ALPHA
        #       INTEGER LDA,LDB,M,N
        #       CHARACTER DIAG,SIDE,TRANSA,UPLO
        # *     .. Array Arguments ..
        #       DOUBLE PRECISION A(LDA,*),B(LDB,*)
        function trsm!(side::Char, uplo::Char, transa::Char, diag::Char,
                       alpha::$elty, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            m, n = size(B)
            k = checksquare(A)
            if k != (side == 'L' ? m : n)
                throw(DimensionMismatch("size of A is ($k,$k), size of B is ($m,$n), side is $side, and transa='$transa'"))
            end
            chkstride1(A)
            chkstride1(B)
            ccall((@blasfunc($smname), libblas), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
                 Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}),
                 side, uplo, transa, diag,
                 m, n, alpha, A,
                 max(1,stride(A,2)), B, max(1,stride(B,2)))
            B
        end
        function trsm(side::Char, uplo::Char, transa::Char, diag::Char, alpha::$elty, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty})
            trsm!(side, uplo, transa, diag, alpha, A, copy(B))
        end
    end
end

end # module

function copyto!(dest::Array{T}, rdest::Union{UnitRange{Ti},AbstractRange{Ti}},
                 src::Array{T}, rsrc::Union{UnitRange{Ti},AbstractRange{Ti}}) where {T<:BlasFloat,Ti<:Integer}
    if minimum(rdest) < 1 || maximum(rdest) > length(dest)
        throw(ArgumentError("range out of bounds for dest, of length $(length(dest))"))
    end
    if minimum(rsrc) < 1 || maximum(rsrc) > length(src)
        throw(ArgumentError("range out of bounds for src, of length $(length(src))"))
    end
    if length(rdest) != length(rsrc)
        throw(DimensionMismatch("ranges must be of the same length"))
    end
    GC.@preserve src dest BLAS.blascopy!(length(rsrc),
                                              pointer(src) + (first(rsrc) - 1) * sizeof(T),
                                              step(rsrc),
                                              pointer(dest) + (first(rdest) - 1) * sizeof(T),
                                              step(rdest))
    dest
end

ext: jl

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
### Fetch required packages and reaction networks ###
using Catalyst, DiffEqJump, Random, Statistics, Test
include("test_networks.jl")

### Compares to the manually calcualted function ###
identical_networks = Vector{Pair}()

rate_1_1(u,p,t) = p[1]
rate_1_2(u,p,t) = p[2]*u[1]
rate_1_3(u,p,t) = p[3]*u[2]
rate_1_4(u,p,t) = p[4]*u[2]
rate_1_5(u,p,t) = p[5]*u[3]
rate_1_6(u,p,t) = p[6]*u[3]
rate_1_7(u,p,t) = p[7]*u[4]
rate_1_8(u,p,t) = p[8]*u[4]
affect_1_1!(integrator) = (integrator.u[1] += 1;)
affect_1_2!(integrator) = (integrator.u[1] -= 1; integrator.u[2] += 1;)
affect_1_3!(integrator) = (integrator.u[2] -= 1; integrator.u[1] += 1)
affect_1_4!(integrator) = (integrator.u[2] -= 1; integrator.u[3] += 1)
affect_1_5!(integrator) = (integrator.u[3] -= 1; integrator.u[2] += 1)
affect_1_6!(integrator) = (integrator.u[3] -= 1; integrator.u[4] += 1)
affect_1_7!(integrator) = (integrator.u[4] -= 1; integrator.u[3] += 1)
affect_1_8!(integrator) = (integrator.u[4] -= 1;)
jump_1_1 = ConstantRateJump(rate_1_1,affect_1_1!)
jump_1_2 = ConstantRateJump(rate_1_2,affect_1_2!)
jump_1_3 = ConstantRateJump(rate_1_3,affect_1_3!)
jump_1_4 = ConstantRateJump(rate_1_4,affect_1_4!)
jump_1_5 = ConstantRateJump(rate_1_5,affect_1_5!)
jump_1_6 = ConstantRateJump(rate_1_6,affect_1_6!)
jump_1_7 = ConstantRateJump(rate_1_7,affect_1_7!)
jump_1_8 = ConstantRateJump(rate_1_8,affect_1_8!)
jumps_1 = (jump_1_1,jump_1_2,jump_1_3,jump_1_4,jump_1_5,jump_1_6,jump_1_7,jump_1_8)
push!(identical_networks, reaction_networks_standard[6] => jumps_1)

rate_2_1(u,p,t) = p[1]/10 + u[1]^p[3]/(u[1]^p[3] + p[2]^p[3])
rate_2_2(u,p,t) = p[4]*u[1]*u[2]
rate_2_3(u,p,t) = p[5]*u[3]
rate_2_4(u,p,t) = p[6]*u[3]
rate_2_5(u,p,t) = p[7]*u[1]
rate_2_6(u,p,t) = p[7]*u[2]
rate_2_7(u,p,t) = p[7]*u[3]
affect_2_1!(integrator) = (integrator.u[1] += 1; integrator.u[2] += 1;)
affect_2_2!(integrator) = (integrator.u[1] -= 1; integrator.u[2] -= 1; integrator.u[3] += 1;)
affect_2_3!(integrator) = (integrator.u[1] += 1; integrator.u[2] += 1; integrator.u[3] -= 1;)
affect_2_4!(integrator) = (integrator.u[3] -= 1; integrator.u[1] += 1;)
affect_2_5!(integrator) = (integrator.u[1] -= 1;)
affect_2_6!(integrator) = (integrator.u[2] -= 1;)
affect_2_7!(integrator) = (integrator.u[3] -= 1;)
jump_2_1 = ConstantRateJump(rate_2_1,affect_2_1!)
jump_2_2 = ConstantRateJump(rate_2_2,affect_2_2!)
jump_2_3 = ConstantRateJump(rate_2_3,affect_2_3!)
jump_2_4 = ConstantRateJump(rate_2_4,affect_2_4!)
jump_2_5 = ConstantRateJump(rate_2_5,affect_2_5!)
jump_2_6 = ConstantRateJump(rate_2_6,affect_2_6!)
jump_2_7 = ConstantRateJump(rate_2_7,affect_2_7!)
jumps_2 = (jump_2_1,jump_2_2,jump_2_3,jump_2_4,jump_2_5,jump_2_6,jump_2_7)
push!(identical_networks, reaction_networks_hill[7] => jumps_2)

rate_3_1(u,p,t) = p[1]*binomial(u[1],1)
rate_3_2(u,p,t) = p[2]*binomial(u[2],2)
rate_3_3(u,p,t) = p[3]*binomial(u[2],2)
rate_3_4(u,p,t) = p[4]*binomial(u[3],3)
rate_3_5(u,p,t) = p[5]*binomial(u[3],3)
rate_3_6(u,p,t) = p[6]*binomial(u[4],4)
affect_3_1!(integrator) = (integrator.u[1] -=1; integrator.u[2] +=2;)
affect_3_2!(integrator) = (integrator.u[2] -=2; integrator.u[1] +=1;)
affect_3_3!(integrator) = (integrator.u[2] -=2; integrator.u[3] +=3;)
affect_3_4!(integrator) = (integrator.u[3] -=3; integrator.u[2] +=2;)
affect_3_5!(integrator) = (integrator.u[3] -=3; integrator.u[4] +=4;)
affect_3_6!(integrator) = (integrator.u[4] -=4; integrator.u[3] +=3;)
jump_3_1 = ConstantRateJump(rate_3_1,affect_3_1!)
jump_3_2 = ConstantRateJump(rate_3_2,affect_3_2!)
jump_3_3 = ConstantRateJump(rate_3_3,affect_3_3!)
jump_3_4 = ConstantRateJump(rate_3_4,affect_3_4!)
jump_3_5 = ConstantRateJump(rate_3_5,affect_3_5!)
jump_3_6 = ConstantRateJump(rate_3_6,affect_3_6!)
jumps_3 = (jump_3_1,jump_3_2,jump_3_3,jump_3_4,jump_3_5,jump_3_6)
push!(identical_networks, reaction_networks_constraint[5] => jumps_3)

for (i,networks) in enumerate(identical_networks)
    for factor in [1e-2, 1e-1, 1e0, 1e1], repeat = 1:3
        (i==3) && (factor > 1e-1) && continue   # Large numbers seems to crash it.
        u0 = rand(1:Int64(factor*100),length(networks[1].states))
        p = factor*rand(length(networks[1].ps))
        prob1 = JumpProblem(networks[1],DiscreteProblem(networks[1],u0,(0.,1000.),p),Direct())
        sol1 = solve(prob1,SSAStepper())
        prob2 = JumpProblem(DiscreteProblem(u0,(0.,1000.),p),Direct(),networks[2]...)
        sol2 = solve(prob2,SSAStepper())
        for i = 1:length(u0)
            vals1 = getindex.(sol1.u,i);
            vals2 = getindex.(sol1.u,i);
            (mean(vals2)>0.001) && @test 0.8 < mean(vals1)/mean(vals2) < 1.25
            (std(vals2)>0.001) && @test 0.8 < std(vals1)/std(vals2) < 1.25
        end
    end
end

### Tries solving a large number of problem, ensuring there are no errors. ###
for (i,network) in enumerate(reaction_networks_all)
    (i%5 == 0) && println("Iteration "*string(i)*" at line 102 in file solve_jumps.jl")
    for factor in [1e-1, 1e0, 1e1]
        u0 = rand(1:Int64(factor*100),length(network.states))
        p = factor*rand(length(network.ps))
        prob = JumpProblem(network,DiscreteProblem(network,u0,(0.,1.),p),Direct())
        @test solve(prob,SSAStepper()).retcode == :Default
    end
end

no_param_network = @reaction_network begin
    (1.2,5), X1 ↔ X2
end
for factor in [1e1, 1e2]
    u0 = rand(1:Int64(factor*100),length(no_param_network.states))
    prob = JumpProblem(no_param_network,DiscreteProblem(no_param_network,u0,(0.,1000.)),Direct())
    sol = solve(prob,SSAStepper())
    vals1 = getindex.(sol.u[1:end],1)
    vals2 = getindex.(sol.u[1:end],2)
    @test mean(vals1) > mean(vals2)
end


text:
module PackageEvaluator

# Evaluate the package itself
include("package.jl")

# Evaluate the metadata entry
include("metadata.jl")

# Score it (scorePkg)
include("score.jl")

# Output it (featuresToJSON, getDetailsString)
include("output.jl")

###############################################################################
# evalPkg
# Performs all tests on a single package.
# Input:
#   pkg                 Name of the package (without .jl)
#   addremove=true      If true, adds then removes the package.
#                       If false, the package needs to be in .julia
#                       and in .julia/METADATA/
# Output:
#   features            A dictionary of test results
export evalPkg
function evalPkg(pkg, addremove=true)
  # Need to add Pkg first
  if addremove
    Pkg.add(pkg)
  end

  features = Dict{Symbol,Any}()

  # Package
  pkg_path = Pkg.dir(pkg)
  checkREQUIRE(features, pkg_path)
  checkLicense(features, pkg_path)
  checkTesting(features, pkg_path, pkg)
  getInfo(features, pkg_path)

  # Metadata
  metadata_path = joinpath(ENV["HOME"],".julia","METADATA",pkg)
  checkURL(features, metadata_path)
  checkDesc(features, metadata_path)
  checkRequire(features, metadata_path)

  # Remove Pkg if necessary
  if addremove
    Pkg.rm(pkg)
  end

  # Return a dictionary of test results
  return features
end

end


text:
"""
    standardMutation(genotype::BinaryGenotype, rng::Random.AbstractRNG,
                     probability::AbstractFloat = 0.1)::BinaryGenotype

Mutates a **binary individual**, modifying each gene with a `probability`.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `genotype::BinaryGenotype`: genotype that will be mutated.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
- `probability::AbstractFloat`: lower limit for random generation.
"""
function standardMutation(genotype::BinaryGenotype, rng::Random.AbstractRNG,
                          probability::AbstractFloat)::BinaryGenotype

    indRep = genotype._representation
    genLen = length(indRep)
    mutatedIndRep = Array{Bool}(undef, genLen)

    for i=eachindex(mutatedIndRep)
        random = rand(rng)
        if random < probability
            mutatedIndRep[i] = !indRep[i]
        else
            mutatedIndRep[i] = indRep[i]
        end
    end
    return BinaryGenotype(mutatedIndRep)
end # function

"""
    standardMutation(genotype::BinaryGenotype, rng::Random.AbstractRNG,
                     nGens::Integer=1)::BinaryGenotype

Mutates a **binary individual**, modifying a fixed set of genes.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `genotype::BinaryGenotype`: genotype that will be mutated.
- `rng::Random.AbstractRNG`: random number generator for random number consistency along an experiment.

# User Arguments
- `nGens::Integer=1`: number of genes of the genotype that will be mutated.
"""
function standardMutation(genotype::BinaryGenotype, rng::Random.AbstractRNG,
                          nGens::Integer=1)::BinaryGenotype

    indRep = genotype._representation
    genLen = length(indRep)
    mutatedIndRep = Array{Bool}(undef, genLen)

    if nGens>50
        indexes = randomIndexSelection2(genLen, nGens, rng)
    else
        indexes = randomIndexSelection(genLen, nGens, rng)
    end

    mutatedIndRep = copy(indRep)

    for i=1:nGens
        mutatedIndRep[i] = !indRep[indexes[i]]
    end

    return BinaryGenotype(mutatedIndRep)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using GaussianProcesses
using StaticArrays
using ConstrainedDynamics

mutable struct MDCache
    key::AbstractArray
    data::AbstractArray

    MDCache() = new(Vector{Float64}(), Vector{Float64}())
end

struct MeanDynamics <: GaussianProcesses.Mean
    mechanism::Mechanism
    getμ::Function
    μID::Int
    cache::MDCache
    xtransform::Function

    function MeanDynamics(mechanism::Mechanism, getμ::Function, μID::Int, cache::MDCache; xtransform=(x, _) -> x)
        new(mechanism, getμ, μID, cache, xtransform)
    end

    function MeanDynamics(mechanism::Mechanism, getμ::Function, μID::Int; xtransform=(x, _) -> x)
        new(mechanism, getμ, μID, MDCache(), xtransform)
    end
end

GaussianProcesses.num_params(::MeanDynamics) = 0
GaussianProcesses.grad_mean(::MeanDynamics, ::AbstractVector) = Float64[]
GaussianProcesses.get_params(::MeanDynamics) = Float64[]
GaussianProcesses.get_param_names(::MeanDynamics) = Symbol[]

function GaussianProcesses.set_params!(::MeanDynamics, hyp::AbstractVector)
    length(hyp) == 0 || throw(ArgumentError("Mean dynamics function has no parameters"))
end

function GaussianProcesses.mean(mDynamics::MeanDynamics, x::AbstractVector)
    if mDynamics.cache.key != x  # Cache is invalid
        mDynamics.cache.key = x  # Set cache to current input state
        mechanism = mDynamics.mechanism
        oldstates = getStates(mechanism)
        x = mDynamics.xtransform(x, mechanism)
        setstates!(mechanism, CState(x))
        ConstrainedDynamics.newton!(mechanism)
        mDynamics.cache.data = mDynamics.getμ(mechanism)  # Set cache data to result
        for (id, state) in enumerate(oldstates)
             mechanism.bodies[id].state = state  # Reset mechanism to default values
        end
    end
    return mDynamics.cache.data[mDynamics.μID]
end

function getμ(ids::AbstractArray{Int,1})
    _getμ(mech) = return CState(mech, usesolution=true)[ids]
    return _getμ
end


text:
examples = ["UTURUniChannelCSK_example"]

println("Running tests:")

for t in examples
    test_fn = joinpath("examples", "$t.jl")
    println("  * $test_fn *")
    include(test_fn)
end


text:
struct SimpleLink{T<:Integer} <: AbstractLink{T}
    upn::T # upnode
    dwn::T # downnode
    props::Union{Nothing,Dict}
end

function SimpleLink{T}(data::Dict; upnkey, dwnkey) where T <: Integer
    (upn = data[upnkey]; delete!(data, upnkey))
    (dwn = data[dwnkey]; delete!(data, dwnkey))
    SimpleLink{T}(upn, dwn, data)
end

SimpleLink{T}(dfrow::DataFrameRow, args...; kwargs...) where T <: Integer = SimpleLink{T}(Dict(pairs(dfrow)); upnkey=kwargs[:upnkey], dwnkey=kwargs[:dwnkey])

# Accessors
upn(e::SimpleLink) = e.upn
dwn(e::SimpleLink) = e.dwn
props(e::SimpleLink) = e.props


text:
module VulkanApp

using QML

function main(qml_file)
    load(qml_file)
    exec()
end

export main

end


text:
# ---
# title: 743. Network Delay Time
# id: problem743
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Heap, Depth-first Search, Breadth-first Search, Graph
# link: <https://leetcode.com/problems/network-delay-time/description/>
# hidden: true
# ---
#
# There are `N` network nodes, labelled `1` to `N`.
#
# Given `times`, a list of travel times as **directed** edges `times[i] = (u, v,
# w)`, where `u` is the source node, `v` is the target node, and `w` is the time
# it takes for a signal to travel from source to target.
#
# Now, we send a signal from a certain node `K`. How long will it take for all
# nodes to receive the signal? If it is impossible, return `-1`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/05/23/931_example_1.png)
#
#
#
#     Input: times = [[2,1,1],[2,3,1],[3,4,1]], N = 4, K = 2
#     Output: 2
#
#
#
#
# **Note:**
#
#   1. `N` will be in the range `[1, 100]`.
#   2. `K` will be in the range `[1, N]`.
#   3. The length of `times` will be in the range `[1, 6000]`.
#   4. All edges `times[i] = (u, v, w)` will have `1 <= u, v <= N` and `0 <= w <= 100`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
importall Bukdu
importall Bukdu.Octo
importall Bukdu.Tag

type UserController <: ApplicationController
end

type User
    name::String
    age::Int
    description::String
    happiness::Float64
end

user = User("tom", 20, "", 0.5)

function post_result(c::UserController)
    change(c, user)
end

function test_form(changes::Assoc)
    form = Changeset(user, changes)
    form_for(form, action=post_result, method=post) do f
"""

    Name: $(text_input(f, :name))

    Age: $(select(f, :age, 18:20))

    Happiness: $(text_input(f, :happiness))

$(submit("Submit"))
"""
    end
end

layout(::Layout, body) = "
$body
"
function index(::UserController)
    contents = test_form(Assoc(name="foo bar"))
    render(HTML/Layout, contents)
end

Router() do
    get("/", UserController, index)
    post("/post_result", UserController, post_result)
end

using Base.Test

contents = test_form(Assoc(name="foo bar"))

@test """
<form action="/post_result" method="post" accept-charset="utf-8">

    Name: <input id="user_name" name="user[name]" type="text" value="foo bar" />

    Age: <select id="user_age" name="user[age]">
    <option value="18">18</option>
    <option value="19">19</option>
    <option value="20" selected>20</option>
</select>

    Happiness: <input id="user_happiness" name="user[happiness]" type="text" value="0.5" />

<input type="submit" value="Submit" />
</form>""" == contents

conn = (Router)(get, "/")
@test "
$contents
" == conn.resp_body

tom = User("tom",20,"",0.5)

conn = (Router)(post, "/post_result", user_name="jack")
@test Changeset(tom,Assoc(name="jack")) == conn.resp_body

conn = (Router)(post, "/post_result", user_age="20")
@test Changeset(tom,Assoc()) == conn.resp_body

conn = (Router)(post, "/post_result", user_age="19")
@test Changeset(tom,Assoc(age=19)) == conn.resp_body

conn = (Router)(post, "/post_result", user_undefined="undefined")
@test Changeset(tom,Assoc()) == conn.resp_body

form = change(default(User), name="jack")
@test_throws NoRouteError form_for(()->"", form, action=post_result)
@test """
<form method="post" action="/post_result" accept-charset="utf-8">
</form>""" == form_for((f)->"", form, method=post, action=post_result)
@test """
<form id="ex" method="post" action="/post_result" accept-charset="utf-8">
</form>""" == form_for((f)->"", form, id="ex", method=post, action=post_result)
@test """
<form class="ex" method="post" action="/test" accept-charset="utf-8">
</form>""" == form_for((f)->"", form, class="ex", method=post, action="/test")
@test """
<form class="ex" action="/test" method="get" accept-charset="utf-8">
</form>""" == form_for((f)->"", form, class="ex", action="/test")

import Requests: statuscode, text

Bukdu.start(8082)
resp1 = Requests.get("http://localhost:8082/")
@test 200 == statuscode(resp1)
@test "
$contents
" == text(resp1)

resp2 = Requests.post("http://localhost:8082/post_result", data=Dict("user[name]"=>"foo bar"))
@test 200 == statuscode(resp2)
@test """Bukdu.Octo.Changeset(User("tom",20,"",0.5),Bukdu.Octo.Assoc(Tuple{Symbol,Any}[(:name,"foo bar")]))""" == text(resp2)

sleep(0.1)
Bukdu.stop()


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "FXdiv"
version = v"0.0.20181117"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/Maratyszcza/FXdiv.git", "b742d1143724d646cd0f914646f1240eacf5bd73"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/FXdiv
install -Dvm 644 include/fxdiv.h "${includedir}/fxdiv.h"
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("include/fxdiv.h", :fxdiv_h)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    julia_compat="1.6")


text:
abstract type AbstractFeasibility end
abstract type AbstractFeasibilityHandler end

struct IgnoreFeasibility <: AbstractFeasibility end

handle_feasibility(::IgnoreFeasibility) = false
num_cuts(::IgnoreFeasibility) = 0

struct HandleFeasibility{T <: AbstractFloat} <: AbstractFeasibility
    cuts::Vector{SparseFeasibilityCut{T}}

    function HandleFeasibility(::Type{T}) where T <: AbstractFloat
        return new{T}(Vector{SparseFeasibilityCut{T}}())
    end
end

handle_feasibility(::HandleFeasibility) = true
num_cuts(feasibility::HandleFeasibility) = length(feasibility.cuts)