function sol=linearsol(mat, b)
% LINEARSOL(MAT,B) は線形一次方程式:
%             mat x = b
% の解ｘをＱＵ分解法を用いて計算するプログラムである．
% ＱＵ分解にはハウスホールダー変換を用いている．
% 注意!!
% 'mat' は正則な正方マトリックスでなければならない．
% 'b' のサイズは'mat'の列のサイズと同じでなければならない.
%
% Copyright (C) 2007 清水 信行, 曽我部 潔, 小金澤 鋼一 All Rights Reserved.

% =================================================================================
% プログラムの使用条件
% 1. 本プログラムおよびデータを商用で，頒布してはならない．
% 2. 本プログラムおよびデータを著作者の許可なく営利目的で利用してはならない．
% 3. 本プログラムおよびデータの改変は，営利目的でない限り自由である．
% 4. 本プログラムおよびデータを使用することにより生じた損害などについては，著者，
%    プロクラム著作者，(社)日本機械学会，および(株)コロナ社は一切の責任を負わない．
% 5. 本プログラムの著作権は，本書編著者ならびにプログラム著作者に帰す．
% =================================================================================

nm=size(mat);
nr=nm(1);
nc=nm(2);
bs=size(b);
if bs(1)==1 bc=b'  ;   %  もしベクトルｂが横ベクトルならば縦ベクトルにする
else bc=b ;   
end
bs=size(bc);
nb=bs(1);

% 行列が正方でなければエラーを返す．
if nr~=nc
    error(' Specified matrix is not square')
end

% ベクトルと行列の]サイズが一致していなければエラーを返す．
if nr~=nb
     error(' Specified matrix and vector is not corresponding in size')
end

n=nr;                     % 行列の行数
idm=eye(n);            % n×ｎの単位行列
q=idm;                   % 直交行列の初期設定
x=zeros(n,1);          %  解の初期化
utb=[mat bc];         %  ｎ×（ｎ＋１）の行列をつくる

% Householder transformation
for i=1: n-1
    mcol=utb(:,i);                    % ｉ列をとる
    mcol=mcol(i:n);                 % i列のｉ行目からｎ行目までをとる
    cs=sqrt(mcol' * mcol);       % 二乗和平方根
    if abs(cs) < 0.000000001
        error(' Specified matrix is singular')
        exit;
    end
    if mcol(1) >0 cs =-cs;
    end
    alpha=cs/(cs-mcol(1));
    w=-mcol/cs;
    w(1)=w(1)+1;
    w=[zeros(i-1,1) ;w];
    h=idm-alpha*(w*w');
    utb=h*utb;
    end

% Solving x recurrently from the last to the first:
for i=n : -1 : 1
   res=0;
    for j=i+1: n
        res=res+utb(i, j)*x(j);
    end
    x(i)=(utb(i, n+1)-res)/utb(i,i);
end
sol=x;