If I understand your question correctly, I believe there are two ways to tackle this problem:
If the issue is just in the generation of the constraints $Ax \leq a$ and $Bx=b$, then these forms can be forced out using automatic differentiation. Specifically, let $g(x) = Ax -a$ and $h(x) = Bx - b$. To generate, $a$ and $b$, we evaluate $g(0)$ and $h(0)$, respectively. In order to generate $A$ and $B$, take an arbitrary point $x$ and evaluate $g^\prime(x)$ and $h^\prime(x)$ using the Jacobian code from an available automatic differentiation toolbox to generate $A$ and $B$, respectively. Even with a linear operator, it's sometimes easier to just use AD to generate them in the case that we have a matrix-free implementation of the operator. Note, it's likely that the AD code will want to fix the number of variables, so this may not be the appropriate solution. However, it depends on the code and the setup that you have.
Matrix-free implementations of an interior point method do not require the number of variables to be known. Typically, they just require certain vector operations to be defined. There may be a different list, but I typically work with
init(x)
- Initialize memory
copy(x)
- copy a vector
scal(alpha,x)
- $\alpha x$
axpy(alpha,x,y)
- $\alpha x + y$
innr(x,y)
- $\langle x,y\rangle$
prod(x,y)
- pointwise product
linv(x,y)
- inverse of the pointwise product
srch(x,y)
- how far we can go in the direction y form x till we violate nonnegativity
barr(x)
- log-barrier function
Technically, we can generalize this to work with general symmetric cone constraints, but we should ignore that for now. Anyway, none of these operations require the size of the vector to be defined as long as we can compute them. Then, based on these operations, we have enough information to implement an interior point method to take care of the inequality constraints and an SQP method to take care of the equality constraints. Now, candidly, if the number of variables is large, it's going to perform poorly. Basically, there's a number of linear systems that need to be solved. However, if the number of equality constraints is bigger than about 50 the Krylov methods used to solve these linear systems will start to deteriorate badly. The way to get around this is to use a preconditioner, but without a matrix you're going to have to be very clever. That'll depend on the problem structure.
In terms of existing codes and algorithms, Optizelle can do this. However, again, you really, really will need an effective preconditioner for the equality constraints in order to have a shot at making this work well. Specifically, a preconditioner for $BB^T$.
As a final note, in general, for a linear program, you typically want to use a linear programming specific algorithm like a specialized primal-dual interior point method or the simplex method. Typically, for the LP specific primal-dual interior point methods, you'll need to be able to take the Choleski factorization of a very particular system. If you don't know the number of variables, this won't work. That pushes you back to more general algorithms like those that Optizelle uses. That said, this code is not specialized for linear programs, so the performance will suffer.