46

Your question is a bit like asking for which screwdriver to choose depending on the drive (slot, Phillips, Torx, ...): Besides there being too many, the choice also depends on whether you want to just tighten one screw or assemble a whole set of library shelves. Nevertheless, in partial answer to your question, here are some of the issues you should keep in ...


33

The important thing when choosing iterative solvers is the spectrum of the operator, see this paper. However, there are so many negative results, see this paper where no iterative solver wins for all problems and this paper in which they prove they can get any convergence curve for GMRES for any spectrum. Thus, it seems impossible to predict the behavior of ...


29

There are many more out there, all with different goals and views of the problems. It really depends on what you are trying to solve. Here is an incomplete list of packages out there. Feel free to add more details. Large Distributed Iterative Solver Packages PETSc — packages focused around Krylov subspace methods and easy switching between linear ...


19

Introduce the vector $y:=-A^{-1}Gx$ and solve the large coupled system $Ay+Gx=0$, $G^Ty=-b$ for $(y,x)$ simultaneously, using an iterative method. If $A$ is symmetric (as seems likely though you don't state it explicitly) then the system is symmetric (but indefinite, though quasidefinite if $A$ is positive definite), which might help you to choose an ...


18

The choice between direct and iterative methods is dependent on goals and problem at hand. For Direct methods, we can note: The coefficient matrix of the linear system changes over the course of computation and may for sparse systems exhaust memory requirements and increase work effort due to fill-in Must complete to give useful results Factorization can ...


16

Both of them are direct solver to solve linear systems (opposing to iterative solver). mldivide does perform the tests for $A$ in solving $Ax = b$. Please see Allan's answer in this thread for more information. Also see MATLAB's help on mldivide algorithm here. mldivide for square matrices: If A is symmetric and has real, positive diagonal elements, ...


15

Have a look at Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods (Barrett et al.). You can find it here. Here's why I'm recommending this over other references: the "flowchart of iterative methods" in appendix D (last page) covers both linear solvers and preconditioners, it is short (100 pages or so), does not go into too ...


14

There are two choices which, if you use reasonable data structures, can solve the problem with $n>10^6$ on a laptop and $n \approx 10^{12}$ on a supercomputer. Note that for efficiency, you should use multigrid to solve with $\Delta$. The cost in either case will be a small factor more expensive than just solving with $\Delta$. The two approaches are ...


14

The most obvious thing you can do is to precompute [L,U] = lu(A) ~ O(n^3) Then you just compute x = U \ (L \ b) ~ O(2 n^2) This would reduce the cost enormously and make it faster. Accuracy would be the same.


14

Interesting that this question came yesterday, since I just finished an implementation yesterday that does this. My Background Just to start of, let me know that while my education background is from scientific computing, all work I have done since graduating, including my current Ph.D. work, has been in computational electromagnetics. So, I guess our ...


14

You can't beat an explicit formula. You can write down the formulas for the solution $x=A^{-1}b$ on a piece of paper. Let the compiler optimize things for you. Any other method will almost inevitably have if statements or for loops (e.g., for iterative methods) that will make your code slower than any straight line code.


13

When is a matrix ill conditioned? It depends on the accuracy of the solution you are looking for, as much as "beauty is in the eye of the beholder"... May be your question should better rephrased as are there cheap and robust condition number estimators based on the $LU$ factorization? Assuming you are interested in the real general (dense, non symmetric) ...


13

When you use ZGELSS to sovle this problem, you're using the truncated singular value decomposition to regularize this extremely ill-conditioned problem. it's important to understand that this library routine is not attempting to find a least squares solution to $Ax=b$, but rather it is attempting to balance finding a solution that minimizes $\| x \|$ ...


13

Ill conditioning is a feature of the system of the equations, not of the algorithm used to solve the system of equations. If your systems are that badly conditioned ($10^{15}$), then you can expect the solution to the system to be extremely sensitive to any perturbation of the problem data, even if the solution is done in extremely high precision (e.g. 500 ...


12

Your matrix $A$ isn't a circulant matrix- it's just Toeplitz. The method that you're trying to use fundamentally only works for circulant systems. Furthermore, your $a$ vector doesn't have the "-1" in it anywhere, so you clearly don't have sufficient information. A method that involves embedding the $n$ by $n$ Toeplitz matrix in a double-sized ...


12

I believe comparing an iterative method (multigrid) to a direct/exact method (Thomas) in terms of exact operation count isn't really meaningful. IIRC, Thomas operation count is $8N$ for any tridiagonal system. The only time I can imagine multigrid conceivably beating that is for a trivial case of having a linear solution, and even then the cost of evaluating ...


11

In general, there is no shortcut other than completely re-factoring the matrix. There have been a few similar questions on this SE that cover the topic in more depth than I can: Can diagonal plus fixed symmetric linear systems be solved in quadratic time after precomputation? LU Decom of PSD Matrix + Diagonal Matrix Perturbation of Cholesky decomposition ...


11

The short answer is that the Thomas algorithm will be faster than any iterative scheme for almost all cases. The exception would perhaps be applying a single iteration of a very simple iterative scheme such as Gauss-Seidel, but this is highly unlikely to give an acceptable solution. Also, this is ignoring parallel processing concerns. Multigrid is an ...


11

Practical experience shows that trying to get good initial iterates has little value. For example, in the context of solving partial differential equations, if you take the solution from one mesh, interpolate it onto a finer mesh, and use that as the starting guess for something like a CG iteration to solve the same problem on the finer mesh, then it turns ...


11

Defining the auxiliary variable $y=Bx$ yields the following algebraically equivalent expanded system, $$\underbrace{\begin{bmatrix} 0 & A \\ B & -I \end{bmatrix}}_{K} \underbrace{\begin{bmatrix} x \\ y \end{bmatrix}}_{u} = \underbrace{\begin{bmatrix} b \\ 0 \end{bmatrix}}_{f},$$ which you could solve with GMRES or another nonsymmetric Krylov method. ...


10

All direct solvers supported by PETSc are available in Python under a common interface via petsc4py. Supported sparse direct solver packages include the PETSc native direct solvers, MUMPS, PasTiX, SuperLU, SuperLU_DIST, Umfpack, CHOLMOD, Spooles, LUSOL, MATLAB, and ESSL. See the MATSOLVER* man pages here.


10

Without taking sides the discussion about whether to use direct or iterative solvers, I just want to add two points: There exist Krylov methods for systems with multiple right-hand sides (called block Krylov methods). As an added bonus, these often have faster convergence than standard Krylov methods since the Krylov space is built from a larger collection ...


10

Trivial answer for square $A$: use dgesvx which solves also for $A^T x = b$ when TRANS = 'T'. Please note that with BLAS or LAPACK you hardly have to transpose (swapping elements in memory) a matrix: most of the subroutines have a TRANS argument to accommodate for operation on the transpose matrix or on a matrix stored with a different memory layout. (...


10

This is called "structurally symmetric". It simplifies graph traversal, such as occurs when setting up aggregates in algebraic multigrid, but doesn't offer much structure to improve convergence rates. Note that all common PDE discretizations have this property so this is still a huge class of matrices including many instances for which no truly good ...


10

In some cases, (F)MG provides an algorithm with optimal properties. For instance, properly tuned FMG can solve some elliptic problems in a small number of "work units", where a work unit is defined to be the computational effort required to express the problem itself - in this case the operations to form the residual $b-Ax$ on the finest grid. This is such ...


9

MKL does not do distributed parallelism (e.g. MPI), and the support for sparse solvers is rudimentary, definitely not at the level of the other two. Currently, there is only one meaningful benchmark: scalable performance of Sparse Matrix-Vector product (SpMV). Since this is memory bandwidth limited, you can only screw it up. Both PETSc and Trilinos do fine ...


9

A preconditioner, say M, is an approximation on the system matrix, say A that changes the problem into another problem with improved eigenvalue spectrum. A perfect preconditioner would be inverse of A i.e inv(M) = A. Unfortunately, this inverse is normally not avaiable, too complicated to compute, requires more space to store because of the fill-in's ...


9

There is typically a trade-off between the amount of work you put into constructing a good preconditioner for an iterative solver and the work you save by using a good preconditioner when actually solving the linear systems. In your case, the case is pretty clear: put as much work as you can into constructing a good preconditioner because you have to solve ...


9

The MUMPS sparse direct solver can handle symmetric indefinite systems and is freely available (http://graal.ens-lyon.fr/MUMPS/). Ian Duff was one of the authors of both MUMPS and MA57 so the algorithms have many similarities. MUMPS was designed for distributed-memory parallel computers but it also works well on single-processor machines. If you link it ...


9

Jed Brown has already pointed this out in the comments to the question, but there is really not very much you can do in usual double precision if your condition number is large: in most cases, you will likely not get a single digit of accuracy in your solution and, worse, you can't even tell because you can't accurately evaluate the residual corresponding to ...


Only top voted, non community-wiki answers of a minimum length are eligible