text:
include("RBFBase.jl");
import .RBF
using LinearAlgebra: qr, I, norm, givens, Hermitian, cholesky, inv
using LinearAlgebra: diag # remove

@with_kw mutable struct RbfModel <: SurrogateModel
    model :: Union{Nothing,RBF.RBFModel} = nothing;

    # matrices for the sampling and improvement algorithm
    Y :: RMat = Matrix{Real}(undef, 0, 0);
    Z :: RMat = Matrix{Real}(undef, 0, 0);
    # matrices for the sampling and improvement algorithm
    Y₂ :: RMat = Matrix{Real}(undef, 0, 0);
    Z₂ :: RMat = Matrix{Real}(undef, 0, 0);

    fully_linear :: Bool = false;

    # for gathering training data in form of `AbstractResults`
    # dict key is the "construction round" for metadata (plotting)
    # dict values are vectors of training data
    tdata :: Dict{Int,Vector{<:AbstractResult}} = Dict( i => AbstractResult[] for i = 0:4);
end

function RbfModel( model :: RBF.RBFModel, n_vars :: Int )
    Y = Matrix{Real}(undef, n_vars, 0);
    Z = Matrix{Int}(I(n_vars));
    return RbfModel(;
        model = model,
        Y = Y,
        Y₂ = Y,
        Z = Z,
        Z₂ = Z
    )
end

tdata_length( rbf :: RbfModel ) :: Int = sum( length.(Base.values(rbf.tdata) ) )
db_indices( rbf :: RbfModel, i :: Int ) = convert( Vector{NothInt}, get_id.( rbf.tdata[i]) ) # WHY DO I HAVE TO CONVERT???
db_indices( rbf :: RbfModel ) = vcat( [ db_indices(rbf, i) for i ∈ Base.keys(rbf.tdata) ]... );
tdata_results(rbf :: RbfModel ) = vcat( (rbf.tdata[i] for i ∈ keys(rbf.tdata) )... );

"""
    RbfConfig(; kwarg1 = val1, … )

Configuration type for local RBF surrogate models.

To choose a kernel, use the kwarg `kernel` and a value of either
`:cubic` (default), `:multiquadric`, `:exp` or `:thin_plate_spline`.
The kwarg `shape_parameter` takes a constant number or a string
that defines a calculation on `Δ`, e.g, "Δ/10".

To see other configuration parameters use `fieldnames(Morbit.RbfConfig)`.
They have individual docstrings attached.
"""
@with_kw mutable struct RbfConfig <: SurrogateConfig
    "(default `:cubic`) RBF kernel (Symbol), either `:cubic`, `:multiquadric`, `:exp` or `:thin_plate_spline`."
    kernel :: Symbol = :cubic;
    "(default `1`) RBF shape paremeter, either a number or a string containing `Δ`."
    shape_parameter :: Union{AbstractString, Real} = 1;
    "(default `1`) Degree of polynomial attached to RBF. `-1` means no polynomial."
    polynomial_degree :: Int64 = 1;

    "(default `2`) Local enlargment factor of trust region for sampling."
    θ_enlarge_1 :: Real = 2;
    "(default `5`) Maximum enlargment factor of maximum trust region for sampling."
    θ_enlarge_2 :: Real = 2;  # reset
    "(default `1/(2*θ_enlarge_1)` Sampling parameter to generate Λ-poised set. The higher, the more poised."
    θ_pivot :: Real = 1 / (2 * θ_enlarge_1);
    "(default `1e-7`) Parameter for 2nd sampling algorithm to ensure boundedness of Cholesky factors."
    θ_pivot_cholesky :: Real = 1e-7;

    "(default `false`) Require models to be fully linear in each iteration."
    require_linear :: Bool = false;

    "(default `-1`) Maximum number of training sites. `-1` is reset to `2n+1`."
    max_model_points :: Int64 = -1; # is probably reset in the algorithm
    "(default `false`) Sample new sites to always use the maximum number of points."
    use_max_points :: Bool = false;

    "(default `:orthogonal`) Algorithm to use for finding affinely independent set."
    sampling_algorithm :: Symbol = :orthogonal # :orthogonal or :monte_carlo
    "(default `:standard_rand`) Algorithm to use if additional points are required."
    sampling_algorithm2 :: Symbol = :standard_rand

    "(default `typemax(Int64)`) Maximum number of objective evaluations."
    max_evals :: Int64 = typemax(Int64);

    @assert sampling_algorithm ∈ [:orthogonal, :monte_carlo] "Sampling algorithm must be either `:orthogonal` or `:monte_carlo`."
    @assert kernel ∈ Symbol.(["exp", "multiquadric", "cubic", "thin_plate_spline"]) "Kernel '$kernel' not supported yet."
    @assert kernel != :thin_plate_spline || shape_parameter isa Int && shape_parameter >= 1
    #@assert θ_enlarge_1 >=1 && θ_enlarge_2 >=1 "θ's must be >= 1."
end

# meta data object to be used during *sophisticated* sampling
@with_kw mutable struct RbfMeta <: SurrogateMeta
    center_index :: Int = 1;
    round1_indices :: Vector{NothInt} = [];
    round2_indices :: Vector{NothInt} = [];
    round3_indices :: Vector{NothInt} = [];
    round4_indices :: Vector{NothInt} = [];
    fully_linear :: Bool = false;
end

max_evals( cfg :: RbfConfig ) :: Int = cfg.max_evals;

fully_linear( rbf :: RbfModel ) :: Bool = rbf.fully_linear;

combinable( cfg :: RbfConfig ) :: Bool = true;

combine(cfg1 :: RbfConfig, :: RbfConfig) :: RbfConfig = cfg1;

eval_models( rbf :: RbfModel, x̂ :: Vec ) = RBF.output( rbf.model, x̂ );
eval_models( rbf :: RbfModel, x̂ :: RVec, ℓ :: Int ) = RBF.output( rbf.model, ℓ, x̂ );
get_gradient( rbf :: RbfModel, x̂ :: RVec, ℓ :: Int  ) = RBF.grad( rbf.model, ℓ, x̂ );
get_jacobian( rbf :: RbfModel, x̂ :: RVec) = RBF.jac( rbf.model, x̂ );

function _init_model( cfg ::RbfConfig, objf:: AbstractObjective, mop:: AbstractMOP,
    id :: AbstractIterData, ac :: AbstractConfig ) :: Tuple{ RbfModel, RbfMeta }

    n_vars = num_vars(mop);

    # let `update_model(…)` do all the work
    rbf = RbfModel();
    rmeta = RbfMeta();

    return update_model( rbf, objf, rmeta, mop, id, ac; ensure_fully_linear = true)
end

function update_model( rbf :: RbfModel, objf :: AbstractObjective, rmeta:: RbfMeta, mop :: AbstractMOP,
    id :: AbstractIterData, ac :: AbstractConfig;
    ensure_fully_linear :: Bool = false ) :: Tuple{ RbfModel, RbfMeta }

    @logmsg loglevel3 "Updating RBF with indices $(output_indices(objf, mop))."
    cfg = model_cfg(objf);
    DB = db(id);
    Δ = Δᵗ(id);
    x = xᵗ(id);
    n_vars = length(x);

    MAX_EVALS = min( max_evals(ac), max_evals(cfg) ) - 1 ;
    # we completly start from scratch here, nothing to recycle
    # this way, tdata is empty and properly initialized
    inner_model = RBF.RBFModel(;
        n_in = n_vars,
        kernel = cfg.kernel,
        polynomial_degree = cfg.polynomial_degree,
        shape_parameter = _get_shape_param(cfg, id),
    );

    rbf =  RbfModel(inner_model, n_vars);
    push!(rbf.tdata[0], get_result( id, xᵗ_index(id) ));

    rmeta = RbfMeta();

    # find affinely independent training points in radius Δ₁
    find_box_independent_points1!( rbf, cfg, id, mop, ac );

    n_missing = n_vars - length(rbf.tdata[1]);
    if n_missing > 0 && !(ensure_fully_linear || cfg.require_linear)
        # the model can be not fully linear
        # Search for more points in trust region of maximum possible radius
       find_box_independent_points2!(rbf, cfg, id, mop, ac)
    end

    n_missing -= length(rbf.tdata[2]);
    n_new = max(0, min( n_missing, MAX_EVALS) );

    add_new_sites!(rbf, cfg, mop, id; max_new = n_new);

    if length(rbf.tdata[3]) < n_new
        return rebuild_model( rbf, objf, rmeta, mop, id, ac);
    end

    rbf.fully_linear = isempty(rbf.tdata[2]) && length(rbf.tdata[1]) + length(rbf.tdata[3]) == n_vars
    @logmsg loglevel4 "RBF Model is $(rbf.fully_linear ? "" : "not ")fully linear."

    n_evals_left = MAX_EVALS - n_new;
    add_old_sites!( rbf, cfg, mop, id; max_new = n_evals_left);

    #Evaluate resuts...
    _eval_new_sites!( rbf, rmeta, objf, mop, id );
    RBF.train!(rbf.model)

    # TODO REMOVE!!
    #=
    for (s,v) ∈ zip(rbf.model.training_sites, rbf.model.training_values)
        @show eval_models(rbf, s ) .- v
    end
    =#

    return rbf, rmeta
end

function improve_model( rbf:: RbfModel, objf::AbstractObjective, rmeta :: RbfMeta,
    mop ::AbstractMOP, id :: AbstractIterData, ac:: AbstractConfig;
    ensure_fully_linear :: Bool = false ) :: Tuple{RbfModel, RbfMeta}
    if size(rbf.Z, 2) > 0
        @logmsg loglevel3 "Performing an improvement step for RBF with indices $(output_indices(objf, mop))."
        cfg = model_cfg(objf);

        x = xᵗ(id);
        Δ = Δᵗ(id);

        dir = copy( rbf.Z[:,1] );
        Δ₁ = cfg.θ_enlarge_1 * Δ;
        pivot = cfg.θ_pivot * Δ;

        len = intersect_bounds( mop, x, Δ₁, dir; return_vals = :absmax)
        if abs(len) > pivot
            offset = len .* dir;
            push!(rbf.tdata[3], init_res( Result, x .+ offset));

            rbf.Y = hcat( rbf.Y, offset );
            rbf.Z = rbf.Z[:, 2:end];

            rbf.fully_linear = size(rbf.Z, 2) == 0;
            _eval_new_sites!(rbf, rmeta, objf, mop, id);
            RBF.train!(rbf.model);
        else
            return rebuild_model( rbf, objf, rmeta, mop, id, ac);
        end
    end
    return rbf, rmeta;
end

function rebuild_model( rbf :: RbfModel, objf :: AbstractObjective, rmeta :: RbfMeta,
        mop :: AbstractMOP, id :: AbstractIterData, ac :: AbstractConfig ) :: Tuple{RbfModel, RbfMeta}
    @logmsg loglevel3 "Rebuild model along coordinates..."

    cfg = model_cfg(objf);
    x = xᵗ(id);
    n_vars = length(x);
    Δ = Δᵗ(id);
    Δ₁ = cfg.θ_enlarge_1 * Δ;

    MAX_EVALS = min( max_evals(ac), max_evals(cfg) ) - 1;
    max_new = min(MAX_EVALS, n_vars);

    rbf = RbfModel(;
        model = rbf.model,
        Y = Matrix{Real}(undef, n_vars, 0),
        Z = Matrix{Int}(I( n_vars )),
    );

    push!(rbf.tdata[0], get_result( id, xᵗ_index(id) ));

    for i = 1 : max_new
        dir = zeros(Int, n_vars);
        dir[i] = 1;
        len = intersect_bounds( mop, x, Δ₁, dir; return_vals = :absmax)
        offset = len .* dir;
        push!(rbf.tdata[3], init_res( Result, x .+ offset));
        rbf.Y = hcat( rbf.Y, offset );
        rbf.Z = rbf.Z[:, 2:end];
        if abs(len) <= cfg.θ_pivot * Δ
            @warn("Cannot make model fully linear!")
        end
    end#for
    rbf.Y₂ = rbf.Y;
    rbf.Z₂ = rbf.Z;

    n_evals_left = MAX_EVALS - max_new;
    add_old_sites!( rbf, cfg, mop, id; max_new = n_evals_left);

    _eval_new_sites!(rbf, rmeta, objf, mop, id );
    RBF.train!(rbf.model);
    return rbf, rmeta
end

# optional helper
num_outputs( rbf :: RbfModel ) :: Int = Rbf.num_outputs( rbf.model );

function _transfer_training_sites!( rbf :: RbfModel ) :: Nothing
    basis_results = tdata_results(rbf);
    empty!(rbf.model.training_sites)
    push!(rbf.model.training_sites, [ get_site(res) for res ∈ basis_results ]... ) ;
    nothing
end

function _transfer_training_values!( rbf :: RbfModel, out_indices :: Vector{Int}) :: Nothing
    basis_results = tdata_results(rbf);
    empty!(rbf.model.training_values)
    push!(rbf.model.training_values, [ get_value(res)[out_indices] for res ∈ basis_results ]... );
    nothing
end

# evaluate at new sites
# set `training_sites` and `training_values` for inner model
# set meta data
function _eval_new_sites!( rbf :: RbfModel, rmeta :: RbfMeta, objf :: AbstractObjective,
    mop :: AbstractMOP, id :: AbstractIterData )  :: Nothing

    basis_results = vcat( (rbf.tdata[i] for i = 1:4 )... )
    _eval_and_store_new_results!(id, basis_results, mop);

    NotXint = Union{Nothing,Int};
    rmeta.center_index = convert( NotXint, xᵗ_index(id) );
    rmeta.round1_indices = [ convert( NotXint, get_id( res ) ) for res ∈ rbf.tdata[1] ] ;
    rmeta.round2_indices = [ convert( NotXint, get_id( res ) ) for res ∈ rbf.tdata[2] ] ;
    rmeta.round3_indices = [ convert( NotXint, get_id( res ) ) for res ∈ rbf.tdata[3] ] ;
    rmeta.round4_indices = [ convert( NotXint, get_id( res ) ) for res ∈ rbf.tdata[4] ] ;

    # "transfer" training data to inner model
    _transfer_training_sites!(rbf);
    _transfer_training_values!(rbf, output_indices( objf, mop));
    nothing
end

function parse_shape_param_string( Δ, expr_str)
    ex = Meta.parse(expr_str)
    return @eval begin
        let Δ=$Δ
            $ex
        end
    end
end

"Get real-valued shape parameter for RBF model from current iter data.
`cfg` allows for a string expression which would be evaluated here."
function _get_shape_param( cfg :: RbfConfig, id :: AbstractIterData ) :: Real
    if isa(cfg.shape_parameter, AbstractString)
        return parse_shape_param_string( Δᵗ(id), cfg.shape_parameter )
    elseif isa(cfg.shape_parameter, Real)
        return cfg.shape_parameter
    end
end

## Construction

function _orthogonal_complement_matrix( Y :: RMat )
    Q, _ = qr(Y);
    Z = Q[:, size(Y,2) + 1 : end];
    if size(Z,2) > 0
        Z ./= norm.( eachcol(Z), Inf )';
    end
    return Z
end

"""
Find affinely independent results in database box of radius `Δ` around `x`
AbstractResults are saved in `rbf.tdata[tdata_index]`.
Both `rbf.Y` and `rbf.Z` are changed.
"""
function find_box_independent_points1!( rbf :: RbfModel, cfg :: RbfConfig, id :: AbstractIterData,
    mop :: AbstractMOP, :: AbstractConfig ) :: Nothing

    res, Y, Z = _find_box_independent_points(
        rbf, id, mop, xᵗ(id), Δᵗ(id) * cfg.θ_enlarge_1;
        θ_pivot = cfg.θ_pivot
    );
    push!(rbf.tdata[1], res...);
    rbf.Y₂ = rbf.Y = Y;
    rbf.Z₂ = rbf.Z = Z;
    nothing
end

function find_box_independent_points2!( rbf :: RbfModel, cfg :: RbfConfig, id :: AbstractIterData,
    mop :: AbstractMOP, ac :: AbstractConfig )

    Δ₂ = Δᵘ(ac) * cfg.θ_enlarge_1;

    res, Y, Z = _find_box_independent_points(
        rbf, id, mop, xᵗ(id), Δ₂;
        θ_pivot = (Δᵗ(id) * cfg.θ_enlarge_1 / Δ₂) * cfg.θ_pivot,
        exclude_indices = [ xᵗ_index(id); db_indices(rbf, 1) ]
    );

    push!(rbf.tdata[2], res...);
    rbf.Y₂ = Y;
    rbf.Z₂ = Z;
    nothing
end

# if we had more than two rounds we would have to pass Y & Z here too
function _find_box_independent_points( rbf :: RbfModel,
    id :: AbstractIterData, mop :: AbstractMOP, x :: RVec, Δ :: Union{Real,RVec};
    θ_pivot = 1e-3, exclude_indices :: Vector{<:NothInt} = NothInt[] ) :: Tuple{Vector{<:AbstractResult},RMat,RMat}

    lb, ub = local_bounds(mop, x, Δ );
    box_indices = find_points_in_box( id, lb, ub; exclude_indices );

    box_results = get_result.( id, box_indices );

    pivot = Δ * θ_pivot;

    affinely_independent_box_indices, Y, Z = _affinely_independent_points(
        box_results,
        x,
        rbf.Y,
        rbf.Z,
        piv_val = pivot
    )

    @logmsg loglevel4 "Found $(length(affinely_independent_box_indices)) points in box of radius $(Δ) with pivot $(pivot)."
    return (box_results[affinely_independent_box_indices], Y, Z);
end

function _affinely_independent_points( list_of_results :: Vector{<:AbstractResult}, x₀ :: RVec,
    Y :: RMat = Matrix{Real}(undef, 0, 0), Z :: RMat = Matrix{Real}(undef, 0, 0);
    piv_val :: Real = 1e-3 ) :: Tuple{Vector{<:Int},RMat,RMat}
    return _affinely_independent_points( get_site.(list_of_results), x₀, Y, Z; piv_val )
end

function _affinely_independent_points( list_of_points :: RVecArr, x₀ :: RVec,
    Ŷ :: RMat = Matrix{Real}(undef, 0, 0), Ẑ :: RMat = Matrix{Real}(undef, 0, 0);
    piv_val :: Real = 1e-3 ) :: Tuple{Vector{<:Int},RMat,RMat}

    n_vars = length(x₀);
    Y = copy(Ŷ);    # because we are possible modifying the matrices
    Z = copy(Ẑ);

    if isempty( list_of_points )
        return Int[], Matrix{Real}(undef, n_vars, 0), Matrix{Int}(I(n_vars))
    else
        num_points = length(list_of_points);

        rY, cY = size(Y);
        if rY == 0
            Y = Matrix{Real}(undef, n_vars, 0)
            Z = Matrix{Int}(I(n_vars));
        end

        # TODO remove this
        @assert size(Y,1) == n_vars;
        @assert size(Y,2) + size(Z,2) == n_vars;
        @assert all( isapprox.(norm.(eachcol(Z), Inf),1 ) )

        num_missing = n_vars - size(Y,2);

        accepted_indices = Int[];
        candidate_indices = collect(eachindex(list_of_points));
        for i = 1 : num_points
            # mapping to translate a vector χ and project it
            # onto the orthogonal complement of Y, spanned by columns of Z
            proj =  χ -> norm(Z*(Z'*(χ.-x₀)),Inf)
            # find maximizer of projections
            # NOTE: this inner loop is not mandatory
            best_val, best_index = findmax( proj.( list_of_points[ candidate_indices ] ) );

            if best_val > piv_val
                true_index = candidate_indices[best_index];
                push!(accepted_indices, true_index );
                deleteat!(candidate_indices, best_index)
                Y = hcat( Y, list_of_points[ true_index ] .- x₀ );
                Z = _orthogonal_complement_matrix(Y);
            end
            if length(accepted_indices) == num_missing
                break;
            end
        end# outer for
        return accepted_indices, Y, Z
    end# if list_of_points empty
end

function add_new_sites!(rbf :: RbfModel, cfg :: RbfConfig, mop :: AbstractMOP,
    id :: AbstractIterData; max_new :: Int = typemax(Int) ) :: Nothing
    return add_new_sites!(rbf, cfg, mop, id, Val( cfg.sampling_algorithm ); max_new )
end

function add_new_sites!(rbf :: RbfModel, cfg :: RbfConfig, mop :: AbstractMOP,
    id :: AbstractIterData, ::Val{:orthogonal}; max_new :: Int ) :: Nothing
    if max_new > 0
        x = xᵗ(id);
        Δ = Δᵗ(id);
        n_vars = length(x);

        Δ₁ = Δ * cfg.θ_enlarge_1;
        for dir = eachcol( rbf.Z )
            if length(rbf.tdata[3]) >= max_new
                break;
            end
            @assert norm(dir, Inf) ≈ 1;
            len = intersect_bounds( mop, x, Δ₁, copy(dir); return_vals = :absmax)
            if abs(len) > cfg.θ_pivot * Δ
                # pivot big enough, accept new site
                offset = len .* dir;
                push!(rbf.tdata[3], init_res( Result, x .+ offset));
                rbf.Y = hcat( rbf.Y, offset );
                rbf.Z = rbf.Z[:, 2:end];
            else
                # could not sample far enough from current iterate, break
                break;
            end

        end
    end
    return nothing;
end

function add_old_sites!(rbf :: RbfModel, cfg :: RbfConfig, mop :: AbstractMOP, id :: AbstractIterData; kwargs... ) :: Nothing
    return add_old_sites!(rbf, cfg, mop, id, Val( cfg.sampling_algorithm2); kwargs...)
end

function add_old_sites!(rbf :: RbfModel, cfg :: RbfConfig, mop :: AbstractMOP,
    id :: AbstractIterData, ::Val{:standard_rand};
    max_new :: Int = typemax(Int) ):: Nothing
    x = xᵗ(id);
    Δ₂ = cfg.θ_enlarge_2 * Δᵗ(id);
    n_vars = length(x);

    chol_pivot = cfg.θ_pivot_cholesky;

    N₀ = N = tdata_length( rbf )
    MAX_POINTS = cfg.max_model_points <= 0 ? 2*n_vars+1 : cfg.max_model_points;

    if N < MAX_POINTS
        # make sure, training sites are set
        _transfer_training_sites!(rbf);

        # prepare matrices as by Wild, R has to be augmented by rows of zeros
        Φ = RBF.get_Φ( rbf.model );
        Π = RBF.get_Π( rbf.model );
        Q, R = qr( transpose(Π) );
        R = [
            Matrix(R);
            zeros( size(Q,1) - size(R,1), size(R,2) )
        ];
        Z = Q[:, RBF.min_num_sites(rbf.model) + 1 : end ];

        ZΦZ = Hermitian(Z'Φ*Z);
        L = cholesky( ZΦZ ).L     # should also be empty at this point
        Lⁱ = inv(L);

        φ₀ = Φ[1,1]; # constant value

        # TEST FOR ADDITIONAL POINTS
        # we reverse the order of database indices because in our
        # implementation this favors recently added points, i.e.,
        # points that are near the current iterate. other heuristics
        # are very much possible
        DB = db(id);
        all_ids = sort( collect( eachindex( DB ) ); rev = true );
        training_ids = db_indices( rbf )

        lb, ub = local_bounds(mop, x, Δ₂);

        for id ∈ all_ids
            if N >= MAX_POINTS
                break;
            end
            if id ∉ training_ids
                old_res = get_result(DB, id);
                x̂ = get_site( old_res );
                if _point_in_box( x̂, lb, ub)
                    test_retval = _test_new_cholesky_site( rbf.model, x̂, φ₀, Q, R, Z, L, Lⁱ, Φ, Π, chol_pivot );
                    if !isnothing(test_retval)
                        Q, R, Z, L, Lⁱ, Φ, Π = test_retval;
                        push!(rbf.tdata[4], old_res );
                        push!(rbf.model.training_sites, x̂ ); # so that the _test… works properly
                        N +=1 ; # increase point counter
                    end
                end
            end
        end
        @logmsg loglevel4 "Found $(N - N₀) additional sites to add to the model."
        if max_new > 0 && cfg.use_max_points && N < MAX_POINTS
            n_new = min(MAX_POINTS, max_new);
            max_tries = 10*n_new;
            n_tries = 0;
            @logmsg loglevel4 "`use_max_points` … trying to find $(n_new - N) new sites."
            while N < n_new && n_tries < max_tries;
                # simple hit-or-miss, not very effective nor good for model quality
                x̂ = _rand_box_point( lb, ub )
                test_retval = _test_new_cholesky_site( rbf.model, x̂, φ₀, Q, R, Z, L, Lⁱ, Φ, Π, chol_pivot );
                if !isnothing(test_retval)
                    Q, R, Z, L, Lⁱ, Φ, Π = test_retval;
                    push!(rbf.tdata[4], init_res(Result, x̂));
                    push!(rbf.model.training_sites, x̂ ); # so that the _test… works properly
                    N +=1 ; # increase point counter
                end
                n_tries += 1;
            end
        end
    end
    nothing
end

function _test_new_cholesky_site(m :: RBF.RBFModel, y :: RVec, φ₀ :: Real,
    Q :: AbstractArray, R :: AbstractArray,  Z :: AbstractArray, L :: AbstractArray,
    Lⁱ :: AbstractArray, Φ :: AbstractArray, Π :: AbstractArray,  θ_pivot_cholesky :: Real )

    # evaluate RBF basis at new site
    φy = RBF.φ(m, y)
    Φy = [
        [Φ φy];
        [φy' φ₀]
    ]

    # get polynomial basis column and stack transpose below R
    πy = RBF.Π_col( m, y );
    Ry = [
        R ;
        πy'
    ]

    # perform givens rotations to turn last row in Ry to zeros
    row_index = size( Ry, 1)
    G = Matrix(I, row_index, row_index) # whole orthogonal matrix
    for j = 1 : size(R,2)
        # in each column, take the diagonal as pivot to turn last elem to zero
        g = givens( Ry[j,j], Ry[row_index, j], j, row_index )[1];
        Ry = g*Ry;
        G = g*G;
    end
    # now, from G we can update the other matrices
    Gᵀ = transpose(G);
    g̃ = Gᵀ[1 : end-1, end];
    ĝ = Gᵀ[end, end];

    Qg = Q*g̃;
    v_y = Z'*( Φ*Qg + φy .* ĝ );
    σ_y = Qg'*Φ*Qg + (2*ĝ) * φy'*Qg + ĝ^2*φ₀;

    τ_y² = σ_y - norm( Lⁱ * v_y, 2 )^2
    # τ_y (and hence τ_y^2) must be bounded away from zero
    # for the model to remain fully linear
    if τ_y² > θ_pivot_cholesky
        τ_y = sqrt(τ_y²)

        Qy = [
            [ Q zeros( size(Q,1), 1) ];
            [ zeros(1, size(Q,2)) 1 ]
        ] * Gᵀ

        z = [
            Q * g̃;
            ĝ
        ]

        Zy = [
            [
                Z;
                zeros(1, size(Z,2))
            ] z
        ]

        Lyⁱ = [
            [Lⁱ zeros(size(Lⁱ,1),1)];
            [ -(v_y'Lⁱ'Lⁱ)./τ_y 1/τ_y ]
        ];

        Ly = [
            [ L zeros(size(L,1), 1) ];
            [ v_y'Lⁱ' τ_y ]
        ]

        Πy = [ Π πy ];

        @assert diag(Ly * Lyⁱ) ≈ ones(size(Ly,1))

        return Qy, Ry, Zy, Ly, Lyⁱ, Φy, Πy
    else
        return nothing
    end
end


text:
# Julia wrapper for header: /usr/include/scip/cons_bounddisjunction.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrBounddisjunction(scip)
    ccall((:SCIPincludeConshdlrBounddisjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsBounddisjunction(scip, cons, name, nvars, vars, boundtypes, bounds, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsBounddisjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{SCIP_BOUNDTYPE}, Ptr{Cdouble}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, boundtypes, bounds, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicBounddisjunction(scip, cons, name, nvars, vars, boundtypes, bounds)
    ccall((:SCIPcreateConsBasicBounddisjunction, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{SCIP_BOUNDTYPE}, Ptr{Cdouble}), scip, cons, name, nvars, vars, boundtypes, bounds)
end

function SCIPgetNVarsBounddisjunction(scip, cons)
    ccall((:SCIPgetNVarsBounddisjunction, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsBounddisjunction(scip, cons)
    ccall((:SCIPgetVarsBounddisjunction, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetBoundtypesBounddisjunction(scip, cons)
    ccall((:SCIPgetBoundtypesBounddisjunction, libscip), Ptr{SCIP_BOUNDTYPE}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetBoundsBounddisjunction(scip, cons)
    ccall((:SCIPgetBoundsBounddisjunction, libscip), Ptr{Cdouble}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

ext: jl


text:
mutable struct DAQConfig
    "Name for referening the device"
    devname::String
    "IP address where the device is located"
    ip::String
    "Model of device"
    model::String
    "Serial number of the device"
    sn::String
    "Storage tag of the device"
    tag::String
    "Integer configuration parameters"
    ipars::Dict{String,Int}
    "Float64  configuration parameters"
    fpars::Dict{String,Float64}
    "String configuration parameters"
    spars::Dict{String,String}
end

"""
`DAQConfig(;devname="", ip="", model="", sn="",tag="")`

Create data structure to store device configuration.

Hopefuly, with the information stored in this struct, reproducible
data configuration of the devices can be achieved.

Besides device characterization information (`ip`, `model`, `sn`, `tag`), the configuration
is stored in three dictionaries:

 * `ipars`: integer parameters
 * `fpars`: floating point parameters
 * `spars`: string parameters.

"""
function DAQConfig(;devname="", ip="", model="", sn="",tag="")

    fpars = Dict{String,Float64}()
    ipars = Dict{String,Int}()
    spars = Dict{String,String}()
    return DAQConfig(devname, ip, model, sn, tag, ipars, fpars, spars)
end

function DAQConfig(ipars, fpars, spars;devname="", ip="", model="", sn="",tag="")

    return DAQConfig(devname, ip, model, sn, tag, ipars, fpars, spars)
end

"Retrieve integer configuration parameter"
iparameters(dconf::DAQConfig, param) = dconf.ipars[param]
"Retrieve string configuration parameter"
sparameters(dconf::DAQConfig, param) = dconf.spars[param]
"Retrieve float configuration parameter"
fparameters(dconf::DAQConfig, param) = dconf.fpars[param]

iparameters(dconf::DAQConfig) = dconf.ipars
sparameters(dconf::DAQConfig) = dconf.spars
fparameters(dconf::DAQConfig) = dconf.fpars

"Retrieve device name"
devname(dconf::DAQConfig) = dconf.devname
"Retrieve device IP address"
daqdevip(dconf::DAQConfig) = dconf.ip
"Retrieve device model"
daqdevmodel(dconf::DAQConfig) = dconf.model
"Retrieve device serial number"
daqdevserialnum(dconf::DAQConfig) = dconf.sn
"Retrieve device tag"
daqdevtag(dconf::DAQConfig) = dconf.tag

iparameters(dev::AbstractDAQ, param) = dev.conf.ipars[param]
sparameters(dev::AbstractDAQ, param) = dev.conf.spars[param]
fparameters(dev::AbstractDAQ, param) = dev.conf.fpars[param]

iparameters(dev::AbstractDAQ) = dev.conf.ipars
sparameters(dev::AbstractDAQ) = dev.conf.spars
fparameters(dev::AbstractDAQ) = dev.conf.fpars

devname(dev::AbstractDAQ) = dev.devname
daqdevip(dev::AbstractDAQ) = dev.conf.ip
daqdevmodel(dev::AbstractDAQ) = dev.conf.model
daqdevserialnum(dev::AbstractDAQ) = dev.conf.sn
daqdevtag(dev::AbstractDAQ) = dev.conf.tag


text:
function get_window_bounds(img::GenericGrayImage, pixel::CartesianIndex{2}, w::Integer)
    row, col = pixel.I
    min_row, min_col = first.(axes(img))
    max_row, max_col = last.(axes(img))

    row₀ = max(row - w, min_row)
    col₀ = max(col - w, min_col)
    row₁ = min(row + w, max_row)
    col₁ = min(col + w, max_col)

    return row₀, col₀, row₁, col₁
end

function μ_in_window(I::AbstractArray{T,2}, row₀::Integer, col₀::Integer, row₁::Integer, col₁::Integer) where T <: Real
    μ = boxdiff(I, row₀:row₁, col₀:col₁) / ((row₁-row₀+1)*(col₁-col₀+1))
    return μ < eps() ? zero(eltype(I)) : μ
end

function σ²_in_window(I²::AbstractArray{T,2}, μ::Real, row₀::Integer, col₀::Integer, row₁::Integer, col₁::Integer) where T <: Real
    σ² = μ_in_window(I², row₀, col₀, row₁, col₁) - μ^2
    return σ² < eps() ? zero(eltype(I²)) : σ²
end

function σ_in_window(I²::AbstractArray{T,2}, μ::Real, row₀::Integer,  col₀::Integer, row₁::Integer, col₁::Integer) where T <: Real
    σ² = σ²_in_window(I², μ, row₀, col₀, row₁, col₁)
    return σ² < eps() ? zero(eltype(I²)) : sqrt(σ²)
end


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
using MAT
np = pyimport("numpy")

λ = 2.0
μ = 0.5
η = 1/12

β = 1/4; γ = 1/2
a = b = 0.1
m = 10
n = 5
h = 0.01
NT = 20
Δt = 1/NT
bdedge = []
for j = 1:n
  push!(bdedge, [(j-1)*(m+1)+m+1 j*(m+1)+m+1])
end
bdedge = vcat(bdedge...)

bdnode = Int64[]
for j = 1:n+1
  push!(bdnode, (j-1)*(m+1)+1)
end

G = [1/Δt+μ/η -μ/3η 0.0
  -μ/3η 1/Δt+μ/η-μ/3η 0.0
  0.0 0.0 1/Δt+μ/η]
S = [2μ/Δt+λ/Δt λ/Δt 0.0
    λ/Δt 2μ/Δt+λ/Δt 0.0
    0.0 0.0 μ/Δt]
invG = inv(G)
H = invG*S

M = compute_fem_mass_matrix1(m, n, h)
Zero = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = [M Zero;Zero M]

K = compute_fem_stiffness_matrix(H, m, n, h)
C = a*M + b*K # damping matrix

L = M + γ*Δt*C + β*Δt^2*K
L, Lbd = fem_impose_Dirichlet_boundary_condition(L, bdnode, m, n, h)

Forces = zeros(NT, 2(m+1)*(n+1))
for i = 1:NT
  T = eval_f_on_boundary_edge((x,y)->0.1, bdedge, m, n, h)
  T = [T zeros(length(T))]
  rhs = compute_fem_traction_term(T, bdedge, m, n, h)

  if i*Δt>0.5
    rhs = zero(rhs)
  end
  Forces[i, :] = rhs
end

a = zeros(2(m+1)*(n+1))
v = zeros(2(m+1)*(n+1))
d = zeros(2(m+1)*(n+1))
U = zeros(2(m+1)*(n+1),NT+1)
Sigma = zeros(NT+1, 4m*n, 3)
Varepsilon = zeros(NT+1, 4m*n, 3)
for i = 1:NT
  @info i
    global a, v, d
    F = compute_strain_energy_term(Sigma[i,:,:]*invG/Δt, m, n, h) - K * U[:,i]
    # @show norm(compute_strain_energy_term(Sigma[i,:,:]*invG/Δt, m, n, h)), norm(K * U[:,i])
    rhs = Forces[i,:] - Δt^2 * F

    td = d + Δt*v + Δt^2/2*(1-2β)*a
    tv = v + (1-γ)*Δt*a
    rhs = rhs - C*tv - K*td
    rhs[[bdnode; bdnode.+(m+1)*(n+1)]] .= 0.0

    a = L\rhs
    d = td + β*Δt^2*a
    v = tv + γ*Δt*a
    U[:,i+1] = d

    Varepsilon[i+1,:,:] = eval_strain_on_gauss_pts(U[:,i+1], m, n, h)
    Sigma[i+1,:,:] = Sigma[i,:,:]*invG/Δt +  (Varepsilon[i+1,:,:]-Varepsilon[i,:,:])*(invG*S)
end

matwrite("U.mat", Dict("U"=>U'|>Array))

# visualize_displacement(U, m, n, h; name = "_eta$η", xlim_=[-0.01,0.5], ylim_=[-0.05,0.22])
# visualize_displacement(U, m, n, h;  name = "_viscoelasticity")
# visualize_stress(H, U, m, n, h;  name = "_viscoelasticity")

close("all")
figure(figsize=(15,5))
subplot(1,3,1)
idx = div(n,2)*(m+1) + m+1
plot((0:NT)*Δt, U[idx,:])
xlabel("time")
ylabel("\$u_x\$")

subplot(1,3,2)
idx = 4*(div(n,2)*m + m)
plot((0:NT)*Δt, Sigma[:,idx,1])
xlabel("time")
ylabel("\$\\sigma_{xx}\$")

subplot(1,3,3)
idx = 4*(div(n,2)*m + m)
plot((0:NT)*Δt, Varepsilon[:,idx,1])
xlabel("time")
ylabel("\$\\varepsilon_{xx}\$")
savefig("visco_eta$η.png")


text:
# AbstractMatrix wrapper around Petsc Mat
export Mat, petscview, SubMat

"""
  A Petsc matrix.

  Unlike Vecs, the Petsc implementation keeps track of the local assembly
  state, so the Julia type does not have to.
  `verify_assembled`: if true, verify all processes are assembled, if false,
                      only local process
  `insertmode`: C.InsertMode used by `setindex!`
"""
abstract type PetscMat{T} <: AbstractSparseMatrix{T, PetscInt} end
mutable struct Mat{T} <: PetscMat{T}
  p::C.Mat{T}
  verify_assembled::Bool # check all processes assembled state or just current
  insertmode::C.InsertMode # current mode for setindex!
  data::Any # keep a reference to anything needed for the Mat
            # -- needed if the Mat is a wrapper around a Julia object,
            #    to prevent the object from being garbage collected.
  function Mat{T}(p::C.Mat{T}, data=nothing; first_instance::Bool=true, verify_assembled=true) where {T}
    A = new{T}(p, verify_assembled, C.INSERT_VALUES, data)
    if first_instance  # if the pointer p has not been put into a Mat before
      if (!deactivate_finalizers)
        finalizer(PetscDestroy, A)
      end
    end

    return A
  end
end

mutable struct SubMat{T} <: PetscMat{T}
  p::C.Mat{T}
  verify_assembled::Bool
  insertmode::C.InsertMode # current mode for setindex!
  data::Any # keep a reference to anything needed for the Mat
            # -- needed if the Mat is a wrapper around a Julia object,
            #    to prevent the object from being garbage collected.
            # in general, we *must* keep a copy of the parent matrix, because
            # creating a submatrix does not increase the reference count
            # in Petsc (for some unknown reason)
  function SubMat(p::C.Mat{T}, data=nothing; verify_assembled=true) where {T}
    A = new{T}(p, verify_assembled, C.INSERT_VALUES, data)
    if (!deactivate_finalizers)
      finalizer(SubMatRestore, A)
    end
    return A
  end
end

"""
  Get the communicator for the object
"""
comm(a::PetscMat{T}) where {T} = C.PetscObjectComm(T, a.p.pobj)

"""
  Create an empty, unsized matrix
"""
function Mat(::Type{T}, mtype::C.MatType; comm::MPI.Comm=MPI.COMM_WORLD) where {T}
  p = Ref{C.Mat{T}}()
  chk(C.MatCreate(comm, p))
  chk(C.MatSetType(p[], mtype))
  Mat(p[])
end

"""
  Create a matrix of a particular size, optionally specifying the pre-allocation.
  If pre-allocation is not specified, no preallocation is done
"""
function Mat(::Type{T}, m::Integer=C.PETSC_DECIDE, n::Integer=C.PETSC_DECIDE;
                mlocal::Integer=C.PETSC_DECIDE, nlocal::Integer=C.PETSC_DECIDE,
                bs=1, nz::Integer=0, nnz::AbstractVector=PetscInt[],
                onz::Integer=0, onnz::AbstractVector=PetscInt[],
                comm::MPI.Comm=MPI.COMM_WORLD,
                mtype::Symbol=C.MATMPIAIJ) where {T}

  mat = Mat(T, mtype, comm=comm)
  resize!(mat, m, n, mlocal=mlocal, nlocal=nlocal)

  # don't preallocate unless the user specified to do so, because
  # after pre-allocation it is an error to try to change the sparsity pattern
  # Before preallocation Petsc will dynamically allocate memory as needed
  # slow but flexible.
  if nz==0 && onz == 0  && nnz == PetscInt[] && onnz == PetscInt[]
    if bs != 1
      set_block_size(mat, bs)
    end
    chk(C.MatSetUp(mat.p))
  else  # preallocate
    setpreallocation!(mat, nz=nz, nnz=nnz, onz=onz, onnz=onnz, bs=bs)
  end
  setoption!(mat, C.MAT_ROW_ORIENTED, false)  # julia data is column major

  return mat
end

"""
  Make a MATSEQ Petsc matrix for a SparseMatrixCSC.  This preserve the
  sparsity pattern of the matrix
"""
function Mat(A::SparseMatrixCSC{T}) where {T}
  m, n = size(A)

  # count number of non-zeros in each row
  nz = zeros(PetscInt, m)
  for i=1:n  # loop over columns
    idx_start = A.colptr[i]
    idx_end = A.colptr[i+1]-1
    for j=idx_start:idx_end
      rownum = A.rowval[j]
      nz[rownum] += 1
    end
  end

  # create the matrix
  PA = Mat(T, m, n, nnz=nz, mtype=C.MATSEQAIJ)

  # copy values
  # create array large enough to hold all values in a column
  maxcol = 0
  for i=1:n
    nvals = A.colptr[i+1] - A.colptr[i]
    if nvals > maxcol
      maxcol = nvals
    end
  end
  idx = Array{PetscInt,1}(undef,maxcol)
  idy = PetscInt[0]

  for i=1:n
    idy[1] = i-1
    idx_start = A.colptr[i]
    idx_end = A.colptr[i+1]-1
    # convert indices to PetscInts, zero-based
    pos=1
    for j=idx_start:idx_end
      idx[pos] = A.rowval[j] - 1
      pos += 1
    end

    # get subarray of only the needed entries
    idx_view = view(idx, 1:(pos-1))
    vals = view(A.nzval, idx_start:idx_end)
    set_values!(PA, idx_view, idy, vals)
  end

  return PA
end

"""
  Construct at MATSEQAIJ from an AbstractArray.  The argument droptol is
  used to determine what size entry is considered non-zero
"""
function Mat(A::AbstractArray{T}; droptol=0.0) where {T}

  m, n = size(A)

  # count non-zeros
  nz = zeros(PetscInt, m)
  for i=1:n
    for j=1:m
      val = A[j, i]
      if abs(val) > droptol
        nz[j] += 1
      end
    end
  end

  # create matrix
  PA = Mat(T, m, n, nnz=nz, mtype=C.MATSEQAIJ)

  # copy values
  maxrow = maximum(nz)
  idx = PetscInt[0]
  idy = zeros(PetscInt, maxrow)
  vals = zeros(T, maxrow)
  for i=1:m
    idx[1] = i-1
    pos = 1
    for j=1:n
      val = A[i, j]
      if abs(val) > droptol
        idy[pos] = j - 1
        vals[pos] = val
        pos += 1
      end
    end
    idy_view = view(idy, 1:(pos-1))
    set_values!(PA, idx, idy_view, vals)
  end

  return PA
end

"""
  Gets the a submatrix that references the entries in the original matrix.
  isrow and iscol contain the *local* indicies of the rows and columns to get.
  The matrix must have a LocalToGlobalMapping for this to work, therefore a
  default one is created if the matrix does not already have one registered.
  The default mapping assumes the matrix is divided up into contiguous block
  of rows.  This is true of AIJ matrices but may not be for other matrix types.
"""
function SubMat(mat::Mat{T}, isrow::IS{T}, iscol::IS{T}) where {T}

  # create the local to global mapping for mat first
  # this only needs to be done the first time
  if !has_local_to_global_mapping(mat)
    rmap, cmap = local_to_global_mapping(mat)
    set_local_to_global_mapping(mat, rmap, cmap)
  end

  # now we can actually create the submatrix
  submat = Ref{C.Mat{T}}()
  chk(C.MatGetLocalSubMatrix(mat.p, isrow.p, iscol.p, submat))
  # keep the data needed for the finalizer
  return SubMat(submat[], (mat, isrow, iscol), verify_assembled=mat.verify_assembled)
end

export SubMatRestore
function SubMatRestore(smat::SubMat{T}) where {T}

  C.MatRestoreLocalSubMatrix(smat.data[1].p, smat.data[2].p, smat.data[3].p, Ref(smat.p))
end

##### MatShell functions #####
export MatShell, setop!, getcontext
"""
  Create a high level matrix from an already created matrix pointer
"""
function Mat(ptr::C.Mat{T}) where {T}
# this is not type stable. Grr
  return Mat{T}(ptr, nothing, first_instance=false)
end

"""
  Create a shell matrix with specified size.  The ctx tuple contains can be
  accessed by any callback function.
"""
# must rename this because of (silent) method ambiguity
function MatShell(::Type{T}, mlocal::Integer, nlocal::Integer, ctx::Tuple=();  m::Integer=C.PETSC_DECIDE, n::Integer=C.PETSC_DECIDE, comm=MPI.COMM_WORLD) where {T}

  mat_ptr = Ref{C.Mat{T}}()
  #TODO: fix this
  #ctx_ptr = pointer_from_objref(ctx)
  #chk(C.MatCreateShell(comm, mlocal, nlocal, m, n, ctx_ptr, mat_ptr))
  return Mat{T}(mat_ptr[], ctx)  # protect ctx from gc
end

"""
  Provide a callback function for a particular matrix operation.  op is a
  Petsc enum value inidcating the operation, and func is a void pointer
  (obtained from cfunction() ) that performs the operation.

  The function should take the low level Petsc objects (defined in the C module)
  rather than the high level ones defined in this file.  There are constructors
  to create a high level object from a low level one

"""
function setop!(mat::Mat{T}, op::C.MatOperation, func::Ptr{Cvoid}) where {T}
  if gettype(mat) != C.MATSHELL
    throw(ArgumentError("Mat must be a MatShell"))
  end

  chk(C.MatShellSetOperation(mat.p, op, func))
end

"""
  Get the tuple of user provided data passed in when the shell matrix was
  created.
"""
function getcontext(mat::Mat)
  if gettype(mat) != C.MATSHELL
    throw(ArgumentError("Mat must be a MatShell"))
  end

  ctx_ptr = C.MatShellGetContext(mat.p)
  return unsafe_pointer_to_objref(ctx_ptr)
end

"""
  Destroy a Mat object and the underlying data structure, if the object
  has not already been finalized
"""
function PetscDestroy(mat::Mat{T}) where {T}
  if !PetscFinalized(T)
    C.MatDestroy(Ref(mat.p))
    mat.p = C.Mat{T}(C_NULL)  # indicate the vector is finalized
  end
end

"""
  Check if PetscDestroy has been called on this object already
"""
function isfinalized(mat::PetscMat)
  return isfinalized(mat.p)
end

function isfinalized(mat::C.Mat)
  return mat.pobj == C_NULL
end

"""
  Print a Petsc matrix to STDOUT
"""
function petscview(mat::PetscMat{T}) where {T}
  viewer = C.PetscViewer{T}(C_NULL)
  chk(C.MatView(mat.p, viewer))
end

"""
  Print a Petsc matrix to a named file, in text format
"""
function petscwrite(mat::PetscMat{T}, fname) where {T}
  viewer_ref = Ref{C.PetscViewer{T}}()
  chk(C.PetscViewerASCIIOpen(comm(mat), fname, viewer_ref))
  chk(C.MatView(mat.p, viewer_ref[]))
  chk(C.PetscViewerDestroy(viewer_ref))
end

function set_block_size(A::Mat{T}, bs::Integer) where {T<:Scalar}
  chk(C.MatSetBlockSize(A.p, PetscInt(bs)))
end

function get_blocksize(A::Mat{T}) where {T<:Scalar}
  bs = Ref{PetscInt}()
  chk(C.MatGetBlockSize(A.p, bs))
  return Int(bs[])
end

export setoption!, gettype

"""
  Pass values to the Petsc function MatSetOption.  Note that the handful of
  options that can be passed here should not be confused with those for the
  global options database
"""
function setoption!(m::Mat, option::C.MatOption, val::Bool)
  chk(C.MatSetOption(m.p, option, PetscBool(val)))
  m
end

"""
  Get the format of the matrix.
"""
function gettype(a::PetscMat)
  sym_arr = Array{C.MatType}(undef,1)
  chk(C.MatGetType(a.p, sym_arr))
  return sym_arr[1]
end

function Base.resize!(a::Mat, m::Integer=C.PETSC_DECIDE, n::Integer=C.PETSC_DECIDE;
  mlocal::Integer=C.PETSC_DECIDE, nlocal::Integer=C.PETSC_DECIDE)
  if m == mlocal == C.PETSC_DECIDE
    throw(ArgumentError("either the global (m) or local (mlocal) #rows must be specified"))
  end
  if n == nlocal == C.PETSC_DECIDE
    throw(ArgumentError("either the global (n) or local (nlocal) #cols must be specified"))
  end
  chk(C.MatSetSizes(a.p, PetscInt(mlocal), PetscInt(nlocal), PetscInt(m), PetscInt(n)))
  a
end

"""
  Preallocates the sparsity pattern for (B)AIJ matrices.
"""
function setpreallocation!(a::Mat{T};
  nz::Integer=16, nnz::AbstractVector=PetscInt[],
  onz::Integer=0, onnz::AbstractVector=PetscInt[],
  bs::Integer=1) where {T}
  MType = gettype(a)
  if MType == C.MATSEQAIJ
    pnnz = if isempty(nnz)
      Ptr{PetscInt}(0)
    else
      if length(nnz) != size(a,1)
        throw(ArgumentError("length(nnz) must be # rows"))
      end
      isa(nnz,Vector{PetscInt}) ? nnz : PetscInt[ i for i in nnz ]
    end
    chk(C.MatSeqAIJSetPreallocation(a.p, PetscInt(nz), pnnz))
  elseif MType == C.MATMPIAIJ
    mlocal = sizelocal(a,1)
    pnnz = if isempty(nnz)
      Ptr{PetscInt}(0)
    else
      if length(nnz) != mlocal
        throw(ArgumentError("length(nnz) must be # local rows"))
      end
      isa(nnz,Vector{PetscInt}) ? nnz : PetscInt[ i for i in nnz ]
    end
    ponnz = if isempty(onnz)
      Ptr{PetscInt}(0)
    else
      if length(onnz) != mlocal
        throw(ArgumentError("length(onnz) must be # local rows"))
      end
      isa(onnz,Vector{PetscInt}) ? onnz : PetscInt[ i for i in onnz ]
    end
    chk(C.MatMPIAIJSetPreallocation(a.p, PetscInt(nz), pnnz, PetscInt(onz), ponnz))

  elseif MType == C.MATMPIBAIJ
    mlocal = sizelocal(a,1)
    pnnz = if isempty(nnz)
      Ptr{PetscInt}(0)
    else
      if length(nnz) != mlocal
        throw(ArgumentError("length(nnz) must be # local rows"))
      end
      isa(nnz,Vector{PetscInt}) ? nnz : PetscInt[ i for i in nnz ]
    end
    ponnz = if isempty(onnz)
      Ptr{PetscInt}(0)
    else
      if length(onnz) != mlocal
        throw(ArgumentError("length(onnz) must be # local rows"))
      end
      isa(onnz,Vector{PetscInt}) ? onnz : PetscInt[ i for i in onnz ]
    end
    chk(C.MatMPIBAIJSetPreallocation(a.p, bs, nz, pnnz, onz, ponnz))
  elseif MType == C.MATMPISBAIJ
    mlocal = sizelocal(a,1)
    pnnz = if isempty(nnz)
      Ptr{PetscInt}(0)
    else
      if length(nnz) != mlocal
        throw(ArgumentError("length(nnz) must be # local rows"))
      end
      isa(nnz,Vector{PetscInt}) ? nnz : PetscInt[ i for i in nnz ]
    end
    ponnz = if isempty(onnz)
      Ptr{PetscInt}(0)
    else
      if length(onnz) != mlocal
        throw(ArgumentError("length(onnz) must be # local rows"))
      end
      isa(onnz,Vector{PetscInt}) ? onnz : PetscInt[ i for i in onnz ]
    end
    chk(C.MatMPISBAIJSetPreallocation(a.p, bs, nz, pnnz, onz, ponnz))
  elseif MType == C.MATBLOCKMAT
    pnnz = if isempty(nnz)
      Ptr{PetscInt}(0)
    else
      if length(nnz) != size(a,1)
        throw(ArgumentError("length(nnz) must be # rows"))
      end
      isa(nnz,Vector{PetscInt}) ? nnz : PetscInt[ i for i in nnz ]
    end
    chk(C.MatBlockMatSetPreallocation(a.p, bs, nz, pnnz))
  elseif MType == C.MATSEQBAIJ
    pnnz = if isempty(nnz)
      Ptr{PetscInt}(0)
    else
      if length(nnz) != size(a,1)
        throw(ArgumentError("length(nnz) must be # rows"))
      end
      isa(nnz,Vector{PetscInt}) ? nnz : PetscInt[ i for i in nnz ]
    end
    chk(C.MatSeqBAIJSetPreallocation(a.p, bs, nz, pnnz))
  elseif MType == C.MATSEQSBAIJ
    pnnz = if isempty(nnz)
      Ptr{PetscInt}(0)
    else
      if length(nnz) != size(a,1)
        throw(ArgumentError("length(nnz) must be # rows"))
      end
      isa(nnz,Vector{PetscInt}) ? nnz : PetscInt[ i for i in nnz ]
    end
    chk(C.MatSeqSBAIJSetPreallocation(a.p, bs, nz, pnnz))
  else # TODO
    throw(ArgumentError("unsupported matrix type $T"))
  end
  a
end

"""
  Maps Matrix formats to the corresponding vector format
"""
# construct Vec for multiplication by a::Mat or transpose(a::Mat)
const mat2vec = Dict{C.MatType, C.MatType}( :mpiaij => :aij, :seqaij => :seq )

"""
  Construct a vector suitable for multiplying by the given matrix
"""
Vec{T2}(a::PetscMat{T2}, transposed=false) where {T2} =
  transposed ? Vec(T2, size(a,1), comm=comm(a), T=mat2vec[gettype(a)],
  mlocal=sizelocal(a,1)) :
  Vec(T2, size(a,2), comm=comm(a), T=mat2vec[gettype(a)],
  mlocal=sizelocal(a,2))

#############################################################################
Base.convert(::Type{C.Mat}, a::PetscMat) = a.p

export sizelocal, localranges, lengthlocal

"""
  Returns the global dimensions of the matrix
"""
function Base.size(a::PetscMat)
  m = Ref{PetscInt}()
  n = Ref{PetscInt}()
  chk(C.MatGetSize(a.p, m, n))
  (Int(m[]), Int(n[]))
end

"""
  Returns the local dimensions of the matrix
"""
function sizelocal(a::PetscMat)
  m = Ref{PetscInt}()
  n = Ref{PetscInt}()
  chk(C.MatGetLocalSize(a.p, m, n))
  (Int(m[]), Int(n[]))
end

"""
  This function returns two Range object corresponding the global indices
  of the rows and columns of the matrix A.

  The function has the same limiations as Petsc's MatGetOwnershipRange,
  in that it assumes the rows of the matrix are divided up contigously.
"""
function localranges(a::PetscMat)
  start_ref = Ref{PetscInt}()
  end_1_ref = Ref{PetscInt}()
  chk(C.MatGetOwnershipRange(a.p, start_ref, end_1_ref))
  start_idx = start_ref[] + 1  # convert to 1 based index
  end_idx = end_1_ref[]  # subtract 1 because petsc supplies 1 past the end,
                         # then add 1 to convert to 1 based index

  return start_idx:end_idx, 1:size(a, 2)
end

"""
  Similar to localpart, but returns the range of block indices
"""
function localpart_block(A::Mat)
  low = Ref{PetscInt}()
  high = Ref{PetscInt}()
  chk(C.MatGetOwnershipRange(A.p, low, high))
  bs = get_blocksize(A)
  low_b = div(low[], bs); high_b = div(high[]-1, bs)
  rows = (low_b+1):(high_b+1)
  cols = 1:div(size(A, 2), bs)

  return rows, cols
end

"""
  Constructs 2 index sets that map from the local row and columns to the
  global rows and columns
"""
function localIS(A::PetscMat{T}) where {T}

  rows, cols = localranges(A)
  rowis = IS(T, rows, comm=comm(A))
  colis = IS(T, cols, comm=comm(A))
  return rowis, colis
end

"""
  Like localIS, but returns a block index IS
"""
function localIS_block(A::Mat{T}) where {T}
  rows, cols = localpart_block(A)
  bs = get_blocksize(A)
  rowis = ISBlock(T, bs, rows, comm=comm(A))
  colis = ISBlock(T, bs, cols, comm=comm(A))
#  set_blocksize(rowis, get_blocksize(A))
  return rowis, colis
end

"""
  Construct ISLocalToGlobalMappings for the the rows and columns of the matrix
"""
function local_to_global_mapping(A::PetscMat)

  # localIS creates strided index sets, which require only constant
  # memory
  if get_blocksize(A) == 1
    rowis, colis = localIS(A)
  else
    rowis, colis = localIS_block(A)
  end
  row_ltog = ISLocalToGlobalMapping(rowis)
  col_ltog = ISLocalToGlobalMapping(colis)
  return row_ltog, col_ltog
end

# need a better name
"""
  Registers the ISLocalToGlobalMappings with the matrix
"""
function set_local_to_global_mapping(A::PetscMat{T}, rmap::ISLocalToGlobalMapping{T}, cmap::ISLocalToGlobalMapping{T}) where {T}

  chk(C.MatSetLocalToGlobalMapping(A.p, rmap.p, cmap.p))
end

"""
  Check if the local to global mappings have been registered
"""
function has_local_to_global_mapping(A::PetscMat{T}) where {T}

  rmap_ref = Ref{C.ISLocalToGlobalMapping{T}}()
  cmap_ref = Ref{C.ISLocalToGlobalMapping{T}}()
  chk(C.MatGetLocalToGlobalMapping(A.p, rmap_ref, cmap_ref))

  rmap = rmap_ref[]
  cmap = cmap_ref[]

  return rmap.pobj != C_NULL && cmap.pobj != C_NULL
end

"""
  prod(sizelocal))
"""
lengthlocal(a::PetscMat) = prod(sizelocal(a))

# this causes the assembly state of the underlying petsc matrix to be copied
function Base.similar(a::PetscMat{T}) where {T}
  p = Ref{C.Mat{T}}()
  chk(C.MatDuplicate(a.p, C.MAT_DO_NOT_COPY_VALUES, p))
  Mat{T}(p[])
end

Base.similar(a::PetscMat{T}, ::Type{T}) where {T} = similar(a)

# TODO: make T2 -> Type{T2}
function Base.similar(a::PetscMat{T}, T2::Type) where {T}
  MType = gettype(a)
  return Mat(T2, size(a)..., comm=comm(a), mtype=MType)
end

function Base.similar(a::PetscMat{T}, T2::Type, m::Integer, n::Integer) where {T}
  MType = gettype(a)
  (m,n) == size(a) && T2==T ? similar(a) : Mat(T2, m,n, comm=comm(a), mtype=MType)
end

Base.similar(a::PetscMat{T}, m::Integer, n::Integer) where {T} = similar(a, T, m, n)
#Base.similar(a::PetscMat, T::Type, d::Dims{2}) = similar(a, T, d...)
Base.similar(a::PetscMat, T::Type, d::Dims{2}) = similar(a, T, d...)
Base.similar(a::PetscMat{T}, d::Dims{2}) where {T}  = similar(a, T, d)

function Base.copy(a::PetscMat{T}) where {T}
  p = Ref{C.Mat{T}}()
  chk(C.MatDuplicate(a.p, C.MAT_COPY_VALUES, p))
  A=Mat{T}(p[])
  A.verify_assembled = a.verify_assembled
  A.data             = a.data
  A.insertmode       = a.insertmode
  A
end

"""
  Get the MatInfo struct for the matrix
"""
function getinfo(m::Mat, infotype::Integer=C.MAT_GLOBAL_SUM)
  info = Ref{C.MatInfo}()
  chk(C.MatGetInfo(m.p, C.MatInfoType(infotype), info))
  info[]
end

"""
  Number of non-zero entries that have been assigned to
"""
SparseArrays.nnz(m::Mat) = Int(getinfo(m).nz_used)

#############################################################################

# for efficient matrix assembly, put all calls to A[...] = ... inside
# assemble(A) do ... end
"""
  Start assembling the matrix (the implmentations probably post
  non-blocking sends and received)

"""
function AssemblyBegin(x::PetscMat, t::C.MatAssemblyType=C.MAT_FLUSH_ASSEMBLY)
  chk(C.MatAssemblyBegin(x.p, t))
end

"""
  Finish assembling the matrix
"""
function AssemblyEnd(x::PetscMat, t::C.MatAssemblyType=C.MAT_FLUSH_ASSEMBLY)
  chk(C.MatAssemblyEnd(x.p, t))
end

"""
  Check if the matrix is assembled or not
"""
function isassembled(p::C.Mat)
  b = Ref{PetscBool}()
  chk(C.MatAssembled(p, b))
  return b[] != 0
end

"""
  Check if the matrix is assembled.  Whether all processes assembly state
  is checked or only the local process is determined by `x.verify_assembled`.

  `local_only` forces only the local process to be checked, regardless of
  `x.verify_assembled`.
"""
function isassembled(x::PetscMat; local_only=false)
  val = isassembled(x.p)
  if !local_only
    if x.verify_assembled
      val = MPI.Allreduce(Int8(val), MPI.LAND, comm(x))
    end
  end

  return Bool(val)
end

"""
  This function provides a mechanism for efficiently inserting values into
  and then assembling Petsc matrices and vectors.  The function f must be a
  zero argument function.

  This function can be used with the do block syntax.
"""
function assemble(f::Function, x::Union{Vec,PetscMat},
  insertmode=x.insertmode,
  assemblytype::C.MatAssemblyType=C.MAT_FINAL_ASSEMBLY)
  if x.insertmode != insertmode && !isassembled(x, local_only=true)
    error("nested assemble with different insertmodes not allowed")
  end
  old_insertmode = x.insertmode
  try  # what is the purpose of the try - finally?
    x.insertmode = insertmode
    result = f()
    return result
  finally
    AssemblyBegin(x, assemblytype)
    yield() # do async computations while messages are in transit
    AssemblyEnd(x, assemblytype)
    x.insertmode = old_insertmode
  end
end

# force assembly even if it might not be necessary
"""
  Assemble the Petsc object
"""
function assemble(x::AbstractArray, t::C.MatAssemblyType=C.MAT_FINAL_ASSEMBLY)
  AssemblyBegin(x, t)
  AssemblyEnd(x,t)
end

# in ksp solve we need to finalize assembly from raw pointer:
"""
  Low level assemble function
"""
function assemble(p::C.Mat, t::C.MatAssemblyType=C.MAT_FINAL_ASSEMBLY)
  if !isassembled(p)
    chk(C.MatAssemblyBegin(p, t))
    chk(C.MatAssemblyEnd(p, t))
  end
  return nothing
end

# intermediate assembly, before it is finally compressed for use
"""
  Perform a flush assembly (take stashed values and put them into the matrix,
  but don't squeeze out any preallocated space that has not been used yet
"""
iassemble(x::Union{Vec,PetscMat}) = assemble(() -> nothing, x, x.insertmode, C.MAT_FLUSH_ASSEMBLY)
iassemble(f::Function, x::PetscMat, insertmode=x.insertmode) =
  assemble(f, x, insertmode, C.MAT_FLUSH_ASSEMBLY)

iassemble(x::AbstractArray) = assemble(x, C.MAT_FLUSH_ASSEMBLY)

# like x[i,j] = v, but requires i,j to be 0-based indices for Petsc
function setindex0!(x::Mat{T}, v::Array{T},
                       i::Array{PetscInt}, j::Array{PetscInt}) where {T}
  ni = length(i)
  nj = length(j)
  if length(v) != ni*nj
    throw(ArgumentError("length(values) != length(indices)"))
  end
  chk(C.MatSetValues(x.p, PetscInt(ni), i, PetscInt(nj), j, v, x.insertmode))
  x
end

# like setindex0 above, but for submatrices.  the indices i and j must
# be *local* indices
function setindex0!(x::SubMat{T}, v::Array{T}, i::Array{PetscInt}, j::Array{PetscInt}) where {T}
  ni = length(i)
  nj = length(j)
  if length(v) != ni*nj
    throw(ArgumentError("length(values) != length(indices)"))
  end
  chk(C.MatSetValuesLocal(x.p, PetscInt(ni), i, PetscInt(nj), j, v, x.insertmode))
  x
end

import Base: setindex!

function setindex!(x::PetscMat{T}, v::Number, i::Integer, j::Integer) where {T}
  # can't call MatSetValue since that is a static inline function
  setindex0!(x, T[ v ],
  PetscInt[ i - 1 ],
  PetscInt[ j - 1 ])
  v
end
function setindex!(x::PetscMat{T3}, v::Array{T3},
  I::AbstractArray{T1},
  J::AbstractArray{T2}) where {T3, T1<:Integer, T2<:Integer}
  I0 = PetscInt[ i-1 for i in I ]
  J0 = PetscInt[ j-1 for j in J ]
  setindex0!(x, v, I0, J0)
end
function setindex!(x::PetscMat{T2}, v::Array{T2},
  i::Integer, J::AbstractArray{T}) where {T2, T<:Integer}
  I0 = PetscInt[ i-1 ]
  J0 = PetscInt[ j-1 for j in J ]
  setindex0!(x, v, I0, J0)
end
function setindex!(x::PetscMat{T2}, v::Array{T2},
  I::AbstractArray{T}, j::Integer) where {T2, T<:Real}
  I0 = PetscInt[ i-1 for i in I ]
  J0 = PetscInt[ j-1 ]
  setindex0!(x, v, I0, J0)
end

setindex!(x::PetscMat, v::Number,
I::AbstractArray{T1},
J::AbstractArray{T2}) where {T1<:Integer, T2<:Integer} = iassemble(x) do
  for i in I
    for j in J
      x[i,j] = v
    end
  end
  x
end
setindex!(x::PetscMat, v::Number,
i::Integer, J::AbstractArray{T}) where {T<:Integer} = iassemble(x) do
  for j in J
    x[i,j] = v
  end
  x
end
setindex!(x::PetscMat, v::Number,
I::AbstractArray{T}, j::Integer) where {T<:Integer} = iassemble(x) do
  for i in I
    x[i,j] = v
  end
  x
end

function setindex!(x::PetscMat, v::AbstractArray{T0},
  I::AbstractArray{T1},
  J::AbstractArray{T2}) where {T0<:Number, T1<:Integer, T2<:Integer}
  if length(v) != length(I)*length(J)
    throw(ArgumentError("length(values) != length(indices)"))
  end
  v0 = eltype(Mat)[ z for z in v ]
  setindex!(x, v0, I, J)
end

# fill! is not a very sensible function to call for sparse matrices,
# but we might as well have it, especially for the v=0 case
"""
  Fill the matrix with the specified values.

  Currently, this function either destroys the sparsity pattern or
  gives an error, unless v = 0, in which case it zeros out the non-zero
  entries without changing the sparsity pattern
"""
function Base.fill!(x::PetscMat, v::Number)
  if v == 0
    chk(C.MatZeroEntries(x.p))
  else
    # FIXME: don't write to non-allocated entries
    iassemble(x) do
      rows, cols = localranges(x)
      for i in rows
        for j in cols
          x[i,j] = v
        end
      end
    end
  end
  return x
end

#############################################################################
import Base.getindex

# like getindex but for 0-based indices i and j
function getindex0(x::Mat{T}, i::Vector{PetscInt}, j::Vector{PetscInt}) where {T}
  ni = length(i)
  nj = length(j)
  v = Array{T}(undef, nj, ni) # row-major!
  chk(C.MatGetValues(x.p, PetscInt(ni), i, PetscInt(nj), j, v))
  ni <= 1 || nj <= 1 ? reshape(v, ni, nj) : transpose(v)
end

getindex(a::PetscMat, i0::Integer, i1::Integer) =
  getindex0(a, PetscInt[ i0-1 ], PetscInt[ i1-1 ])[1]

getindex(a::PetscMat, I0::AbstractArray{T0}, I1::AbstractArray{T1}) where {T0<:Integer,T1<:Integer} =
  getindex0(a, PetscInt[ i0-1 for i0 in I0 ], PetscInt[ i1-1 for i1 in I1 ])

getindex(a::PetscMat, I0::AbstractArray{T0}, i1::Integer) where {T0<:Integer} =
  reshape(getindex0(a, PetscInt[ i0-1 for i0 in I0 ], PetscInt[ i1-1 ]), length(I0))

getindex(a::PetscMat, i0::Integer, I1::AbstractArray{T1}) where {T1<:Integer} =
  getindex0(a, PetscInt[ i0-1 ], PetscInt[ i1-1 for i1 in I1 ])

#############################################################################
# zero based indexing
#TODO: in 0.5, use boundscheck macro to verify stride=1

# global, non-block
function set_values!(x::Mat{T}, idxm::StridedVecOrMat{PetscInt}, idxn::StridedVecOrMat{PetscInt}, v::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar}
  # v should be m x n
  chk(C.MatSetValues(x.p, PetscInt(length(idxm)), idxm, PetscInt(length(idxn)), idxn, v, o))
end

function set_values!(x::Mat{T}, idxm::StridedVecOrMat{I1}, idxn::StridedVecOrMat{I2}, v::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar, I1 <: Integer, I2 <: Integer}

  _idxm = PetscInt[ i for i in idxm]
  _idxn = PetscInt[ i for i in idxn]
  set_values!(x, _idxm, _idxn, v, o)
end

function set_values!(x::Matrix, idxm::AbstractArray, idxn::AbstractArray, v::AbstractArray, o::C.InsertMode=C.INSERT_VALUES)
  if o == C.INSERT_VALUES
    for col = 1:length(idxn)
      colidx = idxn[col]
      for row = 1:length(idxm)
        rowidx = idxm[row]
        x[rowidx+1, colidx+1] = v[row, col]
      end
    end
  else
    for col = 1:length(idxn)
      colidx = idxn[col]
      for row = 1:length(idxm)
        rowidx = idxm[row]
        x[rowidx+1, colidx+1] += v[row, col]
      end
    end
  end

end

# global, block
function set_values_blocked!(x::Mat{T}, idxm::StridedVecOrMat{PetscInt}, idxn::StridedVecOrMat{PetscInt}, v::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar}

  # vals should be m*bs x n*bs
  chk(C.MatSetValuesBlocked(x.p, PetscInt(length(idxm)), idxm, PetscInt(length(idxn)), idxn, v, o))
end

function set_values_blocked!(x::Mat{T}, idxm::StridedVecOrMat{I1}, idxn::StridedVecOrMat{I2}, v::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar, I1 <: Integer, I2 <: Integer}

  _idxm = PetscInt[ i for i in idxm]
  _idxn = PetscInt[ i for i in idxn]
  set_values_blocked!(x, _idxm, _idxn, v, o)
end

# local, non-block
function set_values_local!(x::Mat{T}, idxm::StridedVecOrMat{PetscInt}, idxn::StridedVecOrMat{PetscInt}, v::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar}

  chk(C.MatSetValuesLocal(x.p, PetscInt(length(idxm)), idxm, PetscInt(length(idxn)), idxn, v, o))
end

function set_values_local!(x::Mat{T}, idxm::StridedVecOrMat{I1}, idxn::StridedVecOrMat{I2}, v::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar, I1 <: Integer, I2 <: Integer}

  _idxm = PetscInt[ i for i in idxm]
  _idxn = PetscInt[ i for i in idxn]
  set_values_local!(x, _idxm, _idxn, v, o)
end

function set_values_local!(x::Matrix, idxm::AbstractArray, idxn::AbstractArray, v::AbstractArray, o::C.InsertMode=C.INSERT_VALUES)

  set_values!(x, idxm, idxn, v, o)
end

# local, block
function set_values_blocked_local!(x::Mat{T}, idxm::StridedVecOrMat{PetscInt}, idxn::StridedVecOrMat{PetscInt}, v::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar}

  chk(C.MatSetValuesBlockedLocal(x.p, PetscInt(length(idxm)), idxm, PetscInt(length(idxn)), idxn, v, o))
end

function set_values_blocked_local(x::Mat{T}, idxm::StridedVecOrMat{I1}, idxn::StridedVecOrMat{I2}, v::StridedVecOrMat{T}, o::C.InsertMode=x.insertmode) where {T <: Scalar, I1 <: Integer, I2 <: Integer}

  _idxm = PetscInt[ i for i in idxm]
  _idxn = PetscInt[ i for i in idxn]
  set_values_blocked_local!(x, _idxm, _idxn, v, o)
end

###############################################################################
# transposition etc.

export MatTranspose

"""
  Create a dense Julia matrix for a Petsc sparse matrix.  This only works
  for SEQ matrices
"""
function Array(a::PetscMat)
  Array{eltype(PetscMat),2}(a)
end

function Array{T,2}(a::PetscMat) where {T}
  m,n = size(a)
  convert(Array{T,2},a[1:m, 1:n])
end

#=
"""
 create a new matrix wrapping around A for matrix-vector multiplication
 but which does not actually require new storage
"""
=#
for (f,pf) in ((:MatTranspose,:MatCreateTranspose), # acts like A.'
  (:MatNormal, :MatCreateNormal))      # acts like A'*A
  pfe = Expr(:quote, pf)
  @eval function $f(a::PetscMat{T}) where {T}
    p = Ref{C.Mat{T}}()
    chk(C.$pf(a.p, p))
    Mat(p[], a)
  end
end

#for (f,pf) in ((:transpose,:MatTranspose),(:ctranspose,:MatHermitianTranspose))
for (f,pf) in ((:transpose,:MatTranspose),)
  fb = Symbol(string(f,"!"))
  pfe = Expr(:quote, pf)
  @eval begin
    function LinearAlgebra.$fb(a::PetscMat)
      pa = [a.p]
      chk(C.$pf(a.p, C.MAT_REUSE_MATRIX, pa))
      a
    end

    function LinearAlgebra.$f(a::PetscMat{T}) where {T}
      p = Ref{C.Mat{T}}()
      chk(C.$pf(a.p, C.MAT_INITIAL_MATRIX, p))
      Mat(p[])
    end
  end
end

function Base.conj!(a::PetscMat)
  chk(C.MatConjugate(a.p))
  a
end
Base.conj(a::PetscMat) = conj!(copy(a))

#############################################################################
# simple math operations

# skip for now
#=
function chop!(x::Mat, tol::Real)
chk(ccall((:MatChop, petsc), PetscErrorCode, (pMat, PetscReal), x, tol))
x
end
=#

#import Base: .*, ./, .\, *, +, -, ==
#import LinearAlgebra: At_mul_B, At_mul_B!, Ac_mul_B, Ac_mul_B!, A_mul_Bt, A_mul_Bt!

function LinearAlgebra.tr(A::PetscMat{T}) where {T}
  t = Ref{T}()
  chk(C.MatGetTrace(A.p,t))
  return t[]
end

function Base.real(A::PetscMat{T}) where {T<:Complex}
  N = copy(A)
  chk(C.MatRealPart(N.p))
  return N
end
Base.real(A::PetscMat{T}) where {T<:Real} = A

function Base.imag(A::PetscMat{T}) where {T<:Complex}
  N = copy(A)
  chk(C.MatImaginaryPart(N.p))
  return N
end

#function Base.LinAlg.ishermitian(A::PetscMat{T}, tol::Real=eps(real(float(one(T))))) where {T}
#  bool_arr = Ref{PetscBool}()
#  chk(C.MatIsHermitian(A.p, tol, bool_arr))
#  return bool_arr[] != 0
#end

#function Base.LinAlg.issymmetric(A::PetscMat{T}, tol::Real=eps(real(float(one(T))))) where {T}
#  bool_arr = Ref{PetscBool}()
#  chk(C.MatIsSymmetric(A.p, tol, bool_arr))
#  return bool_arr[] != 0
#end

#currently ONLY gets the main diagonal
#function Base.diag(A::PetscMat{T},vtype::C.VecType=C.VECSEQ) where {T}
#  m = size(A, 1)
#  b = Vec(T, m, vtype=vtype, comm=comm(A), mlocal=sizelocal(A,1))
#  chk(C.MatGetDiagonal(A.p,b.p))
#  return b
#end

function (*)(A::PetscMat{T}, x::Vec{T}) where {T}
  m = size(A, 1)
  MType = gettype(x)
  b = Vec(T, m, vtype=MType, comm=comm(A), mlocal=sizelocal(A,1))
  chk(C.MatMult(A.p, x.p, b.p))
  return b
end

function (*)(A::PetscMat{T}, x::Vec{T}, b::Vec{T}) where {T}
  chk(C.MatMult(A.p, x.p, b.p))
  return b
end

#function (.*)(A::PetscMat{T}, x::Number) where {T}
#  Y = copy(A)
#  chk(C.MatScale(Y.p, T(x)))
#  return Y
#end
#(.*)(x::Number, A::PetscMat{T}) where {T} = A .* x
#(./)(A::PetscMat{T}, x::Number) where {T} = A .* inv(x)
#(.\)(x::Number, A::PetscMat{T}) where {T} = A .* inv(x)

function (*)(A::PetscMat{T}, B::PetscMat{T}) where {T}
  p = Ptr{Float64}(0)
  p_arr = Ref{C.Mat{T}}()
  chk(C.MatMatMult(A.p, B.p, C.MAT_INITIAL_MATRIX, real(T)(C.PETSC_DEFAULT), p_arr))
  new_mat = Mat(p_arr[])
  return new_mat
end

#these two only work for SEQAIJ
function At_mul_B(A::PetscMat{T}, B::PetscMat{T}) where {T}
  p = Ptr{Float64}(0)
  p_arr = Ref{C.Mat{T}}()
  chk(C.MatTransposeMatMult(A.p, B.p, C.MAT_INITIAL_MATRIX, real(T)(C.PETSC_DEFAULT), p_arr))
  new_mat = Mat(p_arr[])
  return new_mat
end

function A_mul_Bt(A::PetscMat{T}, B::PetscMat{T}) where {T}
  p = Ptr{Float64}(0)
  p_arr = Ref{C.Mat{T}}()
  chk(C.MatMatTransposeMult(A.p, B.p, C.MAT_INITIAL_MATRIX, real(T)(C.PETSC_DEFAULT), p_arr))
  new_mat = Mat(p_arr[])
  return new_mat
end

function At_mul_B!(A::PetscMat{T}, x::Vec{T}, y::Vec{T}) where {T}
  chk(C.MatMultTranspose(A.p, x.p, y.p))
  return y
end

function At_mul_B(A::PetscMat{T}, x::Vec{T}) where {T}
  m = size(A, 1)
  MType = gettype(x)
  b = Vec(T, m, vtype=MType, comm=comm(A), mlocal=sizelocal(A,1))
  chk(C.MatMultTranspose(A.p, x.p, b.p))
  return b
end

function Ac_mul_B!(A::PetscMat{T}, x::Vec{T}, y::Vec{T}) where {T}
  chk(C.MatMultHermitianTranspose(A.p, x.p, y.p))
  return y
end

function Ac_mul_B(A::PetscMat{T}, x::Vec{T}) where {T}
  m = size(A, 1)
  MType = gettype(x)
  b = Vec(T, m, vtype=MType, comm=comm(A), mlocal=sizelocal(A,1))
  chk(C.MatMultHermitianTranspose(A.p, x.p, b.p))
  return b
end

for (f,s) in ((:+,1), (:-,-1))
  @eval function ($f)(A::PetscMat{T}, B::PetscMat{T}) where {T}
    Y = copy(A)
    chk(C.MatAXPY(Y.p, T($s), B.p, C.DIFFERENT_NONZERO_PATTERN))
    return Y
  end
end

(-)(A::PetscMat{T}) where {T} = A * (-1)

# there don't appear to be PETSc functions for pointwise
# operations on matrices

function (==)(A::PetscMat{T}, b::PetscMat{T}) where {T}
  bool_arr = Ref{PetscBool}()
  chk(C.MatEqual(A.p, b.p, bool_arr))
  return bool_arr[] != 0
end
#=
# needed for disambiguation
function (==)(A::PetscMat{T}, B::PetscMat{T})
  bool_arr = Ref{PetscBool}()
  chk(C.MatEqual(A.p, b.p, bool_arr))
  return bool_arr[] != 0
end
=#
"""
  Equality test for AbstractArray and PetscMat, SEQ only
"""
function (==)(A::PetscMat{T}, B::AbstractArray) where {T}
  if MPI.Comm_size(comm(A)) != 1
    throw(ArgumentError("Mat must reside on a single MPI rank"))
  end
  if size(A) != size(B)
    throw(ArgumentError("Matrices must be same size"))
  end

  m, n = size(B)

  isame = true
  for i=1:n
    for j=1:m
      isame = isame && A[i, j] == B[i,j]
    end
  end

  return isame
end

###############################################################################
# MatRow: accessing the structure of each row of a matrix
"""
  Object that enables access to a single row of a sparse matrix.

  Users *must* call restore when done with a MatRow, before attempting to
  create another one.
"""
struct MatRow{T}
  mat::C.Mat{T}  # the matrix to which the rows belong
  row::Int
  ncols::Int
  cols_ptr::Ptr{PetscInt}
  vals_ptr::Ptr{T}
end

"""
  Preferred constructor for a MatRow
"""
function MatRow(A::Mat{T}, row::Integer) where {T}

  ref_ncols = Ref{PetscInt}()
  ref_cols = Ref{Ptr{PetscInt}}()
  ref_vals = Ref{Ptr{T}}()
  chk(C.MatGetRow(A.p, PetscInt(row-1), ref_ncols, ref_cols, ref_vals))
  return MatRow(A.p, row, Int(ref_ncols[]), ref_cols[], ref_vals[])
end

"""
  Count the number of non-zeros in a row of the matrix.  This temporarily
  creates a MatRow object, so it cannot be used if one already exists
"""
function count_row_nz(A::Mat{T}, row::Integer) where {T}
  ref_ncols = Ref{PetscInt}()
  ref_cols = Ref{Ptr{PetscInt}}(C_NULL)
  ref_vals = Ref{Ptr{T}}(C_NULL)
  chk(C.MatGetRow(A.p, PetscInt(row-1), ref_ncols, ref_cols, ref_vals))
  ncols = ref_ncols[]
  chk(C.MatRestoreRow(A.p, PetscInt(row-1), ref_ncols, ref_cols, ref_vals))

  return ncols
end
function restore(row::MatRow{T}) where {T}
#  if !PetscFinalized(T) && !isfinalized(row.mat)
    C.MatRestoreRow(row.mat, PetscInt(row.row-1), Ref(PetscInt(row.ncols)), Ref{Ptr{PetscInt}}(C_NULL), Ref{Ptr{T}}(C_NULL))

#    return nothing  # return type stability
#  end
end

### indexing on a MatRow ###
import Base: length, size
length(A::MatRow) = A.ncols
size(A::MatRow) = (A.ncols,)
getcol(A::MatRow, i) = unsafe_load(A.cols_ptr, i) + 1
getval(A::MatRow, i) = unsafe_load(A.vals_ptr, i)

import Base.kron
"""
  Kronecker product for SEQ matrices only.  The output is a non-block matrix
  even if the inputs are block
"""
function kron(A::Mat{T}, B::Mat{T}) where {T}
  if (A.p == B.p)
    throw(ArgumentError("A and B cannot be same matrix"))
  end

  if MPI.Comm_size(comm(A)) != 1 || MPI.Comm_size(comm(B)) != 1
    throw(ArgumentError("A and B must reside on a single MPI rank"))
  end

  Am, An = size(A)
  Bm, Bn = size(B)
  # step 1: figure out size, sparsity pattern of result
  A_nz = zeros(Int, Am)
  B_nz = zeros(Int, Bm)
  for i=1:Am
    A_nz[i] = count_row_nz(A, i)
  end
  for i=1:Bm
    B_nz[i] = count_row_nz(B, i)
  end

  Dm = Am*Bm
  Dn = An*Bn
  D_nz = getDnzs(A_nz, B_nz)
  # create matrix
  # can't use C becaue that is the module name
  D = Mat(T, Dm, Dn, nnz=D_nz, mtype=C.MATSEQAIJ)

  # step 2: populate C, one row at a time
  max_entries = maximum(D_nz)
  D_colidx = zeros(PetscInt, max_entries)
  D_vals = zeros(T, max_entries)
  D_rowidx = zeros(PetscInt, 1)

  for i=1:Am
    rowA = MatRow(A, i)
    for j=1:Bm
      rowB = MatRow(B, j)
      rowidx = (i-1)*Bm + j

      pos = 1  # position in C_vals, C_colidx
      for k=1:length(rowA)
        Aval = getval(rowA, k)
        Aidx = getcol(rowA, k)
        for p=1:length(rowB)
          Bval = getval(rowB, p)
          Bidx = getcol(rowB, p)
          D_colidx[pos] = (Aidx - 1)*Bn + Bidx - 1
          D_vals[pos] = Aval*Bval
          pos += 1
        end
      end

      D_rowidx[1] = rowidx - 1
      idx_extract = view(D_colidx, 1:(pos-1))
      vals_extract = view(D_vals, 1:(pos-1))
      set_values!(D, D_rowidx, idx_extract, vals_extract)
      restore(rowB)
    end
    restore(rowA)
  end

  return D
end

function getDnzs(A_nz, B_nz)
# this function computes the number of non-zeros in D = kron(A, B), where A_nz and
# B_nz are the number of non-zeros in each row of A and B, respectively

  Am = length(A_nz)
  Bm = length(B_nz)

  Dm = length(A_nz)*length(B_nz)
  D_nz = zeros(PetscInt, Dm)
  for i=1:Am
    A_nz_i = A_nz[i]
    for j=1:Bm
      pos = (i-1)*Bm + j
      D_nz[pos] = A_nz_i*B_nz[j]
    end
  end

  return D_nz
end

function getmax_nz_col(A::SparseMatrixCSC)
  Am, An = size(A)
  A_maxnz = 0
  for i=1:An
    val = A.colptr[i+1] - A.colptr[i]
    if val > A_maxnz
      A_maxnz = val
    end
  end

  return A_maxnz
end

"""
  Kronecker product of A and B where the result is a Petsc Mat
"""
function PetscKron(A::SparseMatrixCSC{T}, B::SparseMatrixCSC{T}) where {T <: Scalar}

  Am = size(A, 1); An = size(A, 2)
  Bm = size(B, 1); Bn = size(B, 2)

  A_nz = zeros(Int, Am)
  B_nz = zeros(Int, Bm)

  for i in A.rowval
    A_nz[i] += 1
  end
  for i in B.rowval
    B_nz[i] += 1
  end

  Dm = Am*Bm
  Dn = An*Bn
  D_nz = getDnzs(A_nz, B_nz)
  # create matrix
  # can't use C becaue that is the module name
  D = Mat(T, Dm, Dn, nnz=D_nz, mtype=C.MATSEQAIJ)
  # because the Mat constructor is type unstable, use a function barrier
  PetscKron(A, B, D)

  return D
end

"""
  Kronecker product of A and B, storing the result in D.  D should already be
  pre-allocated with the right sparsity pattern
"""
@noinline function PetscKron(A::SparseMatrixCSC{T}, B::SparseMatrixCSC{T}, D::Mat{T}) where {T <: Scalar}

  Am = size(A, 1); An = size(A, 2)
  Bm = size(B, 1); Bn = size(B, 2)

  # now figure out the maximum number of non-zeros in any column of D
  A_maxnz = getmax_nz_col(A)
  B_maxnz = getmax_nz_col(B)
  max_entries = A_maxnz * B_maxnz

  D_rowidx = zeros(PetscInt, max_entries)
  D_vals = zeros(T, max_entries)
  D_colidx = zeros(PetscInt, 1)

  # loop over columns of A and B
  for i=1:An
    col_start = A.colptr[i]
    col_end = A.colptr[i+1] - 1
    A_rowidx = view(A.rowval, col_start:col_end)
    A_vals = view(A.nzval, col_start:col_end)
    for j=1:Bn
      col_start = B.colptr[j]
      col_end = B.colptr[j+1] - 1
      B_rowidx = view(B.rowval, col_start:col_end)
      B_vals = view(B.nzval, col_start:col_end)

      pos = 1
      for k=1:length(A_vals)
        Aval = A_vals[k]
        Aidx = A_rowidx[k]
        for p=1:length(B_vals)
          Bval = B_vals[p]
          Bidx = B_rowidx[p]

          D_vals[pos] = Aval*Bval
          D_rowidx[pos] = (Aidx-1)*Bm + Bidx - 1
          pos += 1
        end
      end

      D_colidx[1] = (i-1)*Bn + j - 1
      rowidx_extract = view(D_rowidx, 1:(pos-1))
      vals_extract = view(D_vals, 1:(pos-1))

#      chk(C.MatSetValues(D.p, length(rowidx_extract), rowidx_extract, length(D_colidx), D_colidx, vals_extract, C.INSERT_VALUES))
      set_values!(D, rowidx_extract, D_colidx, vals_extract)
    end
#    assemble(D, C.MAT_FLUSH_ASSEMBLY)
  end

  return D
end


text:
using ItemGraphs
using Test

abstract type Item end

struct A <: Item end
struct B <: Item end
struct C <: Item end
struct D <: Item end

@testset "ItemGraphs" begin
    g = ItemGraph{Item}(lazy=true)
    add_edge!(g, A(), B())
    add_edge!(g, B(), C())
    add_edge!(g, B(), D())

    @test ItemGraphs.getid(g, A()) == 1
    @test ItemGraphs.getid(g, B()) == 2
    @test ItemGraphs.getid(g, C()) == 3
    @test ItemGraphs.getid(g, D()) == 4
    @test ItemGraphs.getitem(g, 1) == A()
    @test ItemGraphs.getitem(g, 2) == B()
    @test ItemGraphs.getitem(g, 3) == C()
    @test ItemGraphs.getitem(g, 4) == D()

    @test items(g, A(), C()) == [A(), B(), C()]
    @test items(g, A(), D()) == [A(), B(), D()]

    g = ItemGraph{Int}()
    add_edge!(g, 101, 202)
    add_edge!(g, 202, 303)
    add_edge!(g, 202, 404)

    @test ItemGraphs.getid(g, 101) == 1
    @test ItemGraphs.getid(g, 202) == 2
    @test ItemGraphs.getid(g, 303) == 3
    @test ItemGraphs.getid(g, 404) == 4
    @test ItemGraphs.getitem(g, 1) == 101
    @test ItemGraphs.getitem(g, 2) == 202
    @test ItemGraphs.getitem(g, 3) == 303
    @test ItemGraphs.getitem(g, 4) == 404

    @test items(g, 101, 303) == [101, 202, 303]
    @test items(g, 101, 404) == [101, 202, 404]
    @test edgeitems(g, 101, 404) == [0.0, 0.0]

    @test_throws ItemGraphException items(g, 101, 102)
    @test_throws ItemGraphException items(g, 102, 101)

    add_edge!(g, 505, 606)
    @test_throws ItemGraphException items(g, 101, 505)

    g = ItemGraph{Int,Int}()
    add_edge!(g, 101, 202, 102)
    add_edge!(g, 202, 303, 203)
    add_edge!(g, 202, 404, 204)
    @test edgeitem(g, 101, 202) == 102
    @test edgeitems(g, 101, 303) == [102, 203]
    @test edgeitems(g, 101, 404) == [102, 204]
end


text:
const ArbMatIdx0 = zero(Int32)

mutable struct ArbComplexMatrix{P} <: AbstractArbMatrix{P, ArbComplex}
    eachcell::Ptr{ArbComplex{P}}
    rowcount::Int
    colcount::Int
    eachrow::Ptr{Ptr{ArbComplex{P}}}

    function ArbComplexMatrix{P}(rowcount::Int, colcount::Int) where {P}
        z = new{P}() # z = new{P}(Ptr{ArbComplex{P}}(0), 0, 0, Ptr{Ptr{ArbComplex{P}}}(0))
        acb_mat_init(z, rowcount, colcount)
        finalizer(acb_mat_clear, z)
        return z
    end
end

#ArbComplexMatrix(x::ArbComplexMatrix) = x
ArbComplexMatrix(x::ArbComplexMatrix{P}) where {P} = x
ArbComplexMatrix{P}(x::ArbComplexMatrix{P}) where {P} = x

@inline function acb_mat_clear(x::ArbComplexMatrix{P}) where {P}
    ccall(@libarb(acb_mat_clear), Cvoid, (Ref{ArbComplexMatrix}, ), x)
end

@inline function acb_mat_init(x::ArbComplexMatrix{P}, rowcount::I, colcount::I) where {P, I<:Signed}
    ccall(@libarb(acb_mat_init), Cvoid, (Ref{ArbComplexMatrix}, Cint, Cint),
                                         x, rowcount, colcount)
end

@inline function ArbComplexMatrix(nrows::Int, ncols::Int)
    P = workingprecision(ArbComplex)
	return ArbComplexMatrix{P}(nrows, ncols)
end

function ArbComplexMatrix(fpm::Array{ArbComplex, 2})
    P = workingprecision(ArbComplex)
    nrows, ncols = size(fpm)
    arm = ArbComplexMatrix{P}(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm[r,c] = fpm[r,c]
	    end
	end
	return arm
end

function ArbComplexMatrix(fpm::Array{ArbComplex{P},2}) where {P}
    nrows, ncols = size(fpm)
    arm = ArbComplexMatrix{P}(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm[r,c] = fpm[r,c]
	    end
	end
	return arm
end

function ArbComplexMatrix{P}(fpm::Array{ArbComplex{P},2}) where {P}
    nrows, ncols = size(fpm)
    arm = ArbComplexMatrix{P}(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm[r,c] = fpm[r,c]
	    end
	end
	return arm
end

function ArbComplexMatrix{P}(fpm::Array{ArbComplex,2}) where {P}
    nrows, ncols = size(fpm)
    arm = ArbComplexMatrix{P}(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm[r,c] = ArbComplex{P}(fpm[r,c])
	    end
	end
	return arm
end

function ArbComplexMatrix(fpm::Array{ArbFloat{P},2}) where {P}
    nrows, ncols = size(fpm)
    arm = ArbComplexMatrix(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm[r,c] = ArbComplex{P}(fpm[r,c])
	    end
	end
	return arm
end

function ArbComplexMatrix(fpm::Array{F,2}) where {F<:Union{Integer,AbstractFloat}}
    P = workingprecision(ArbComplex)
    nrows, ncols = size(fpm)
    arm = ArbComplexMatrix{P}(nrows, ncols)
    for r = 1:nrows
		for c = 1:ncols
			arm[r,c] = ArbComplex{P}(fpm[r,c])
	    end
	end
	return arm
end

function Matrix(arm::ArbComplexMatrix{P}) where {P}
    nrows, ncols = arm.rowcount, arm.colcount
    m = zeros(ArbComplex{P}, nrows, ncols)
    for r = 1:nrows
        for c = 1:ncols
            m[r,c] = arm[r,c]
        end
    end
    return m
end

Matrix(x::Array{ArbComplex, 2}) = x
Matrix(x::Array{ArbComplex{P}, 2}) where {P} = x

@inline rowcount(x::ArbComplexMatrix{P}) where {P} = getfield(x, :rowcount)
@inline colcount(x::ArbComplexMatrix{P}) where {P} = getfield(x, :colcount)
@inline eachcell(x::ArbComplexMatrix{P}) where {P} = getfield(x, :eachcell)
@inline eachrow(x::ArbComplexMatrix{P}) where {P} = getfield(x, :eachrow)

@inline cellvalue(x::ArbComplexMatrix{P}, row::Int, col::Int) where {P} = eachrow(x)[row][col-1]

arbzeros(::Type{ArbComplex{P}},rowcount::SI, colcount::SI) where {P, SI<:Signed} =
    ArbComplexMatrix{P}(rowcount, colcount)

arbzeros(::Type{ArbComplex},rowcount::SI, colcount::SI) where {P, SI<:Signed} =
    ArbComplexMatrix(rowcount, colcount)

@inline Base.isempty(x::ArbComplexMatrix{P}) where {P} =
    rowcount(x) === ArbMatIdx0 || colcount(x) === ArbMatIdx0

@inline Base.size(x::ArbComplexMatrix{P}) where {P} = (rowcount(x), colcount(x))

@inline function issquare(x::ArbComplexMatrix{P}) where {P}
	return rowcount(x) === colcount(x)
end

@inline function Base.getindex(x::ArbComplexMatrix{P}, linearidx::Int) where {P}
    rowidx, colidx = linear_to_cartesian(rowcount(x), linearidx)
    return getindex(x, rowidx, colidx)
end

@inline function Base.getindex(x::ArbComplexMatrix{P}, rowidx::Int, colidx::Int) where {P}
    checkbounds(x, rowidx, colidx)
    return getindexˌ(x, rowidx, colidx)
end

@inline function getindexˌ(x::ArbComplexMatrix{P}, rowidx::Int, colidx::Int) where {P}
    z = ArbComplex{P}()
    GC.@preserve x begin
        v = ccall(@libarb(acb_mat_entry_ptr), Ptr{ArbComplex},
                  (Ref{ArbComplexMatrix}, Int, Int), x, rowidx - 1, colidx - 1)
        ccall(@libarb(acb_set), Cvoid, (Ref{ArbComplex}, Ptr{ArbComplex}), z, v)
    end
    return z
end

function Base.setindex!(x::ArbComplexMatrix{P}, z::ArbComplex{P}, linearidx::Int) where {P}
    rowidx, colidx = linear_to_cartesian(rowcount(x), linearidx)
    return setindex!(x, z, rowidx, colidx)
end

Base.setindex!(x::ArbComplexMatrix{P1}, z::ArbComplex{P2}, linearidx::Int) where {P1,P2} =
    setindex!(x, ArbComplex{P1}(z), linearidx)

Base.setindex!(x::ArbComplexMatrix{P}, z::ArbFloat{P}, linearidx::Int) where {P} =
    setindex!(x, ArbComplex{P}(z), linearidx)

Base.setindex!(x::ArbComplexMatrix{P1}, z::ArbFloat{P2}, linearidx::Int) where {P1,P2} =
    setindex!(x, ArbComplex{P1}(ArbFloat{P1}(z)), linearidx)

Base.setindex!(x::ArbComplexMatrix{P}, z::F, linearidx::Int) where {P, F<:Union{Signed,IEEEFloat}} =
    setindex!(x, ArbComplex{P}(z), linearidx)

Base.setindex!(x::ArbComplexMatrix{P}, z::R, linearidx::Int) where {P, R<:Real} =
    setindex!(x, ArbComplex{P}(BigFloat(z)), linearidx)

function Base.setindex!(x::ArbComplexMatrix{P}, z::ArbComplex{P},
                        rowidx::Int, colidx::Int) where {P}
    checkbounds(x, rowidx, colidx)
    return setindexˌ!(x, z, rowidx, colidx)
end

Base.setindex!(x::ArbComplexMatrix{P1}, z::ArbComplex{P2}, rowidx::Int, colidx::Int) where {P1,P2} =
    setindex!(x, ArbComplex{P1}(z), rowidx, colidx)

Base.setindex!(x::ArbComplexMatrix{P}, z::ArbFloat{P}, rowidx::Int, colidx::Int) where {P} =
    setindex!(x, ArbComplex{P}(z), rowidx, colidx)

Base.setindex!(x::ArbComplexMatrix{P1}, z::ArbFloat{P2}, rowidx::Int, colidx::Int) where {P1,P2} =
    setindex!(x, ArbComplex{P1}(ArbFloat{P1}(z)), rowidx, colidx)

Base.setindex!(x::ArbComplexMatrix{P}, z::F, rowidx::Int, colidx::Int) where {P, F<:Union{Signed,IEEEFloat}} =
    setindex!(x, ArbComplex{P}(z), rowidx, colidx)

Base.setindex!(x::ArbComplexMatrix{P}, z::R, rowidx::Int, colidx::Int) where {P, R<:Real} =
    setindex!(x, ArbComplex{P}(BigFloat(z)), rowidx, colidx)

@inline function setindexˌ!(x::ArbComplexMatrix{P}, z::ArbComplex{P},
                            rowidx::Int, colidx::Int) where {P}
    GC.@preserve x begin
        ptr = ccall(@libarb(acb_mat_entry_ptr), Ptr{ArbComplex},
        	        (Ref{ArbComplexMatrix}, Cint, Cint), x, rowidx-1, colidx-1)
        ccall(@libarb(acb_set), Cvoid, (Ptr{ArbComplex}, Ref{ArbComplex}), ptr, z)
    end
    return z
end

# matrix multiply

function mul!(z::ArbComplexMatrix{P}, x::ArbComplexMatrix{P}, y::ArbComplexMatrix{P}) where {P}
    ccall(@libarb(acb_mat_mul), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, x, y, P)
    return nothing
end

function matmul(x::ArbComplexMatrix{P}, y::ArbComplexMatrix{P}) where {P}
    z = ArbComplexMatrix{P}(rowcount(x), colcount(y))
    ccall(@libarb(acb_mat_mul), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, x, y, P)
    return Matrix(z)
end

function matmul(x::Array{ArbComplex, 2}, y::Array{ArbComplex, 2})
    xx = ArbComplexMatrix(x)
    yy = ArbComplexMatrix(x)
    xy = *(xx, yy)
    return Matrix(xy)
end

function matmul(x::Array{ArbComplex{P}, 2}, y::Array{ArbComplex{P},2}) where {P}
    xx = ArbComplexMatrix{P}(x)
    yy = ArbComplexMatrix{P}(x)
    xy = *(xx, yy)
    return Matrix(xy)
end

function Base.:(*)(x::ArbComplexMatrix{P}, y::ArbComplexMatrix{P}) where {P}
    checkmulable(x, y)
    return matmul(x, y)
end

@inline function Base.:(*)(x::Array{ArbComplex{P},2}, y::Array{ArbComplex{P},2}) where {P}
    checkmulable(x, y)
    return matmul(ArbComplexMatrix{P}(x), ArbComplexMatrix{P}(y))
end

@inline function Base.:(*)(x::Array{ArbComplex,2}, y::Array{ArbComplex,2})
    checkmulable(x, y)
    P = workingprecision(ArbComplex)
    return matmul(ArbComplexMatrix{P}(x), ArbComplexMatrix{P}(y))
end

# checks for validity

@inline function checkbounds(x::ArbComplexMatrix{P}, r::Int, c::Int) where {P}
    withinbounds = 0 < r <= rowcount(x) && 0 < c <= colcount(x)
    withinbounds && return nothing
    throw(BoundsError("($r, $c) not in 1:$(rowcount(x)), 1:$(colcount(x))"))
end

@inline function checkbounds(x::ArbComplexMatrix{P}, linear_rc::Int) where {P}
    withinbounds = 0 < linear_rc <= rowcount(x) * colcount(x)
    withinbounds && return nothing
    throw(BoundsError("($rc) not in 1:$(rowcount(x) * colcount(x))"))
end

@inline function checksquare(x::ArbComplexMatrix{P}) where {P}
    issquare(x) && return nothing
    throw(DomainError("matrix is not square ($rowcount(x), $colcount(x))"))
end

@inline function checkcompat(x::ArbComplexMatrix{P}, y::ArbComplexMatrix{P}) where {P}
    compat = (rowcount(x) === colcount(y) && rowcount(y) === colcount(x))
    compat && return nothing
    throw(DomainError("incompatible matrices: ($rowcount(x), $colcount(x)), ($rowcount(y), $colcount(y))"))
end

@inline function checkmulable(x::ArbComplexMatrix{P}, y::ArbComplexMatrix{P}) where {P}
    mulable = colcount(x) === rowcount(y)
    mulable && return nothing
    throw(ErrorException("Dimension Mismatach: ($rowcount(x), $colcount(x)), ($rowcount(y), $colcount(y))"))
end

#=
(::Type{Array{ArbComplex{P},2}})(::UndefInitializer, nrows::Int, ncols::Int) where {P} =
ArbComplexMatrix{P}(nrows, ncols)
Base.Matrix{ArbComplex{P}}(undef::UndefInitializer, r::I, c::I) where {P, I<:Integer} = ArbComplexMatrix{P}(m, n)
Base.zeros(::Type{ArbComplex{P}}, r::I, c::I) where {P, I<:Integer} = ArbComplexMatrix{P}(m,n)
=#

# operators over a matrix

function conj(src::ArbComplexMatrix{P}) where {P}
    dest = ArbComplexMatrix(rowcount(src), colcount(src))
    ccall(@libarb(acb_mat_conjugate), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}), dest, src)
	return dest
end

function conj(src::Array{ArbComplex{P},2}) where {P}
    return conj(ArbComplexMatrix{P}(src))
end

function transpose(src::ArbComplexMatrix{P}) where {P}
    if issquare(src)
    	dest = copy(src)
    else
    	dest = ArbComplexMatrix{P}(colcount(src), rowcount(src))
    end

    ccall(@libarb(acb_mat_transpose), Cvoid,
          (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}), dest, src)
    return dest
end

function transpose(src::Array{ArbComplex{P}, 2}) where {P}
    result = transpose(ArbComplexMatrix{P}(src))
    return Matrix(result)
end

function transpose!(x::ArbComplexMatrix{P}) where {P}
    checksquare(x)
    y = copy(x)
    ccall(@libarb(acb_mat_transpose), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}), x, y)
    return x
end

function transpose!(src::Array{ArbComplex{P}, 2}) where {P}
    result = transpose(ArbComplexMatrix{P}(src))
    src[:] = Matrix(result)[:]
    return src
end

function transpose!(dest::ArbComplexMatrix{P}, src::ArbComplexMatrix{P}) where {P}
    checkcompat(dest, src)
    ccall(@libarb(acb_mat_transpose), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}), dest, src)
    return dest
end

function transpose!(dest::Array{ArbComplex{P},2}, src::Array{ArbComplex{P},2}) where {P}
    result = transpose!(ArbComplexMatrix{P}(dest), ArbComplexMatrix{P}(src))
    dest[:] = Matrix(result)[:]
    return(dest)
end

function norm(m::ArbComplexMatrix{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_mat_frobenius_norm), Cvoid,
    	  (Ref{ArbComplex}, Ref{ArbComplexMatrix}, Cint), z, m, P)
    return z
end

function norm(m::Array{ArbComplex{P},2}) where {P}
    return Matrix(norm(ArbComplexMatrix{P}(m)))
end

function tr(x::ArbComplexMatrix{P}) where {P}
    checksquare(x)
    return trˌ(x)
end

function trˌ(x::ArbComplexMatrix{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_mat_trace), Cvoid,
    	  (Ref{ArbComplex}, Ref{ArbComplexMatrix}, Cint), z, x, P)
    return z
end

function tr(m::Array{ArbComplex{P},2}) where {P}
    return Matrix(tr(ArbComplexMatrix{P}(m)))
end

function det(x::ArbComplexMatrix{P}) where {P}
    checksquare(x)
    return detˌ(x)
end

function detˌ(x::ArbComplexMatrix{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_mat_det), Cvoid,
          (Ref{ArbComplex}, Ref{ArbComplexMatrix}, Cint), z, x, P)
    return z
end

function det(m::Array{ArbComplex{P},2}) where {P}
    return Matrix(det(ArbComplexMatrix{P}(m)))
end

"""
    determinant(<:ArbMatrix)
Computes the determinant using a more stable, albeit slower
algorithm than that used for `det`.
"""
function determinant(x::ArbComplexMatrix{P}) where {P}
    checksquare(x)
    P2 = P + P>>1
    z = ArbComplex{P2}()
    ccall(@libarb(acb_mat_det_precond), Cvoid,
          (Ref{ArbComplex}, Ref{ArbComplexMatrix}, Cint), z, x, P2)
    return ArbComplex{P}(z)
end

determinant(x::Array{ArbComplex{P},2}) where {P} = Matrix(determinant(ArbComplexMatrix{P}(x)))

function inv(x::ArbComplexMatrix{P}) where {P}
    checksquare(x)
    return invˌ(x)
end

function invˌ(x::ArbComplexMatrix{P}) where {P}
    z = ArbComplexMatrix{P}(rowcount(x), colcount(x))
    ok = ccall(@libarb(acb_mat_inv), Cint, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, x, P)
    ok !== Cint(0) && return z
    throw(ErrorException("cannot invert $(x)"))
end

inv(x::Array{ArbComplex{P},2}) where {P} = Matrix(inv(ArbComplexMatrix{P}(x)))

"""
    inverse(ArbMatrix)WARNING: Method definition transpose(ArbNumerics.ArbRealMatrix{P} where P) in module ArbNumerics at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbRealMatrix.jl:326 overwritten at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbRealMatrix.jl:338.
WARNING: Method definition rowcount(Array{T, 2}) where {T} in module ArbNumerics at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbRealMatrix.jl:138 overwritten at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbComplexMatrix.jl:127.
WARNING: Method definition colcount(Array{T, 2}) where {T} in module ArbNumerics at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbRealMatrix.jl:139 overwritten at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbComplexMatrix.jl:128.
WARNING: Method definition checkmulable(Array{T, 2}, Array{T, 2}) where {T} in module ArbNumerics at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbRealMatrix.jl:308 overwritten at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbComplexMatrix.jl:297.
WARNING: Method definition transpose(ArbNumerics.ArbComplexMatrix{P} where P) in module ArbNumerics at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbComplexMatrix.jl:314 overwritten at /home/jas/.julia/packages/ArbNumerics/HCeSA/src/libarb/ArbComplexMatrix.jl:326.

A version of `inv` with greater accuracy.
"""
function inverse(x::ArbComplexMatrix{P}) where {P}
    checksquare(x)
    P2 = P + P>>1
    z = ArbComplexMatrix{P2}(rowcount(x), colcount(x))
    ok = ccall(@libarb(acb_mat_inv), Cint, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, x, P2)
    ok !== Cint(0) && return z
    throw(ErrorException("cannot invert $(x)"))
end

inverse(x::Array{ArbComplex{P},2}) where {P} = Matrix(inverse(ArbComplexMatrix{P}(x)))

# matrix functions: exp

function exp(x::Array{ArbComplex, 2})
    checksquare(x)
    P = workingprecision(ArbComplex)
    y = ArbComplexMatrix(x)
    z = ArbComplexMatrix{P}(rowcount(x), colcount(x))
    ccall(@libarb(acb_mat_exp), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, y, P)
    return Matrix(z)
end

function exp(x::Array{ArbComplex{P}, 2}) where {P}
    checksquare(x)
    y = ArbComplexMatrix(x)
    z = ArbComplexMatrix{P}(rowcount(x), colcount(x))
    ccall(@libarb(acb_mat_exp), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, y, P)
    return Matrix(z)
end

function exp(x::ArbComplexMatrix{P}) where {P}
    checksquare(x)
    z = ArbComplexMatrix{P}(rowcount(x), colcount(x))
    ccall(@libarb(acb_mat_exp), Cvoid, (Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Cint), z, x, P)
    return Matrix(z)
end

## eigenvalues and eigenvectors

# sort complex values
complex_lt = (a,b)->(real(a) < real(b)) || (real(a) == real(b) && imag(a) < imag(b))

#=
int acb_mat_approx_eig_qr(acb_ptr E, acb_mat_t L, acb_mat_t R, const acb_mat_t A, const mag_t tol, slong maxiter, slong prec)

            acb_mat_t C, R;
            acb_ptr E;

            acb_mat_init(R, n, n);
            acb_mat_init(C, n, n);
            E = _acb_vec_init(n);

            acb_mat_set_arb_mat(C, A);
            acb_mat_approx_eig_qr(E, NULL, R, C, NULL, 0, prec)

Computes floating-point approximations of all the n eigenvalues (and optionally eigenvectors) of the given n by n matrix A.
The approximations of the eigenvalues are written to the vector E, in no particular order.
If L is not NULL, approximations of the corresponding left eigenvectors are written to the rows of L.
If R is not NULL, approximations of the corresponding right eigenvectors are written to the columns of R.

The parameters tol and maxiter can be used to control the target numerical error and the maximum number of iterations allowed
before giving up. Passing NULL and 0 respectively results in default values being used.

Uses the implicitly shifted QR algorithm with reduction to Hessenberg form.
No guarantees are made about the accuracy of the output.
A nonzero return value indicates that the QR iteration converged numerically,
but this is only a heuristic termination test and does not imply any statement whatsoever about error bounds.

The output may also be accurate even if this function returns zero.
=#
#=

function LinearAlgebra.eigvals(m::ArbComplexMatrix{P}) where {P}
    checksquare(m)
    n = rowcount(m)
    eigvalues = ccall(@libarb(_acb_vec_init), Ptr{ArbComplex}, (Cint,), n)
    eigvectors = ArbComplexMatrix(rowcount(m), colcount(m))
    eigvectors2 = copy(m)
    tol = Base.C_NULL
    maxiter = 0
    result = ccall(@libarb(acb_mat_approx_eig_qr), Cint,
		  (Ptr{ArbComplex}, Ptr{Nothing}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ptr{Nothing}, Clong, Clong),
		  eigvalues, Base.C_NULL, eigvectors, eigvectors2, Base.C_NULL, maxiter, P)
    eigvalues = sort(eigvalues, lt=complex_lt)
    return eigvalues
end

function LinearAlgebra.eigvecs(m::ArbComplexMatrix{P}) where {P}
    checksquare(m)
    eigvalues = zeros(ArbComplex, rowcount(m))
    eigvectors = ArbComplexMatrix(rowcount(m), colcount(m))
    tol = Base.C_NULL
    maxiter = 0
    result = ccall(@libarb(acb_mat_approx_eig_qr), Cint,
		  (Ref{Vector{ArbComplex}}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ptr{Nothing}, Clong, Clong),
		  eigvalues, Base.C_NULL, eigvectors, m, tol, maxiter, P)
    return eigvectors
end

function LinearAlgebra.eigen(m::ArbComplexMatrix{P}) where {P}
    checksquare(m)
    eigvalues = zeros(ArbComplex, rowcount(m))
    eigvectors = ArbComplexMatrix(rowcount(m), colcount(m))
    tol = Base.C_NULL
    maxiter = 0
    result = ccall(@libarb(acb_mat_approx_eig_qr), Cint,
		  (Ref{Vector{ArbComplex}}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ref{ArbComplexMatrix}, Ptr{Nothing}, Clong, Clong),
		  eigvalues, Base.C_NULL, eigvectors, m, tol, maxiter, P)
    eigvalues = sort(eigvalues, lt=complex_lt)
    return eigvalues, eigvectors
end

LinearAlgebra.eigvals(m::Array{ArbComplex{P},2}) where {P} = eigvals(ArbComplexMatrix(m))
LinearAlgebra.eigvecs(m::Array{ArbComplex{P},2}) where {P} = eigvecs(ArbComplexMatrix(m))
LinearAlgebra.eigen(m::Array{ArbComplex{P},2}) where {P} = eigs(ArbComplexMatrix(m))
=#