I am solving a 3D heat transfer equation with variable boundaries (insulated, convective, radiative or free) using a F.D.M. technique. My geometry of choice is a cube. The purpose of my work is to get more familiar with such problems -- out of interest -- and their parallelizations using both distributed (MPI) and shared memory (OpenMP/Pthreads).

In order to make this question as neat as possible, I will avoid tediously long expressions -- as is usually the case in 3D -- and go straight to the point. My stencil is a 7-point kernel. I am solving on my cube by traversing the discretization of the inner points using an SOR approach (i.e. excluding boundary surfaces, vertices and corners). Hence I find myself confused about how one usually solves such equations via Finite Differences.

To set-up a (mass) matrix, instead of the naive point-by-point stencil approach is greatly perplexing. I say this because per each row of the mass matrix I end up having only 7 non-zero elements.

In other words, if my cube is discretized equally in x, y and z directions while having 100 points in each respectively, then my overall mass matrix dimensions would become: $(100x100x100)^2 = 1,000,000-by-1,000,000$ with only 7 nnz elements per row. This is too sparse and inefficient from the perspective of memory and certain algorithm applications.

This makes me wonder how inexperienced I currently am when trying to formulate a solution for such problems. That being said, how does one usually formulate such a huge sparse matrix to tackle such problems ? Surely my naive approach isn't the best nor the only solution out there (duh). Any explanation and/or examples would be very welcome.


1 Answer 1


I am not certain I understand your question but it appears to have to do with storage and manipulation of sparse matrices. If so, then I suggest you take a a look at the lecture notes from a course Professor Saad teaches on sparse matrix computations: http://www-users.cselabs.umn.edu/classes/Spring-2014/csci8314/ As you will see, lecture 2 discusses different ways of representing sparse matrices. Lecture 3 specifically discusses the issue you are concerned about-- finite difference equations. Later lectures discuss solution techniques for linear equations when the matrices are sparse.

  • $\begingroup$ The documentation of SuiteSparse might be helpful here as well. $\endgroup$ Jul 7, 2016 at 14:37
  • $\begingroup$ @Deathbreath Thank you for your post. I like the graphical representations - even though they are pretty abstract. However this is mostly, if not explicitly, a GPU implementation ? very useful, just not for someone that has no NVIDEA card for CUDA support. $\endgroup$ Jul 8, 2016 at 14:55
  • $\begingroup$ What is a gpu implementation? SuiteSparse is not. $\endgroup$ Jul 8, 2016 at 14:59
  • $\begingroup$ Thanks Bill, I will vote your answer as valid. But to recap, from what I understand, in order to solve such a conventional 3D FDM system of equations via iterative solvers (gauss-seidel, SOR, conjugate gradient...) one does not explicitly store such mass matrices (due to extreme sparsity). Instead, it has a special sparse technique for not only setting up but solving as well ? That actually sounds very complicated to implement on a compiled language such as C -- needless to mention parallelization... $\endgroup$ Jul 8, 2016 at 15:02
  • 1
    $\begingroup$ Sparse storage formats do a good job of compactly storing even very sparse matrices. The easiest way to exploit black-box iterative sparse solvers is by first storing the matrix in one of the standard formats that Saad describes. The Cholmod direct sparse solver, that is part of SuiteSparse, does support GPUs but also runs very well on modern CPUs. $\endgroup$ Jul 8, 2016 at 15:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.