Another possibility, at least for relatively small matrices, is to take the determinant (strictly speaking it is the permanent that is required, I suppose).
For example, for an $11 \times 11$ matrix (o=5), I find there are 7 solutions.
primePositions5 = Position[With[{o = 5}, Table[If[PrimeQ[n], 1, 0], {m, 0, Prime[o]^2 - Prime[o], Prime[o]}, {n, m + 1, m + Prime[o]}]], 1]; mylist = List @@ (Det@ SparseArray[## -> Subscript[a, ##] & /@ primePositions5] /. {-x_ -> x})
gives the following:

Matrix plots of all seven solutions:
MatrixPlot[ Normal@SparseArray[(List @@ #) /. Subscript[a, {x_, y_}] -> {x, y} -> 1], Mesh -> All, ImageSize -> 200] & /@ mylist
I'll give then as a grid:

I reckon it needs to be emphasized that the Mathematica's Det command is slow.
With o=7 which gives a $17 \times 17$ matrix, I obtain 2144 solutions. For 0 =8 ($19 \times 19$), the figure is 2641. I could not go beyond this with the computer I am using (with Mathematica 7, as it so happens).
For o=4 ($7 \times 7$), I get two solutions:

Update for Mathematica 11
In Mma 11, we can use the Permanent function
myListAlt = List @@ (SparseArray[## -> Subscript[a, ##] & /@ primePositions5] // Permanent // Expand)
The behaviour of Det seems to have changed somewhat since this question was posted.
I now need to Expand the result of the Det function:
mylist = List @@ (Expand@ Det@SparseArray[## -> Subscript[a, ##] & /@ primePositions5] /. {-x_ -> x})
and
mylist == myListAlt
True