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

Question

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.

Answer 1

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.