How to efficiently populate a symmetric matrix?

One way to create your data structure is to realize that your list of integers is closely related to all the Tuples taken in groups. For the case I=J=2, the permutations are of 1 and 2 taken 4 at a time:

mat = Partition[Tuples[{1, 2}, 4], 4];
mask = UpperTriangularize[ConstantArray[1, {4, 4}]];
outMat = mat mask + Transpose[mat mask] - mat IdentityMatrix[4]

outMat//MatrixForm

Your $m-1, n-1$ and $p-1, q-1$ are the two digits of $k-1$ and $l-1$, respectively, in base $J$. This should be pretty fast:

dim = i*j;
A = ConstantArray[0, {dim, dim}];
Do[
   A[[k, l]] = A[[l, k]] = cf[Nfunc, xi, yi, wix, wiy, Quotient[k - 1, j] + 1, Mod[k - 1, j] + 1, Quotient[l - 1, j] + 1, Mod[l - 1, j]+ 1]
, {k, 1, dim}, {l, k, dim}
]

If $IJ$ is huge and your cf is very very fast, then you can probably shave off a little bit more time by making an outer Do over k where you compute the quotient and mod once for that k, followed by an inner Do over l.