Assuming A, B, C, and D are square, which seems reasonable since you could otherwise have an over or under determined system, the usual thing to do would be to use Kronecker product properties to rewrite this as a matrix equation
where x=\mbox{vec}(X) and v = \mbox{vec}(V). Here \mbox{vec} stacks the columns of a matrix into a vector. This is a generalized Sylvester equation. There are then solvers that compute generalized Schur decompositions of (B,D) and (A,C) to allow for the efficient direct column-wise solution of X using a back-substitution procedure. If every matrix is n\times n, it’s an O(n^3) algorithm to compute the n^2 elements of X. I believe that MatrixEquations.jl has this implemented with the function gsylv
.
Even if your matrices are sparse there is limited benefit in trying to use an iterative method to exploit sparsity, since X will not in general be sparse. The direct method should be the way to go.