Brandt,<br><br><div class="gmail_quote">On Tue, Nov 15, 2011 at 11:12 AM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov">jedbrown@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div class="gmail_quote"><div class="im">Jack Poulson (Cc'd) may have some advice because he has been doing this for high frequency Helmholtz.</div></div></blockquote><div><br></div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">
> The matrix is block tri/penta-diagonal, depending on the stencil, and the blocks are also </div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">
> tri/penta-diagonal. Correct me if I'm wrong, but I believe these types of matrix equations can be </div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">
> solved directly and cheaply on one node. </div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); "><br></div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">
These solves should be cheap, but since the bandwidth is not small, it makes more sense to use a sparse-direct solver (e.g., MUMPS or SuperLU). The complexity of banded factorizations for matrices of size N x N with bandwidth of size b is O(b^2 N), and the memory and solve complexity is O(bN). For 2d sparse-direct the factorization complexity is O(N^{3/2}) and the memory and solve complexities are O(N log(N)). Thus, when the bandwidth is larger than O(N^{1/4}), it makes sense to consider sparse direct in order to make the factorization cheaper. I suspect that in your case, b=sqrt(N).</div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); "><br></div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">
> The two options I'm comparing are:</div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> 1. Distribute the data in z and solve x-y plane matrices with LAPACK or other serial or shared-</div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> memory libraries. Then do an MPI all-to-all to distribute the data in x and/or y, and do all </div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> computations in z. This method allows all calculations to be done with all necessary data </div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> available to one node, so serial or shared-memory algorithms can be used. The disadvantages </div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> are that the MPI all-to-all can be expensive and the number of nodes is limited by the number of</div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> points in the z direction.</div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">
> </div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> 2. Distribute the data in x-y only and use petsc to do matrix solves in the x-y plane across </div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> nodes. The data would always be contiguous in z. The possible disadvantage is that the x-y </div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> plane matrix solves could be slower. However, there is no need for an all-to-all and the number</div>
<div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> of nodes is roughly limited by nx*ny instead of nz.</div><div style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">
> </div><div><span class="Apple-style-span" style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> The size of the grid will be about 1000 x 100 x 500 in x, y, and z, so matrices would be about </span></div>
<div><span class="Apple-style-span" style="color: rgb(34, 34, 34); font-family: arial, sans-serif; font-size: 13px; background-color: rgba(255, 255, 255, 0.917969); ">> 100,000 x 100,000, but the grid size could vary. </span> </div>
<div><br></div><div>I would look into a modification of approach 1, where, if you have p processes and an nx x ny x nz grid, you use p/nz processes for each xy plane solve simultaneously, and then to perform the AllToAll communication to rearrange the solutions. I think that you are overestimating the cost of the AllToAll communication.</div>
<div><br></div><div>Jack</div></div>