一、串行LU分解算法（详细见MIT线性代数）

1. LU分解

矩阵分解	LU分解
分解形式	L（下三角矩阵）、U（上三角矩阵）
目的	提高计算效率
前提	（1）矩阵A为方阵；（2）矩阵可逆（满秩矩阵）；（3）消元过程中没有0主元出现，也就是消元过程中不能出现行交换的初等变换

LU分解其实就是将线性方程组：
Ax=bAx = bAx=b
分解为：
LUx=bLUx = bLUx=b
这样一来就会有：
{Ly=bUx=y\begin{cases}Ly = b \\ Ux = y \end{cases}{Ly=bUx=y
先由下三角矩阵求解出y，再通过上三角矩阵完成对x的求解。

2. Gussian Elimination

高斯迭代求解其实就是我们在小学时候学习的二元一次方程组的解法。将一个方程中的未知数，用含有另一未知数的代数式来表示。比如：
{x1−2x2+x3=02x2−8x3=85x1−5x3=10\left\{\begin{array}{l}x_{1}-2 x_{2}+x_{3}=0 \\ 2 x_{2}-8 x_{3}=8 \\ 5 x_{1}-5 x_{3}=10\end{array}\right.⎩⎨⎧x1−2x2+x3=02x2−8x3=85x1−5x3=10
写成增广形式则有：
[1−21002−8850−510]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 & 10\end{array}\right]⎣⎡105−2201−8−50810⎦⎤
使用高斯迭代求解，先使用式1消去式3中的x1x_1x1，则有：
[1−21002−88010−1010]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 10 & -10 & 10\end{array}\right]⎣⎡100−22101−8−100810⎦⎤
再使用式2消去式3中的x2x_2x2，则有：
[1−21002−880030−30]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30\end{array}\right]⎣⎡100−2201−83008−30⎦⎤
最后通过这个上三角矩阵，可以很轻松的求解出x：
{x1−2x2+x3=0x2−4x3=4x3=−1\left\{\begin{array}{l}x_{1}-2 x_{2}+x_{3}=0 \\ x_{2}-4 x_{3}=4 \\ x_{3}=-1\end{array}\right.⎩⎨⎧x1−2x2+x3=0x2−4x3=4x3=−1
简单进行回代即可。

void Guass(float *A,float *x,float *b,int n)
{float *At = (float *)malloc(sizeof(float)*n*n);memcpy(At,A,sizeof(float)*n*n);float *bt = (float *)malloc(sizeof(float)*n);memcpy(bt,b,sizeof(float)*n);for (int i=0;i<n;i++){float pivot = At[i*n+i];   //对角线元素for (int j=i+1;j<n;j++){float scale=-At[j*n+i]/pivot;   //计算消元时所乘的系数for (int k=0;k<n;k++){At[j*n+k]+=scale*At[i*n+k];   //消去元素，一行每一个元素都要相加}bt[j]+=scale*bt[i];   //增广形式，向量也需要操作}}printf("At=\n");printMat(At,n);printf("bt=\n");printVec(bt,n);//回代解方程for (int i=n-1;i>=0;i--){float sum=0;for (int j=i+1;j<n;j++){sum+=At[i*n+j]*x[j];  //此处从最底层开始计算时不会进入此循环}x[i]=(bt[i]-sum)/At[i*n+i];}printf("x=\n");printVec(x,n);free(At);free(bt);free(x);
}

执行结果如下图：

3. 二者的联系

可能大家已经看出来了最后消除完得到了一个上三角矩阵，其实高斯迭代最后所得的即为：Ux=yUx=yUx=y。

但是有个疑问，下三角矩阵L又去哪了？
这个地方就涉及到线代中的矩阵初等变换了，如图所示：

我们在高斯迭代时对于矩阵做了两次初等行变换，也就相当于在矩阵左边乘了两个初等矩阵：
[1000100−51][100010−501][1−21002−8850−510]=[1−21002−880030−30]\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -5 & 1\end{array}\right]\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & 0 & 1\end{array}\right]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 & 10\end{array}\right]=\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30\end{array}\right]⎣⎡10001−5001⎦⎤⎣⎡10−5010001⎦⎤⎣⎡105−2201−8−50810⎦⎤=⎣⎡100−2201−83008−30⎦⎤
将消元矩阵合并，用E表示则为：
EA=[100010−5−51][1−21002−8850−510]=[1−21002−880030−30]=UE A=\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ -5 & -5 & 1\end{array}\right]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 & 10\end{array}\right]=\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30\end{array}\right]=UEA=⎣⎡10−501−5001⎦⎤⎣⎡105−2201−8−50810⎦⎤=⎣⎡100−2201−83008−30⎦⎤=U
两边在左乘E−1E^{-1}E−1则有：
A=[1−21002−8850−510]=E−1U=[100010551][1−21002−880030−30]=LUA=\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 5 & 0 & -5 & 10\end{array}\right]=E^{-1} U=\left[\begin{array}{ccc}1 & 0 & 0 \\ 0 & 1 & 0 \\ 5 & 5 & 1\end{array}\right]\left[\begin{array}{ccc|c}1 & -2 & 1 & 0 \\ 0 & 2 & -8 & 8 \\ 0 & 0 & 30 & -30\end{array}\right]=L UA=⎣⎡105−2201−8−50810⎦⎤=E−1U=⎣⎡105015001⎦⎤⎣⎡100−2201−83008−30⎦⎤=LU
这样就成功将矩阵A分解为了一个下三角矩阵（L）和一个上三角矩阵（U）的乘积。但是需要注意的是并不是所有矩阵都可以分解，要求参见表1的前提！

void Lu_doolittle(float *A,float *x,float *b,int n)
{float *At = (float *)malloc(sizeof(float)*n*n);memcpy(At,A,sizeof(float)*n*n);float *bt = (float *)malloc(sizeof(float)*n);memcpy(bt,b,sizeof(float)*n);//为三角型矩阵L和U分配内存float *L = (float *)malloc(sizeof(float)*n*n);float *U = (float *)malloc(sizeof(float)*n*n);memset(L,0,sizeof(float)*n*n);memset(U,0,sizeof(float)*n*n);//一边高斯消元一边记录相应的值，构成L和Ufor (int i=0;i<n;i++){float pivot =At[i*n+i];for (int j=i+1;j<n;j++){float scale=-At[j*n+i]/pivot;At[j*n+i]=-scale;for (int k=i+1;k<n;k++){At[j*n+k]+=scale*At[i*n+k];}}}//给L矩阵赋值for (int i=0;i<n;i++){L[i*n+i]=1.00;    //对角线置为1for (int j=0;j<i;j++)L[i*n+j]=At[i*n+j];//给下方小三角赋值即可}//给U矩阵赋值for (int i=0;i<n;i++)for (int j=i;j<n;j++)U[i*n+j]=At[i*n+j];printf("L=\n");printMat(L,n);printf("U=\n");printMat(U,n);float *y=(float *)malloc(sizeof(float)*n);//从上向下解出向量yfor (int i=0;i<n;i++){float sum=0;for (int j=0;j<i;j++)sum+=L[i*n+j]*y[j];y[i]=(bt[i]-sum)/L[i*n+i];}//从下向上解出xfor (int i=n-1;i>=0;i--){float sum=0;for (int j=i+1;j<n;j++)sum+=U[i*n+j]*x[j];x[i]=(y[i]-sum)/U[i*n+i];}printf("x=\n");printVec(x,n);
}

4. 部分选主元的LU分解算法

4.1 为什么要选主元？

在进行LU分解时，涉及到除法（计算pilot），如果分母太大则会导致不稳定的情况，如下所示：
假设三位十进制浮点数下计算：
(0.0011.001.002.00)(x1x2)=(1.003.00)\left(\begin{array}{ll}0.001 & 1.00 \\ 1.00 & 2.00\end{array}\right)\left(\begin{array}{l}x_{1} \\ x_{2}\end{array}\right)=\left(\begin{array}{l}1.00 \\ 3.00\end{array}\right)(0.0011.001.002.00)(x1x2)=(1.003.00)
则有：l21=a21a11=1.000.001=1000l_{21}=\frac{a_{21}}{a_{11}}=\frac{1.00}{0.001}=1000l21=a11a21=0.0011.00=1000
LA=(10−10001)(0.0011.001.002.00)=(0.0011.000−1000)=UL A=\left(\begin{array}{cc}1 & 0 \\ -1000 & 1\end{array}\right)\left(\begin{array}{cc}0.001 & 1.00 \\ 1.00 & 2.00\end{array}\right)=\left(\begin{array}{cc}0.001 & 1.00 \\ 0 & -1000\end{array}\right)=ULA=(1−100001)(0.0011.001.002.00)=(0.00101.00−1000)=U
此处的-1000在计算时不够精确，
−1000×1.00+1×2.00=−0.100×104+0.200×10=−0.100×104+0.0002×104=−0.100×104=−1000-1000 \times 1.00+1 \times 2.00=-0.100 \times 10^{4}+0.200 \times 10\quad=-0.100 \times 10^{4}+0.0002 \times 10^{4}=-0.100 \times 10^{4}=-1000−1000×1.00+1×2.00=−0.100×104+0.200×10=−0.100×104+0.0002×104=−0.100×104=−1000
此处的0.0002被舍去了，最后结果为：x=(0,1)Tx=(\mathbf{0}, \mathbf{1})^{\mathrm{T}}x=(0,1)T，与精确值差距很大。

4.2 列主元的LU分解算法（部分选主元）

选主元其实就是为了避免分母很小影响精度的情况出现，因此选取消元后的低阶矩阵内绝对值最大的元素作为主元。此处仅接受列主元的LU分解，列主元只做行变换，相比全主元的方式，可以减少选择主元时的计算量，也可以避免记录交换信息。

列选主元的LU分解算法其实就是在进行高斯消元的过程中需要先选主元在消元。在矩阵变换中相当于先左乘了一个初等矩阵进行行交换，然后又乘了一个矩阵进行消元。
如下：
A=(−32610−705−15)A=\left(\begin{array}{ccc}-3 & 2 & 6 \\ 10& -7 & 0 \\ 5 & -1 & 5\end{array}\right)A=⎝⎛−31052−7−1605⎠⎞
选主元10，交换第一行和第二行（等同于左乘初等矩阵P），然后进行高斯消元再左乘矩阵F，则有：
P1=(010100001),F1=(13101−5101)P_{1}=\left(\begin{array}{lll}0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1\end{array}\right), \quad F_{1}=\left(\begin{array}{ccc}1 & \\ \frac{3}{10} & 1 & \\ -\frac{5}{10} & & 1\end{array}\right)P1=⎝⎛010100001⎠⎞,F1=⎝⎛1103−10511⎠⎞
此时：
A(1)=F1P1A=(10−700−11060525)A^{(1)}=F_{1} P_{1} A=\left(\begin{array}{ccc}10 & -7 & 0 \\ 0 & -\frac{1}{10} & 6 \\ 0 & \frac{5}{2} & 5\end{array}\right)A(1)=F1P1A=⎝⎛1000−7−10125065⎠⎞
再次进行选主元，选择52\frac{5}{2}25，交换二三行，再进行高斯消元得：
A(2)=F2P2A(1)=(10−70525315)\boldsymbol{A}^{(2)}=\boldsymbol{F}_{2} \boldsymbol{P}_{2} \boldsymbol{A}^{(1)}=\left(\begin{array}{rrr}10 & -7 & 0 \\ & \frac{5}{2} & 5 \\ & & \frac{31}{5}\end{array}\right)A(2)=F2P2A(1)=⎝⎛10−72505531⎠⎞

记A(2)A^{(2)}A(2)为U，则有F2P2F1P1A=UF_{2} P_{2} F_{1} P_{1} A=UF2P2F1P1A=U，令L−1L^{-1}L−1为F2P2F1P1F_{2} P_{2} F_{1} P_{1}F2P2F1P1可得：
A=LU=(−310−12511001210)(10−70525315)A=LU=\left(\begin{array}{ccc}-\frac{3}{10} & -\frac{1}{25} & 1 \\ 1 & 0 & 0 \\ \frac{1}{2} & 1 & 0\end{array}\right)\left(\begin{array}{rrr}10 & -7 & 0 \\ & \frac{5}{2} & 5 \\ & \frac{31}{5}\end{array}\right)A=LU=⎝⎛−103121−25101100⎠⎞⎝⎛10−72553105⎠⎞
伪代码如下所示：

二、列主元并行LU分解算法

每次选主元时，只涉及到同一列，因此没有数据通信。
确定主元后，计算 L 的第 k 列
将 l 和 L 的第 k 列传给其他结点
各结点更新自己的数据

伪代码如下所示：
每个节点处理一列，确定主元后，计算消元所需的值，将其和主元位置进行广播。其他节点接收后先进行主元的交换，然后进行消元。

LU分解算法（串行、并行）相关推荐

串行并行程序在效率上的简单比较
串行&并行程序在效率上的简单比较分类: Multi-X 2010-10-15 10:33 1198人阅读评论(0) 收藏举报 parallel工作程序开发作业语言测试开头: 这 ...
基于串行并行ADMM算法的主从配电网分布式优化控制研究
基于串行并行ADMM算法的主从配电网分布式优化控制研究关键词:ADMM 串行并行算法主动配电网无功优化分布式优化参考文档:非复现,仅参考部分模型: 1)<主动配电网分布式无功优化控制方 ...
【转载】串行并行工序混合的生产线数学模型
串行并行工序混合的生产线数学模型串行并行工序混合的生产线数学模型生产线简介示意图目标函数求解结果生产线简介很多生产线工序并不是简单的串行或并行关系,而实两种同时存在的混合并发关系,本文主 ...
【51单片机】串行口连接74LS164进行串行/并行转换，输出到一个七段数码管。数码管循环显示0-9。采用串行通信方式0，定时间隔1秒。
实验内容:51单片机的串行口连接74LS164进行串行/并行转换,然后输出到一个七段数码管.数码管循环显示0-9这10个数字.要求采用串行通信方式0,定时间隔1秒. 工具:proteus+keil # ...
串行并行同步异步通信
终端与其他设备(例如其他终端.计算机和外部设备)通过数据传输进行通信.数据传输可以通过两种方式进行:并行通信和串行通信. 1.串行通信是指使用一条数据线,将数据一位一位地依次传输,每一位数据占据一个 ...
CRC校验串行并行长除移位查表矩阵
CRC校验的几种类型: 长除法,也叫直接计算法移位寄存器,也叫线性移位查表法并行算法一些有用的网页: CRC并行推导 https://blog.csdn.net/Old_Street/arti ...
计算机网络之物理层：1、接口特性、同步异步、串行并行、双工
物理层:1.物理层相关概念思维导图:(学习任务) 接口特性: 数据通信基础: 设计数据通信系统要考虑的三个问题: 三种通信方式: 串行.并行传输: 同步.异步传输: 思维导图:(学习任务) 接口特性 ...
L9825_用于电阻和电感负载的八通道低侧驱动器，具有串行/并行输入控制、输出保护和诊断
描述 L9825是一种八进制低压侧驱动电路,专用于汽车应用.当感应负载被驱动时,输出电压箝位用于反激电流再循环.芯片选择和串行外围接口,用于输出控制和诊断数据传输.两个输出的并联控制输入. 所有功能 ...
【Java8新特性串行/并行流-Optional容器类-时间格式化线程安全等】
Java8新特性二一.并行流与顺序流 1.概念 2.Fork/Join框架 3. Fork/Join框架代码示例: 二.Optional类 1. 什么是Optional对象 2. Optional类 ...
【计算机网络】物理层 : 数据通信 ( 数据通信模型 | 信源 | 信宿 | 信道 | 通信方式 | 单工 | 半双工 | 全双工 | 数据传输方式 | 串行 | 并行 )
文章目录一.数据通信模型示例二.数据通信模型三.数据通信模型分类四.数据通信术语五.三种通信方式六.数据传输方式一.数据通信模型示例数据通信模型示例 : ① 通信场景 : 两台计 ...

LU分解算法（串行、并行）