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

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("

");
    }
    printf("///
");
    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("

");
    }
}

每一行的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("

");
    }
    printf("///
");
    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("

");
    }

    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("

");
    }
    else if (m <= k && k < 2 * m)//printf("///
");    //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("

");
    }

    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("

");
    }
    else if (m <= k && k < 2 * m)//printf("///
");    //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("

");
    }

    return 0;
}


int main()
{
    int n = 40;
    all_sweep_order_symmetric_matrix(n);
    printf("


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

");

    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("
");

    }
    printf("
");

    return 0;
}