How to get the dimension of Atkin-Lehner eigenspace or do you have any data already obtained? I need the dimensions of the Atkin-Lehner eigenspace for the paper I'm writing.
As is well known, the cuspidal space $S_{k}(\Gamma_{0}(N))$ can be decomposed by Atkin Lehner involution. For example, when $N=21=3\times7$, we have the decomposition
$$
S_{k}(\Gamma_{0}(N))=S^{(++)} \oplus S^{(+-)} \oplus S^{(-+)} \oplus S^{(--)}
$$
where $S^{(++)}$, $S^{(+-)}$, $S^{(-+)}$, and $S^{(--)}$ are the subspaces of $S_{k}(\Gamma_{0}(N))$ for which pairs of eigenvalues for $W_3$, $W_7$ are (+1,+1), (+1,-1), (-1,+1), (-1,-1) respectively.
In addition, dimensions of each summand spaces are 0,0,1,0.
This is from Table 5 of Antwerp IV.
For k=2, the dimension of each Atkin-Lehner eigenspace has already been given by Table 5 of Antwerp IV and David Kohel.
And for any weight k, I found the following documentation.
However, using this, it is difficult for me to find the dimensions of the eigenspaces.
So, is there any idea or data?
 A: Exact formulas for dimensions of Atkin-Lehner eigenspaces follow from trace formulas of Yamauchi and Skoruppa-Zagier.  Skoruppa-Zagier corrected some clerical errors in Yamauchi's paper.  See:

Nils-Peter Skoruppa and Don Zagier, Jacobi forms and a certain space of modular forms, Invent. Math. 94 (1988), no. 1, 113–146.

In the case of squarefree level, I worked things out explicitly in this paper:

Kimball Martin, Refined dimensions of cusp forms, and equidistribution and bias of signs, J. Number Theory 188 (2018), 1–17.

The main focus is dimensions of new parts of Atkin-Lehner eigenspaces, but dimensions which include oldspaces can be computed similarly.  Either you can take the newspace dimensions, and add in the old form contribution, or just use the trace formulas on the full spaces $S_k(N)$ to get a formula mimicking what I did for the newspaces.
Code to compute dimensions (and a link to my paper) in the squarefree level case is available here:
https://math.ou.edu/~kmartin/data/
Using these formulas should be faster than the direct calculations in David Loeffler's answer.  You can also code up the non-squarefree level case using Skoruppa-Zagier's trace formula without too much trouble.
Also, in case you're not familiar with it, the LMFDB modular forms page has tabulated a lot of data for newforms, including Atkin-Lehner signs.  You can use this check dimensions in many cases.
A: You can compute these dimensions using modular symbols (an auxiliary space which has the same Hecke action as modular forms, but is easier to compute). Here's a Sage example for weight 4 cusp forms of level Gamma0(17):
sage: S=ModularSymbols(Gamma0(17), weight=4, sign=1).cuspidal_submodule()
sage: S.atkin_lehner_operator().charpoly().factor()
(x + 17) * (x - 17)^3

Sage normalises the Atkin–Lehner operator so that its square is multiplication by $(-N)^{k-2}$, since that normalisation works better for odd weights (it avoids introducing square roots). If you want to normalise so that the operator becomes an involution, as is usual for $\Gamma_0$ levels, you need to scale by $N^{(k - 2) / 2}$. So, in my example, the scaling factor is 17, and the output means that the +1 eigenspace has dimension 3, and the -1 eigenspace has dimension 1.
You can compute the local AL operators for each prime $p \mid N$ similarly, and get the simultaneous eigenspaces by taking intersections of kernels – here's a code snippet which does this:
sage: def AL_eigenspace_dims(N, k):
....:     facs = [p^r for (p, r) in N.factor()]
....:     S = ModularSymbols(Gamma0(N), sign=1, weight=k).cuspidal_submodule()
....:     Wmats = [S.atkin_lehner_operator(f).matrix() for f in facs]
....:     for s in cartesian_product([ [-1, 1] for f in facs ]):
....:         Vs = [(Wmats[i] - s[i]*facs[i]**((k-2)/2)).kernel() for i in range(len(facs))]
....:         V = reduce(lambda u,v: u.intersection(v), Vs)
....:         print(s, V.dimension())
....: 
sage: AL_eigenspace_dims(21, 2)
(-1, -1) 0
(-1, 1) 1
(1, -1) 0
(1, 1) 0
sage: AL_eigenspace_dims(21, 12)
(-1, -1) 7
(-1, 1) 6
(1, -1) 7
(1, 1) 8

This modular symbol method will be slow if k or N is large, because it involves linear algebra with matrices having about $kN$ rows and columns. (For $k N$ around 500, it takes roughly ten seconds on my machine). In contrast, the algorithm in Lloyd Kilford's draft code from your link uses the Riemann–Roch theorem, and the running time is dominated by factorising N, so it should be practical for N up to about 100 digits; but his implementation only covers the case of weight 2. In principle it should be possible to extend the Riemann–Roch method to higher weights, but I don't know if anyone has implemented this in Sage.
