Sparse diagonal matrix solver

Question

I want to solve, in MatLab, a linear system (corresponding to a PDE system of two equations written in finite difference scheme). The action of the system matrix (corresponding to one of the diffusive terms of the PDE system) reads, symbolically (u is one of the unknown fields, n is the time step, j is the grid point):

and fully:

The above matrix has to be intended as A, where A*U^n+1 = B is the system. U contains the 'u' and the 'v' (second unknown field of the PDE system) alternatively: U = [u_1,v_1,u_2,v_2,...,u_J,v_J]. So far I have been filling this matrix using spdiags and diag in the following expensive way:

    E=zeros(2*J,1);

    E(1:2:2*J) = 1;
    E(2:2:2*J) = 0;

    Dvec=zeros(2*J,1);

        for i=3:2:2*J-3
                 Dvec(i)=D_11((i+1)/2);    
        end

        for i=4:2:2*J-2
                 Dvec(i)=D_21(i/2);
        end

    A = diag(Dvec)*spdiags([-E,-E,2*E,2*E,-E,-E],[-3,-2,-1,0,1,2],2*J,2*J)/(dx^2);`

and for the solution

[L,U]=lu(A);
 y = L\B; 
 U(:) =U\y; 

where B is the right hand side vector.

This is obviously unreasonably expensive because it needs to build a JxJ matrix, do a JxJ matrix multiplication, etc.

Then comes my question: is there a way to solve the system without passing MatLab a matrix, e.g., by passing the vector Dvec or alternatively directly D_11 and D_22? This would spare me a lot of memory and processing time!

Answer 1

Matlab doesn't store sparse matrices as JxJ arrays but as lists of size O(J). See http://au.mathworks.com/help/matlab/math/constructing-sparse-matrices.html Since you are using the spdiags function to construct A, Matlab should already recognize A as sparse and you should indeed see such a list if you display A in console view.

For a tridiagonal matrix like yours, the L and U matrices should already be sparse.

So you just need to ensure that the \ operator uses the appropriate sparse algorithm according to the rules in http://au.mathworks.com/help/matlab/ref/mldivide.html. It's not clear whether the vector B will already be considered sparse, but you could recast it as a diagonal matrix which should certainly be considered sparse.