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?
2
$\begingroup$
$\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$
$\endgroup$
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.
