I have to solve a system $Ax^{(n)} = b^{(n)}$ many times, $A$ being a sparse (pentadiagonal in most part of its structure), unsymmetric, constant matrix.Currently, I am performing the LU factorization with rows and columns permutations, i.e., $PAQ = LU$, and then I am using this factorization to compute $x^{(k)}$ in the following way:$x^{(n)} = Q*(U \backslash L \backslash (P*b^{(n)}))$I have read recently about iterative methods for sparse systems, but I am not very familiar with them. Could you please tell if there is any iterative method that could perform faster and if it is already implemented in Matlab?Thanks, Manu
Fastest way to solve a sparse unsymmetric system many times
sparse;iterative method;linear system;matrix factorization
Generally speaking, for many right-hand side (RHS) problems, a direct solver is a more feasible solution for several reasons:Major computations are performed during the factorization step (which is done only once for all RHS), and the solution find (for each RHS) is much cheaper.Direct solvers do not suffer from poor conditioning of the matrix or, in other words, you don't have to look for a nice preconditioner (required for iterative solvers) ensuring $\mathcal O(1)$ iterations for every RHS (relatively easy) and every matrix (much harder) resulting from a numerical simulation.NB. Sure, if the matrix is extremely ill-conditioned, the direct solver would also suffer, however, direct solvers extend the dynamic range of problems that can be solved without re-formulating a problem.There are a couple of things that I should mention:Usually, it is ill-advised to use an inverse based on LU-decomposition to find a solution. You might want to use a back-substitution LU subroutine that would be more efficient. From your formula, it is unclear if you are doing it in the right way.Since your matrix is sparse and has a stable structure, I would look into sparse direct solvers as opposed to full direct solver. It should speed up the computations significantly.In your case, I would consider doing partial pivoting ($PA=LU$) only, as opposed to full-pivoting ($PAQ=LU$). That might speed up the computations.Connecting a high-performance LAPACK library to MATLAB is also an option.Regarding sparse iterative solvers:You can try different iterative solvers, including GMRES, BiCGStab, TFQMR, etc. Your task would be to provide a fast sparse matrix-vector product to them. Most of them are already included in MATLAB. And if you are using sparse matrix storage there, using them should not be too hard.Repeating an item from above about advantages of a direct solver: you would have to implement a decent preconditioner for all those methods that would reduce the number of iterations they require to converge. One choice would be using incomplete LU (iLU) as a preconditioner. There are more sophisticated and easier options available.Using iterative solvers is a bit tricky. You have to figure out how well is you system conditioned, what is a reasonable choice for a preconditioner, what is a good tolerance to use with an iterative solver (accuracy of the obtained solution), have a feel what kind of problems (in the subset of problems you are solving using numerical simulations) cause non- or slow-convergence of your iterative algorithm.