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