an Order of a sweep of parallel two-side Jacobi method evd

原创已于 2023-05-14 00:59:14 修改 · 131 阅读

0 ·

本内容遵循CC 4.0 BY-SA版权协议

于 2023-05-13 01:22:54 首次发布

这段代码展示了如何使用并行批处理对对称矩阵执行雅可比方法。每个批次处理m对(p,q)元素，目标是通过2m-1次迭代消除所有非对角线元素。算法分为两部分，分别处理不同区域的元素。

evd or svd a symmetric matrix A(nxn),

m = (n+1)/2;

In a sweep, there are 2m-1 parallel orthogonal transformations； m elements in upper triangle，is eliminated to be zero in each orthogonal transformation;

// Ordering_twoside_Jacobi_method.cpp
// ordering_evd_Parallel_Jacobi_method.cpp
//

#include <iostream>

int main()
{
    int n = 16;
    int m = (n + 1) / 2;// 8
    int p, q;

    for (int k = 1; k <= m - 1; k++)//k-th sweep
    {
        for (int batch = 1; batch <= n - m; batch++)
        {
            int branch;
            branch = 0;
            p = 9999;
            q = 9999;
            q = m - k + batch;//  m-k+1, m-k+2, ... , n-k

            if (m - k + 1 <= q && q <= 2 * m - 2 * k)
            {
                branch = 1;
                p = (2 * m - 2 * k + 1) - q;
            }
            else if (2 * m - 2 * k < q && q <= 2 * m - k - 1)
            {
                branch = 2;
                p = (4 * m - 2 * k) - q;
            }
            else if (2 * m - k - 1 < q)
            {
                branch = 3;
                p = n;
            }
            printf("(%d, %d) ", p, q);
        }
        printf("\n\n");
    }
    printf("///////////////////////////////////////////////////////////////\n");
    for (int k = m; k < 2 * m; k++)
    {
        for (int batch = 0; batch <= n - m - 1; batch++)
        {
            p = 9999;
            q = 9999;
            q = 4*m - n - k + batch;
            if(q<2*m-k+1)
            {
                p = n;
            }
            else if (2 * m - k + 1 <= q && q <= 4 * m - 2 * k - 1)
            {
                p = 4 * m - 2 * k - q;
            }
            else if (4 * m - 2 * k - 1 < q)
            {
                p = 6 * m - 2 * k - 1 - q;
            }
            printf("(%d, %d) ", p, q);
        }
        printf("\n\n");
    }
}

每一行的m个(p, q) pair 是一个并行batch，这里 m(=8)个pairs 是一个并行batch；有8个元素被清零；

总共有 2m-1 = 2*8-1=15个batch；这15个batch，会清零全部offdiagonal 元素一次，即，一个sweep；

on jacobi and jacobi-like algorithms for a parallel computer

// Ordering_twoside_Jacobi_method.cpp
// ordering_evd_Parallel_Jacobi_method.cpp
//

#include <iostream>





int all_sweep_order_symmetric_matrix(int n)//
{
    //int n = 16;
    int m = (n + 1) / 2;// 8
    int p, q;

    for (int k = 1; k <= m - 1; k++)//k-th sweep
    {
        for (int ele = 1; ele <= n - m; ele++)
        {
            p = 9999;
            q = 9999;
            q = m - k + ele;//  m-k+1, m-k+2, ... , n-k

            if (m - k + 1 <= q && q <= 2 * m - 2 * k)
            {
                p = (2 * m - 2 * k + 1) - q;
            }
            else if (2 * m - 2 * k < q && q <= 2 * m - k - 1)
            {
                p = (4 * m - 2 * k) - q;
            }
            else if (2 * m - k - 1 < q)
            {
                p = n;
            }
            printf("(%d, %d) ", p, q);
        }
        printf("\n\n");
    }
    printf("///\n");
    for (int k = m; k < 2 * m; k++)
    {
        for (int ele = 1; ele <= n - m; ele++)
        {
            p = 9999;
            q = 9999;
            q = 4 * m - n - k + ele - 1;
            if (q < 2 * m - k + 1)
            {
                p = n;
            }
            else if (2 * m - k + 1 <= q && q <= 4 * m - 2 * k - 1)
            {
                p = 4 * m - 2 * k - q;
            }
            else if (4 * m - 2 * k - 1 < q)
            {
                p = 6 * m - 2 * k - 1 - q;
            }
            printf("(%d, %d) ", p, q);
        }
        printf("\n\n");
    }

    return 0;
}

int kth_batch_of_all_sweep_order_symmetric_matrix(int n, int k)//
{
    //int n = 16;
    int m = (n + 1) / 2;// 8
    int p, q;

    //for (int k = 1; k <= m - 1; k++)//k-th sweep    {
    if (1 <= k && k <= m - 1) {
        for (int ele = 1; ele <= n - m; ele++)
        {
            p = 9999;
            q = 9999;
            q = m - k + ele;//  m-k+1, m-k+2, ... , n-k

            if (m - k + 1 <= q && q <= 2 * m - 2 * k)
            {
                p = (2 * m - 2 * k + 1) - q;
            }
            else if (2 * m - 2 * k < q && q <= 2 * m - k - 1)
            {
                p = (4 * m - 2 * k) - q;
            }
            else if (2 * m - k - 1 < q)
            {
                p = n;
            }
            printf("(%d, %d) ", p, q);
        }
        printf("\n");
    }
    else if (m <= k && k < 2 * m)//printf("///\n");    //for (int k = m; k < 2 * m; k++){
    {
        for (int ele = 1; ele <= n - m; ele++)
        {
            p = 9999;
            q = 9999;
            q = 4 * m - n - k + ele - 1;
            if (q < 2 * m - k + 1)
            {
                p = n;
            }
            else if (2 * m - k + 1 <= q && q <= 4 * m - 2 * k - 1)
            {
                p = 4 * m - 2 * k - q;
            }
            else if (4 * m - 2 * k - 1 < q)
            {
                p = 6 * m - 2 * k - 1 - q;
            }
            printf("(%d, %d) ", p, q);
        }
        printf("\n");
    }

    return 0;
}


int ith_ele_of_kth_batch_of_all_sweep_order_symmetric_matrix(int n, int k, int tid)//
{
    //int n = 16;
    int m = (n + 1) / 2;// 8
    int p, q;
    int ele = tid;
    //for (int k = 1; k <= m - 1; k++)//k-th sweep    {
    if (1 <= k && k <= m - 1) {
        //for (int ele = 1; ele <= n - m; ele++)
        if (1 <= ele && ele <= n - m)
        {
            p = 9999;
            q = 9999;
            q = m - k + ele;//  m-k+1, m-k+2, ... , n-k

            if (m - k + 1 <= q && q <= 2 * m - 2 * k)
            {
                p = (2 * m - 2 * k + 1) - q;
            }
            else if (2 * m - 2 * k < q && q <= 2 * m - k - 1)
            {
                p = (4 * m - 2 * k) - q;
            }
            else if (2 * m - k - 1 < q)
            {
                p = n;
            }
            printf("(%d, %d) ", p, q);
        }
        //printf("\n\n");
    }
    else if (m <= k && k < 2 * m)//printf("///\n");    //for (int k = m; k < 2 * m; k++){
    {
        //for (int ele = 1; ele <= n - m; ele++)
        if (1 <= ele && ele <= n - m)
        {
            p = 9999;
            q = 9999;
            q = 4 * m - n - k + ele - 1;
            if (q < 2 * m - k + 1)
            {
                p = n;
            }
            else if (2 * m - k + 1 <= q && q <= 4 * m - 2 * k - 1)
            {
                p = 4 * m - 2 * k - q;
            }
            else if (4 * m - 2 * k - 1 < q)
            {
                p = 6 * m - 2 * k - 1 - q;
            }
            printf("(%d, %d) ", p, q);
        }
        //printf("\n\n");
    }

    return 0;
}


int main()
{
    int n = 40;
    //all_sweep_order_symmetric_matrix(n);
    printf("\n\n\n");
    for (int k = 1; k <= 2 * ((n + 1) / 2) - 1; k++)
    {
        kth_batch_of_all_sweep_order_symmetric_matrix(n, k);
    }
    printf("\n\n");

    int batch = 2 * ((n + 1) / 2) - 1;//last batch
    int tid = (n + 1) / 2;
    //ith_ele_of_kth_batch_of_all_sweep_order_symmetric_matrix(n, 39, 20);

    int m = (n + 1) / 2;
    for (int batch = 1; batch <= 2 * m - 1; batch++) {
        for (int ele = 1; ele <= m; ele++) {
            ith_ele_of_kth_batch_of_all_sweep_order_symmetric_matrix(n, batch, ele);
        }
        printf("\n");

    }
    printf("\n");

    return 0;
}

标签

#算法 #EVD #JacobiMethod