function solution=linsolver(A,b)
%This is an introduction of this M file. You need to give inputs A and b.
%A should be an n by n matrix and b should be an n by 1 vector. If A is nonsingular, you will
%get one output which is the solution of Ay=b. If A is singular, the answer
%will be showed as The matrix is singular. When you use this one, please
%make sure that this is the only file in Matlab directory in case there is
%some overlap of M files.
n=length(A);
check=1;
xindex=1:n;
for k=1:n
    big=0;
    for i=k:n
        for j=k:n
            if abs(A(i,j))>big
                big=abs(A(i,j));
                irow=i;
                icol=j;
            end
        end
    end
    if big<=10^-15 %singular case
        sprintf('The matrix is singular.')
        check=0;
        break
    end
    q=b(irow);
    b(irow)=b(k);
    b(k)=q;
    for i=1:n;
        p=A(irow,i);
        A(irow,i)=A(k,i);
        A(k,i)=p;
    end
    for i=1:n;
        p=A(i,icol);
        A(i,icol)=A(i,k);
        A(i,k)=p;
    end
    r=xindex(icol);
    xindex(icol)=xindex(k);
    xindex(k)=r;
    for i=(k+1):n
        l=A(i,k);
        b(i)=b(i)-l*b(k)/A(k,k);
        for j=k:n
        A(i,j)=A(i,j)-l*A(k,j)/A(k,k);
        end
    end
end
x=b;
if check==1
    x(n)=x(n)/A(n,n);
    for i=n-1:-1:1
        sum=0;
        for j=n:-1:(i+1)
            sum=sum+A(i,j)*x(j);
        end
        x(i)=(x(i)-sum)/A(i,i);
    end
    solution=zeros(n,1);
    for i=1:n
        solution(xindex(i))=x(i);
    end
end