skyline.m

skyline method. for solving [K]{x}={F}

function[U]=Skyline(k,R)
% ************ Reference Bathe Sec 8.2 *****************
[M,N]=size(k);
l=zeros(M,N);   %l matrix Renrence Bathe page 709
g=zeros(M,N);   %g matrix Renrence Bathe page 709
m=zeros(N,1);   %band for skyline
V=zeros(N,1);   %force matrix Renrence Bathe page 712
U=zeros(N,1);   %U matrix Renrence Bathe page 712
a=0;
b=0;
c=0;
V(1)=R(1);
for i=1:N
    bond=find(k(:,i));
    m(i)=bond(1);
end
l(1,1)=k(1,1);
for j=2:N
    g(m(j),j)=k(m(j),j);
    for i=m(j)+1:j-1
        mm=max(m(i),m(j));
        for r=mm:i-1
            b=l(r,i)*g(r,j)+b;
        end
        g(i,j)=k(i,j)-b;
        b=0;
    end
    for i=m(j):(j-1)
        l(i,j)=g(i,j)/l(i,i);
        c=l(i,j)*V(i)+c;
    end
    V(j)=R(j)-c;
    c=0;
    for r=m(j):j-1
        a=l(r,j)*g(r,j)+a;
    end
    l(j,j)=k(j,j)-a;
    a=0;
end
%% split D form l matrix
% [7  -  +  *  &]                 \            [1  -  +  *  &]         [7]
% |-  8  $  *  #| .................\           |-  1  $  *  #|         |8|
% |-  -  3  ^  @| ................./        l= |-  -  1  ^  @|   &&  D=|3|
% [-  -  -  4  !]                 /            [-  -  -  1  !]         [4]


D=diag(l);
for i=1:N
    l(i,i)=1;
end

%% Backward Substitution
% Reference Bathe page 713
Vbar=V./D;
for i=N:-1:1
    U(i)=Vbar(i);
    for r=m(i):i-1
        Vbar(r)=Vbar(r)-l(r,i)*U(i);
    end

end