text:
using CUDA.CUSPARSE

using LinearAlgebra
using SparseArrays

@test CUSPARSE.version() isa VersionNumber

m = 25
n = 35
k = 10
blockdim = 5

@testset "util" begin
    x = sprand(m,0.2)
    d_x = CuSparseVector(x)
    @test length(d_x) == m
    @test size(d_x)   == (m, 1)
    @test size(d_x,1) == m
    @test size(d_x,2) == 1
    @test ndims(d_x)  == 1
    CUDA.@allowscalar begin
        @test Array(d_x[:])        == x[:]
        @test d_x[firstindex(d_x)] == x[firstindex(x)]
        @test d_x[div(end, 2)]     == x[div(end, 2)]
        @test d_x[end]             == x[end]
        @test Array(d_x[firstindex(d_x):end]) == x[firstindex(x):end]
    end
    @test_throws BoundsError d_x[firstindex(d_x) - 1]
    @test_throws BoundsError d_x[end + 1]
    @test nnz(d_x)    == nnz(x)
    @test Array(nonzeros(d_x)) == nonzeros(x)
    @test Array(SparseArrays.nonzeroinds(d_x)) == SparseArrays.nonzeroinds(x)
    @test nnz(d_x)    == length(nonzeros(d_x))
    x = sprand(m,n,0.2)
    d_x = CuSparseMatrixCSC(x)
    @test length(d_x) == m*n
    @test size(d_x)   == (m,n)
    @test size(d_x,1) == m
    @test size(d_x,2) == n
    @test size(d_x,3) == 1
    @test ndims(d_x)  == 2
    CUDA.@allowscalar begin
        @test Array(d_x[:])        == x[:]
        @test d_x[firstindex(d_x)] == x[firstindex(x)]
        @test d_x[div(end, 2)]     == x[div(end, 2)]
        @test d_x[end]             == x[end]
        @test d_x[firstindex(d_x), firstindex(d_x)] == x[firstindex(x), firstindex(x)]
        @test d_x[div(end, 2), div(end, 2)]         == x[div(end, 2), div(end, 2)]
        @test d_x[end, end]        == x[end, end]
        @test Array(d_x[firstindex(d_x):end, firstindex(d_x):end]) == x[:, :]
    end
    @test_throws BoundsError d_x[firstindex(d_x) - 1]
    @test_throws BoundsError d_x[end + 1]
    @test_throws BoundsError d_x[firstindex(d_x) - 1, firstindex(d_x) - 1]
    @test_throws BoundsError d_x[end + 1, end + 1]
    @test_throws BoundsError d_x[firstindex(d_x) - 1:end + 1, :]
    @test_throws BoundsError d_x[firstindex(d_x) - 1, :]
    @test_throws BoundsError d_x[end + 1, :]
    @test_throws BoundsError d_x[:, firstindex(d_x) - 1:end + 1]
    @test_throws BoundsError d_x[:, firstindex(d_x) - 1]
    @test_throws BoundsError d_x[:, end + 1]
    @test nnz(d_x)    == nnz(x)
    @test Array(nonzeros(d_x)) == nonzeros(x)
    @test nnz(d_x)    == length(nonzeros(d_x))
    @test !issymmetric(d_x)
    @test !ishermitian(d_x)
    @test_throws ArgumentError size(d_x,0)
    @test_throws ArgumentError CUSPARSE.CuSparseVector(x)
    y = sprand(k,n,0.2)
    d_y = CuSparseMatrixCSC(y)
    @test_throws ArgumentError copyto!(d_y,d_x)
    d_y = CUSPARSE.switch2csr(d_y)
    d_x = CUSPARSE.switch2csr(d_x)
    @test_throws ArgumentError copyto!(d_y,d_x)
    d_y = CUSPARSE.switch2bsr(d_y,convert(Cint,blockdim))
    d_x = CUSPARSE.switch2bsr(d_x,convert(Cint,blockdim))
    @test_throws ArgumentError copyto!(d_y,d_x)
    x = sprand(m,0.2)
    d_x = CuSparseVector(x)
    @test size(d_x, 1) == m
    @test size(d_x, 2) == 1
    @test_throws ArgumentError size(d_x, 0)
    y = sprand(n,0.2)
    d_y = CuSparseVector(y)
    @test_throws ArgumentError copyto!(d_y,d_x)
    x = sprand(m,m,0.2)
    d_x = Symmetric(CuSparseMatrixCSC(x + transpose(x)))
    @test issymmetric(d_x)
    x = sprand(ComplexF64, m, m, 0.2)
    d_x = Hermitian(CuSparseMatrixCSC(x + x'))
    @test ishermitian(d_x)
    x = sprand(m,m,0.2)
    d_x = UpperTriangular(CuSparseMatrixCSC(x))
    @test istriu(d_x)
    @test !istril(d_x)
    d_x = LowerTriangular(CuSparseMatrixCSC(x))
    @test !istriu(d_x)
    @test istril(d_x)
end

@testset "char" begin
    @test_throws ArgumentError CUSPARSE.cusparseop('Z')
    @test_throws ArgumentError CUSPARSE.cusparsetype('Z')
    @test_throws ArgumentError CUSPARSE.cusparsefill('Z')
    @test_throws ArgumentError CUSPARSE.cusparsediag('Z')
    @test_throws ArgumentError CUSPARSE.cusparsedir('Z')
    @test_throws ArgumentError CUSPARSE.cusparseindex('A')

    @test CUSPARSE.cusparseop('N') == CUSPARSE.CUSPARSE_OPERATION_NON_TRANSPOSE
    @test CUSPARSE.cusparseop('C') == CUSPARSE.CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE
    @test CUSPARSE.cusparseop('T') == CUSPARSE.CUSPARSE_OPERATION_TRANSPOSE
    @test CUSPARSE.cusparsefill('U') == CUSPARSE.CUSPARSE_FILL_MODE_UPPER
    @test CUSPARSE.cusparsefill('L') == CUSPARSE.CUSPARSE_FILL_MODE_LOWER
    @test CUSPARSE.cusparsediag('U') == CUSPARSE.CUSPARSE_DIAG_TYPE_UNIT
    @test CUSPARSE.cusparsediag('N') == CUSPARSE.CUSPARSE_DIAG_TYPE_NON_UNIT
    @test CUSPARSE.cusparseindex('Z') == CUSPARSE.CUSPARSE_INDEX_BASE_ZERO
    @test CUSPARSE.cusparseindex('O') == CUSPARSE.CUSPARSE_INDEX_BASE_ONE
end

@testset "conversion" begin
    @testset for elty in [Float32, Float64, ComplexF32, ComplexF64]
        @testset "make_csc" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSC(x)
            h_x = collect(d_x)
            @test h_x == x
            @test eltype(d_x) == elty
        end

        @testset "make_csr" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSR(x)
            h_x = collect(d_x)
            @test h_x == x
        end

        @testset "convert_r2c" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSR(x)
            d_x = CUSPARSE.switch2csc(d_x)
            h_x = collect(d_x)
            @test h_x.rowval == x.rowval
            @test h_x.nzval ≈ x.nzval
        end

        @testset "convert_r2b" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSR(x)
            d_x = CUSPARSE.switch2bsr(d_x,convert(Cint,blockdim))
            d_x = CUSPARSE.switch2csr(d_x)
            h_x = collect(d_x)
            @test h_x ≈ x
        end

        @testset "convert_c2b" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSC(x)
            d_x = CUSPARSE.switch2bsr(d_x,convert(Cint,blockdim))
            d_x = CUSPARSE.switch2csc(d_x)
            h_x = collect(d_x)
            @test h_x ≈ x
        end

        @testset "convert_c2h" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSC(x)
            d_x = CUSPARSE.switch2hyb(d_x)
            d_y = CUSPARSE.switch2csc(d_x)
            CUSPARSE.cusparseDestroyHybMat(d_x.Mat)
            h_x = collect(d_y)
            @test h_x.rowval == x.rowval
            @test h_x.nzval ≈ x.nzval
        end

        @testset "convert_r2h" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSR(x)
            d_x = CUSPARSE.switch2hyb(d_x)
            d_y = CUSPARSE.switch2csr(d_x)
            CUSPARSE.cusparseDestroyHybMat(d_x.Mat)
            h_x = collect(d_y)
            @test h_x.rowval == x.rowval
            @test h_x.nzval ≈ x.nzval
        end

        @testset "convert_d2h" begin
            x = rand(elty,m,n)
            d_x = CuArray(x)
            d_x = CUSPARSE.sparse(d_x,'H')
            d_y = Array(d_x)
            CUSPARSE.cusparseDestroyHybMat(d_x.Mat)
            h_x = collect(d_y)
            @test h_x ≈ x
        end

        @testset "convert_d2b" begin
            x = rand(elty,m,n)
            d_x = CuArray(x)
            d_x = CUSPARSE.sparse(d_x,'B')
            d_y = Array(d_x)
            h_x = collect(d_y)
            @test h_x ≈ x
        end

        @testset "convert_c2r" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSC(x)
            d_x = CUSPARSE.switch2csr(d_x)
            h_x = collect(d_x)
            @test h_x.rowval == x.rowval
            @test h_x.nzval ≈ x.nzval
        end

        @testset "convert_r2d" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSR(x)
            d_x = Array(d_x)
            h_x = collect(d_x)
            @test h_x ≈ Array(x)
        end

        @testset "convert_c2d" begin
            x = sparse(rand(elty,m,n))
            d_x = CuSparseMatrixCSC(x)
            d_x = Array(d_x)
            h_x = collect(d_x)
            @test h_x ≈ Array(x)
        end

        @testset "convert_d2c" begin
            x = rand(elty,m,n)
            d_x = CuArray(x)
            d_x = CUSPARSE.sparse(d_x,'C')
            h_x = collect(d_x)
            @test h_x ≈ sparse(x)
        end

        @testset "convert_d2r" begin
            x = rand(elty,m,n)
            d_x = CuArray(x)
            d_x = CUSPARSE.sparse(d_x)
            h_x = collect(d_x)
            @test h_x ≈ sparse(x)
        end
    end
end
@testset "bsric02" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = rand(elty, m, m)
        A += adjoint(A)
        A += m * Diagonal{elty}(I, m)

        @testset "bsric02!" begin
            d_A = CuSparseMatrixCSR(sparse(tril(A)))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_A = CUSPARSE.ic02!(d_A,'O')
            h_A = collect(CUSPARSE.switch2csr(d_A))
            Ac = sparse(Array(cholesky(Hermitian(A))))
            h_A = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
            d_A = CuSparseMatrixCSR(sparse(tril(rand(elty,m,n))))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.ic02!(d_A,'O')
        end

        @testset "bsric02" begin
            d_A = CuSparseMatrixCSR(sparse(tril(A)))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_B = CUSPARSE.ic02(d_A,'O')
            h_A = collect(CUSPARSE.switch2csr(d_B))
            Ac = sparse(Array(cholesky(Hermitian(A))))
            h_A = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
end

@testset "bsrilu02" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = rand(elty,m,m)
        A += transpose(A)
        A += m * Diagonal{elty}(I, m)
        @testset "bsrilu02!" begin
            d_A = CuSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_A = CUSPARSE.ilu02!(d_A,'O')
            h_A = collect(CUSPARSE.switch2csr(d_A))
            Alu = lu(Array(A), Val(false))
            Ac = sparse(Alu.L*Alu.U)
            h_A = adjoint(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
            d_A = CuSparseMatrixCSR(sparse(rand(elty,m,n)))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.ilu02!(d_A,'O')
        end

        @testset "bsrilu02" begin
            d_A = CuSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_B = CUSPARSE.ilu02(d_A,'O')
            h_A = collect(CUSPARSE.switch2csr(d_B))
            Alu = lu(Array(A),Val(false))
            Ac = sparse(Alu.L*Alu.U)
            h_A = adjoint(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
end

@testset "bsrsm2" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "bsrsm2!" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CuArray(X)
            d_A = CuSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_X = CUSPARSE.bsrsm2!('N','N',alpha,d_A,d_X,'O')
            h_Y = collect(d_X)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_X = CUDA.rand(elty,n,n)
            @test_throws DimensionMismatch CUSPARSE.bsrsm2!('N','N',alpha,d_A,d_X,'O')
            @test_throws DimensionMismatch CUSPARSE.bsrsm2!('N','T',alpha,d_A,d_X,'O')
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSR(A)
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.bsrsm2!('N','N',alpha,d_A,d_X,'O')
        end

        @testset "bsrsm2" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CuArray(X)
            d_A = CuSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_Y = CUSPARSE.bsrsm2('N','N',alpha,d_A,d_X,'O')
            h_Y = collect(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSR(A)
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.bsrsm2('N','N',alpha,d_A,d_X,'O')
        end
    end
end
@testset "bsrsv2" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "bsrsv2!" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CuArray(X)
            d_A = CuSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_X = CUSPARSE.sv2!('N','U',alpha,d_A,d_X,'O')
            h_Y = collect(d_X)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_X = CUDA.rand(elty,n)
            @test_throws DimensionMismatch CUSPARSE.sv2!('N','U',alpha,d_A,d_X,'O')
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSR(A)
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.sv2!('N','U',alpha,d_A,d_X,'O')
        end

        @testset "bsrsv2" begin
            A = rand(elty,m,m)
            A = triu(A)
            Al = tril(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CuArray(X)
            d_A = CuSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_Y = CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
            h_Y = collect(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_Y = UpperTriangular(d_A)\d_X
            h_Y = collect(d_Y)
            @test h_Y ≈ A\X
            #=d_Y = UpperTriangular(d_A)'\d_X
            h_Y = collect(d_Y)
            @test h_Y ≈ A'\X=#
            d_Y = transpose(UpperTriangular(d_A))\d_X
            h_Y = collect(d_Y)
            @test h_Y ≈ transpose(A)\X
            d_Y = LowerTriangular(d_A)\d_X
            h_Y = collect(d_Y)
            @test h_Y ≈ Al\X
            #=d_Y = LowerTriangular(d_A)'\d_X
            h_Y = collect(d_Y)
            @test h_Y ≈ A'\X=#
            d_Y = transpose(LowerTriangular(d_A))\d_X
            h_Y = collect(d_Y)
            @test h_Y ≈ transpose(Al)\X
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSR(A)
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
        end
    end
end

@testset "ilu0" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = rand(elty,m,m)
        @testset "csr" begin
            d_A = CuSparseMatrixCSR(sparse(A))
            info = CUSPARSE.sv_analysis('N','G','U',d_A,'O')
            d_B = CUSPARSE.ilu0('N',d_A,info,'O')
            h_B = collect(d_B)
            Alu = lu(Array(A),Val(false))
            Ac = sparse(Alu.L*Alu.U)
            h_A = adjoint(h_B) * h_B
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
        @testset "csc" begin
            d_A = CuSparseMatrixCSC(sparse(A))
            info = CUSPARSE.sv_analysis('N','G','U',d_A,'O')
            d_B = CUSPARSE.ilu0('N',d_A,info,'O')
            h_B = collect(d_B)
            Alu = lu(Array(A),Val(false))
            Ac = sparse(Alu.L*Alu.U)
            h_A = adjoint(h_B) * h_B
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
end

@testset "ilu02" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "csr" begin
            A = rand(elty,m,m)
            A += transpose(A)
            A += m * Diagonal{elty}(I, m)
            d_A = CuSparseMatrixCSR(sparse(A))
            d_B = CUSPARSE.ilu02(d_A,'O')
            h_A = collect(d_B)
            Alu = lu(Array(A),Val(false))
            Ac = sparse(Alu.L*Alu.U)
            h_A = adjoint(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
        @testset "csc" begin
            A = rand(elty,m,m)
            A += transpose(A)
            A += m * Diagonal{elty}(I, m)
            d_A = CuSparseMatrixCSC(sparse(A))
            d_B = CUSPARSE.ilu02(d_A,'O')
            h_A = collect(d_B)
            Alu = lu(Array(A),Val(false))
            Ac = sparse(Alu.L*Alu.U)
            h_A = adjoint(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
end

@testset "ic0" begin
    @testset for elty in [Float32,Float64]
        @testset "csr" begin
            A    = rand(elty, m ,m)
            A    = A + transpose(A)
            A   += m * Diagonal{elty}(I, m)
            d_A  = Symmetric(CuSparseMatrixCSR(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N', 'S', 'U', d_A, 'O')
            d_B  = CUSPARSE.ic0('N', 'S', d_A, info, 'O')
            h_A  = collect(d_B)
            Ac   = sparse(Array(cholesky(A)))
            h_A  = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
        @testset "csc" begin
            A    = rand(elty, m ,m)
            A    = A + transpose(A)
            A   += m * Diagonal{elty}(I, m)
            d_A  = Symmetric(CuSparseMatrixCSC(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N', 'S', 'U', d_A, 'O')
            d_B  = CUSPARSE.ic0('N', 'S', d_A, info, 'O')
            h_A  = collect(d_B)
            Ac   = sparse(Array(cholesky(A)))
            h_A  = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
    @testset for elty in [ComplexF32,ComplexF64]
        @testset "csr" begin
            A    = rand(elty,m,m)
            A    = A + adjoint(A)
            A   += m * Diagonal{elty}(I, m)
            d_A  = Hermitian(CuSparseMatrixCSR(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N', 'H', 'U', d_A, 'O')
            d_B  = CUSPARSE.ic0('N', 'H', d_A, info, 'O')
            h_A  = collect(d_B)
            Ac   = sparse(Array(cholesky(A)))
            h_A  = adjoint(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
        @testset "csc" begin
            A    = rand(elty,m,m)
            A    = A + adjoint(A)
            A   += m * Diagonal{elty}(I, m)
            d_A  = Hermitian(CuSparseMatrixCSC(sparse(triu(A))))
            info = CUSPARSE.sv_analysis('N', 'H', 'U', d_A, 'O')
            d_B  = CUSPARSE.ic0('N', 'H', d_A, info, 'O')
            h_A  = collect(d_B)
            Ac   = sparse(Array(cholesky(A)))
            h_A  = adjoint(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
end

@testset "ic2" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "csr" begin
            A   = rand(elty, m, m)
            A  += adjoint(A)
            A  += m * Diagonal{elty}(I, m)
            d_A = CuSparseMatrixCSR(sparse(tril(A)))
            d_B = CUSPARSE.ic02(d_A, 'O')
            h_A = collect(d_B)
            Ac  = sparse(Array(cholesky(Hermitian(A))))
            h_A = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
            A   = rand(elty,m,n)
            d_A = CuSparseMatrixCSR(sparse(tril(A)))
            @test_throws DimensionMismatch CUSPARSE.ic02(d_A, 'O')
        end
        @testset "csc" begin
            A   = rand(elty, m, m)
            A  += adjoint(A)
            A  += m * Diagonal{elty}(I, m)
            d_A = CuSparseMatrixCSC(sparse(tril(A)))
            d_B = CUSPARSE.ic02(d_A, 'O')
            h_A = collect(d_B)
            Ac  = sparse(Array(cholesky(Hermitian(A))))
            h_A = transpose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
            A   = rand(elty,m,n)
            d_A = CuSparseMatrixCSC(sparse(tril(A)))
            @test_throws DimensionMismatch CUSPARSE.ic02(d_A, 'O')
        end
    end
end

@testset "csrsm" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "csr" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CuArray(X)
            d_A = CuSparseMatrixCSR(sparse(A))
            info = CUSPARSE.sm_analysis('N','U',d_A,'O')
            d_Y = CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            h_Y = collect(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_y = UpperTriangular(d_A)\d_X
            h_y = collect(d_y)
            y = A\X
            @test y ≈ h_y
            #=d_y = UpperTriangular(d_A)'\d_X
            h_y = collect(d_y)
            y = A'\X=#
            @test y ≈ h_y
            d_y = transpose(UpperTriangular(d_A))\d_X
            h_y = collect(d_y)
            y = transpose(A)\X
            @test y ≈ h_y
            d_X = CUDA.rand(elty,n,n)
            @test_throws DimensionMismatch CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.sm_analysis('T','U',d_A,'O')
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end
        @testset "csc" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CuArray(X)
            d_A = CuSparseMatrixCSC(sparse(A))
            info = CUSPARSE.sm_analysis('N','U',d_A,'O')
            d_Y = CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            h_Y = collect(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_y = UpperTriangular(d_A)\d_X
            h_y = collect(d_y)
            y = A\X
            @test y ≈ h_y
            if elty <: Real
                d_y = UpperTriangular(d_A)'\d_X
                h_y = collect(d_y)
                y = A'\X
                @test y ≈ h_y
            end
            d_y = transpose(UpperTriangular(d_A))\d_X
            h_y = collect(d_y)
            y = transpose(A)\X
            @test y ≈ h_y
            d_X = CUDA.rand(elty,n,n)
            @test_throws DimensionMismatch CUSPARSE.sm_solve('N','U',alpha,d_A,d_X,info,'O')
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.sm_analysis('T','U',d_A,'O')
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end
    end
end

@testset "cssv" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "csrsv" begin
            A = rand(elty,m,m)
            A = triu(A)
            x = rand(elty,m)
            alpha = rand(elty)
            d_x = CuArray(x)
            d_A = CuSparseMatrixCSR(sparse(A))
            d_y = CUSPARSE.sv('N','T','U',alpha,d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\(alpha * x)
            @test y ≈ h_y
            d_A = CuSparseMatrixCSR(sparse(A))
            d_y = CUSPARSE.sv('N','T','U',d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\x
            @test y ≈ h_y
            d_A = UpperTriangular(CuSparseMatrixCSR(sparse(A)))
            d_y = CUSPARSE.sv('N',d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\x
            @test y ≈ h_y
            d_A = CuSparseMatrixCSR(sparse(A))
            x = rand(elty,n)
            d_x = CuArray(x)
            info = CUSPARSE.sv_analysis('N','T','U',d_A,'O')
            @test_throws DimensionMismatch CUSPARSE.sv_solve('N','U',alpha,d_A,d_x,info,'O')
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.sv_analysis('T','T','U',d_A,'O')
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end

        @testset "cscsv" begin
            A = rand(elty,m,m)
            A = triu(A)
            x = rand(elty,m)
            alpha = rand(elty)
            d_x = CuArray(x)
            d_A = CuSparseMatrixCSC(sparse(A))
            d_y = CUSPARSE.sv('N','T','U',alpha,d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\(alpha * x)
            @test y ≈ h_y
            d_x = CuArray(x)
            d_A = CuSparseMatrixCSC(sparse(A))
            d_y = CUSPARSE.sv('N','T','U',d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\x
            @test y ≈ h_y
            d_x = CuArray(x)
            d_A = UpperTriangular(CuSparseMatrixCSC(sparse(A)))
            d_y = CUSPARSE.sv('N',d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\x
            @test y ≈ h_y
            d_A = CuSparseMatrixCSC(sparse(A))
            x = rand(elty,n)
            d_x = CuArray(x)
            info = CUSPARSE.sv_analysis('N','T','U',d_A,'O')
            @test_throws DimensionMismatch CUSPARSE.sv_solve('N','U',alpha,d_A,d_x,info,'O')
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.sv_analysis('T','T','U',d_A,'O')
            CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
        end

        @testset "csrsv2" begin
            A = rand(elty,m,m)
            A = triu(A)
            Al = tril(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CuArray(X)
            d_A = CuSparseMatrixCSR(sparse(A))
            d_Y = CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
            h_Y = collect(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_y = UpperTriangular(d_A)\d_X
            h_y = collect(d_y)
            y = A\X
            @test y ≈ h_y
            d_y = transpose(UpperTriangular(d_A))\d_X
            h_y = collect(d_y)
            y = transpose(A)\X
            @test y ≈ h_y
            #=d_y = UpperTriangular(d_A)'\d_X
            h_y = collect(d_y)
            y = A'\X
            @test y ≈ h_y=#
            d_y = LowerTriangular(d_A)\d_X
            h_y = collect(d_y)
            y = Al\X
            @test y ≈ h_y
            d_y = transpose(LowerTriangular(d_A))\d_X
            h_y = collect(d_y)
            y = transpose(Al)\X
            @test y ≈ h_y
            #=d_y = LowerTriangular(d_A)'\d_X
            h_y = collect(d_y)
            y = A'\X
            @test y ≈ h_y=#
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
        end

        @testset "cscsv2" begin
            A = rand(elty,m,m)
            A = triu(A)
            Al = tril(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CuArray(X)
            d_A = CuSparseMatrixCSC(sparse(A))
            d_Y = CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
            h_Y = collect(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_Y = CUSPARSE.sv2('T','U',alpha,d_A,d_X,'O')
            h_Y = collect(d_Y)
            Y = transpose(A)\(alpha * X)
            @test Y ≈ h_Y
            d_y = UpperTriangular(d_A)\d_X
            h_y = collect(d_y)
            y = A\X
            @test y ≈ h_y
            d_y = transpose(UpperTriangular(d_A))\d_X
            h_y = collect(d_y)
            y = transpose(A)\X
            @test y ≈ h_y
            d_y = LowerTriangular(d_A)\d_X
            h_y = collect(d_y)
            y = Al\X
            @test y ≈ h_y
            d_y = transpose(LowerTriangular(d_A))\d_X
            h_y = collect(d_y)
            y = transpose(Al)\X
            @test y ≈ h_y
            #=d_y = UpperTriangular(d_A)'\d_X
            h_y = collect(d_y)
            y = A'\X
            @test y ≈ h_y=#
            # shouldn't work for now bc sv2 has no way to do conj...
            A = sparse(rand(elty,m,n))
            d_A = CuSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
        end
    end
end

@testset "cssm" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "csrsm" begin
            A = rand(elty,m, m)
            A = triu(A)
            x = rand(elty,m, n)
            alpha = rand(elty)
            d_x = CuArray(x)
            d_A = CuSparseMatrixCSR(sparse(A))
            d_y = CUSPARSE.sm('N','U',alpha,d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\(alpha * x)
            @test y ≈ h_y
            d_A = CuSparseMatrixCSR(sparse(A))
            d_y = CUSPARSE.sm('N','U',d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\x
            @test y ≈ h_y
            d_A = UpperTriangular(CuSparseMatrixCSR(sparse(A)))
            d_y = CUSPARSE.sm('N',alpha,d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\(alpha*x)
            @test y ≈ h_y
        end
        @testset "cscsm" begin
            A = rand(elty,m, m)
            A = triu(A)
            x = rand(elty,m, n)
            alpha = rand(elty)
            d_x = CuArray(x)
            d_A = CuSparseMatrixCSC(sparse(A))
            d_y = CUSPARSE.sm('N','U',alpha,d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\(alpha * x)
            @test y ≈ h_y
            d_A = CuSparseMatrixCSC(sparse(A))
            d_y = CUSPARSE.sm('N','U',d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\x
            @test y ≈ h_y
            d_A = UpperTriangular(CuSparseMatrixCSC(sparse(A)))
            d_y = CUSPARSE.sm('N',alpha,d_A,d_x,'O')
            h_y = collect(d_y)
            y = A\(alpha*x)
            @test y ≈ h_y
        end
    end
end
@testset "doti" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        x = sparsevec(rand(1:m,k), rand(elty,k), m)
        y = rand(elty,m)
        d_x = CuSparseVector(x)
        d_y = CuArray(y)
        ddot = CUSPARSE.doti(d_x,d_y,'O')
        #compare
        dot = zero(elty)
        for i in 1:length(x.nzval)
            dot += x.nzval[i] * y[x.nzind[i]]
        end
        @test ddot ≈ dot
    end
end

@testset "dotci" begin
    @testset for elty in [ComplexF32,ComplexF64]
        x = sparsevec(rand(1:m,k), rand(elty,k), m)
        y = rand(elty,m)
        d_x = CuSparseVector(x)
        d_y = CuArray(y)
        ddot = CUSPARSE.dotci(d_x,d_y,'O')
        #compare
        dot = zero(elty)
        for i in 1:length(x.nzval)
            dot += conj(x.nzval[i]) * y[x.nzind[i]]
        end
        @test ddot ≈ dot
    end
end

@testset "axpyi" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "axpyi!" begin
            x = sparsevec(rand(1:m,k), rand(elty,k), m)
            y = rand(elty,m)
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            alpha = rand(elty)
            d_y = CUSPARSE.axpyi!(alpha,d_x,d_y,'O')
            #compare
            h_y = collect(d_y)
            y[x.nzind] += alpha * x.nzval
            @test h_y ≈ y
        end

        @testset "axpyi" begin
            x = sparsevec(rand(1:m,k), rand(elty,k), m)
            y = rand(elty,m)
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            alpha = rand(elty)
            d_z = CUSPARSE.axpyi(alpha,d_x,d_y,'O')
            #compare
            h_z = collect(d_z)
            z = copy(y)
            z[x.nzind] += alpha * x.nzval
            @test h_z ≈ z
            d_z = CUSPARSE.axpyi(d_x,d_y,'O')
            #compare
            h_z = collect(d_z)
            z = copy(y)
            z[x.nzind] += x.nzval
            @test h_z ≈ z
        end
    end
end

@testset "geam" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,n))
        B = sparse(rand(elty,m,n))
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_A = CuSparseMatrixCSR(A)
            d_B = CuSparseMatrixCSR(B)
            C = alpha * A + beta * B
            d_C = CUSPARSE.geam(alpha,d_A,beta,d_B,'O','O','O')
            h_C = collect(d_C)
            @test C ≈ h_C
            d_C = CUSPARSE.geam(d_A,beta,d_B,'O','O','O')
            h_C = collect(d_C)
            C = A + beta * B
            @test C ≈ h_C
            d_C = CUSPARSE.geam(alpha,d_A,d_B,'O','O','O')
            h_C = collect(d_C)
            C = alpha * A + B
            @test C ≈ h_C
            d_C = CUSPARSE.geam(d_A,d_B,'O','O','O')
            h_C = collect(d_C)
            C = A + B
            @test C ≈ h_C
            d_C = d_A + d_B
            h_C = collect(d_C)
            C = A + B
            @test C ≈ h_C
            d_C = d_A - d_B
            h_C = collect(d_C)
            C = A - B
            @test C ≈ h_C
            B_ = sparse(rand(elty,k,n))
            d_B = CuSparseMatrixCSR(B_)
            @test_throws DimensionMismatch CUSPARSE.geam(d_B,d_A,'O','O','O')
        end
        @testset "csc" begin
            d_A = CuSparseMatrixCSC(A)
            #=d_B = CuSparseMatrixCSC(B)
            C = alpha * A + beta * B
            d_C = CUSPARSE.geam(alpha,d_A,beta,d_B,'O','O','O')
            h_C = transpose!(collect(d_C))
            @test C ≈ h_C
            d_C = CUSPARSE.geam(d_A,beta,d_B,'O','O','O')
            h_C = transpose!(collect(d_C))
            C = A + beta * B
            @test C ≈ h_C
            d_C = CUSPARSE.geam(alpha,d_A,d_B,'O','O','O')
            h_C = transpose!(collect(d_C))
            C = alpha * A + B
            @test C ≈ h_C
            d_C = CUSPARSE.geam(d_A,d_B,'O','O','O')
            h_C = transpose!(collect(d_C))
            C = A + B
            @test C ≈ h_C
            d_C = d_A + d_B
            h_C = transpose!(collect(d_C))
            C = A + B
            @test C ≈ h_C
            d_C = d_A - d_B
            h_C = transpose!(collect(d_C))
            C = A - B
            @test C ≈ h_C=#
            B_ = sparse(rand(elty,k,n))
            d_B = CuSparseMatrixCSC(B_)
            @test_throws DimensionMismatch CUSPARSE.geam(d_B,d_A,'O','O','O')
        end
        @testset "mixed" begin
            A = spdiagm(0=>rand(elty,m))
            B = spdiagm(0=>rand(elty,m)) + sparse(diagm(1=>rand(elty,m-1)))
            d_A = CuSparseMatrixCSR(A)
            d_B = CuSparseMatrixCSC(B)
            d_C = d_B + d_A
            h_C = collect(d_C)
            @test h_C ≈ A + B
            d_C = d_A + d_B
            h_C = collect(d_C)
            @test h_C ≈ A + B
            d_C = d_A - d_B
            h_C = collect(d_C)
            @test h_C ≈ A - B
            d_C = d_B - d_A
            h_C = collect(d_C)
            @test h_C ≈ B - A
        end
    end
end

@testset "gemm" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        # CSR
        A = sparse(rand(elty,m,k)); d_A = CuSparseMatrixCSR(A)
        B = sparse(rand(elty,k,n)); d_B = CuSparseMatrixCSR(B)
        d_C = CUSPARSE.gemm('N','N',d_A,d_B,'O','O','O')
        @test A * B ≈ collect(d_C)
        #
        @test_throws DimensionMismatch CUSPARSE.gemm('N','T',d_A,d_B,'O','O','O')
        @test_throws DimensionMismatch CUSPARSE.gemm('T','T',d_A,d_B,'O','O','O')
        @test_throws DimensionMismatch CUSPARSE.gemm('T','N',d_A,d_B,'O','O','O')
        @test_throws DimensionMismatch CUSPARSE.gemm('N','N',d_B,d_A,'O','O','O')
        #
        A = sparse(rand(elty,m,k)); d_A = CuSparseMatrixCSR(A)
        B = sparse(rand(elty,n,k)); d_B = CuSparseMatrixCSR(B)
        d_C = CUSPARSE.gemm('N','T',d_A,d_B,'O','O','O')
        @test A * transpose(B) ≈ collect(d_C)

        # CSC
        A = sparse(rand(elty,m,k)); d_A = CuSparseMatrixCSC(A)
        B = sparse(rand(elty,k,n)); d_B = CuSparseMatrixCSC(B)
        d_C = CUSPARSE.gemm('N','N',d_A,d_B,'O','O','O')
        @test A * B ≈ collect(d_C)
        #
        @test_throws(DimensionMismatch,CUSPARSE.gemm('N','T',d_A,d_B,'O','O','O'))
        @test_throws(DimensionMismatch,CUSPARSE.gemm('T','T',d_A,d_B,'O','O','O'))
        @test_throws(DimensionMismatch,CUSPARSE.gemm('T','N',d_A,d_B,'O','O','O'))
        @test_throws(DimensionMismatch,CUSPARSE.gemm('N','N',d_B,d_A,'O','O','O'))
        #
        A = sparse(rand(elty,m,k)); d_A = CuSparseMatrixCSC(A)
        B = sparse(rand(elty,n,k)); d_B = CuSparseMatrixCSC(B)
        d_C = CUSPARSE.gemm('N','T',d_A,d_B,'O','O','O')
        @test A * transpose(B) ≈ collect(d_C)
    end
end

@testset "gthr and gthrz" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        x = sparsevec(rand(1:m,k), rand(elty,k), m)
        y = rand(elty,m)
        @testset "gthr!" begin
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            d_y = CUSPARSE.gthr!(d_x,d_y,'O')
            h_x = collect(d_x)
            @test h_x ≈ SparseVector(m,x.nzind,y[x.nzind])
        end

        @testset "gthr" begin
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            d_z = CUSPARSE.gthr(d_x,d_y,'O')
            h_z = collect(d_z)
            @test h_z ≈ SparseVector(m,x.nzind,y[x.nzind])
        end

        @testset "gthrz!" begin
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            d_x,d_y = CUSPARSE.gthrz!(d_x,d_y,'O')
            h_x = collect(d_x)
            h_y = collect(d_y)
            @test h_x ≈ SparseVector(m,x.nzind,y[x.nzind])
            #y[x.nzind] = zero(elty)
            #@test h_y ≈ y
        end

        @testset "gthrz" begin
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            d_z,d_w = CUSPARSE.gthrz(d_x,d_y,'O')
            h_w = collect(d_w)
            h_z = collect(d_z)
            @test h_z ≈ SparseVector(m,x.nzind,y[x.nzind])
            #y[x.nzind] = zero(elty)
            #@test h_w ≈ y
        end
    end
end

@testset "gtsv" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        @testset "gtsv!" begin
            dl = rand(elty,m-1)
            du = rand(elty,m-1)
            d = rand(elty,m)
            B = rand(elty,m,n)
            d_dl = CuArray(vcat([0],dl))
            d_du = CuArray(vcat(du,[0]))
            d_d = CuArray(d)
            d_B = CuArray(B)
            d_B = CUSPARSE.gtsv!(d_dl,d_d,d_du,d_B)
            C = diagm(0=>d) + diagm(1=>du) + diagm(-1=>dl)
            h_B = collect(d_B)
            @test h_B ≈ C\B
        end

        @testset "gtsv" begin
            dl = rand(elty,m-1)
            du = rand(elty,m-1)
            d = rand(elty,m)
            B = rand(elty,m,n)
            d_dl = CuArray(vcat([0],dl))
            d_du = CuArray(vcat(du,[0]))
            d_d = CuArray(d)
            d_B = CuArray(B)
            d_C = CUSPARSE.gtsv(d_dl,d_d,d_du,d_B)
            C = diagm(0=>d) + diagm(1=>du) + diagm(-1=>dl)
            h_C = collect(d_C)
            @test h_C ≈ C\B
        end

        @testset "gtsv_nopivot!" begin
            dl = rand(elty,m-1)
            du = rand(elty,m-1)
            d = rand(elty,m)
            B = rand(elty,m,n)
            d_dl = CuArray(vcat([0],dl))
            d_du = CuArray(vcat(du,[0]))
            d_d = CuArray(d)
            d_B = CuArray(B)
            d_B = CUSPARSE.gtsv_nopivot!(d_dl,d_d,d_du,d_B)
            C = diagm(0=>d) + diagm(1=>du) + diagm(-1=>dl)
            h_B = collect(d_B)
            @test h_B ≈ C\B
        end

        @testset "gtsv_nopivot" begin
            dl = rand(elty,m-1)
            du = rand(elty,m-1)
            d = rand(elty,m)
            B = rand(elty,m,n)
            d_dl = CuArray(vcat([0],dl))
            d_du = CuArray(vcat(du,[0]))
            d_d = CuArray(d)
            d_B = CuArray(B)
            d_C = CUSPARSE.gtsv_nopivot(d_dl,d_d,d_du,d_B)
            C = diagm(0=>d) + diagm(1=>du) + diagm(-1=>dl)
            h_C = collect(d_C)
            @test h_C ≈ C\B
        end

        @testset "gtsvStridedBatch!" begin
            dla = rand(elty,m-1)
            dua = rand(elty,m-1)
            da = rand(elty,m)
            dlb = rand(elty,m-1)
            dub = rand(elty,m-1)
            db = rand(elty,m)
            xa = rand(elty,m)
            xb = rand(elty,m)
            d_dl = CuArray(vcat([0],dla,[0],dlb))
            d_du = CuArray(vcat(dua,[0],dub,[0]))
            d_d = CuArray(vcat(da,db))
            d_x = CuArray(vcat(xa,xb))
            d_x = CUSPARSE.gtsvStridedBatch!(d_dl,d_d,d_du,d_x,2,m)
            Ca = diagm(0=>da) + diagm(1=>dua) + diagm(-1=>dla)
            Cb = diagm(0=>db) + diagm(1=>dub) + diagm(-1=>dlb)
            h_x = collect(d_x)
            @test h_x[1:m] ≈ Ca\xa
            @test h_x[m+1:2*m] ≈ Cb\xb
        end

        @testset "gtsvStridedBatch" begin
            dla = rand(elty,m-1)
            dua = rand(elty,m-1)
            da = rand(elty,m)
            dlb = rand(elty,m-1)
            dub = rand(elty,m-1)
            db = rand(elty,m)
            xa = rand(elty,m)
            xb = rand(elty,m)
            d_dl = CuArray(vcat([0],dla,[0],dlb))
            d_du = CuArray(vcat(dua,[0],dub,[0]))
            d_d = CuArray(vcat(da,db))
            d_x = CuArray(vcat(xa,xb))
            d_y = CUSPARSE.gtsvStridedBatch(d_dl,d_d,d_du,d_x,2,m)
            Ca = diagm(0=>da) + diagm(1=>dua) + diagm(-1=>dla)
            Cb = diagm(0=>db) + diagm(1=>dub) + diagm(-1=>dlb)
            h_y = collect(d_y)
            @test h_y[1:m] ≈ Ca\xa
            @test h_y[m+1:2*m] ≈ Cb\xb
        end
    end
end

@testset "hybsv" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = rand(elty,m,m)
        A = triu(A)
        x = rand(elty,m)
        alpha = rand(elty)
        beta = rand(elty)
        d_x = CuArray(x)
        d_A = CuSparseMatrixCSR(sparse(A))
        d_A = CUSPARSE.switch2hyb(d_A)
        d_y = CUSPARSE.sv('N','T','U',alpha,d_A,d_x,'O')
        h_y = collect(d_y)
        y = A\(alpha * x)
        @test y ≈ h_y

        #=d_y = UpperTriangular(d_A) \ d_x
        h_y = collect(d_y)
        @test h_y ≈ A\x=#

        d_x = CUDA.rand(elty,n)
        info = CUSPARSE.sv_analysis('N','T','U',d_A,'O')
        @test_throws DimensionMismatch CUSPARSE.sv_solve('N','U',alpha,d_A,d_x,info,'O')
        A = sparse(rand(elty,m,n))
        d_A = CuSparseMatrixCSR(A)
        d_A = CUSPARSE.switch2hyb(d_A)
        @test_throws DimensionMismatch CUSPARSE.sv_analysis('T','T','U',d_A,'O')
        CUSPARSE.cusparseDestroySolveAnalysisInfo(info)
    end
end

@testset "mm" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        Ah = sparse(rand(elty,k,k))
        Ah += Ah'
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = CuSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = collect(d_D)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,beta,d_C,'O')
            h_D = collect(d_D)
            D = A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,d_C,'O')
            h_D = collect(d_D)
            D = A * B + C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,'O')
            h_D = collect(d_D)
            D = alpha * A * B
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,'O')
            h_D = collect(d_D)
            D = A * B
            @test D ≈ h_D
            d_D = d_A*d_B
            h_D = collect(d_D)
            D = A * B
            @test D ≈ h_D
            d_A = Hermitian(CuSparseMatrixCSR(Ah))
            d_D = CUSPARSE.mm('N',d_A,d_B,'O')
            h_D = collect(d_D)
            D = Ah * B
            @test D ≈ h_D
        end
        @testset "csc" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = CuSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = collect(d_D)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,beta,d_C,'O')
            h_D = collect(d_D)
            D = A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,d_C,'O')
            h_D = collect(d_D)
            D = A * B + C
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,'O')
            h_D = collect(d_D)
            D = alpha * A * B
            @test D ≈ h_D
            d_D = CUSPARSE.mm('N',d_A,d_B,'O')
            h_D = collect(d_D)
            D = A * B
            @test D ≈ h_D
            d_D = d_A*d_B
            h_D = collect(d_D)
            D = A * B
            @test D ≈ h_D
        end
    end
end

@testset "mm_symm" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,m))
        A = A + transpose(A)
        B = rand(elty,m,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = Symmetric(CuSparseMatrixCSR(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_C = transpose(d_A) * d_B
            h_C = collect(d_C)
            D = transpose(A) * B
            @test D ≈ h_C
            d_B = CUDA.rand(elty,k,n)
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
        end
        @testset "csc" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = Symmetric(CuSparseMatrixCSC(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_C = transpose(d_A) * d_B
            h_C = collect(d_C)
            D = transpose(A) * B
            @test D ≈ h_C
            d_B = CUDA.rand(elty,k,n)
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
        end
    end
end
@testset "mm_herm" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,m))
        A = A + A'
        B = rand(elty,m,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = Hermitian(CuSparseMatrixCSR(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_B = CUDA.rand(elty,k,n)
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
        end
        @testset "csc" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = Hermitian(CuSparseMatrixCSC(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_B = CUDA.rand(elty,k,n)
            @test_throws DimensionMismatch CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O')
        end
    end
end

@testset "mm2" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = CuSparseMatrixCSR(A)
            @test_throws DimensionMismatch CUSPARSE.mm2!('N','T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('T','N',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('T','T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_B,'O')
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = collect(d_D)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,beta,d_C,'O')
            h_D = collect(d_D)
            D = A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,d_C,'O')
            h_D = collect(d_D)
            D = A * B + C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,'O')
            h_D = collect(d_D)
            D = alpha * A * B
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,'O')
            h_D = collect(d_D)
            D = A * B
            @test D ≈ h_D
        end
        @testset "csc" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = CuSparseMatrixCSC(A)
            @test_throws DimensionMismatch CUSPARSE.mm2!('N','T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('T','N',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('T','T',alpha,d_A,d_B,beta,d_C,'O')
            @test_throws DimensionMismatch CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_B,'O')
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = collect(d_D)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,beta,d_C,'O')
            h_D = collect(d_D)
            D = A * B + beta * C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,d_C,'O')
            h_D = collect(d_D)
            D = A * B + C
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,'O')
            h_D = collect(d_D)
            D = alpha * A * B
            @test D ≈ h_D
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,'O')
            h_D = collect(d_D)
            D = A * B
            @test D ≈ h_D
        end
    end
end
@testset "bsrmm2" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        d_B = CuArray(B)
        d_C = CuArray(C)
        d_A = CuSparseMatrixCSR(A)
        d_A = CUSPARSE.switch2bsr(d_A,convert(Cint,blockdim))
        @test_throws DimensionMismatch CUSPARSE.mm2('N','T',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2('T','N',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2('T','T',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_B,'O')
        d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_C,'O')
        h_D = collect(d_D)
        D = alpha * A * B + beta * C
        @test D ≈ h_D
        d_D = CUSPARSE.mm2('N','N',d_A,d_B,beta,d_C,'O')
        h_D = collect(d_D)
        D = A * B + beta * C
        @test D ≈ h_D
        d_D = CUSPARSE.mm2('N','N',d_A,d_B,d_C,'O')
        h_D = collect(d_D)
        D = A * B + C
        @test D ≈ h_D
        d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,'O')
        h_D = collect(d_D)
        D = alpha * A * B
        @test D ≈ h_D
        d_D = CUSPARSE.mm2('N','N',d_A,d_B,'O')
        h_D = collect(d_D)
        D = A * B
        @test D ≈ h_D
    end
end

@testset "mv" begin
    for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,m))
        x = rand(elty,m)
        y = rand(elty,m)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "mv_symm" begin
            A_s = A + transpose(A)
            @testset "csr" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = Symmetric(CuSparseMatrixCSR(A_s))
                d_y = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
                h_y = collect(d_y)
                z = alpha * A_s * x + beta * y
                @test z ≈ h_y
                x_  = rand(elty,n)
                d_x = CuArray(x_)
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
            end
            @testset "csc" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = Symmetric(CuSparseMatrixCSC(A_s))
                d_y = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
                h_y = collect(d_y)
                z = alpha * A_s * x + beta * y
                @test z ≈ h_y
                x_  = rand(elty,n)
                d_x = CuArray(x_)
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
            end
        end

        @testset "mv_herm" begin
            A_h = A + adjoint(A)
            @testset "csr" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = Hermitian(CuSparseMatrixCSR(A_h))
                d_y = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
                h_y = collect(d_y)
                z = alpha * A_h * x + beta * y
                @test z ≈ h_y
                x_  = rand(elty,n)
                d_x = CuArray(x_)
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
            end
            @testset "csc" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = Hermitian(CuSparseMatrixCSC(A_h))
                d_y = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
                h_y = collect(d_y)
                z = alpha * A_h * x + beta * y
                @test z ≈ h_y
                x_  = rand(elty,n)
                d_x = CuArray(x_)
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
            end
        end
        A = sparse(rand(elty,m,n))
        x = rand(elty,n)
        y = rand(elty,m)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "mv" begin
            @testset "csr" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = CuSparseMatrixCSR(A)
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
                h_z = collect(d_z)
                z = alpha * A * x + beta * y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,beta,d_y,'O')
                h_z = collect(d_z)
                z = A * x + beta * y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,d_y,'O')
                h_z = collect(d_z)
                z = A * x + y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,d_y,'O')
                h_z = collect(d_z)
                z = alpha * A * x + y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,'O')
                h_z = collect(d_z)
                z = alpha * A * x
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,'O')
                h_z = collect(d_z)
                z = A * x
                @test z ≈ h_z
                d_z = d_A*d_x
                h_z = collect(d_z)
                z = A * x
                @test z ≈ h_z
            end
            @testset "csc" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = CuSparseMatrixCSC(A)
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
                h_z = collect(d_z)
                z = alpha * A * x + beta * y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,beta,d_y,'O')
                h_z = collect(d_z)
                z = A * x + beta * y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,d_y,'O')
                h_z = collect(d_z)
                z = A * x + y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,d_y,'O')
                h_z = collect(d_z)
                z = alpha * A * x + y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,'O')
                h_z = collect(d_z)
                z = alpha * A * x
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,'O')
                h_z = collect(d_z)
                z = A * x
                @test z ≈ h_z
                d_z = d_A*d_x
                h_z = collect(d_z)
                z = A * x
                @test z ≈ h_z
            end
            @testset "bsr" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = CuSparseMatrixCSR(A)
                d_A = CUSPARSE.switch2bsr(d_A,convert(Cint,blockdim))
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
                h_z = collect(d_z)
                z = alpha * A * x + beta * y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,beta,d_y,'O')
                h_z = collect(d_z)
                z = A * x + beta * y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,d_y,'O')
                h_z = collect(d_z)
                z = A * x + y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,d_y,'O')
                h_z = collect(d_z)
                z = alpha * A * x + y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,'O')
                h_z = collect(d_z)
                z = alpha * A * x
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,'O')
                h_z = collect(d_z)
                z = A * x
                @test z ≈ h_z
                d_z = d_A*d_x
                h_z = collect(d_z)
                z = A * x
                @test z ≈ h_z
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
            end
            @testset "hyb" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = CuSparseMatrixCSR(A)
                d_A = CUSPARSE.switch2hyb(d_A)
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,beta,d_y,'O')
                h_z = collect(d_z)
                z = alpha * A * x + beta * y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,beta,d_y,'O')
                h_z = collect(d_z)
                z = A * x + beta * y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,d_y,'O')
                h_z = collect(d_z)
                z = A * x + y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,d_y,'O')
                h_z = collect(d_z)
                z = alpha * A * x + y
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',alpha,d_A,d_x,'O')
                h_z = collect(d_z)
                z = alpha * A * x
                @test z ≈ h_z
                d_z = CUSPARSE.mv('N',d_A,d_x,'O')
                h_z = collect(d_z)
                z = A * x
                @test z ≈ h_z
                d_z = d_A*d_x
                h_z = collect(d_z)
                z = A * x
                @test z ≈ h_z
                @test_throws DimensionMismatch CUSPARSE.mv('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv('N',alpha,d_A,d_y,beta,d_x,'O')
            end
        end
    end
end

@testset "mm!" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = CuSparseMatrixCSR(A)
            CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_C = CuArray(C)
            mul!(d_C, d_A, d_B)
            h_C = collect(d_C)
            D = A * B
            @test D ≈ h_C
        end
        @testset "csc" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = CuSparseMatrixCSC(A)
            CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_C = CuArray(C)
            mul!(d_C, d_A, d_B)
            h_C = collect(d_C)
            D = A * B
            @test D ≈ h_C
        end
    end
end
@testset "mm_symm!" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,m))
        A = A + transpose(A)
        B = rand(elty,m,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = Symmetric(CuSparseMatrixCSR(A))
            CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_C = CuArray(C)
            mul!(d_C, d_A, d_B)
            h_C = collect(d_C)
            D = A * B
            @test D ≈ h_C
        end
        @testset "csc" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = Symmetric(CuSparseMatrixCSC(A))
            CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_C = CuArray(C)
            mul!(d_C, d_A, d_B)
            h_C = collect(d_C)
            D = A * B
            @test D ≈ h_C
        end
    end
end
@testset "mm_herm!" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,m))
        A = A + A'
        B = rand(elty,m,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = Hermitian(CuSparseMatrixCSR(A))
            CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_C = CuArray(C)
            mul!(d_C, d_A, d_B)
            h_C = collect(d_C)
            D = A * B
            @test D ≈ h_C
        end
        @testset "csc" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = Hermitian(CuSparseMatrixCSC(A))
            CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_C
            d_C = CuArray(C)
            mul!(d_C, d_A, d_B)
            h_C = collect(d_C)
            D = A * B
            @test D ≈ h_C
        end
    end
end
@testset "mm2!" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        Bᵀ = collect(transpose(B))
        alpha = rand(elty)
        beta = rand(elty)
        @testset "csr" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = CuSparseMatrixCSR(A)
            CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_C = CuArray(C)
            mul!(d_C, d_A, d_B)
            h_C = collect(d_C)
            D = A * B
            @test D ≈ h_C
            d_Bᵀ = CuArray(Bᵀ)
            d_C = CuArray(C)
            mul!(d_C, d_A, transpose(d_Bᵀ))
            h_C = collect(d_C)
            D = A * transpose(Bᵀ)
            @test D ≈ h_C
        end
        @testset "csc" begin
            d_B = CuArray(B)
            d_C = CuArray(C)
            d_A = CuSparseMatrixCSC(A)
            CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = collect(d_C)
            D = alpha * A * B + beta * C
            @test D ≈ h_D
            d_C = CuArray(C)
            mul!(d_C, d_A, d_B)
            h_C = collect(d_C)
            D = A * B
            @test D ≈ h_C
            d_Bᵀ = CuArray(Bᵀ)
            d_C = CuArray(C)
            @test_throws CUSPARSEError mul!(d_C, d_A, transpose(d_Bᵀ))
        end
    end
end
@testset "bsrmm2!" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,k))
        B = rand(elty,k,n)
        C = rand(elty,m,n)
        alpha = rand(elty)
        beta = rand(elty)
        d_B = CuArray(B)
        d_C = CuArray(C)
        d_A = CuSparseMatrixCSR(A)
        d_A = CUSPARSE.switch2bsr(d_A,convert(Cint,blockdim))
        @test_throws DimensionMismatch CUSPARSE.mm2!('N','T',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2!('T','N',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2!('T','T',alpha,d_A,d_B,beta,d_C,'O')
        @test_throws DimensionMismatch CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_B,'O')
        CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_C,'O')
        h_D = collect(d_C)
        D = alpha * A * B + beta * C
        @test D ≈ h_D
        d_C = CuArray(C)
        mul!(d_C, d_A, d_B)
        h_C = collect(d_C)
        D = A * B
        @test D ≈ h_C
    end
end
@testset "mv!" begin
    for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,m))
        x = rand(elty,m)
        y = rand(elty,m)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "mv_symm!" begin
            A_s = A + transpose(A)
            @testset "csr" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = Symmetric(CuSparseMatrixCSR(A_s))
                CUSPARSE.mv!('N',alpha,d_A,d_x,beta,d_y,'O')
                h_y = collect(d_y)
                z = alpha * A_s * x + beta * y
                @test z ≈ h_y
                x_  = rand(elty,n)
                d_x = CuArray(x_)
                @test_throws DimensionMismatch CUSPARSE.mv!('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv!('N',alpha,d_A,d_y,beta,d_x,'O')
                d_x = CuArray(x)
                mul!(d_y, d_A, d_x)
                h_y = collect(d_y)
                z = A_s * x
                @test z ≈ h_y
            end
            @testset "csc" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = Symmetric(CuSparseMatrixCSC(A_s))
                mv!('N',alpha,d_A,d_x,beta,d_y,'O')
                h_y = collect(d_y)
                z = alpha * A_s * x + beta * y
                @test z ≈ h_y
                x_  = rand(elty,n)
                d_x = CuArray(x_)
                @test_throws DimensionMismatch CUSPARSE.mv!('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv!('N',alpha,d_A,d_y,beta,d_x,'O')
                d_x = CuArray(x)
                mul!(d_y, d_A, d_x)
                h_y = collect(d_y)
                z = A_s * x
                @test z ≈ h_y
            end
        end

        @testset "mv_herm!" begin
            A_h = A + adjoint(A)
            @testset "csr" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = Hermitian(CuSparseMatrixCSR(A_h))
                CUSPARSE.mv!('N',alpha,d_A,d_x,beta,d_y,'O')
                h_y = collect(d_y)
                z = alpha * A_h * x + beta * y
                @test z ≈ h_y
                x_  = rand(elty,n)
                d_x = CuArray(x_)
                @test_throws DimensionMismatch CUSPARSE.mv!('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv!('N',alpha,d_A,d_y,beta,d_x,'O')
                d_x = CuArray(x)
                mul!(d_y, d_A, d_x)
                h_y = collect(d_y)
                z = A_h * x
                @test z ≈ h_y
            end
            @testset "csc" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = Hermitian(CuSparseMatrixCSC(A_h))
                CUSPARSE.mv!('N',alpha,d_A,d_x,beta,d_y,'O')
                h_y = collect(d_y)
                z = alpha * A_h * x + beta * y
                @test z ≈ h_y
                x_  = rand(elty,n)
                d_x = CuArray(x_)
                @test_throws DimensionMismatch CUSPARSE.mv!('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv!('N',alpha,d_A,d_y,beta,d_x,'O')
                d_x = CuArray(x)
                mul!(d_y, d_A, d_x)
                h_y = collect(d_y)
                z = A_h * x
                @test z ≈ h_y
            end
        end
        A = sparse(rand(elty,m,n))
        x = rand(elty,n)
        y = rand(elty,m)
        alpha = rand(elty)
        beta = rand(elty)
        @testset "mv!" begin
            @testset "csr" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = CuSparseMatrixCSR(A)
                @test_throws DimensionMismatch CUSPARSE.mv!('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv!('N',alpha,d_A,d_y,beta,d_x,'O')
                CUSPARSE.mv!('N',alpha,d_A,d_x,beta,d_y,'O')
                h_z = collect(d_y)
                z = alpha * A * x + beta * y
                @test z ≈ h_z
                mul!(d_y, d_A, d_x)
                h_y = collect(d_y)
                z = A * x
                @test z ≈ h_y
            end
            @testset "csc" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = CuSparseMatrixCSC(A)
                @test_throws DimensionMismatch CUSPARSE.mv!('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv!('N',alpha,d_A,d_y,beta,d_x,'O')
                CUSPARSE.mv!('N',alpha,d_A,d_x,beta,d_y,'O')
                h_z = collect(d_y)
                z = alpha * A * x + beta * y
                @test z ≈ h_z
                mul!(d_y, d_A, d_x)
                h_y = collect(d_y)
                z = A * x
                @test z ≈ h_y
            end
            @testset "bsr" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = CuSparseMatrixCSR(A)
                d_A = CUSPARSE.switch2bsr(d_A,convert(Cint,blockdim))
                @test_throws DimensionMismatch CUSPARSE.mv!('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv!('N',alpha,d_A,d_y,beta,d_x,'O')
                CUSPARSE.mv!('N',alpha,d_A,d_x,beta,d_y,'O')
                h_z = collect(d_y)
                z = alpha * A * x + beta * y
                @test z ≈ h_z
                mul!(d_y, d_A, d_x)
                h_y = collect(d_y)
                z = A * x
                @test z ≈ h_y
            end
            @testset "hyb" begin
                d_x = CuArray(x)
                d_y = CuArray(y)
                d_A = CuSparseMatrixCSR(A)
                d_A = CUSPARSE.switch2hyb(d_A)
                CUSPARSE.mv!('N',alpha,d_A,d_x,beta,d_y,'O')
                h_z = collect(d_y)
                z = alpha * A * x + beta * y
                @test z ≈ h_z
                @test_throws DimensionMismatch CUSPARSE.mv!('T',alpha,d_A,d_x,beta,d_y,'O')
                @test_throws DimensionMismatch CUSPARSE.mv!('N',alpha,d_A,d_y,beta,d_x,'O')
                mul!(d_y, d_A, d_x)
                h_y = collect(d_y)
                z = A * x
                @test z ≈ h_y
            end
        end
    end
end

@testset "sctr" begin
    @testset for elty in [Float32,Float64,ComplexF32,ComplexF64]
        x = sparsevec(rand(1:m,k), rand(elty,k), m)
        y = zeros(elty,m)
        @testset "sctr!" begin
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            d_y = CUSPARSE.sctr!(d_x,d_y,'O')
            h_y = collect(d_y)
            y[x.nzind]  += x.nzval
            @test h_y ≈ y
        end
        y = zeros(elty,m)

        @testset "sctr" begin
            d_x = CuSparseVector(x)
            d_y = CUSPARSE.sctr(d_x,'O')
            h_y = collect(d_y)
            y = zeros(elty,m)
            y[x.nzind]  += x.nzval
            @test h_y ≈ y
        end
    end
end

@testset "roti" begin
    @testset for elty in [Float32,Float64]
        x = sparsevec(rand(1:m,k), rand(elty,k), m)
        y = rand(elty,m)
        @testset "roti!" begin
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            angle = rand(elty)
            d_x,d_y = CUSPARSE.roti!(d_x,d_y,cos(angle),sin(angle),'O')
            h_x = collect(d_x)
            h_y = collect(d_y)
            z = copy(x)
            w = copy(y)
            y[x.nzind] = cos(angle)*w[z.nzind] - sin(angle)*z.nzval
            @test h_x ≈ SparseVector(m,x.nzind,cos(angle)*z.nzval + sin(angle)*w[z.nzind])
            @test h_y ≈ y
        end

        @testset "roti" begin
            d_x = CuSparseVector(x)
            d_y = CuArray(y)
            angle = rand(elty)
            d_z,d_w = CUSPARSE.roti(d_x,d_y,cos(angle),sin(angle),'O')
            h_w = collect(d_w)
            h_z = collect(d_z)
            z = copy(x)
            w = copy(y)
            w[z.nzind] = cos(angle)*y[x.nzind] - sin(angle)*x.nzval
            @test h_z ≈ SparseVector(m,z.nzind, cos(angle)*x.nzval + sin(angle)*y[x.nzind])
            @test h_w ≈ w
        end
    end
end

@testset "mul!" begin
    for elty in [Float32,Float64,ComplexF32,ComplexF64]
        A = sparse(rand(elty,m,m))
        x = rand(elty,m)
        y = rand(elty,m)
        @testset "csr -- $elty" begin
            d_x  = CuArray(x)
            d_y  = CuArray(y)
            d_A  = CuSparseMatrixCSR(A)
            d_Aᵀ = transpose(d_A)
            d_Aᴴ = adjoint(d_A)
            CUSPARSE.mul!(d_y, d_A, d_x)
            h_y = collect(d_y)
            z = A * x
            @test z ≈ h_y
            CUSPARSE.mul!(d_y, d_Aᵀ, d_x)
            h_y = collect(d_y)
            z = transpose(A) * x
            @test z ≈ h_y
            CUSPARSE.mul!(d_y, d_Aᴴ, d_x)
            h_y = collect(d_y)
            z = adjoint(A) * x
            @test z ≈ h_y
        end
        @testset "csc -- $elty" begin
            d_x = CuArray(x)
            d_y = CuArray(y)
            d_A = CuSparseMatrixCSC(A)
            d_Aᵀ = transpose(d_A)
            d_Aᴴ = adjoint(d_A)
            CUSPARSE.mul!(d_y, d_A, d_x)
            h_y = collect(d_y)
            z = A * x
            @test z ≈ h_y
            CUSPARSE.mul!(d_y, d_Aᵀ, d_x)
            h_y = collect(d_y)
            z = transpose(A) * x
            @test z ≈ h_y
            CUSPARSE.mul!(d_y, d_Aᴴ, d_x)
            h_y = collect(d_y)
            z = adjoint(A) * x
            @test z ≈ h_y
        end
    end
end


text:
using ConeProgramDiff

using MathOptInterface
using SCS
using Random
using LinearAlgebra

const MOI = MathOptInterface

function check_kkt(A, b, c, cone_prod, x, y, s; tol=1e-8)
    try
        @assert isapprox(A*x + s, b,atol=tol)
        @assert isapprox(A'*y + c, zeros(length(c)),atol=tol)
        @assert isapprox(s, project_onto_cone(s, cone_prod),atol=tol)
        @assert isapprox(y, project_onto_cone(y, [MOI.dual_set(c) for c in cone_prod]),atol=tol)
        @assert isapprox(s'*y, 0, atol=tol)
    catch e
        error("KKT conditions failed:
               \nprimal constraint: $(norm(A*x + s - b))
               \ndual constraint: $(norm(A'*y +c))
               \nK error: $(norm(s - project_onto_cone(s, cone_prod)))
               \nK* error: $(norm(y - project_onto_cone(y, [MOI.dual_set(c) for c in cone_prod])))
               \nCS cond: $(norm(s'*y))"
        )
    end
end

function check_adjoint(tol)
    Random.seed!(0)
    dims = (15, 10)
    A, b, c, cone_prod = l1_minimization_program(dims; solve=false)
    x, y, s, pf, pb = solve_and_diff(A, b, c, cone_prod; eps=1e-10, max_iters=10000)
    check_kkt(A, b, c, cone_prod, x, y, s, tol=1e-8)

    # Check adjoint
    pstar = dot(x, c)
    dA, db, dc = pb(c)
    del = 1e-6
    x_pert, _, _, _, _ = solve_and_diff(A+del*dA, b+del*db, c+del*dc, cone_prod;
                                                  eps=1e-10, max_iters=10000)
    pstar_pert = dot(x_pert, c)
    df = pstar_pert - pstar
    dp = del * (sum(dA .* dA) + db'*db + dc'*dc)
    if !isapprox(df, dp, atol=tol)
        error("Adjoint mismatch: $(abs(df-dp))\n\tdf = $df\n\tdp = $dp")
    end
    return true
end

function check_derivative()
    Random.seed!(0)
    dims = (15, 10)
    A, b, c, cone_prod = l1_minimization_program(dims; solve=false)
    x, y, s, pf, pb = solve_and_diff(A, b, c, cone_prod; eps=1e-10, max_iters=10000)
    check_kkt(A, b, c, cone_prod, x, y, s, tol=1e-6)

    # Check derivative
    del = 1e-6
    dA, db, dc = del*randn(size(A)), del*randn(size(b)), del*randn(size(c))
    dx, dy, ds = pf(dA, db, dc)
    x_pert, y_pert, s_pert, _, _ = solve_and_diff(A+dA, b+db, c+dc, cone_prod;
                                                  eps=1e-10, max_iters=10000)

    # TODO: tols are high here -- should investigate
    @assert isapprox(x_pert - x, dx, atol=1e-4)
    @assert isapprox(y_pert - y, dy, atol=1e-4)
    @assert isapprox(s_pert - s, ds, atol=2e-4)
    return true
end

@testset "solve_and_diff Tests" begin
    @test check_adjoint(2e-5)
    @test check_derivative()
end

ext: jl


text:
# ArgMacros
# help.jl

# Extend only for the Argument and Help types used here
# Extended methods will not be exported
import Base.show

#=
Macros have no function themselves
Only used as markers in the AST for generating Help
=#

"""
    @helpusage usage_text::String

Add usage text for the help screen
Automatically prepended with "Usage: "

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    @helpusage "example.jl foo [bar] [-v]"
    ...
end
```
"""
macro helpusage(usage::String) end

"""
    @helpdescription description::String

Add description for the help screen

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @helpdescription "Lorem ipsum dolor sit amet"
    ...
end
```
"""
macro helpdescription(description::String) end

"""
    @helpepilog epilog::String

Add epilog for the help screen, displayed after rest of help

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @helpepilog "Lorem ipsum dolor sit amet"
    ...
end
```
"""
macro helpepilog(epilog::String) end

"""
    @arghelp help_text::String

Add help string for an argument.
Applied to the preceding declared argument

Must be used in `@beginarguments begin ... end` block

# Example
```julia
@beginarguments begin
    ...
    @argumentflag v "-v" "--verbose"
    @arghelp "Display additional output"
    ...
end
```
"""
macro arghelp(helptext::String) end

"Print a message, ask user to try --help, and exit the program"
function _quit_try_help(message::String)
    println(message)
    println("Try the --help option")
    exit(1)
end

#=
Help is generated on-demand, may want to change this for precompiled scripts
=#

"Check for -h/--help flag, end program and print help if needed"
function _help_check(args::Vector{String}, block::Expr)
    if !isnothing(_get_option_idx(args, ["-h", "--help"]))
        _make_help(block) |> print
        exit(0)
    end
end

"Hold information to print an argument"
@kwdef mutable struct Argument
    names::Vector{Union{String, Symbol}} = []
    type::Symbol = :Any
    required::Bool = false
    positional::Bool = false
    default::String = ""
    description::String = ""
end

"Hold information to print entire help screen"
@kwdef mutable struct Help
    usage_text::String = ""
    description_text::String = ""
    epilog_text::String = ""
    positionals::Vector{Argument} = []
    options::Vector{Argument} = []
end

"Assemble the Help struct from the macro contents"
function _make_help(block::Expr)::Help
    help = Help()

    lastpushed::Symbol = :none

    # Loop through all macros in the block
    for arg in _getmacrocalls(block)
        macroname::Symbol = _get_macroname(arg)

        # Modify the help struct as needed by the encountered macro
        if macroname == USAGE_SYMBOL
            help.usage_text = arg.args[3] # 3rd arg for these is the string passed
        elseif macroname == DESCRIPTION_SYMBOL
            help.description_text = arg.args[3]
        elseif macroname == EPILOG_SYMBOL
            help.epilog_text = arg.args[3]
        elseif macroname == ARGHELP_SYMBOL
            # Add description to last pushed element in positionals or options
            argvector::Vector{Argument} = getfield(help, lastpushed)
            argvector[end].description = arg.args[3]
        elseif macroname == ARGUMENT_REQUIRED_SYMBOL
            push!(help.options, Argument(
                arg.args[5:end], arg.args[3], true, false, "", "")
            )
            lastpushed = :options
        elseif macroname == ARGUMENT_DEFAULT_SYMBOL
            push!(help.options, Argument(
                arg.args[6:end], arg.args[3], false, false, string(arg.args[4]), "")
            )
            lastpushed = :options
        elseif macroname == ARGUMENT_OPTIONAL_SYMBOL
            push!(help.options, Argument(
                arg.args[5:end], arg.args[3], false, false, "", "")
            )
            lastpushed = :options
        elseif macroname == ARGUMENT_FLAG_SYMBOL
            push!(help.options, Argument(
                arg.args[4:end], :Flag, false, false, "", "")
            )
            lastpushed = :options
        elseif macroname == ARGUMENT_COUNT_SYMBOL
            push!(help.options, Argument(
                arg.args[4:end], :Count, false, false, "", "")
            )
            lastpushed = :options
        elseif macroname == POSITIONAL_REQUIRED_SYMBOL
            push!(help.positionals, Argument(
                [arg.args[end]], arg.args[3], true, true, "", "")
            )
            lastpushed = :positionals
        elseif macroname == POSITIONAL_DEFAULT_SYMBOL
            push!(help.positionals, Argument(
                [arg.args[end]], arg.args[3], false, true, string(arg.args[4]), "")
            )
            lastpushed = :positionals
        elseif macroname == POSITIONAL_OPTIONAL_SYMBOL
            push!(help.positionals, Argument(
                [arg.args[end]], arg.args[3], false, true, "", "")
            )
            lastpushed = :positionals
        end
    end

    return help
end

function show(io::IO, arg::Argument)
    # Print as "required" or "[optional]"
    # For positional arguments
    if arg.positional && !arg.required
        print(io, "[")
        print(io, join(arg.names, " "))
        print(io, "]")
    else
        print(io, join(arg.names, " "))
    end
    println()

    if !isempty(arg.description)
        println_wrapped(io, arg.description, initial_indent=6, subsequent_indent=6)
    end

    println(io, " "^6 * (
        isempty(arg.default) ?
            arg.required ? "(Type: $(arg.type), Required)" : "(Type: $(arg.type))" :
            "(Type: $(arg.type), Default: $(arg.default))"
    ))
end

function show(io::IO, help::Help)
    if !isempty(help.usage_text)
        println(io, "Usage: $(help.usage_text)")
        println(io)
    end

    if !isempty(help.description_text)
        println(io, help.description_text)
        println(io)
    end

    if !isempty(help.positionals)
        println(io, "Positional Arguments:")
        for arg in help.positionals
            print(io, arg)
        end
        println(io)
    end

    if !isempty(help.options)
        println(io, "Option Arguments:")
        for arg in help.options
            print(io, arg)
        end
        println(io)
    end

    if !isempty(help.epilog_text)
        println(io, help.epilog_text)
        println(io)
    end
end


text:
import Base.Random: rand

export moransel, moransel!

@doc """moransel(p::Population, iters::Int64)

Conduct Moran selection `iters` number of times. See below for more information
about a Moran process. Basically, this means to first select an individual
within the population to "die", then proportionally selection an individual
to replace it.

https://en.wikipedia.org/wiki/Moran_process
"""
function moransel(p::Population, iters::Int64)
  np = Population(p)
  moransel!(np, iters)
  return np
end

function moransel(mp::MetaPopulation, iters::Int64)
  np = MetaPopulation(mp)
  moransel!(np, iters)
  return np
end

@doc """moransel!(p::Population, iters::Int64)

Conduct Moran selection in-place.
"""
function moransel!(p::Population, iters::Int64)
  fs = popfits(p)
  fmax = maximum(fs)
  n = popsize(p)

  for _ = 1:iters
    indkill = rand(1:n)
    indbirth = 0
    while indbirth == 0
      i = rand(1:n)
      w = fs[i] / fmax
      if rand() < w
        indbirth = i
      end
    end

    p.genotypes[indkill] = Genotype(p.genotypes[indbirth])
  end
end

function moransel!(mp::MetaPopulation, iters::Int64)
  for p = mp.populations
    moransel!(p, iters)
  end
end


text:
baremodule UnixIOHeaders

import Base
using Base:@ccall, Cstring, Cint, Sys, Vector, signed, (&), (+), (>>), (>), (==)

const __spawn_action = Cvoid

ioctl(fd, cmd, arg) = @ccall ioctl(fd::Cint, cmd::Cint, arg::Ptr{Cint})::Cint

using CInclude

@cinclude([
    "<errno.h>",
    "<string.h>",
    "<limits.h>",
    "<stdlib.h>",
    "<pthread.h>",
    "<sys/ioctl.h>",
    "<termios.h>",
    "<fcntl.h>",
    "<poll.h>",
    (Sys.islinux() ? ("<sys/epoll.h>",) : ())...,
    "<unistd.h>",
    "<sys/stat.h>",
    "<sys/socket.h>",
    "<signal.h>",
    "<sys/wait.h>",
    "<sys/syscall.h>",
    "<spawn.h>"],

    args = ["-D_GNU_SOURCE"],

    quiet,

    exclude=r"""
        sv_onstack
      | sched_priority | sigcontext_struct
      | ru_first | ru_last
      | NR_OPEN | ARG_MAX | LINK_MAX
      | MACH_MSG_TYPE_INTEGER_T | msgh_reserved | msgh_kind | mach_msg_kind_t
      | ^SIOC | TC[SG]ETS[FW]?2
    """x
    )

# Need Ptr{Ptr{UInt8}} not Ptr{String} for NULL-terminated string vectors:
execv(path, args::Vector{Ptr{UInt8}}) =
    @ccall execv(path::Ptr{UInt8}, args::Ptr{Ptr{UInt8}})::Cint

posix_spawn(pid, path, file_actions, attrp, argv::Vector{Ptr{UInt8}},
                                            envp::Vector{Ptr{UInt8}}) =
    @ccall posix_spawn(pid::Ptr{pid_t},
                       path::Cstring,
                       file_actions::Ptr{posix_spawn_file_actions_t},
                       attrp::Ptr{posix_spawnattr_t},
                       argv::Ptr{Ptr{UInt8}},
                       envp::Ptr{Ptr{UInt8}})::Cint

# Signal wrongly wrapped as 3-args by CInclude.jl.
if Sys.isapple()
    signal(sig, func) = @ccall signal(sig::Cint, func::sig_t)::sig_t
end
const SIG_DFL = sig_t(0)

# Need multiple methods for `open` and `fcntl`:
open(pathname::AbstractString, flags) =
    @ccall open(pathname::Cstring, flags::Cint)::Cint

open(pathname::AbstractString, flags, mode) =
    @ccall open(pathname::Cstring, flags::Cint, mode::Cint)::Cint

fcntl(fd, cmd) = @ccall fcntl(fd::Cint, cmd::Cint)::Cint
fcntl(fd, cmd, arg) = @ccall fcntl(fd::Cint, cmd::Cint, arg::Cint)::Cint

# Need variants for _m struct
tcsetattr_m(fd, action, p) =
    @ccall tcsetattr(fd::Cint, action::Cint, p::Ptr{termios_m})::Cint

tcgetattr_m(fd, p) =
    @ccall tcgetattr(fd::Cint, p::Ptr{termios_m})::Cint

cfsetspeed_m(p, speed) =
    @ccall cfsetspeed(p::Ptr{termios_m}, speed::speed_t)::Cint

# Not yet in glibc.
const SYS_pidfd_open=434 # https://git.io/J4j1A
pifd_open(pid, flags) = @ccall syscall(SYS_pidfd_open::Cint,
                                       pid::pid_t, flags::Cint)::Cint

const P_PIDFD = 3
waitid(idtype, id, infop, options) =
    @ccall syscall(SYS_waitid::Cint,
                   idtype::Cint,
                   id::id_t,
                   infop::Ptr{siginfo_t},
                   options::Cint,
                   Base.C_NULL::Ptr{Cvoid})::Cint

# Function-like mactos not yet wrapped.
WIFSIGNALED(x) = (((x & 0x7f) + 1) >> 1) > 0
WTERMSIG(x) = (x & 0x7f)
WIFEXITED(x) = WTERMSIG(x) == 0
WEXITSTATUS(x) = signed(UInt8((x >> 8) & 0xff))
WIFSTOPPED(x) = (x & 0xff) == 0x7f
WSTOPSIG(x) = WEXITSTATUS(x)
WIFCONTINUED(x) = x == 0xffff
WCOREDUMP(x) = x & 0x80

end # module


text:
function eval_testfunction(out, f)
    theta = deepcopy(out.theta0)
    x = f(theta)
    d = length(x); n = size(out.Theta)[2]
    F = zeros(eltype(x), d, n)
    for i = 1:n
        theta .= out.Theta[:, i]
        F[:,i] .= f(theta)
    end
    F
end

function iact_sokal(X::T) where {FT <: AbstractFloat, T <: AbstractVector{FT}}
    n = length(X)
    # Calculate centred & normalised X
    X_ = (X.-mean(X))/sqrt(var(X))
    # ...and this gives then ACF at specified lag
    C = FT(max(5.0, log10(n))) # Suggested by Sokal according to
                               # http://dfm.io/posts/autocorr/
    tau = one(FT)
    for k in 1:(n-1)
        tau += 2.0*dot(X_, 1:(n-k), X_, (1+k):n)/(n-k)
        if k>C*tau
            break
        end
    end
    FT(max(0.0, tau))
end

# Default IACT is pessimistic (for safety) -- maximum of three estimators:
# 1) Sokal's estimator defined above
# 2) Batch means from MonteCarloMarkovKernels
# 3) Consistent spectral variance also from MonteCarloMarkovKernels
@inline function default_iact(X)
    vX = var(X)
    BM_iact = MonteCarloMarkovKernels.estimateBM(X)/vX
    SV_iact = MonteCarloMarkovKernels.estimateSV(X, :ModifiedBartlett)/vX
    sokal_iact = iact_sokal(X)
    max(SV_iact, BM_iact, sokal_iact)
end

# When correction kernel is simple, we may calculate estimates for all epsilon
# (and then find the suitable epsilon among)
function is_est(X::Matrix{FT}, D::Vector{FT}, phi::ABCCutoff,
    kernel::ABCSimpleCutoff, eps; normal_quantile=FT(1.96),
    iact=nothing) where {FT <: AbstractFloat}

    eps0 = phi.tol

    d, n = size(X)
    I = sortperm(D); D_ = D[I]
    if typeof(eps) == Nothing
        n_eps = n
        eps = D_
    else
        n_eps = length(eps)
    end

    E = Matrix{FT}(undef,d,n)
    S = Matrix{FT}(undef,d,n)
    E_eps = Matrix{FT}(undef,d,n_eps)
    S_eps = Matrix{FT}(undef,d,n_eps)
    ci_L = Matrix{FT}(undef,d,n_eps)
    ci_U = Matrix{FT}(undef,d,n_eps)

    I = sortperm(D); D_ = D[I]
    phi_0 = phi.(D_)
    invalid = phi_0 .== zero(FT)
    xi_1 = one(FT)./phi_0
    xi_1[invalid] .= zero(FT)
    xi_1 /= sum(xi_1)
    for i in 1:d
        X_ = X[i,I]
        xi_f = X_.*xi_1
        E[i,:] = cumsum(xi_f)./cumsum(xi_1)
        S[i,:] = (cumsum(xi_f.^2) - 2E[i,:].*cumsum(xi_f.*xi_1)
                  + cumsum(xi_1.^2).*E[i,:].^2) ./ cumsum(xi_1).^2
    end
    if typeof(eps) == Nothing
        E_eps = E; S_eps = S
    else
        j = 1
        for k in 1:n_eps
            while j<n && D_[j+1] < eps[k]
                j += 1
            end
            E_eps[:,k] = E[:,j]
            S_eps[:,k] = S[:,j]
        end
    end
    for i in 1:d
        dx = normal_quantile*sqrt.(max.(S_eps[i,:],0)*iact(X[i,:]))
        ci_L[i,:] = E_eps[i,:] - dx; ci_U[i,:] = E_eps[i,:] + dx
    end
    (E=E, S=S, eps=eps, ci_L=ci_L, ci_U=ci_U, normal_quantile=normal_quantile)
end

# For general correction, calculate directly from the definition
function is_est(X::Matrix{FT}, D::Vector{FT}, phi::ABCCutoff,
    kernel::ABCGenCutoff, eps; normal_quantile=FT(1.96),
    iact=nothing) where {FT <: AbstractFloat}

    eps isa Nothing || error("Must supply ϵ with general cutoff")
    d, n = size(X)
    n_eps = length(eps)
    E = Matrix{FT}(undef, d, n_eps)
    S = Matrix{FT}(undef, d, n_eps)
    ci_L = Matrix{FT}(undef, d, n_eps)
    ci_U = Matrix{FT}(undef, d, n_eps)
    # Evaluate the cutoff for eps0, and the positive indices
    phi_0 = phi.(D)
    invalid = phi_0 .== zero(FT)
    for k in 1:n_eps
        U = kernel.fun.(D/eps[k]) ./ phi_0
        U[invalid] .= zero(FT)
        W = U/sum(U)
        for i in 1:d
            E[i,k] = sum(W .* X[i,:])
            S[i,k] = sum(W.^2 .* (X[i,:] .- E[i,k]).^2)
        end
    end
    for i in 1:d
        dx = normal_quantile*sqrt.(max.(S[i,:],0)*iact(X[i,:]))
        ci_L[i,:] = E[i,:] - dx
        ci_U[i,:] = E[i,:] + dx
    end
    (E=E, S=S, eps=eps, ci_L=ci_L, ci_U=ci_U, normal_quantile=normal_quantile)
end

"""
    est = abc_postprocess(out, f, [tol]; regress=false, kernel=out.cutoff,
                          normal_quantile=1.96)

Calculate post-correction estimates for ABC-MCMC output.

# Arguments
- `out`: Output of the `abc_mcmc` function.
- `f`: Function of interest.
- `tol`: Vector of tolerance values for which estimates are calculated.
         If not given (or `nothing`), and `out.cutoff` is simple, then
         correction calculated to all tolerances.
- `regress`: Whether regression correction should be used.
- `kernel`: Kernel for post-correction. (NB kernel.tol is ignored!)
- `normal_quantile`: The desired normal quantile (default 1.96 -> 95% coverage;
                     any other quantile may be calculated using Distributions as:
                     `-quantile(Normal(), (1.0-coverage)/2)`
- `est`: Estimates
"""
function abc_postprocess(out, f::Function, tol=nothing;
    regress=false, kernel=out.cutoff, normal_quantile=1.96,
    iact=default_iact)
    if regress
        if out.Summary isa Missing
            error("Summary statistics not recorded; please re-run abc_mcmc with 'record_summaries=true'")
        end
        if tol isa Nothing
            error("Must supply tolerances with regression correction")
        end
        return abc_regress(out, f, tol; kernel=kernel.fun,
                                        normal_quantile=normal_quantile,
                                        iact=iact)
    else
        #F = mapslices(f, out.Theta, dims=1)
        F = eval_testfunction(out, f)
        return is_est(F, out.Dist, out.cutoff, kernel, tol;
                      normal_quantile=normal_quantile, iact=iact)
    end
end