How to solve a sparse matrix efficiently?

Hallo everyone,

I am trying to solve a linear equation system A*x=b that describes a 3D grid and heat transfer between adjecent cells and mass convection by an upwind scheme. As the upwind scheme is somewhat asymmetric, the sparse Matrix that looks like the one below. So far I use x = A\b, but when the problem starts to get bigger this appears to be rather slow. Therefore I am looking for a way to speed up the computation. As the thermal conductivity and so are temperature dependent, I have to solve this equation system a couple of times till a global convergence is achieved.

There are some iterative methods implemented (<http://www.mathworks.de/de/help/matlab/linear-equations-iterative-methods.html>) but I am not sure, which one to chose.

I guess, that lsqr(A,b,tol,maxit,M1,M2,x0) might be a good choice (just maybe ;-) ) but now I am wondering how to chose an appropriate preconditioning M1 and M2. As mentioned before, I have to solve this problem a couple of times, therefore I think that I should be able to provide an good guess on x0 the second time a try to solve this problem - but I am not really sure how to do this... Plus: when I use lsqr or bicq without any changes, the results looks like crap ;-)

You might have figured, that I am more of an engineer than a mathematician, therefore I'm not really sure what to do.

It would be so great, if you could give me decent advice :-)

Best regards Nils

[Patrik Ek]

Have you tried singular value decomposition (svd) or matlabs function pinv? svd requires a little coding, but pinv and methods using svd is specialized on sparse matrices.

[John D'Errico]

Patrik - pinv would be terribly slow compared to backslash.

[Nils]

No, I did not. Can you give me a hint on how svd works? I looked at http://www.mathworks.de/de/help/matlab/ref/svd.html But I'm not sure, what I should do with it.

[Patrik Ek]

You should try wikipedia instead, matlabs function svd, does a singular value decomposition. It does not invert matrices. However, this can be used as a tool for matrix inversion. I know that this is nothing automatical or even simple, but if matlabs functions does not work, this may be necessary. I am not sure it works, but it is worth a try. The thing is that one application of svd is that its relation to pseudo inverses makes it possible to split up the matrix and calculate the inverse of all the non-zero parts.

@John D'Errico, are you sure pinv is still slower if the matrix is sparse?

[Nils]

I think that pinv does now work for sparse matrixes, at least in Matlab 2013.

Answer 1

[Paul]

Did you try defining A and b as sparse? So:

A=sparse(A);
b=sparse(b);
x=A\b;

Answer 2

[Nils]

Hi,

yes I dis this. Really nice speed up :-) But it still starts to get slow as the problem becomes bigger.

[Paul]

See http://www.mathworks.nl/help/matlab/math/sparse-matrix-operations.html#f6-9169 . If your matrix A is symmetric I guess symmlq should be faster.

[Nils]

As the upwind scheme is not symmetric, the matrix is not symmetric. Do you think that a symmetric Matrix would be much faster to solve? What might be the speedup? 10%? 20%?

Thx!

[Paul]

I should've seen it's not symmetric from the figure, my bad. I don't know how much faster it would be. But what's the result if you use lsqr, so:

x = lsqr(A,b,TOL,MAXIT,[],[],x0)

where TOL and MAXIT are the tolerance, and maximum number of iterations respectively. You don't have to specify them (put [] there in that case like I did for M1 and M2). x0 is your initial guess. You can use x as the x0 of the next time by doing:

for i=1:times_you_needtosolve
    blabla; % put your statements that define A and b here
    x{i} = lsqr(A,b,TOL,MAXIT,[],[],x0)
    x0 = x{i};
end

x{i} corresponds to the solution for the ith problem. You can also put all solutions in a matrix by replacing it with x(:,i) and doing the same for the x0 = ... line.

[Nils]

Thank you again for your help! I tried using this method, but although the I reduced the tolerance to 10e-7 and increased iterations to 800 the results looks rather poor. It does not seem to be the right solution. The deviation to the mldevide operator is significant... Seems odd to me...

[Paul]

Well, a solution may not exist. Are you sure you need to solve x=A\b, because afaik an upwind scheme for solving differential equations means that you use values from the previous time step and at the grid points around the point you consider to calculate the new value. What you basically do is :

x_new = A * x_old

Where A is the matrix which defines the points you use for the calculation at each point. So you really don't need to find the inverse since this is just a simple matrix multiplication.

[Nils]

I try to solve a steady temperature distribution. I consider the following phenomena:

• Thermal conduction (central differences)
• Enthalpy flow by mass flux (upwind scheme)
• Boundary conditions
• 3D geometry

As is solve for a steady state problem, I put all the dependencies between the cells into a Matrix and solve for the equilibrium temperature distribution. Of course now the thermal conductivity and so on changes, that is why I have to repeat this step a couple of times. When I solve this problem using mldevide the solution looks physically correct and, by the way, is identical to the solution that I can compute using ANSYS FLUENT. Therefore I think that a solution exists. Sadly the convection term in the equation is rather big compared to the thermal conduction and also the heat capacity behaves quite rough. Therefore I read that an upwind scheme (I use fourth order) should be appropriate. I should be able to use a symmetric central differencing scheme for the convection in large parts of the domain, but I certainly can't use a symmetric scheme on the borders leading to an asymmetric matrix anyway.