LinearSolve and Krylov Method options

Question

I need to solve a large, sparse, linear system. At some point Method -> "Multifrontal" fails and a bit later also "Pardiso" method fails, for insufficient memory I suppose: indeed the Kernel crashes with no messages (maybe there is a better way to handle this situation).

So I want to try with an iterative method. My linear system is naturally not well scaled and not well conditioned so I'll probably need to play with method options ("ResidualNormFunction", "Tolerance", "MaxIterations", "Preconditioner"...) to get acceptable results within a resonable number of iterations.

With few words, my question can be stated as: supposing for example the solution of my linear system represents the displacements (in meters) of some parts and nodes of a mechanical system, how should I set method options if I consider acceptable an error of, says, 1mm (= 0.001 meters)?

I have to admit I'm weak at Krylov iterative methods for linear system. Trying to answer this question I studied the various options and the effects of various setting. Some result contraddict my understanding of the algorithm and/or my expectations. What I did, what I discovered, what surprised me follows.

---

First of all, I noticed a documentation error (I suppose) regarding a sub-method of the Krylov method:

"BiCGStab" iterative method for Hermitian matrices

1. The sub-method, to be recognized by LinearSolve, shoud be spelled "BiCGSTAB".
2. According for example to Wikipedia the sub-method is suitable for nonsymmetric matrices. Indeed, as apparently proved by subsequent experiments, it's the default sub-method.

In numerical linear algebra, the biconjugate gradient stabilized method, often abbreviated as BiCGSTAB, is an iterative method developed by H. A. van der Vorst for the numerical solution of nonsymmetric linear systems.

This code (and other results) shoud suggest that "BiCGSTAB" is the default sub-method, suitable for non-hermitian matrices:

Table[method -> Norm[Table[Module[{
      A = RandomReal[{0, 10.^11}, {20, 20}],
      b = RandomReal[{0, 10.^11}, {20}]},
     Norm[
      LinearSolve[A, b, Method -> "Krylov"] - 
       LinearSolve[A, b, 
        Method -> {"Krylov", "Method" -> method}], \[Infinity]]
     ], {k, 100}], \[Infinity]], {method, {"GMRES", "BiCGSTAB", 
   "ConjugateGradient"}}]

{"GMRES" -> 8.85677, "BiCGSTAB" -> 0., "ConjugateGradient" -> 5.51835*10^17}

With a similar technique I discovered that the default preconditioner is probably "ILUT":

Table[preconditioner -> Norm[Table[Module[{
      A = RandomReal[{0, 10.^11}, {20, 20}],
      b = RandomReal[{0, 10.^11}, {20}]},
     Norm[
      LinearSolve[A, b, Method -> "Krylov"] - 
       LinearSolve[A, b, 
        Method -> {"Krylov", 
          "Preconditioner" -> preconditioner}], \[Infinity]]
     ], {k, 100}], \[Infinity]], {preconditioner, {None, "ILU0", 
   "ILUT", "ILUTP"}}]

{None -> 30.8518, "ILU0" -> 6.6984*10^-6, "ILUT" -> 0., "ILUTP" -> 3.07847*10^-6}

To eventually diagnose convergence problems, it is useful to check the residual. Unfortunately there is not a direct way to do this. An idea to do this, a posteriori, can be to take advantage of "ResidualNormFunction" and Sow. It important to know that the default "ResidualNormFunction" it's probably the 2-Norm.

Table[norm -> Norm[Table[Module[{
      A = RandomReal[{0, 10.^11}, {20, 20}],
      b = RandomReal[{0, 10.^11}, {20}]},
     Norm[
      LinearSolve[A, b, Method -> "Krylov"] - 
       LinearSolve[A, b, 
        Method -> {"Krylov", 
          "ResidualNormFunction" -> (Norm[#, norm] &)}], \[Infinity]]
     ], {k, 100}], \[Infinity]], {norm, {1, 2, 3, 4, 5, \[Infinity]}}]

{1 -> 1.05201*10^-7, 2 -> 0., 3 -> 8.79819*10^-8, 4 -> 1.72088*10^-7, 5 -> 1.07014*10^-6, [Infinity] -> 1.33178*10^-6}

Now setting

A = RandomReal[{0, 10.^11}, {20, 20}];
b = RandomReal[{0, 10.^11}, {20}];

and

{sol, {rlist, rnlist}} = 
  Reap@LinearSolve[A, b, Method -> {"Krylov", "Tolerance" -> 0.01,
      "ResidualNormFunction" -> ((Sow[#, "v"]; Sow[Norm[#, 2]]) &)}];

The first vector submitted to "ResidualNormFunction" is always

rlist[[1]]

{1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}

The second is exactly the residual corresponding to a zero starting vector

A.Table[0, {20}] - b + rlist[[2]]

{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}

But the last is not exactly the returned solution:

A.sol - b + rlist[[-1]] // Norm

0.00431624

The last residuals norm are

Take[rnlist, -10]

{48053.7, 48117.4, 45135.4, 34865.3, 31010.1, 30925.1, 26778.7, 10431.8, 9254.33, 111.13}

Now if I set a Tolerance

{sol, {rlist, rnlist}} = 
  Reap@LinearSolve[A, b, Method -> {"Krylov", "Tolerance" -> 0.1,
      "ResidualNormFunction" -> ((Sow[#, "v"]; Sow[Norm[#, 2]]) &)}];
Take[rnlist, -10]

{8.39886*10^10, 8.65597*10^10, 7.88562*10^10, 4.39645*10^12, 1.56643*10^12, 4.2778*10^10, 3.357*10^10, 1.71505*10^11, 8.4448*10^10, 2.03388*10^10}

so the residual norm is not directly compared to the tolerance as a criteria to stop iterations. Moreover

{sol, {rlist, rnlist}} = 
  Reap@LinearSolve[A, b, Method -> {"Krylov", "MaxIterations" -> 30,
      "ResidualNormFunction" -> ((Sow[#, "v"]; Sow[Norm[#, 2]]) &)}];
Length@rnlist

62

So the ResidualNormFunction is called twice per iteration?

Can someone explain all this behaviors and suggest proper interpretation and settings for these options?