2
$\begingroup$

After realizing that Gauss-Seidel is extremely slow for my simulation, i wanted to try GMRES and luckily found the C++ code here without diving into the theory. The size of the matrix in my case is nxn where n=50,000. However, I even cannot initialize the matrix to zero because the computer freezes while initializing it. Is there a memory friendly version of GMRES?

$\endgroup$
  • $\begingroup$ The example you linked is for dense matrices, which it is no surprise that your matrix takes too much RAM. You are going to want to hunt for a sparse implementation, which should be fairly common (and better performing) than the dense version. $\endgroup$ – Godric Seer May 17 '14 at 23:46
  • $\begingroup$ FYI, initializing a 50000x50000 double precision matrix requires about 18.6 GB of memory. $\endgroup$ – Doug Lipinski May 18 '14 at 16:13
  • $\begingroup$ If your implementation uses a full matrix (as opposed to a sparse matrix) and you only have bad preconditioners, why not use Gauß elimination? $\endgroup$ – Guido Kanschat May 18 '14 at 18:04
2
$\begingroup$

From what I can see, on a normal machine, you would run out of memory if you are using a Dense matrix. The implementation you indicated in the question uses a dense matrix ( storing all the matrix entries ). Generally in most simulations, each element computes fluxes with its adjacent elements and the resulting matrix structure is sparse ( the number of non-zero entries in the matrix scales only linearly with the size ). So, one way to get around this problem is to use a sparse matrix implementation of the GMRES algorithm. There are several open-source implementations available. I would recommend PETSc - http://www.mcs.anl.gov/petsc/ PETSc provides implementations for the data-structures and algorithms you need and above that, you get parallelism for free.

If you prefer a stand-alone code, I can point you to this reference - http://people.sc.fsu.edu/~jburkardt/cpp_src/mgmres/mgmres.html

If you (for some reason) prefer to stick to a dense matrix implementation, you should go for distributed memory (MPI) parallelism. Once again, PETSc provides you with algorithms for the same.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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