There's no need to use the Schur complement here because $A$ and $S$ are already symmetric. The conventional formulation of this problem as an SDP is
$$\min t \quad\text{ subject to}\\
A-S+tI \succeq 0 \\
A-S - tI \preceq 0$$
results in an SDP of the form (I'm specific here because there are so many different "standard forms" for SDP)
$$\min c^{T}x\\
F_{0}+x_{1}F_{1}+...+x_{m}F_{m} \succeq 0$$
where
$$x=[t\;\; S_{1,1} \;\; S_{1,2} \;\; S_{1,3} \; \ldots \; S_{n,n}]^{T}$$
has $1+n(n-1)/2$ elements, and the matrices $F_{i}$ are of sizes $2n \times 2n$.
For your problems with $n=500$, this means that you have about 125,000 variables $x_{j}$, and matrices of size $1000$ by $1000$. Primal-dual interior-point methods (as implemented in SDPA, SDPT3, SeDuMi, CSDP, etc.) would require the solution of a system of 125,000 equations in 125,000 variables in each iteration, and this system would typically be fully dense.
This is not something you can do on a typical desktop machine but would be within reach on a more powerful server with enough RAM (more than 128 gigabytes.) Solution times would be fairly long (10+ hours)
The "splitting cone solver" that Julia implements appears to me to be an implementation of a first-order ADMM method such as the one described in
B. O'Donoghue, E. Chu, N. Parikh, and S. Boyd. Operator Splitting for Conic Optimization via Homogeneous Self-Dual Embedding.
The problem with ADMM for SDP is that it doesn't yield very accurate solutions (you'll typically get two or three digits of accuracy at best) and the method requires lots of tuning to achieve good results. It's not nearly as robust (in the sense that it can consistently solve every problem that you throw at it) as the primal-dual interior-point methods.
I can't vouch for the particular implementation of this method in Julia since I've never worked with it. It might be that some other first-order code for SDP could do a better job than SCSSolver.
There is probably some work in the convex optimization literature on this kind of spectral norm minimization that avoids even formulating an SDP, but I'm no expert on that topic.
Are you sure that you really need the best approximation in the matrix 2-norm sense rather than in the Frobenius norm?
Note added much later: It turns out that the approach that computes this projection in the Frobenius norm also works equally well to compute the projection in the spectral norm. First, perform an orthogonal diagonalization of $A$ as
$A=\sum_{i=1}^{n} \lambda_{i} u^{(i)}u^{(i)^{T}} $
Then
$S=\sum_{i=1}^{n} \lambda_{i}^{+} u^{(i)}u^{(i)^{T}} $
When I wrote this answer, I had known that this worked for the Frobenius norm projection, but not that it also works for the 2-norm projection. For non-symmetric $A$, or if there are additional constraints, this simple formula isn't applicable and the SDP formulation might still be useful.