超级会员免费看
高斯消元法及相关线性求解技术详解
1. 高斯消元法的并行优化策略
在处理线性系统求解时，高斯消元法是一种经典方法。传统的块分配方式可能存在一些效率问题，例如某些处理器可能提前退出计算。为了优化这一过程，可以采用循环或交错分配行（或列）块的方式。
这种方式具有显著优势：
- 所有处理器都能大致承担相同比例（约为 1/P）的工作量，避免了部分处理器过早闲置。
- 由于采用了分块处理，能够有效利用 BLAS2 和 BLAS3 例程，提高计算效率。
还可以进一步结合行和列分块交错存储矩阵 A。具体来说，将矩阵 A 划分为大小为 b 的块，以循环长度为 C 的方式映射到由 P = Pr × Pc 个处理器组成的网状并行计算机上。这种映射可以通过以下函数实现：
Si = F (i, b) 和 Sj = F (j, b)
其中，函数 F 的定义为：
function F (i, b)
floor (i/b) modulo C
return;
这里假设处理的是方阵，且 C² = P。这种方法不仅具有良好的并行性能，还能充分利用 BLAS3 例程的高效性。
2. LU 分解与 BLAS3 的应用
利用 BLAS3 进行 LU 分解时，新的待计算元素位于矩阵 A 右下角的阴影区域。当前要计算的子矩阵 As 大小为 b × b，即行和列的块大小。这种实现方式被包含在用于分布式计算机的 ScaLAPACK 软件中。
3. 并行回代算法
3.1 理想并行回代算法
为了确定并行回代算法的计算复杂度下限，我们先介绍一种理想算法。该算法基于分治策略来求三角矩阵的逆。
考虑一个下三角矩阵 L，将其分解为半大小的子矩阵 L1、L2 和 L3：
L =
[
[L1, 0],
[L2, L3]
]
可以证明：
L⁻¹ =
[
[L1⁻¹, 0],
[-L3⁻¹L2L1⁻¹, L3⁻¹]
]
基于此，可以建立分治算法来求矩阵 L 的逆。以下是主要步骤的伪代码：
Function InvTriangular(L)
if size(L) = 1 return 1/L
else
Set L1 top triangular part of L
Set L2 square part of L
Set L3 bottom triangular part of L
InvL1 = InvTriangular(L1)
InvL3 = InvTriangular(L3)
UpdateL2 = - InvL3 * L2 * InvL1
return L =
[
[InvL1, 0],
[UpdateL2, InvL3]
]
endif
该算法的理想计算成本为 O(log² n)，但在实际中由于硬件限制难以实现。
3.2 实际并行回代算法
对于回代步骤 Ux = y（也是一个三角系统），可以采用以下两种实现方式：
-
按行存储 U
：
for j = n, 1
for j = i + 1, n
bi = bi - uijxj
endfor
xi = bi/uii
endfor
这种方式使用 BLAS1 ddot 例程实现内层循环的点积运算。
-
按列存储 U
：
for j = n, 1
xj = bj/ujj
for i = 1, j - 1
bi = bi - uijxj
endfor
endfor
这种方式被称为左向访问 U，内层循环涉及 daxpy 操作。
4. MPI Bcast 函数
在并行实现高斯消元法之前，我们需要了解一个新的 MPI 函数：MPI Bcast。该函数用于在通信器内的所有进程之间分发相同的数据。
4.1 函数调用语法
int MPI-Bcast(
void* buffer /* in/out */,
int count /* in */,
MPI-Datatype datatype /* in */,
int root /* in */,
MPI-Comm comm /* in */)
4.2 参数解释
-
buffer：发送缓冲区的起始地址。 -
count：发送缓冲区中的元素数量。 -
datatype：发送缓冲区中元素的数据类型。 -
root：广播数据的进程的排名。 -
comm：通信器。
4.3 使用示例
int mynode, totalnodes;
int datasize; // number of data units to be broadcast
int root; // process that is broadcasting its data
MPI-Init(&argc,&argv);
MPI-Comm-size(MPI-COMM-WORLD, &totalnodes);
MPI-Comm-rank(MPI-COMM-WORLD, &mynode);
// Determine datasize and root
double * databuffer = new double[datasize];
// Fill in databuffer array with data to be broadcast
MPI-Bcast(databuffer,datasize,MPI-DOUBLE,root,MPI-COMM-WORLD);
// At this point, every process has received into the
// databuffer array the data from process root
5. 并行高斯消元与回代的实现
以求解 Runge 函数的插值多项式为例，通过基于 Chebyshev 点形成 Vandermonde 矩阵来解决 Ax = b 系统。整个程序分为六个部分：
5.1 第一部分：MPI 初始化与内存分配
#include <iostream.h>
#include <iomanip.h>
#include ”SCmathlib.h”
#include ”SCchapter3.h”
#include<mpi.h>
void ChebyVandermonde(int npts, double *A, int row);
// Global variable to set size of the system
const int size = 10;
int main(int argc, char *argv[]){
int i,j,k,index;
int mynode, totalnodes;
double scaling;
MPI-Status status;
MPI-Init(&argc,&argv);
MPI-Comm-size(MPI-COMM-WORLD, &totalnodes);
MPI-Comm-rank(MPI-COMM-WORLD, &mynode);
int numrows = size/totalnodes;
double **A-local = new double*[numrows];
int * myrows = new int[numrows];
}
注意事项
：
- 使用全局常量变量表示矩阵系统的大小，方便在所有函数中使用。
- 假设矩阵大小能被处理器数量整除，若不能整除则需要特殊处理。
5.2 第二部分：生成矩阵行
double * xpts = new double[size];
ChebyshevPoints(size,xpts);
for(i=0;i<numrows;i++){
A-local[i] = new double[size+1];
index = mynode + totalnodes*i;
myrows[i] = index;
ChebyVandermonde(size,A-local[i],index);
// Set-up right-hand side as the Runge function
A-local[i][size] = 1.0/(1.0+25.0*xpts[index]*xpts[index]);
}
delete[] xpts;
double * tmp = new double[size+1];
double * x = new double[size];
注意事项
：
- 为每行分配 (size + 1) 列，形成增广矩阵，避免额外的通信开销。
- 使用
myrows
数组记录处理器负责的行索引。
5.3 第三部分：高斯消元
int cnt = 0;
for(i=0;i<size-1;i++){
if(i == myrows[cnt]){
MPI-Bcast(A-local[cnt],size+1,MPI-DOUBLE,
mynode,MPI-COMM-WORLD);
for(j=0;j<size+1;j++)
tmp[j] = A-local[cnt][j];
cnt++;
}
else{
MPI-Bcast(tmp,size+1,MPI-DOUBLE,i%totalnodes,
MPI-COMM-WORLD);
}
for(j=cnt;j<numrows;j++){
scaling = A-local[j][i]/tmp[i];
for(k=i;k<size+1;k++)
A-local[j][k] = A-local[j][k] - scaling*tmp[k];
}
}
注意事项
：
- 使用
cnt
变量跟踪每个处理器上已消元的行数。
- 采用循环分布策略，确保所有处理器在大部分阶段都能保持活跃。
5.4 第四部分：回代准备
cnt = 0;
for(i=0;i<size;i++){
if(i==myrows[cnt]){
x[i] = A-local[cnt][size];
cnt++;
}
else
x[i] = 0;
}
注意事项
：将
x
数组初始化为增广矩阵最后一列的值，对于处理器不负责的行初始化为 0。
5.5 第五部分：回代求解
cnt = numrows-1;
for(i=size-1;i>0;i--){
if(cnt>=0){
if(i == myrows[cnt]){
x[i] = x[i]/A-local[cnt][i];
MPI-Bcast(x+i,1,MPI-DOUBLE,mynode,MPI-COMM-WORLD);
cnt--;
}
else
MPI-Bcast(x+i,1,MPI-DOUBLE,i%totalnodes,MPI-COMM-WORLD);
}
else
MPI-Bcast(x+i,1,MPI-DOUBLE,i%totalnodes,MPI-COMM-WORLD);
for(j=0;j<=cnt;j++)
x[myrows[j]] = x[myrows[j]] - A-local[j][i]*x[i];
}
if(mynode==0){
x[0] = x[0]/A-local[cnt][0];
MPI-Bcast(x,1,MPI-DOUBLE,0,MPI-COMM-WORLD);
}
else
MPI-Bcast(x,1,MPI-DOUBLE,0,MPI-COMM-WORLD);
注意事项
：
-
cnt
变量初始化为
numrows - 1
并递减，确保按行回溯求解。
- 检查
cnt >= 0
以避免非法内存访问。
- 在 MPI Bcast 参数中使用指针算术更新传递的地址。
5.6 第六部分：程序结束与清理
if(mynode==0){
for(i=0;i<size;i++)
cout << x[i] << endl;
}
delete[] tmp;
delete[] myrows;
for(i=0;i<numrows;i++)
delete[] A-local[i];
delete[] A-local;
MPI-Finalize();
注意事项 ：仅由第一个处理器输出解，避免所有处理器重复输出。
6. 高斯消元与稀疏系统
许多由偏微分方程离散化产生的线性系统是稀疏的，特别是带状系统。对于对称带状系统，可对 (ijk) 循环进行相应修改以适应稀疏性。例如，kij 循环修改如下：
for k = 1, n - 1
for i = k + 1, min(k + m, n)
ℓik = aik/akk
for j = k + 1, min(k + m, n)
aij = aij - ℓikakj
endfor
endfor
endfor
该算法在串行计算机上的计算复杂度明显低于并行计算机。对于带宽为 m 的 n × n 矩阵，LU 分解的操作计数约为 O(nm²)，回代操作计数为 O(nm)。但在并行计算机上，直接实现该算法可能会导致效率低下，特别是当 m < P 时，只有 m 个处理器能有效工作。
7. 三对角系统的并行循环约简
7.1 循环约简方法概述
循环约简方法是求解三对角线性系统的一种有效并行算法。其主要思想是将未知数分为偶数和奇数编号的条目，类似于黑 - 红高斯 - 赛德尔方法，然后逐步消除奇数编号的条目。
考虑一个三对角系统：
aixi−1 + bixi + cixi+1 = Fi, i = 1, ..., n
假设 n = 2p - 1，如果 n 不满足此条件，则添加额外的平凡方程 xi = 0。
通过线性组合方程来消除奇数编号的未知数，具体步骤如下：
1. 计算 (α2, β2, γ2)、(α4, β4, γ4)、(α6, β6, γ6)。
2. 计算 (ˆb2, ˆc2, ˆF2)、(ˆa4, ˆb4, ˆc4, ˆF4)、(ˆa6, ˆb6, ˆF6)。
3. 计算 (α′4, β′4, γ′4)、(a∗4, F∗4)。
4. 求解 x4, x2, x6, x1, x3, x5, 和 x7。
7.2 串行实现
#include <iostream.h>
#include <iomanip.h>
#include ”SCmathlib.h”
const int size = 15;
void main(){
int i,j,k;
int index1,index2,offset;
double alpha,gamma;
double * x = new double[size];
for(i=0;i<size;i++)
x[i] = 0.0;
double * F = new double[size];
double ** A = new double*[size];
for(i=0;i<size;i++){
A[i] = new double[size];
for(j=0;j<size;j++)
A[i][j] = 0.;
F[i] = (double)i;
}
A[0][0] = -2.0; A[0][1] = 1.0;
A[size-1][size-2] = 1.0; A[size-1][size-1] = -2.0;
for(i=1;i<size-1;i++){
A[i][i] = -2.0;
A[i][i-1] = 1.0;
A[i][i+1] = 1.0;
}
for(i=0;i<log2(size+1)-1;i++){
for(j=pow(2,i+1)-1;j<size;j=j+pow(2,i+1)){
offset = pow(2,i);
index1 = j - offset;
index2 = j + offset;
alpha = A[j][index1]/A[index1][index1];
gamma = A[j][index2]/A[index2][index2];
for(k=0;k<size;k++){
A[j][k] -= (alpha*A[index1][k] + gamma*A[index2][k]);
}
F[j] -= (alpha*F[index1] + gamma*F[index2]);
}
}
int index = (size-1)/2;
x[index] = F[index]/A[index][index];
for(i=log2(size+1)-2;i>=0;i--){
for(j=pow(2,i+1)-1;j<size;j=j+pow(2,i+1)){
offset = pow(2,i);
index1 = j - offset;
index2 = j + offset;
x[index1] = F[index1];
x[index2] = F[index2];
for(k=0;k<size;k++){
if(k!= index1)
x[index1] -= A[index1][k]*x[k];
if(k!= index2)
x[index2] -= A[index2][k]*x[k];
}
x[index1] = x[index1]/A[index1][index1];
x[index2] = x[index2]/A[index2][index2];
}
}
for(i=0;i<size;i++){
cout << x[i] << endl;
}
delete[] x;
delete[] F;
for(i=0;i<size;i++)
delete[] A[i];
delete[] A;
}
注意事项 ：该程序使用了比实际需求更多的内存，因为三对角系统可以减少内存使用。
7.3 并行实现
#include <iostream.h>
#include <iomanip.h>
#include ”SCmathlib.h”
#include<mpi.h>
int main(int argc, char *argv[]){
int i,j,k,size,index;
int index1,index2;
int mynode, totalnodes;
double alpha,gamma;
const int numrows = 5;
MPI-Status status;
MPI-Init(&argc,&argv);
MPI-Comm-size(MPI-COMM-WORLD, &totalnodes);
MPI-Comm-rank(MPI-COMM-WORLD, &mynode);
size = (int) pow(2.0,log2(totalnodes+1)+1)-1;
double ** A = new double*[numrows];
for(i=0;i<numrows;i++){
A[i] = new double[size+1];
for(j=0;j<size+1;j++)
A[i][j] = 0.0;
}
if(mynode==0){
A[0][0] = -2.0; A[0][1] = 1.0;
A[1][0] = 1.0; A[1][1] = -2.0; A[1][2] = 1.0;
A[2][1] = 1.0; A[2][2] = -2.0; A[2][3] = 1.0;
}
else if(mynode==(totalnodes-1)){
index = 2*mynode;
A[0][index-1] = 1.0; A[0][index] = -2.0;
A[0][index+1] = 1.0;
index = 2*mynode+1;
A[1][index-1] = 1.0; A[1][index] = -2.0;
A[1][index+1] = 1.0;
A[2][size-2] = 1.0; A[2][size-1] = -2.0;
}
else{
for(i=0;i<3;i++){
index = i + 2*mynode;
A[i][index-1] = 1.0;
A[i][index] = -2.0;
A[i][index+1] = 1.0;
}
}
for(i=0;i<3;i++)
A[i][size] = 2*mynode+i;
int numactivep = totalnodes;
int * activep = new int[totalnodes];
for(j=0;j<numactivep;j++)
activep[j] = j;
for(j=0;j<size+1;j++){
A[3][j] = A[0][j];
A[4][j] = A[2][j];
}
for(i=0;i<log2(size+1)-1;i++){
for(j=0;j<numactivep;j++){
if(mynode==activep[j]){
index1 = 2*mynode + 1 - pow(2,i);
index2 = 2*mynode + 1 + pow(2,i);
alpha = A[1][index1]/A[3][index1];
gamma = A[1][index2]/A[4][index2];
for(k=0;k<size+1;k++)
A[1][k] -= (alpha*A[3][k] + gamma*A[4][k]);
if(numactivep>1){
if(j==0){
MPI-Send(A[1],size+1,MPI-DOUBLE,activep[1],0,
MPI-COMM-WORLD);
}
else if(j==numactivep-1){
MPI-Send(A[1],size+1,MPI-DOUBLE,activep[numactivep-2],
1,MPI-COMM-WORLD);
}
else if(j%2==0){
MPI-Send(A[1],size+1,MPI-DOUBLE,activep[j-1],
1,MPI-COMM-WORLD);
MPI-Send(A[1],size+1,MPI-DOUBLE,activep[j+1],
0,MPI-COMM-WORLD);
}
else{
MPI-Recv(A[3],size+1,MPI-DOUBLE,activep[j-1],0,
MPI-COMM-WORLD,&status);
MPI-Recv(A[4],size+1,MPI-DOUBLE,activep[j+1],1,
MPI-COMM-WORLD,&status);
}
}
}
}
numactivep = 0;
for(j=activep[1];j<totalnodes;j=j+pow(2,i+1)){
activep[numactivep++]=j;
}
}
double * x = new double[totalnodes];
for(j=0;j<totalnodes;j++)
x[j] = 0.0;
if(mynode==activep[0]){
x[mynode] = A[1][size]/A[1][(size-1)/2];
}
double tmp;
for(i=log2(size+1)-3;i>=0;i--){
tmp = x[mynode];
MPI-Allgather(&tmp,1,MPI-DOUBLE,x,1,MPI-DOUBLE,
MPI-COMM-WORLD);
numactivep = 0;
for(j=activep[0]-pow(2.0,i);j<totalnodes;j=j+pow(2.0,i+1)){
activep[numactivep++]=j;
}
for(j=0;j<numactivep;j++){
if(mynode == activep[j]){
x[mynode] = A[1][size];
for(k=0;k<totalnodes;k++){
if(k!=mynode)
x[mynode] -= A[1][2*k+1]*x[k];
}
x[mynode] = x[mynode]/A[1][2*mynode+1];
}
}
}
tmp = x[mynode];
MPI-Allgather(&tmp,1,MPI-DOUBLE,x,
1,MPI-DOUBLE,MPI-COMM-WORLD);
for(k=0;k<totalnodes;k++){
A[0][size] -= A[0][2*k+1]*x[k];
A[2][size] -= A[2][2*k+1]*x[k];
}
A[0][size] = A[0][size]/A[0][2*mynode];
A[1][size] = x[mynode];
A[2][size] = A[2][size]/A[2][2*mynode+2];
delete[] activep;
for(i=0;i<numrows;i++)
delete[] A[i];
delete[] A;
delete[] x;
MPI-Finalize();
}
注意事项
：
- 与并行高斯消元法类似，通过增广矩阵减少通信开销。
- 通信通过 MPI Send 和 MPI Recv 调用实现，使用
activep
数组跟踪活跃处理器。
- 使用 MPI Allgather 命令确保所有处理器在每个级别都能获取可用的解。
综上所述，通过这些并行算法和技术，可以有效提高线性系统求解的效率，特别是在处理大规模矩阵和并行计算环境中。不同的算法适用于不同的场景，需要根据具体问题选择合适的方法。
高斯消元法及相关线性求解技术详解
8. 并行循环约简的通信与效率分析
在并行循环约简算法中，通信是一个关键环节。当有 $P = \frac{n - 1}{2}$ 个处理器，且每个处理器存储一个三元组时，消除奇数编号未知数的工作可以并行进行。在这个阶段结束后，每个处理器持有一个约简后的方程。接下来，处理器之间需要进行数据交换，这可以通过最近邻通信实现。例如，$P_2$ 会从 $P_1$ 和 $P_3$ 接收数据，$P_4$ 会从 $P_3$ 和 $P_5$ 接收数据，依此类推。这意味着大约一半的处理器（如所有奇数编号的 $P_1$、$P_3$ 等）会处于闲置状态。
经过 $(p - 1)$ 个约简阶段（其中 $n = 2^p - 1$）后，只有一个处理器会求解最终方程。但在回代过程中，所有处理器都会再次参与工作。该算法的串行操作计数为 $O(13n)$ 次乘法，而标准 LU 分解为 $O(5n)$ 次乘法。不过，很多前面的计算可以并行进行，从而提高整体效率。
9. 串行与并行循环约简实现的对比
| 对比项 | 串行实现 | 并行实现 |
|---|---|---|
| 复杂度 | 串行操作计数为 $O(13n)$ 次乘法，通常比标准 LU 分解复杂 | 可并行处理部分计算，提高整体效率，但通信开销需考虑 |
| 内存使用 | 可能使用比实际需求更多的内存，未充分利用三对角系统特性 | 与并行高斯消元法类似，通过增广矩阵减少通信开销 |
| 通信 | 无通信开销 | 需通过 MPI Send、MPI Recv 和 MPI Allgather 等进行数据交换和同步 |
| 适用场景 | 适合理解算法原理，不适合大规模计算 | 适合大规模三对角系统的并行求解 |
10. 并行算法的总结与选择建议
在处理线性系统求解问题时，我们介绍了多种并行算法，包括高斯消元法的并行优化、并行回代算法以及三对角系统的并行循环约简算法。这些算法各有优缺点，适用于不同的场景。
- 高斯消元法的并行优化 ：适用于一般的线性系统求解，通过循环或交错分配行（或列）块以及结合行和列分块交错存储矩阵，可以提高并行性能，充分利用 BLAS 例程。但对于稀疏系统，直接应用可能效率不高。
- 并行回代算法 ：理想算法提供了计算复杂度的下限，但实际中受硬件限制难以实现。实际算法根据矩阵存储方式（按行或按列）有不同的实现方式，需要根据具体情况选择。
- 三对角系统的并行循环约简 ：对于三对角线性系统，该算法是一种有效的并行求解方法。虽然串行操作计数较高，但很多计算可以并行进行，通过合理的通信策略可以提高效率。
选择合适的算法时，需要考虑以下因素：
-
矩阵特性
：矩阵是稠密还是稀疏，是否为三对角矩阵等。
-
计算资源
：可用的处理器数量、内存大小等。
-
问题规模
：大规模问题更适合并行算法，小规模问题串行算法可能更简单高效。
11. 代码优化建议
11.1 高斯消元法代码优化
- 内存管理 ：在处理大规模矩阵时，合理管理内存至关重要。例如，在并行高斯消元法中，可以考虑使用动态内存分配，避免不必要的内存浪费。
- 通信优化 ：减少通信开销是提高并行算法效率的关键。可以通过合并通信操作、减少通信次数等方式优化。例如，在并行循环约简中，使用增广矩阵可以同时传递行和右侧信息，减少通信量。
11.2 三对角系统代码优化
- 内存压缩 ：对于三对角系统，只需要存储三条对角线元素，而不是整个矩阵。在串行和并行实现中，可以修改代码以减少内存使用。
- 通信策略优化 ：在并行循环约简中，合理安排处理器之间的通信顺序和方式，可以减少闲置处理器的数量，提高整体效率。
12. 未来发展趋势
随着计算机技术的不断发展，并行计算领域也在不断进步。未来，可能会出现更高效的并行算法和更强大的并行计算硬件。例如，量子计算技术的发展可能会为线性系统求解带来革命性的变化。同时，人工智能和机器学习领域对大规模线性系统求解的需求也在不断增加，这将推动并行算法的进一步优化和创新。
13. 总结
本文详细介绍了高斯消元法及相关线性求解技术，包括并行优化策略、LU 分解、并行回代算法、稀疏系统处理和三对角系统的并行循环约简。通过对这些算法的原理、实现和性能分析，我们可以看到不同算法适用于不同的场景。在实际应用中，需要根据具体问题的特点和计算资源选择合适的算法，并进行必要的优化，以提高线性系统求解的效率。
以下是并行循环约简算法的流程 mermaid 图：
graph TD
A[初始化矩阵和处理器] --> B[循环约简阶段]
B --> C{是否完成约简?}
C -- 否 --> D[处理器通信与更新方程]
D --> B
C -- 是 --> E[求解最终方程]
E --> F[回代过程]
F --> G[输出解]
总之，线性系统求解是科学计算和工程领域中的重要问题，并行算法的发展为解决大规模问题提供了有力的工具。通过不断研究和优化这些算法，我们可以更好地应对各种复杂的实际问题。
扫一扫
8793
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
