这几天在听麻省理工Gilbert Strang教授的线性代数，其中讲到了高斯消元法等，又加上教授一直提及MATLAB运算思想，所以这个矩阵系列我将用C语言模拟实现MATLAB的内部计算（不一定对）

百科：在线性代数中， LU分解(LU Decomposition)是矩阵分解的一种，可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积（有时是它们和一个置换矩阵的乘积）。LU分解主要应用在数值分析中，用来解线性方程、求反矩阵或计算行列式。实现矩阵的LU分解。矩阵的三角分解又称LU分解，它是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。

LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。这正是所谓的杜尔里特算法：从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。这类算法的复杂度一般在(三分之二的n三次方) 左右

注意：当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU，且当L的对角元全为1时分解唯一

在MATLAB中是这样使用的LU函数使用语句：

Y = lu(A)
[L,U] = lu(A)
[L,U,P] = lu(A)
[L,U,P,Q] = lu(A)
[L,U,P,Q,R] = lu(A)

具体用c怎么实现的呢？当然你也可以用其他语言

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
//

// main.c

// A=LU分解

//

// Created by 杨旭磊 on 16/4/4.

// Copyright © 2016年 杨旭磊. All rights reserved.

//

include <stdio.h>

/第一种方法直接解出方程的根/

//void solve(float l[][100],float u[][100],float b[],float x[],int n)

//{int i,j;

// float t,s1,s2;

// float y[100];

// for(i=1;i<=n;i++) /* 第一次回代过程开始 */

// {s1=0;

// for(j=1;j<i;j++)

// {

// t=-l[i][j];

// s1=s1+t*y[j];

// }

// y[i]=(b[i]+s1)/l[i][i]; }

// for(i=n;i>=1;i--) /* 第二次回代过程开始 */

// {

// s2=0;

// for(j=n;j>i;j--)

// {

// t=-u[i][j];

// s2=s2+t*x[j];

// }

// x[i]=(y[i]+s2)/u[i][i];

// }

//}

//

//void main()

//{float a[100][100],l[100][100],u[100][100],x[100],b[100];

// int i,j,n,r,k;

// float s1,s2;

// for(i=1;i<=99;i++)/将所有的数组置零，同时将L矩阵的对角值设为1/

// for(j=1;j<=99;j++)

// {

// l[i][j]=0,u[i][j]=0;

// if(j==i) l[i][j]=1;

// }

//

// printf ("input n:\n");/输入方程组的个数/

// scanf("%d",&n);

// printf ("input array A:\n");/读取原矩阵A/

// for(i=1;i<=n;i++)

// for(j=1;j<=n;j++)

// scanf("%f",&a[i][j]);

// printf ("input array B:\n");/读取列矩阵B/

// for(i=1;i<=n;i++)

// scanf("%f",&b[i]);

// for(r=1;r<=n;r++)/求解矩阵L和U/

// {

// for(i=r;i<=n;i++)

// {

// s1=0;

// for(k=1;k<=r-1;k++)

// s1=s1+l[r][k]*u[k][i];

// u[r][i]=a[r][i]-s1;

// }

// for(i=r+1;i<=n;i++)

// {s2=0;

// for(k=1;k<=r-1;k++)

// s2=s2+l[i][k]*u[k][r];

// l[i][r]=(a[i][r]-s2)/u[r][r];

// }

// }

// printf("array L:\n");/输出矩阵L/

// for(i=1;i<=n;i++)

// {

// for(j=1;j<=n;j++)

// printf("%7.3f ",l[i][j]);

// printf("\n");

// }

// printf("array U:\n");/输出矩阵U/

// for(i=1;i<=n;i++)

// {

// for(j=1;j<=n;j++)

// printf("%7.3f ",u[i][j]);

// printf("\n");

// }

// solve(l,u,b,x,n);

// printf("解为:\n");

// for(i=1;i<=n;i++)

// printf("x%d=%f\n",i,x[i]);

//}

/第二种方法,只提供思路，部分printf未写/

void main(){

 

int i,j,k,m,n,q,r,t,y;

float a[100][100],c[100][100],d[100][100],l[100][100],u[100][100];

scanf("%d",&n);

/输入矩阵A/

for (i = 1; i < n + 1; i++) {

printf("输入矩阵第%d行元素：%d\n",i);

for (j = 1; j < n +1; j++) {

scanf("%f",&a[i][j]);

}

}

/判断A的各阶顺序主子式（矩阵A下标kk）是否全大于0，若不成立进行LU分解/

for (t = 1, k = d[1][1];t < n + 1;t++) {

if (d[1][1]==0) { //判断矩阵第一行是否为0，若为0则不分解

 

break;

 

}else{

for (i = t + 1; i < n + 1; i++) { //进行第t次行优先的所有元素的消元

 

m = d[i][t]/d[t][t]; //计算行乘数

 

for (j = t + 1; j < n; j++) {

 

d[q][j] = d[q][j] - m*d[t][j]; //计算每次不为0的元素的值

m = 0; //将m归零

 

}

 

n = d[t+1][t+1]; //取主元n

k *= n; //计算各阶的顺序主子式

 

}

}

}

/若可分解，进行LU分解，l为单位下三角矩阵，u为上三角矩阵/

for (r = 2; r < n + 1; r++) { //控制u的行数（L的列数）

 

for (j = r; j < n + 1; j++) { //控制u的第r行的j列变化

 

int x = 0; //初始化中间变量x

for (k = 1; k < r; k++) {

x += l[i][k] * u[k][r]; //求元素之前对应的行和列的元素和

l[i][r] = (c[i][r]-y)/u[r][r]; //计算矩阵l的第r列元素

}

}

}

//可以用两个for打印你想看到的l或u矩阵，在WordPress写代码是在太压抑了，下面就不写了，参考第一种

}