%%%----　サンプルコード　  ----%%%

%%%----------------------------------%%%
%%%  問2.1　ガウスの消去法 　　　　 %%%
%%%----------------------------------%%%

clear all

n=4;
A=[2 1 -2 3;2 0 4 1;3 0 5 2; 1 -1 2 1];
b=[11 -5 -5 -4]';
x=zeros(n,1);

%% Aを上三角行列まで変換する。
for i=1:(n-1)
    if( abs(A(i,i)) < eps ) % eps = 2.2204e-16 
        printf('error');
    end
        for j=(i+1):n
            lamda=-A(j,i)/A(i,i);
            A(j,i+1:n)=A(j,i+1:n)+lamda*A(i,i+1:n);
            b(j)=b(j)+lamda*b(i);
        end
end
%% 後退代入
for i=n:-1:1
    if i<n
        b(i)=b(i)-A(i,i+1:n)*x(i+1:n);
    end
    x(i)=b(i)/A(i,i);
end
printf('Solution x based on Gaussian elimination method:')
x

%%%------------------------------------------%%%
%%%  問2.2　AをLU分解してから、解を計算する  %%%
%%%------------------------------------------%%%

clear all

n=4;
A=[2 1 -2 3;2 0 4 1;3 0 5 2; 1 -1 2 1];
b=[11 -5 -5 -4]';
L=eye(n);
U=zeros(n,n);
x=zeros(n,1);
y=zeros(n,1);

for i=1:(n-1)
    if( abs(A(i,i)) < eps ) % eps = 2.2204e-16 
        printf('error');
    end
    for j=(i+1):n
        lamda = -A(j,i)/A(i,i);
        L(j,i) = -lamda;
        A(j,i+1:n) = A(j,i+1:n)+lamda*A(i,i+1:n);
    end 
end

for i = 1:n
    U(i,i:n) = A(i,i:n);   
end

printf('LU decomposition of A:\n')
printf('L=');
L
printf('U=');
U

%前進代入
for i=1:n
    if i>1
        y(i)=b(i)-L(i,1:i-1)*y(1:i-1);
    else 
        y(i)=b(i);
    end
end

%後退代入
for i=n:-1:1
    if i<n
        y(i)=y(i)-U(i,i+1:n)*x(i+1:n);
    end
    x(i)=y(i)/U(i,i);
end

printf('Solution x=')
x

%%%--------------------------------------------%%%
%%% 問2.3　AをLU分解し、L,Uを元々のAに格納する。  %%%
%%%--------------------------------------------%%%

clear all

n=4;
A=[2 1 -2 3;2 0 4 1;3 0 5 2; 1 -1 2 1];
b=[11 -5 -5 -4]';

for i=1:(n-1)
    if( abs(A(i,i)) < eps ) % eps = 2.2204e-16 
        printf('error');
    end
    for j=(i+1):n
        lamda=-A(j,i)/A(i,i);
        A(j,i)=-lamda;
        A(j,i+1:n)=A(j,i+1:n)+lamda*A(i,i+1:n);
    end  
end

printf('Decomposition L, U saved in A:\n')
A
 


