How to make a directed graph symmetric?

Say I have a directed graph given as an adjacency matrix $A$ in CSR format represented by the arrays ia (row indexes) and ja (column indexes).
In my application the graph serves as the underlying spatial connectivity for an RBF-FD discretization of a PDE. The graph originates from finding the $N$-nearest neighbor stencil for each node $boldsymbol{x}_i$ of a $d$-dimensional point cloud.

For some purposes it is benefical to use symmetric stencils $S$, meaning that $boldsymbol{x}_j in S(boldsymbol{x}_i)$ implies $boldsymbol{x}_i in S(boldsymbol{x}_j)$ for all $i,j$ (or at least for interior nodes). This corresponds to making the directed graph symmetric (undirected). One example of where this is needed is in the METIS partitioning library, which expects graphs to be symmetric. Similarly, the Reverse Cuthill McKee ordering algorithm also expects a symmetric matrix.

Question: Given an adjacency matrix in CSR format as arrays ia and ja, how can I find the symmetric graph adjacency matrix arrays ias, jas?

I have noted in Scipy, the way they achieve this is by forming the matrix $A + A^T$ (see Scipy source here). Note that I don’t want to compute the actual matrix values, but only the structure of the resulting matrix.

Is there any other way to achieve this, apart from the naive way of of initializing $M$ empty adjacency lists (for all $M$ nodes of my point cloud), iterating through all the nodes, and pushing each connection $(i,j)$ into the right list? It seems kind of unneccesary to process each connection twice.