# How to assemble and solve a matrix system in parallel from values generated in different processors?

I am solving a multiscale problem using the Heterogeneous Multiscale Method (HMM). Essentially, my particular procedure uses the following iterative process:

1. Solve many local matrix systems.
2. Compute a value of interest from the solutions of the local systems.
3. Assemble a global matrix system from the local "values of interest"
4. Solve the global matrix system
5. Use the solution of the global matrix system to form new local matrix systems.

Repeat until some convergence criteria is met.

Since there are many local (independent) linear systems of equations and multiple systems can fit into local RAM memory, I figure it is best to load multiple "local" systems into each processor and solve each system sequentially (see this posted question).

My question regards the best strategy to assemble and solve the global matrix system. In my particular case, the global matrix system is small enough that it could fit entirely on any processor's RAM memory. Furthermore, the local & global matrices do not change size between iterations. So, I foresee one of three possible strategies:

1. Gather the "values of interest" onto a single processor, and assemble/solve the global matrix system sequentially on one processor.
2. Copy the values of interest onto every processor, and assemble/solve the same global matrix system sequentially on every processor.
3. Assuming that each processor possesses the "values of interest" necessary to produce contiguous blocks of the global matrix, then we can assemble partitions of the global matrix locally, then solve them together in parallel.

I can see some advantages/disadvantages to each method. In Method 1, no communication is necessary in the solving phase, but communication to and from the root processor may become a bottleneck (especially at scale). Method 2 may require more interprocessor communications to assemble the global matrix than the first method, but no communication is needed in the solving phase or in the local matrix assembly stage that follows. Method 3 requires no interprocessor communication for assembly of the local or global matrices, but requires it in the solving phase.

Suppose that each local system is on the order of $10^3$x$10^3$ and there are $10^3$x$10^3$ local matrix systems. Let's further suppose that the global matrix system has size $10^3$x$10^3$. Under these assumptions, which one the three aforementioned strategies will likely lead to a faster solution of the global system? Are there other mapping strategies for the global matrix that might work faster per iteration?

• Very interesting question. I hope someone has good answers. – Inquest Feb 12 '12 at 17:45
• Do you have an idea about how big the global system is in relation to the local systems? I.e., if there are $n$ local systems to be solved, is the global system $kn \times kn$ for some $k$? Do you have an idea for how big $n$ is? The answers to your questions are likely to depend heavily on the sizes. – Bill Barth Feb 13 '12 at 1:11
• @BillBarth: Let's suppose that n is on the order of $10^6$, and we want k to get increasingly larger. – Paul Feb 13 '12 at 1:20
• So the answer to my first question is "yes"? And how large do you want $k$ to get? I.e., are you eventually going to extract a million parameters from the local systems, or will it stay relatively small compared to $n$? How big are the local systems? Finally, are all the systems denser or sparse? – Bill Barth Feb 13 '12 at 4:06
• @BillBarth: For now, let's say that $k<100$ and the global matrix will extract only one parameter from each of the linear systems. The size of the local systems may vary from are $O(n)$ where n is the size of the global matrix, and all linear systems (local and global) are sparse, symmetric, positive definite and diagonally dominant. – Paul Feb 13 '12 at 4:13

I do not think there is any case where you want to solve on rank 0. Redundant solve is almost always better since, for small things, allreduce is as efficient as reduce, and redundant computation only has one instead of two.

However, whether to compute redundantly on all nodes, or on a subset, or redundant subsets depends on the hardware and system size. Thus, you should have a system that can do any of them. The PCREDUNDANT in PETSc can solve redundantly on all processes, some processes, or subsets of processes in parallel.

But if the global problem is of size $10^6$ as you claim in the comments, it is large enough to benefit significantly from a parallel solve. Parallel assembly is very much the standard and recommended scenario.

• I tried your suggestion of using a redundant global solve (using MPI_Allgatherv), and compared it to a rank-0-only global solve (using MPI_Gatherv and MPI_Bcast), and tested it up to problem sizes of $N=4096$ unknowns in the global system. It seems that the redundant solve has been consistently a bit slower than the rank-0-only solve. I suspect either the MPI implementation or the network hardware to be a likely reason. Does this seem plausible/likely? – Paul Apr 4 '12 at 5:02