gg_jacobi.m

雅克比过关法算法求矩阵的奇异值分解

function [a,v,t]=gg_jacobi(n,a,Eps) %一次跌代需要（n*12+11）次乘法
                                        %3次开方3个除法(设所需乘法次数为p次)
                                        %若使用n个dsp乘法器，跌代40000次，乘法用的时间为：
                                        %t=(n*12+11+p)*40000/(n*500*1000000)约等于0.001秒
                                        %比较次数为迭代次数的10倍左右的数量级
t(1)=0;        %记录跌代次数
t(2)=0;
for i=1:n
    v(i,i)=1.0;
    for j=1:n
        if i~=j 
            v(i,j)=0.0;
        end
    end
end

v=single(v);
ff=single(0.0);
for i=1:n
    for j=i+1:n
       ff=ff+a(i,j)^2;
    end
end
ff=sqrt(2.0*ff)/n;

while(1)
    flag=0;
    for i=1:n
        for j=i+1:n
            cn=abs(a(i,j));
            t(2)=t(2)+1;
            if cn>ff
                t(1)=t(1)+1;
                flag=1;
                [a,v]=jcb(a,v,n,i,j);
            end
        end
    end
    
    if flag==0
        if ff<Eps
            return;
        else
            ff=ff/(1.0*n);
        end        
    end
    
end
end

function [a,v]=jcb(a,v,n,p,q)

    sn=-a(p,q)*2;
    cn=a(q,q)-a(p,p);
    omega=sn/sqrt(sn*sn+cn*cn);
    if cn<0
        omega=-omega;
    end
    sn=1+sqrt(1-omega*omega);
    sn=omega/sqrt(2*sn);
    cn=sqrt(1-sn*sn);
    
    fm=a(p,p);
    a(p,p)=fm*cn*cn+a(q,q)*sn*sn+a(p,q)*omega;
    a(q,q)=fm*sn*sn+a(q,q)*cn*cn-a(p,q)*omega;
    a(p,q)=0;
    a(q,p)=0;
    
    for j=1:n
        if j~=p && j~=q
            fm=a(p,j);
            a(p,j)=fm*cn+a(q,j)*sn;
            a(q,j)=-fm*sn+a(q,j)*cn;
        end
    end
    
    for i=1:n
        if i~=p && i~=q
            %fm=a(i,p);
            %a(i,p)=fm*cn+a(i,q)*sn;
            %a(i,q)=-fm*sn+a(i,q)*cn;
            a(i,p)=a(p,i);
            a(i,q)=a(q,i);
        end
    end
    
    for i=1:n
        fm=v(i,p);
        v(i,p)=fm*cn+v(i,q)*sn;
        v(i,q)=-fm*sn+v(i,q)*cn;
    end
end