Given a weighted digraph, I can check whether a given vertex belongs to a negative cycle in $O(|V|\cdot|E|)$ using Bellman-Ford. But what if I need to find all vertices on negative cycles? Is there a way to do it faster than Floyd-Warshall's $O(|V|^3)$?
If you don't constrain yourself to simple cycles, you can actually use the Bellman-Ford algorithm to find all the relevant vertices.
Start by running a DFS on the graph to find its strongly connected components. Let them be $G_1,G_2,...G_k$. For each $G_i$, run BF on the subgraph. If the BF algorithms detects a negative cycle (which it can after $|V_i|+1$ iterations where $V_i$ is the vertices in $G_i$), add $V_i$ to your list of vertices that belong to a negative (not necessarily simple) cycle.
If one of the strongly connected components contain a negative cycle, then because $G_i$ is finite and so are the weights, you can loop through the negative cycle as much as you need, and then add go in a cycle through all the vertices in $G_i$. You only need to loop an amount that will guarantee a total negative cycle.
It is easy to show that there are no cycles between any $G_i,G_j$ and therefore there are no negative cycles other than inside the strongly connected components.
Since you run BF on each subgraph seperately, you keep the time complexity of BF which is $O(|V|\cdot|E|)$.
I was not able to turn up any better algorithm in my research.
In practice, you could improve the running time for many graphs by first decomposing the graph into strongly connected components, then running Floyd-Warshall on each strongly connected component. This does not improve the worst-case complexity (since in the worst case the entire graph could form one large strongly connected component), but it might help on many graphs you're likely to run into in practice.
I just realized that my answer below may be incorrect. The intended idea was to perform |V|-1 waves of "relaxation", as in Bellman-Ford (see CLRS). At the end, instead of returning FALSE to denote "negative weight cycle found, so Bellman-Ford failed to find a shortest path tree", we just print the cycle, then continue looking for more cycles to print.
However, it is doubtful (or at least unclear) whether this will actually find all negative weight cycles (or at least simple cycles). For example, if the |V|-1 waves of relaxation only recorded tree edges for extremely negative cycles, but not for moderately negative cycles, then the latter may not be picked up by our print statements, since our print statements only detect cycles that are formed by tree edges.
Getting the negative weight cycle is possible with a minor modification to Bellman-Ford.
The implementation of Bellman-Ford in CLRS instantly returns FALSE whenever a negative-weight cycle is found:
BELLMAN-FORD(G, w, s)
for i = 1 to |G.V| - 1
for each edge(u,v) in G.E
for each edge(u,v) in G.E
if v.d > u.d + w(u,v)
Instead of returning FALSE, we look at the condition
v.d > u.d + w(u,v), which indicates even after |V|-1 passes through all edges,
v.d can still be improved. This can only mean that there's a path from
u through tree edges found during the RELAX loop, for which completing the path with edge
(u,v) would cause a negative weight cycle. Hence, to print out the cycle, just follow
v and its "parent" pointers all the way up to
u. That path plus
(u,v) is your negative weight cycle.