function sprandsym(n, density, rc=1e12, kind=1)

    # Create random sparse symmetric matrix
    A = sprand(n, n, density)
    A = (A + A') / 2

    # Get current eigenvalues
    λ_current = eigvals(Matrix(A))

    # Create target eigenvalues
    if kind == 1  # Positive definite
        λ_target = exp.(range(log(1.0), log(rc), length=n))
    elseif kind == 0  # Indefinite
        n_pos = div(n, 2)
        n_neg = n - n_pos
        λ_pos = exp.(range(log(1.0), log(rc), length=n_pos))
        λ_neg = -exp.(range(log(1.0), log(rc), length=n_neg))
        λ_target = vcat(λ_pos, λ_neg)
    else  # Negative definite
        λ_target = -exp.(range(log(1.0), log(rc), length=n))
    end

    shuffle!(λ_target)

    # Get eigenvectors of sparse matrix
    F = eigen(Matrix(A))
    V = F.vectors

    # Reconstruct with target eigenvalues
    A_new = V * Diagonal(λ_target) * V'

    # Sparsify by keeping only the largest elements in magnitude
    # to approximately match the original sparsity pattern
    target_nnz = nnz(A)
    A_vec = vec(A_new)
    threshold = sort(abs.(A_vec), rev=true)[min(target_nnz, length(A_vec))]

    A_new[abs.(A_new) .< threshold] .= 0

    return sparse(Symmetric(A_new+A_new'/2))

end

function MatrixGenerator1(N)
    n = N^2
    # Main diagonal : -4
    main_diag = fill(4.0, n)
    # Off - diagonals at ±N: -1
    off_diag_N = fill(-1.0, N * (N - 1))
    # Off - diagonals at ±1 within each row

    off_diag_1 = fill(-1.0, n - 1)
    off_diag_1[N:N:end] .= 0
    # Create sparse matrix
    A = spdiagm(-1 => off_diag_1 ,
                -N => off_diag_N ,
                0 => main_diag ,
                N => off_diag_N ,
                1 => off_diag_1 )
    return A
end

function MatrixGenerator2(N)
    n = N^2
    # this typically has smaller condition number than the one in MATLAB
    A = sprandsym(N^2,0.2,0.1) + 2*Diagonal(ones(N^2))
end