8

I faintly recall from my early "numerics" lectures that iterative linear solvers for $Ax=b$ often require that when $A$ is decomposed as

$$A=D + M$$

where D is a diagonal matrix and $M$ has zero diagonal, the elements of $D$ should be dominant over the entries in $M$ for iterative solvers to perform well.

What if that's not the case and the entries of $D$ become really small?

Should I use a direct solver then?

To be more specific, the linear system I want to solve involves a matrix $$A(\omega) = D(\omega) + M$$ where the non-diagonal part is constant but the diagonal part depends on a parameter $\omega$ in some non-trivial way. So far, I don't see a way around solving $A(\omega) x = b$ for each $\omega$ anew.

The diagonal entries are of the form $A_{jj} = \omega + z_j + i\eta$ where $z_j$ is some real number that depends on the row we are in whereas $\eta$ is a very small convergence factor and $i$ is the imaginary unit. Could this lead to numerical instabilities when $\omega + z \approx 0$?

EDIT: Well, maybe one more thing about the nature of $A(\omega)$: If one sets $\eta$ to $0$ exactly, then $A(\omega)$ is guaranteed to have poles. This is because ultimately I use this matrix to compute (many-body) Green's functions in frequency domain, and those need a convergence factor $\eta$ to move their poles off the real axis. The sum of absolute values of off-diagonal matrix elements in each row is $10$ at most, but the diagonal will always have some entries whose real part is very close or equal to zero.

Lagerbaer
  • 487
  • 4
  • 9
  • @JackPoulson Stokes is perhaps the canonical example of a PDE saddle point problem. Incompressible Navier-Stokes is sometimes called a generalized saddle point problem because it is non-symmetric, but has the same constraint structure. – Jed Brown Jan 19 '12 at 13:35
  • @JedBrown: Fair enough, it's now obvious that I haven't worked on Stokes! I always think of mixed methods for Darcy when I think of saddle point. – Jack Poulson Jan 19 '12 at 16:00

3 Answers3

4

Although it places some restrictions on which methods will work, lack of diagonal dominance or symmetry is not inherently catastrophic. However, these properties are commonly associated with the more difficult problem of non-local influence and difficulty of coarsening, and many "black-box" solvers will not work. To answer your question in a more substantive way, we would need to know the details of this particular system (physics, discretization, parameter regime).

My practical suggestion is to start with direct solvers and delve into iterative solvers only if necessary. You can spend a career developing robust iterative solvers for a particular difficult problem.

Jed Brown
  • 25,650
  • 3
  • 72
  • 130
3

A diagonally dominant matrix is guaranteed to have either all positive (if the entries of the diagonal are all positive) or all negative (if the entries are all negative) eigenvalues, by Gershgorin's theorem. Most iterative methods only work if the eigenvalues of the iteration matrix are in a particular region of the complex plane, so diagonal dominance ensures that all of the eigenvalues have either a stricly positive or strictly negative real part (or that all the eigenvalues lie within a particular radius of some number).

If your iteration matrix has eigenvalues which lie outside the prescribed region for a particular method, it usually won't work correctly. There are piles of methods out there, but I can't pick a specific method out without knowing more about the spectrum of $A(\omega)$.

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
Dan
  • 3,355
  • 3
  • 21
  • 47
2

How large/sparse is your system? Do you really need to solve this iteratively?

I would suggest trying to solve it in Matlab or Octave using the sparse solver (just initialize the matrix using "sparse" and then use backslash). If that works for you, then use UMFPACK directly, which is what Matlab and Octave use internally for the sparse solves.

Pedro
  • 9,573
  • 1
  • 36
  • 45
  • On the order of $10^5$ to $10^6$ rows with at most 8 entries per row. – Lagerbaer Jan 19 '12 at 17:18
  • 1
    Sparsity does not tell the whole story. If small vertex separators exist, then direct solvers should find them and perform reasonably. If the problem is 2D with that sparsity, then a direct solver will be fine up to $10^6$ degrees of freedom. For 3D problems, that problem size will need quite a lot of memory and time. – Jed Brown Jan 20 '12 at 13:03
  • Thus my suggestion to try it out in Matlab or Octave first -- you'll see if the solver can handle it efficiently or not. The efficiency strongly depends on the structure of the problem so it won't be possible to give any definite recommendations. – Pedro Jan 21 '12 at 14:14