function fact=fsf(mat)

% FSF(MAT) はマトリックスを相似分解するプログラムである．
% 引数として与える行列’mat’は正方行列でなければならない．
% 行列が対称行列であれば，対角行列’d’ と直交行列’q’ が算出される．
% 行列が非対称行列で，全ての固有値が実数ならば，
% 上三角行列’d’ と直交行列’q’ が算出される．
% 行列が非対称行列で，複素数の固有値を持つ場合，エラーが表示される．
% この関数は２つの行列’d’と’q’を，[d q]のフォームで出力する．
% 
% Copyright (C) 2007 清水 信行, 曽我部 潔, 小金澤 鋼一 All Rights Reserved.

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

ci=1;
d=mat;
nm=size(mat);
n=nm(1);                %   行列の行数
idm=eye(n);            %    n×ｎの単位行列
q=idm;                   %  直交行列の初期設定
h=idm;
ite=0;
epsi=0.00001;          %   ゼロに収束したかどうかの判定値

while   ci > 0
    ci=0;
    for i=1:nm
        for j=nm : -1 : i+1
            dji=d(j,i);
            if abs(dji) > epsi
                c=d(i,i)-d(j,j);
                ci=ci+1;
                ite=ite+1;
                g=c^2;
                a=(1.+g/(g+(2*d(i,j)+dji)*dji))/4;
                b=(a*dji)^2;
                r=g+b;
                h(i,i)=(g-b)/r;
                h(j,j)=-h(i,i);
                p=2*a*c*dji/r;
                h(i,j)=p;
                h(j,i)=p;
                d=h*d*h;
                q=q*h;
                h=idm;
             end
         end
    end
 
 if ite > 1000
            error('収束しません.行列は共役複素数の固有値を持っています')
 end
end
ite
fact=[d,q];