Which method should I choose to solve a large (~20 000 variables) dense symmetric positive-definite, possibly ill-conditioned, system of linear equations? The system will be solved for two vectors. I'm interested both in performance and numerical accuracy, so I'm looking for possible solutions to choose the best compromise. I want to use it as a part of a special purpose SDP solver under OCTAVE/MATLAB, so most probably it will be implemented as a mex file, with a single thread.
EDIT: Thank you for your responses. I mean an optimization solver, with primal-dual interior point method. Dual solution is always feasible in my case, and dense (I mean the dual is not block-diagonal).
EDIT2: The reason why I'm writing an SDP solver is that I want to try to get a more efficient tool for solving large NPA problems (http://arxiv.org/abs/0803.4290). This case allows for some simplifications, the interior is always non-empty, I can find a good initial dual feasible solution (and thus keep this feasibility), and the structure allows for easier creation of the equation for dy. The problem is that this equation is large for these problems.
I tried to use SeDuMi and SDPT3, with the latter being able to run much larger programs. For large I got with SDPT3 the error saying "Schur complement matrix not positive definite" meaning that Cholesky factorization failed.
EDIT3: I observe that with this warm start it takes about 3-4 iterations less to solve the program. The problem is the failure of Cholesky for the large Schur complement matrix. In this case it is fully dense, but near the last iteration it becomes very bad conditioned. This is why I'm looking for a method which will be able to solve it.
EDIT4: Thank you for your response. I'll try pivoting with Cholesky. Do you recomend any package which works under OCTAVE for this task?