Thanks for the replies.<div>Sorry, I was mistaken about the dt/Re coefficient. I use an implicit time step on the linear part of incompressible Navier-Stokes, so roughly discretizing du/dt = 1/Re * Lu with implicit Euler for simplicity gives:</div>
<div><br></div><div>(u^{n+1} - u^n) / dt = 1/Re L u^{n+1}</div><div>rearranging:</div><div>(I - dt/Re L) u^{n+1} = u^n<br><br>I incorrectly had dt/Re on I instead of L before.</div><div><br></div><div>For tri-diagonal matrices, I believe the direct solve is O(N), and similar methods exist for block tri-diagonal matrices. I know multigrid is also O(N), but I've heard it tends to be slower.</div>
<div><br></div><div>I think the bandwidth (b) of the matrices is small. For example for a 5-point stencil (center, left, right, below, above) the bandwidth is 5 and for a 9-point stencil it is 9.</div><div><br></div><div>
The other limiting factor of doing option 1, with the all-to-all, is the number of nodes is limited by nz. That might be ok, but I would like the code to scale to larger grids in the future. </div><div><br></div><div>Thanks,</div>
<div>Brandt</div><div><br><br><div class="gmail_quote">On Tue, Nov 15, 2011 at 1:12 PM, Jack Poulson <span dir="ltr"><<a href="mailto:jack.poulson@gmail.com">jack.poulson@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Brandt,<br><br><div class="gmail_quote"><div class="im">On Tue, Nov 15, 2011 at 11:12 AM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">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>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><div class="im"><div style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px">
> 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">
> 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">
> solved directly and cheaply on one node. </div><div style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px"><br></div></div><div style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px">
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 class="im">
<div style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px"><br></div><div style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px">
> The two options I'm comparing are:</div><div style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px">> 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">> 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">> 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">> 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">> 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">> points in the z direction.</div><div style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px">
> </div><div style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px">> 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">> 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">> 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">> 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">
> </div><div><span style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px">> 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 style="color:rgb(34, 34, 34);font-family:arial, sans-serif;font-size:13px">> 100,000 x 100,000, but the grid size could vary. </span> </div>
<div><br></div></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><font color="#888888"><div>Jack</div></font></div>
</blockquote></div><br></div>