15 votes

"Cookbook" about iterative linear solvers and preconditioners

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 ...
GoHokies's user avatar
  • 2,176
12 votes
Accepted

role of initial guess for iterative linear solver

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, ...
Wolfgang Bangerth's user avatar
12 votes
Accepted

Solving linear system of the form $ABx=b$

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{...
Nick Alger's user avatar
  • 3,143
12 votes

Iteration counts of AMG solver changes in parallel

This is something that can happen in almost any numerical algorithm running in parallel. It's important to know that floating-point addition is not associative due to round-off errors. Thus you can't ...
Brian Borchers's user avatar
11 votes
Accepted

What happens when I use a conjugate gradient solver with a symmetric positive semi-definite matrix?

You first have to make sure that your system is solvable: this happens iff the right-hand side $b$ is orthogonal to the kernel of $A$. If $A$ has a dimension-1 kernel spanned by $v$, you need to have $...
Federico Poloni's user avatar
11 votes

Is using iterative methods to solve a linear system always superior to inversing the matrix?

First off, there are basically no scenarios where one would ever actually compute and store $A^{-1}$ in memory, even for small problems. An LU factorization offers both superior efficiency and ...
whpowell96's user avatar
  • 2,054
10 votes
Accepted

How can a CG solver solve a non positive definite sparse matrix

I highly recommend the following read: J.R. Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" In short, if the matrix is non-positive definite, there is no ...
Anton Menshov's user avatar
  • 8,602
10 votes

Which C++ linear algebra library is probably the fastest on solving huge sparse [square matrix] linear system?

Eigen 3 is a nice C++ template library some of whose routines are parallelized. c.f. Eigen documentation The parallelization is OMP only, so if you intend to parallelise using MPI (and OMP) it is ...
Nox's user avatar
  • 341
10 votes
Accepted

Is there an iterative solver for dense matrices with possible zero diagonal entries?

Iterative Krylov-subspace solvers generally only require matrix-vector products and don't care whether or where there are zeros in the matrix. In your case, unless you have other information about the ...
Wolfgang Bangerth's user avatar
9 votes
Accepted

numerical solution of an under-determined linear equation in high dimensions

You want to minimize $\min \| Ax -y \|_{2}^{2} + x^{T}B^{T}Bx=\| Ax -y \|_{2}^{2} + \| Bx \|_{2}^{2}$ Recall that $\| u \|_{2}^{2} + \| v \|_{2}^{2}= \left\| \left[ \begin{array}{c} u \\ v \end{...
Brian Borchers's user avatar
8 votes
Accepted

Finding the matrix inverse given a solver for the matrix equation $Ax=b$

Your two ideas make it much too complicated. If $X$ is the inverse of $A$, $$ AX=I, $$ and $x_i$ is the $i$-th column of $X$ and $e_i$ is the $i$-th column of the identity matrix $I$ ($e_i$ is a ...
Kirill's user avatar
  • 11.4k
8 votes
Accepted

Writing a non-square linear system in standard form $A\cdot{x}=b$

The matrix math in this paper is terribly hard to follow, but here's what equation (19) should look like $\left( \begin{array}{c|c|c|c} c_{11} \mathbf{I} - \mathbf{A}^T & c_{12} \mathbf{I} &...
LedHead's user avatar
  • 1,223
8 votes
Accepted

Solving Poisson equation with current BC using FEM

It is standard procedure for general case like that. You can enforce condition for electrode by Lagrange multipliers, $$ L(u,\lambda) = \int_\Omega \sigma u_{,i} u_{,i} \, \textrm{d}V - \int_{\...
likask's user avatar
  • 906
8 votes

Lost on Matrix Inversion

Several points I want to mention (with an encouragement to other CompSci users that are more familiar with Java specifics to give additional, more Java related answers): The solution of a system of ...
Anton Menshov's user avatar
  • 8,602
8 votes
Accepted

Using matrix exponential to solve linear system

You are effectively asking how to compute $y=(\log M )^{-1}b$, where $M=e^A$ is the given matrix. There are several methods for computing $f(M)b$ without forming $f(M)$, and they are reviewed here. ...
Amit Hochman's user avatar
  • 1,081
7 votes

Linear Systems with Multiple Right Hand sides

It's possible that your performance is limited by the memory bandwidth of your system. It's not at all uncommon to have a situation where two or three of your cores can use all of the available ...
Brian Borchers's user avatar
7 votes
Accepted

Method to solve linear, first order ODE of generalized matrix matrix form

There are a few ways you can do this. Even though RK4 fails unless stepsizes are really small, one thing that can happen a lot is that the model can start out more stiff than it really is in the full ...
Chris Rackauckas's user avatar
7 votes
Accepted

Lost on Matrix Inversion

Normally when you invert a sparse matrix the inverse is dense. This imply to have enough memory to store the inverse, in your case the matrix is not so big for nowdays computers. In double precision (...
Mauro Vanzetto's user avatar
7 votes

fastest linear system solve for small square matrices (10x10)

Using an Eigen matrix type where the number of rows and columns is encoded into the type at compile time gives you an edge over LAPACK, where the matrix size is known only at runtime. This extra ...
Daniel Shapero's user avatar
7 votes
Accepted

Sparse matrix inversion

For a matrix that small, you're probably not going to do better than using dense methods. I wrote up a quick test in C++ for an 18x18 matrix with your sparse structure and randomly generated values ...
LedHead's user avatar
  • 1,223
7 votes

Linear solver recommendation(s) for small problems

For problems this small, sparse direct solvers are faster than most iterative solvers even if you include the cost of factorization. As a result, I don't believe that you will be able to find a ...
Wolfgang Bangerth's user avatar
7 votes
Accepted

Getting to know about various BLAS implementations

I'm the primary author of many Julia libraries geared toward "architecture-specific optimizations", including LoopVectorization.jl and Octavian.jl. For BLAS-like operations, one of the most ...
Chris Elrod's user avatar
7 votes

BiCGSTAB convergence

First of all, size 40 is pretty much microscopic for the sort of purposes that BiCGstab was invented. I assure you, this method is great for matrices of sizes in the millions and beyond. Then: even ...
Victor Eijkhout's user avatar
7 votes

Correctness of direct numerical solution of ill-conditioned linear system

One thing you can do to test it is computing the residual $b - A \tilde{x}$ in higher precision. If the residual is small (and this will always happen with an accurate solution), then you can testify ...
Federico Poloni's user avatar
6 votes

Does this partial eigen-expansion have a name?

I don't know for sure whether there is an existing name for this method, but @jessechan's suggestion of "truncated eigenvector expansion" sounds perfectly fine to me (and most people would understand ...
Wolfgang Bangerth's user avatar
6 votes

Does this partial eigen-expansion have a name?

This is something like a truncated SVD or eigenvector expansion of your solution. If you take $$x_m = \sum_{j=1}^m \frac{q_j\cdot b}{\lambda_j}q_j$$ with $m=n$, this is the exact solution to $Ax=b$...
Jesse Chan's user avatar
  • 3,112
6 votes
Accepted

Thomas algorithm for 3D finite difference

The problem is that 2d and 3d discretizations are not block tridiagonal. The only tridiagonal decomposition would be a 2x2 decomposition that encompasses the entire matrix. Of course, this fact is a ...
Wolfgang Bangerth's user avatar
6 votes
Accepted

Fast c++ library to solve very big sparse systems

I second the idea of using Eigen, which is pretty efficient, but also very simple to include. If you need a lot more performance, you could try to use PETSc or Trilinos. They are very powerful ...
BlaB's user avatar
  • 1,147
6 votes
Accepted

Why is the speed of the parts of the LU-decomposition so different?

First, don't forget to also time the LU decomposition in a loop! Otherwise it's not really a fair comparison. If I do that, I get the following timings: ...
Christian Clason's user avatar
6 votes
Accepted

Efficiency of parallel direct linear solver

Your problem is too small. You have to consider that to get good efficiency, each processor has to have enough work to offset the cost of communication. In other words, there is a threshold how many ...
Wolfgang Bangerth's user avatar

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