0
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ The issue, of course, is that if you want to symmetrize the sparsity pattern, you don't yet know how many entries you need in each row. So you can't allocate the ias and jas arrays up front because you don't know how many entries they are going to have. Some kind of temporary storage is going to be necessary. $\endgroup$ – Wolfgang Bangerth Nov 24 '20 at 21:37
  • $\begingroup$ Assuming I find the transpose of $A^T$ and store it in arrays in ib and jb, I can find the size of jc (equal to number of non-zero elements in $C = A+A^T$) by passing through all of the rows and counting the number of common elements in ja and jb using a switch array, effectively populating the ic array. However, it might be better to use some estimate for the size of jc (knowing that the degrees of my nodes are roughly equal) and also store the column indices with one pass. $\endgroup$ – IPribec Nov 25 '20 at 11:14
  • $\begingroup$ I have looked into the work of Davis. The sparse matrix addition is actually recast as a multiplication $C = \alpha A + \beta B = [A \; B][\alpha I \; \beta I]^T$. The matrix C is allocated with storage of size anz + bnz, and then the extra workspace is freed upon finishing. $\endgroup$ – IPribec Nov 25 '20 at 11:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.