Hi All.I am trying to solve a PDE using an iterative method, and I am trying to store the 3D space into a sparse matrix. I am having trouble understanding one concept though. The 3D space I am working on can have several different (or all the same) boundary conditions at each of the 6 faces of the space. For example, Z-min might have Dirichlet, while all others have Nuemann.In my 2D representation this would mean the first nx*ny values of the matrix would want to be set to a constant, whereas the other values would want to have the 7 values of the stencil (thus 7 columns for that row) factored in the calculation.So - how do I pick a sparse matrix storage format that can take into account that my matrix *may* be completely diagonally dominant, partially, or some other combination.Thanks for any insight!Derek

If I have understood correct, you are solving a 3d elliptic PDE some of whose boundary conditions are Neumann?Since the FD scheme holds at mesh pt j = 1 and we have NBC at j = 0 we can eliminate the value of U at j = 0 and then just solve the FD scheme using the usual iterative schemes?My feeling is that we don't need a sparse matrix as such but a multiarray data structire is good enough. I think not all methods will converge especially if the eigenvaues are near 1 in absolute value.

you could try a hybrid format eg., ELL/COO, I'm doing this with reservoir simulation code. Check out bell/garland's paper even if you're not doing this in GPU, they discuss optimal kernels/formats, or anything written by Schenk. He wrote pardiso and imho is a leading sparse guy, as well as a genuinely nice chap. I also blog a little on my reservoir experiences, I've been moving away from CSR to drag more performance out of the solver(s)

The boundary conditions could be a mix of Neumann, Reflective, Dirichlet. I guess that is one way of doing it though, assume they are all Neumann, do the solve as if it is symettric, then go back in and clear the values that are on the boundary, and then recompute those to be dirichlet or whatever.I think I know what ELL is - however, I couldnt figure out what I would need to use for padding to not change the answer. Any insight on how to do a matrix vector multiplication when you don't have a constant number of non-zero entries per row? (meaning for example dirichlet would have zero non-zeros for those rows, right?)

btw - to be clear, i dont care about space usage for this problem, just solving the cache-locality problem

I was suggesting a simple decomposition like (ELL+COO) * v, where COO is changing and ELL is constant. Is this a moving boundary problem?

It is not a moving boundary problem. It is just a 7 point stencil problem. Each of the boundary conditions is set at runtime for xmin, xmax, ymin, ymax, zmin, zmaxIf I were going to try to force this into a scheme, should I just go with COO only? That paper you cited was saying that the hybrid would require you knowing a max for ELL, but I am not sure what that would be. With no boundary conditions, it would of course be 7, but I can actually end up with more than that in a row.

if you could post a pic of your matrix using spy, that might help?but basically, speaking from the reservoir experience, I have a list of nodes from a very large cartesian grid which are also of the 7 point stencil type. Throughout the calculation, at least once, nodes are trimmed and connections lost due to (for instance) computationally trapped cells. So I have matrix which changes in structure I guess like yours. However the vast majority of nodes/stencils don't and this I keep in ELL, promoting locality; the much smaller variable non-zeros are in COO

I will see what I can do - I dont have matlab on this system at work, so it would be hard to generate a picture (unless you know of an online tool?).I think i see what you are saying though - when possible, store in ELL, when not, store in COO.This is kind of what I envisioned as well, I just have to figure out how to determine if a row in the matrix is one of the ELL rows, or one of the COO rows. It wont be the same for every run in my case. This is for a PIC algorithm, and depending on the scenario being run, the boundary conditions will change. Once they are set for the run, it will be the same matrix every time (the values change of course, but the index values wont change).

I think I just figured something out. The only real boundary condition that messes me up is Dirichlet. Setting the boundary to a constant is the only one that I am not sure how to represent in the sparse matrix. I suppose I could assume that ALL rows in the matrix would have 7 non zero values, and then go back after the calculation and set the values in the solution vector that correspond to the rows of the boundary condition back to the constant. Does this sound reasonable?

so you would solve A.x = b and then reset some of the x? is this physical? sorry I don't know much about PIC, my colleague does the EM around hereI recommend octave, the free version of matlab, it has much of the same functionality including a nice little selection of sparse routines eg., spy

Be reminded that it is allowed to actually store zero values as non-zeros in a sparse matrix!For the Dirichlet nodes you can do either of the following:1. Eliminate Dirichlet values from the linear system. I.e. They are not a degree of freedom anymore, but you need to change the values of the right hand side associated with the neighbouring nodes accordingly.2. Represent them as constants (Diagonal value is one, and all offdiagonal values are zero). BTW that does not mean that your row only has one non-zero column stored, but all non-zero off-diagonal columns actually store zero! Please note that in this case the matrix is not symmetric any more, which might cause problems for some solvers.The linear system for the 1D Laplacian using approach 2. looks like 1 0 0 0 0 0-1 2 -1 0 0 0 0 -1 2 -1 0 0 0 0 -1 2 -1 0 0 0 0 -1 2 -1 0 0 0 0 0 1Note that Dirichlet entries will stay constant (and have the value of the right hand side). There is no need to reset anything after the discretization.

I think I have that part figured out. I am now onto a different problem with this solver, related to option 2(treat dirichlet as constants).In the solver I am working on, the sparse matrices are stored in coordinate form. Basically for every row there is a set of coordinates corresponding to the column in the matrix. During the creation of this sparse representation, things are sorted so that for row n, the column access in J is always ascending. That is, for node 441, the J columns corresponding might be 0, 420, 440, 441, 442, 461, 880 or something like that. For some b.c., the case might be that the absolute index of the Z-positive neighbor node might be less than the Y-positive neighbor node value. That is, the column coordinate might *really* be like this : 0, 420, 440, 441, 442, 461, 450but when storing into the sparse format, there is a sort going on, so that for the node, the column cooridnates would be stored0, 420, 440, 441, 442, 450, 461Why would this need to be the case? Why does it matter in the rest of the PCG solver. Seems to me the matrix-vector multiply, and subtraction shouldnt be affected. Yet when I take out the sort, so that the columns are stored the other way - the solver doesnt converge and I get a lof of NaN output.

very good question. your solve(r) depends critically on your matrix structure or lack thereof, the sort obviously changes this. Look around a little, but Conjugate Gradient(?) while easy isn't terribly general and you might want instead gmres with a suitable preconditioner; your eigenspectrum will dictate the exact choice. CG is only good for positive definite I think.I also recommend direct methods like gauss elimination based solvers, unless your non-zeros are > 1e7 in which case memory quickly becomes an issue due to fill-in rates during preprocessing. Direct methods are non-iterative and thus less sensitive to your matrix type.

You said you are using PCG (preconditioned conjugate gradient). Therefore the preconditioner might want to use the diagonal entries. Due to the storage format I would suspect that referencing the diagonal triggers a binary search. This fails if the row is not sorted in the right way.But you should probably ask the support for your solver.