I have a class of matrices $A$ which are created by a domain decomposition method. Each matrix represents several subproblems of equal size, and I know that for some permutation matrix $P$, $PAP^T$ will be a block diagonal system, i.e. $$ A = \begin{bmatrix} A_1 & 0 & 0 \\ 0 & A_2 & 0 \\ 0 & 0 & A_3 \\ \end{bmatrix} $$ Solving this system for some right hand side is obviously just a matter of inverting each block and permuting back to he original ordering. The question is if there are any methods for such systems when the ordering required to permute the system into block diagonal form is not known.

  1. One approach would be to interpret the matrix system as a graph, and try to find cycles. The system can then be permuted by ordering according to each cycle. The question is then what algorithms are fast and suitable for the purpose?
  2. The condition number of such a system is often much lower than that of an similar system where the blocks are connected because the eigenvalues are that of each subblock. A question of general interest could be what algorithms are good at solving such systems without actually permuting the system - CG would probably converge faster for a more decoupled system because of the lowered condition number. Are some algorithms better for this type of problem?
  3. Is there a known nomenclature for such systems? When searching for inversion/permutation of block diagonal matrices, what I usually find is systems which are block diagonal except for the last row - a different problem, really.

While I have no problems solving each system, I thought this question was interesting. The submatrices can be thought of as similar to the five point stencil for the heat equation in terms of matrix properties - they arise from a similar mass conservation system. (I hope this is of general interest, so focusing on the structure of the submatrices is not very important)

  • 2
    $\begingroup$ Is it common not to know enough about your underlying problem not to be able to block-structure it from the start? It seems to me if you know it's got this structure, there's a way to reorganize it from the beginning to generate the matrices in block form. This should generally be true regardless of whether or not one can explicitly write down the permutation matrices easily. $\endgroup$ – Bill Barth Apr 26 '12 at 13:13

Q1 To find the individual blocks in the matrix (and hence recover the block permutation) you could find the connected components in the underlying graph. This can be done efficiently ( in O(|E|) for a graph with E edges) via repeated breadth-first-searches from unvisited vertices:

for (i = all vertices in graph G)
    if (vertex i is unvisited)
        // start a BFS from vertex i
        // all of the vertices visited during this BFS will be 
        // part of the same diagonal block

Note that this only identifies the block structure in the matrix - there are still a whole family of perumtations that can be applied locally within each block. Depending on your choice of linear solver you would then choose an appropriate ordering scheme (AMD, RCM, etc).

Q2 I don't understand what you're getting at here - since the diagonal blocks are fully decoupled from each other I don't see why you wouldn't just solve them individually - in parallel if you're interested in speed.

Q3 Block diagonal is the correct name - as far as I know.


On Q1: In the FEM/FDM language, you would apply the Cuthill-McKee algorithm to number all degrees of freedom that are connected to each other. Once the algorithm stops without having numbered all degrees of freedom, you know that there's another as yet unnumbered component of the matrix, and you would continue with one of its degrees of freedom, with the next available index. This way, you end up with the block diagonal structure of the matrix -- in other words, the Cuthill-McKee algorithm will compute the permutation for you.

On Q2: All Krylov-space methods have convergence properties independent of the numbering of degrees of freedom (i.e., they are invariant under permutations) as long as you have a preconditioner that is invariant.


Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.