Sparse Martices

I am trying to write a program that will create an nxn matrix A with 3 on the main diagonal, -1 on the super and subdiagonal, and 1/2 in the (i,N + 1 ? i) position for all i = 1, · · · , n, except for i = n/2 and n/2 + 1. Then I make a vector b that has 2.5 at each end, n ? 4 repetitions of 1.5 and 2 repetitions of 1.0 in the middle so that b for n=12 looks like (2.5, 1.5, 1.5, 1.5, 1.5, 1.0, 1.0, 1.5, 1.5, 1.5, 1.5, 2.5). So that's A and b. Then I wrote code for the Jacobian, Gauss-Seidel, and Successive Over-Relaxation methods to solve the system.

All of that I can do. Now I am trying to figure out how to run this on a very large (n=100,000) matrix. I know this is considered a sparse matrix, but I have no idea how to alter my code to accommodate it.

Here's the main program and the Jacobian method

clear
n = input( 'Input size desired: ');
A=eye(n,n);
for i = 1:n
X(i,1)=.5;
end
for i=1:n
    A(i,n+1-i)=1/2;
    A(i,i)=3;
end
for i=1:n-1
    A(i,i+1)=-1;
    A(i+1,i)=-1;
end
b=zeros(n,1);
b(1,1)=2.5;
b(n,1)=2.5;
r=(n-4)/2;
for i=1:r
    b(i+1,1)=1.5;
    b(n-i,1)=1.5;
end
b(r+2,1)=1;
b(r+3,1)=1;
jacobi(n,A,b,X,10^-5,1000);

function [] = jacobi(n,A,b,X,tol,N)
x=zeros(n,1);
mu=zeros(n,1);
%----------Step 1
k=1;
%----------Step 2
num=0;
while (k<=N)
%----------Step 3
    for i=1:n
        sum=0;
        for j=1:n
            if j ~= i
                sum=sum+A(i,j)*X(j,1);
            end
        end
        num=-sum+b(i,1);
        x(i,1)= num/A(i,i);
    end
%---------Step 4
    if norm(x-X,inf)<tol
        display 'Success! Iterations for Jacobian Method:'
        disp(k)
        return
    end
%--------Step 5
    k = k+1;
%--------Step 6
    for i=1:n
        X(i,1)=x(i,1);
    end
end
%--------Step 7
display ' Failure! Iterations exceeded max!'
return
end

Answer 1

[Bruno Luong]

You should take a look at SPARSE command, and more specifically SPDIAGS command that can help to populate the diagonals of sparse matrix.

Answer 2

[veronica]

Thanks for the quick reply. I have looked at SPARSE and SPDIAG, but I don't know how to use those in my code.

[Bruno Luong]

You have tried in 20 minutes at most then gave up. I think you should try longer.

If you want to check what SPARSE or SPDIAGS do, use FULL command

n = 10

e = ones(n,1);

A = spdiags([e -2*e e], -1:1, n, n);

disp(full(A))