Jacobi并行拆解

作者：桂。

时间：2018-04-23  21:12:02

链接：http://www.cnblogs.com/xingshansi/p/8921815.html

---

前言

一、并行拆解思路

　　回顾前文，对于8X8的实数矩阵：

仿真：

clc;clear all;

X = rand(8);

R = X+X';

Iteration = 20;

[D,U] = Jac_sweep(R,Iteration);

[u,s,v] = svd(R);

[sort(diag(s)),sort(abs(diag(D)))]

　　Jac_sweep并行代码：

function [D,U] = Jac_sweep(A,Iteration)

iter = 0;

n = size(A,1);

U = eye(n);

A_c = [1,2;3,4;5,6;7,8;...

    2,3;4,5;6,7;1,8;...

    1,3;2,4;5,7;6,8;...

    3,5;4,6;1,7;2,8;...

    1,4;3,6;5,8;2,7;...

    2,5;4,7;1,6;3,8;...

    1,5;2,6;3,7;4,8];

while iter <Iteration

    iter = iter+1;

    for flag = 1:size(A_c,1)/4

        T = eye(n);

        for t = 1:4

            p = A_c((flag-1)*4+t,1);q = A_c((flag-1)*4+t,2);

            y = 2*A(p,q);

            x = A(p,p)-A(q,q);

            phi = atan2(y,x)/2;

            T(p,p) = cos(phi);

            T(q,q) = cos(phi);

            T(p,q) = -sin(phi);

            T(q,p) = sin(phi);

        end

        D = T'*A*T;

        U = U*T;

        A = D;

    end

end

　　

打印结果与SVD分解的结果一致，并行思路可行。其他维度依次类推。对于维度N的矩阵，

• N为偶数，可并行N/2路；
• N为奇数，可并行[N-1]/2路；

二、改进思路

　　每一次sweep，需要1次Cordic：phi = atan2(y,x)/2，和两次Cordic（两次可并行）：cos(phi) / sin(phi)，二者串行。对于atan2操作，可借鉴复数相位近似估计一文的思路，即对于atan的计算，考虑到CORDIC耗时较长，内存资源充足的情况下，1）直接查表；若内存相对紧张，2）多项式逼近。二者较CORDIC均减少运算时间。

　　另外，关于定点仿真可调用fi工具包，其中CORDIC对应指令：cordicatan2、cordiccos......