How to restore diagonal-symmetric matrix that has been shuffled?
First, note that the matrix can only be restored up to a permutation of the diagonal.
To illustrate on your example, when we swap two bottom rows and two left columns, the resulting matrix is also diagonal-symmetric:

133 199 101 121                             199 133 101 121
142 133 199 101        transforms to        133 142 199 101
156 142 133 199                             156 108 142 133
108 156 142 133                             142 156 133 199


So, the good news is, instead of searching for a row permutation and a column permutation, we can fix the order of the rows and permute only columns.
Indeed, if we took the original matrix and permuted the rows in any fixed order,
by permuting the columns in the exact same order, we get a diagonal-symmetric matrix.



Next, note that, before the permutations, the multiset of numbers in row $1$ was exactly the same as the multiset of numbers in column $1$.
The same holds for row $2$ and column $2$, and so on.

After permutations, say row $1$ landed at row $i$, and column $1$ landed at column $j$.
Then, the multiset of numbers at row $i$ is still the same as the multiset of numbers at column $j$.

Importantly, the problem seems to be about some large and more-or-less real data, so it is likely that the converse also holds: for any two rows, the multisets of numbers in them are different.
If not, perhaps what follows can be tricked to fail by constructing the appropriate data.
Anyway, the check is straightforward.



Having noted the above, I propose the following simple algorithm:


Write down the multiset of numbers $S_1$ in the bottom row.
(A sorted array will suffice for storage and comparison.)
Search for the column which contains the exact same multiset of numbers $S_1$.
Swap this column with the leftmost one.
Proceed with the next-to-bottom row and find its multiset $S_2$.
Search for the column except the first one which contains the exact same multiset of numbers $S_2$.
Swap this column with the second-to-leftmost one.
Et cetera.




The algorithm can run in $O (n^2 \log n + n^2)$, with some hashing.
First, for each of the $n$ rows and each of the $n$ columns, calculate a hash of the multiset of numbers in it (for example, make a sorted array in $O (n \log n)$ and get its polynomial hash).
Then, we can consider hashes instead of the rows and columns of the real matrix.
After we find the needed permutation of rows, apply it once in $O (n^2)$.