Linear equation system: Direct solver works, iterative solver does not

I have to solve for x in b = A*x, where a is sparse. This works fine with Matlab's mldivide: x = A \ b. Since I will have to use an iterative algorithm for very large A, I'm currently testing Matlab's iterative algorithms, bicgstab in particular:

[L, U] = ilu(A, struct('type', 'ilutp', 'droptol', 0.0000001, 'udiag', 0));
[delta2, flag, relres, iter] = bicgstab(A, b, 0.0000001, 2000, L, U);


Basically, this works with a general vector b generated with with randn(...). But it does not work with the vector b from my problem. Not work in this case means: It does one or only a few iterations and returns with flag 1 which says that it did not converge in the allowed number of iterations. The returned vector is then wrong compared to the solution with x = A \b.

The vector b, as well as the matrix A have values over a large magnitude range (about 1e10), matrix A is a general sparse matrix with dimension about 2000x2000 and has condition about 1e14

The only idea which I currently have would be to change physical units of the system behind the matrix to make large values smaller and small values larger.

I know this is very little information about this problem, but I would be glad for any hint, especially how to find out more about the problem.

• What is the value of relres that it returns? A condition number of $10^{14}$ is going to give you lots of problems. – Bill Barth May 8 '14 at 11:53
• @BillBarth: relres is normally equal to 1. So you would definitely recommend to choose units such that the condition number will decrease significantly? – Michael May 8 '14 at 12:20
• I don't think you can change just the units to change the condition number. You have to change the ratio of the biggest to smallest scales in your underlying problem, not just the units. That, typically, means you're solving a different problem from the one you originally set out to solve, and that problem typically has a different answer. Maybe you could say a bit more about where the linear systems comes from? – Bill Barth May 8 '14 at 13:12
• @BillBarth: I opened a new question adressing a specific reason (which I suspect) for the problem: scicomp.stackexchange.com/questions/11581/… – Michael May 9 '14 at 10:56

The problem was in the underlying coupled PDE FEM model. This model includes a potential which is not fixed in any point, so any result potential = phi + C would be a solution, where C is an integration constant. I don't know why, but the direct solver (Matlab's \ using UMFPACK) finds a solution, where the iterative solver gets in massive troubles.

I first checked all the sub matrices on the diagonal - the mapping of each variable in my coupled PDE system to itself. So I found that the matrix for this potential has very poor condition. Second, I found that eliminating any of the rows + corresponding column in that matrix would make condition fine, which is a hint that one of the rows/columns are linearly dependent which tells me, that there's one degree of freedom too much. I suspected that it's the integration constant which proved to be the right guess after fixing one of the points.

• Be wary of Matlab's backslash: it applies quite a few heuristics to select from a variety of direct solvers -- among them a least squares solver. So it might be the case that the direct solver isn't solving the same problem. (You can find out what's going on behind the curtains by setting spparms('spumoni',1).) – Christian Clason May 12 '14 at 11:29
• As mentioned, it tells me that it will use UMFPACK for my equation. – Michael May 12 '14 at 11:52