Write a Gauss-Seidel preconditioner for ML

I am having trouble writing a Gauss-Seidel preconditioner for ML. It says "out of memory," ML is actually a vector of 44229 entries on the first column so its a (44229,1) size. A Gauss-Siedel preconditioner is the diagonal plus the lower, but here we only have a vector (ML). So now the G-S preconditioner is the values from ML along the diagonal surrounded by zeros. PLEASE HELP!

%Lumped mass matrix

for i=1:3

for j=1:3

% \int_E phi_i * phi_i dE = A/6 and % \int_E phi_i * phi_j dE = A/12

if(i==j)

ML(index(n,j)) = ML(index(n,j))+areas(n)/6;

else

ML(index(n,j)) = ML(index(n,j))+areas(n)/12;

end

end

end

% %%%%%%%%%%%%%%%%%%%%%%%%%%% M PRECONDITIONER

[g h] = size(ML);

D = diag(ML);

% b = zeros(g);

% M = D + b;

ML = D^-1.*ML;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Answer 1

[John D'Errico]

Learn what a sparse matrix is. And how to build them, how to use them.

help sparse

    SPARSE Create sparse matrix.
    S = SPARSE(X) converts a sparse or full matrix to sparse form by squeezing
    out any zero elements.
 
    S = sparse(i,j,s,m,n,nzmax) uses vectors i, j, and s to generate an
    m-by-n sparse matrix such that S(i(k),j(k)) = s(k), with space
    allocated for nzmax nonzeros.  Vectors i, j, and s are all the same
    length.  Any elements of s that are zero are ignored, along with the
    corresponding values of i and j.  Any elements of s that have duplicate
    values of i and j are added together.  The argument s and one of the
    arguments i or j may be scalars, in which case the scalars are expanded
    so that the first three arguments all have the same length.
 
    S = SPARSE(i,j,s,m,n) where nzmax = length(s).
 
    S = SPARSE(i,j,s) where m = max(i) and n = max(j).
 
    S = SPARSE(m,n) abbreviates SPARSE([],[],[],m,n,0).  This
    generates the ultimate sparse matrix, an m-by-n all zero matrix.
 
    For example, this dissects and then reassembles a sparse matrix:
 
               [i,j,s] = find(S);
               [m,n] = size(S);
               S = sparse(i,j,s,m,n);
 
    So does this, if the last row and column have nonzero entries:
 
               [i,j,s] = find(S);
               S = sparse(i,j,s);
 
    All of MATLAB's built-in arithmetic, logical and indexing operations
    can be applied to sparse matrices, or to mixtures of sparse and
    full matrices.  Operations on sparse matrices return sparse matrices
    and operations on full matrices return full matrices.  In most cases,
    operations on mixtures of sparse and full matrices return full
    matrices.  The exceptions include situations where the result of
    a mixed operation is structurally sparse, eg.  A .* S is at least
    as sparse as S.
 
    See also ISSPARSE, SPALLOC, SPONES, SPEYE, SPCONVERT, FULL, FIND.

    Documentation for sparse
       doc sparse

    Other uses of sparse

       codistributed/sparse      codistributor2dbc/sparse
       codistributor1d/sparse    gpuArray/sparse

help spdiags

 SPDIAGS Sparse matrix formed from diagonals.
    SPDIAGS, which generalizes the function "diag", deals with three
    matrices, in various combinations, as both input and output.
 
    [B,d] = SPDIAGS(A) extracts all nonzero diagonals from the m-by-n
    matrix A. B is a min(m,n)-by-p matrix whose columns are the p nonzero
    diagonals of A. d is a vector of length p whose integer components
    specify the diagonals in A.
 
    B = SPDIAGS(A,d) extracts the diagonals specified by d.
 
    A = SPDIAGS(B,d,A) replaces the diagonals of A specified by d with
        the columns of B. The output is sparse.
 
    A = SPDIAGS(B,d,m,n) creates an m-by-n sparse matrix from the
        columns of B and places them along the diagonals specified by d.
 
    Roughly, A, B and d are related by
        for k = 1:p
            B(:,k) = diag(A,d(k))
        end
 
    Example: These commands generate a sparse tridiagonal representation of
    the classic second difference operator on n points.
        e = ones(n,1);
        A = spdiags([e -2*e e], -1:1, n, n)
 
    Some elements of B, corresponding to positions "outside" of A, are not
    actually used. They are not referenced when B is an input and are set
    to zero when B is an output. See the documentation for an illustration
    of this behavior.
 
    See also DIAG, SPEYE.

    Documentation for spdiags
       doc spdiags

    Other uses of spdiags

       codistributed/spdiags    gpuArray/spdiags

The point being, you don't want to create a FULL matrix of that size. Instead, a diagonal matrix is better built as a sparse matrix.