# 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. 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. 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? 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

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?