def gauss(matrix_a,matrix_b):if not np.all(matrix_b==0):augmented=np.concatenate([matrix_a,matrix_b],1)  # 拼接增广矩阵else:augmented=matrix_arow=augmented.shape[0]  # 行数col=augmented.shape[1]  # 列数for j in range(col):if j>=row:breakif augmented[j,j]==0:for t in range(j+1,row):if augmented[t,j] !=0:exchange2row(augmented[t,:], augmented[j,:])break# todo 如果除已处理的行，全为零怎么办？此时秩一定少于未知数个数if augmented[j, j] == 0:  # 检验交换是否成功，如果这一列全为零，跳过print("###第"+j+"列全零，无唯一解")breakfor i in range(j+1,row):if augmented[i,j]==0:continuemul=augmented[i,j]/augmented[j,j]  # 倍数augmented[i,:]=augmented[i,:]-mul*augmented[j,:]  # 转变成上三角阵judge(augmented)

# 化简成行最简求解 回代 对角线元素不能为零 m为增广矩阵
def solve(m,rank):row = m.shape[0]  # 行数col = m.shape[1]  # 列数# 收集结果ans=[]for i in range(rank-1,-1,-1):m[i, :]=m[i,:]/m[i,i]  # 归一化ans.append(-m[i,col-1])for j in range(i-1,-1,-1):if m[j,i]==0:continuemul=m[j,i]/m[i,i]m[j,:]=m[j,:]-mul*m[i,:]print("归一化后的增广矩阵：n",m)print("唯一解为：n",ans[::-1])

mport numpy as np
from mpi4py import MPI# 交换矩阵两行 输入矩阵两行
def exchange2row(a,b):t=a.copy()a=bb=t# 判断是否有解
def judge(m):row = m.shape[0]  # 行数col = m.shape[1]  # 列数h=min(col,row)rank=0for i in range(h):if np.all(m[i,:]==0):  # 全零行后置for j in range(h-1,i,-1):if not np.all(m[j, :] == 0):exchange2row(m[i,:],m[j,:])rank+=1else:rank += 1print("矩阵的秩为：",rank)if m[rank-1,col-1]==0:  # rank_a = rank_abif rank==col-1:print("有唯一解")solve(m, rank)if rank<col-1:print("有无穷多解")else:if np.all(m[rank-1, 0:col-1] == 0):print("无解")elif rank==col-1:print("有唯一解")solve(m,rank)elif rank<col-1:print("有无穷多解")# 化简成行最简求解 回代 对角线元素不能为零 m为增广矩阵
def solve(m,rank):row = m.shape[0]  # 行数col = m.shape[1]  # 列数# 收集结果ans=[]for i in range(rank-1,-1,-1):m[i, :]=m[i,:]/m[i,i]  # 归一化ans.append(-m[i,col-1])for j in range(i-1,-1,-1):if m[j,i]==0:continuemul=m[j,i]/m[i,i]m[j,:]=m[j,:]-mul*m[i,:]print("归一化后的增广矩阵：n",m)print("唯一解为：n",ans[::-1])# 高斯消元法
def mpi_gauss(augmented):row=augmented.shape[0]  # 行数col=augmented.shape[1]  # 列数for j in range(col):if j>=row:breakif augmented[j,j]==0:print("主元为0，换行",j)for t in range(j+1,row):if augmented[t,j] !=0:exchange2row(augmented[t,:], augmented[j,:])break# todo 如果除已处理的行，全为零怎么办？此时秩一定少于未知数个数if augmented[j, j] == 0:  # 检验交换是否成功，如果这一列全为零，跳过print("###第"+j+"列全零，无唯一解")breakcomm.bcast(j, root=0)comm.bcast(augmented, root=0)for p in range(1,comm_size):data = comm.recv(source=p)length = n / (comm_size - 2)start = int((p - 1) * length)end = int(min(start + length, n))start1 = int(max(j + 1, start))for i in range(start1,end):augmented[i, :]=data[i,:]judge(augmented)# 生成指定规模问题
def pre_ques(n):np.random.seed(n)a=np.random.random(size=(n, n))b=np.random.random(size=(n, 1))print("测试数据a=n",a,"n测试数据b=n",b)return np.concatenate([a,b],1)# 初始化mpi
comm = MPI.COMM_WORLD
comm_rank = comm.Get_rank()
comm_size = comm.Get_size()n=4  # 指定问题规模
if comm_rank==0:data0 = pre_ques(n)  # 构建增广矩阵print("rank=", comm_rank,"size=",comm_size)mpi_gauss(data0)
else:print("rank=", comm_rank)col=0length=n/(comm_size-2)  # 0进程不处理，留一个进程处理余数start=int((comm_rank-1)*length)end=int(min(start+length,n))while col<n-1:# 获取当前处理列col = comm.bcast(None, root=0)data0 = comm.bcast(None, root=0)start=int(max(col+1,start))for i in range(start,end):mul = data0[i, col] / data0[col, col]  # 倍数data0[i, :] = data0[i, :] - mul * data0[col, :]  # 转变成上三角阵comm.send(data0, dest=0)