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;
}