void SolutionSlau lu std vector std vector double std vector std vecto

 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
void SolutionSlau::lu(std::vector<std::vector<double> > &L,
std::vector<std::vector<double> > &U, std::vector<double> &x)
{
int size = A.size();
for (int i = 0; i < size; ++i)
{
U[0][i] = A[0][i];
}
for (int i = 0; i < size; ++i)
{
L[i][0] = A[i][0] / U[0][0];
}
for (int i = 1; i < size; ++i)
{
for (int j = i; j < size; ++j)
{
double sum = 0;
for (int k = 0; k < i; ++k)
{
sum += L[i][k] * U[k][j];
}
U[i][j] = A[i][j] - sum;
if (i > j)
L[j][i] = 0;
else
{
double sum = 0;
for (int k = 0; k < i; ++k)
{
sum += L[j][k] * U[k][i];
}
L[j][i] = (A[j][i] - sum) / U[i][i];
}
}
}
std::vector<double> y(size);
for (int i = 0; i < size; ++i)
{
for (int j = 0; j < i; ++j)
{
b[i] -= L[i][j] * y[j];
}
y[i] = b[i] / L[i][i];
}
for (int i = size - 1; i >= 0; --i)
{
for (int j = i + 1; j < size; ++j)
{
y[i] -= U[i][j] * x[j];
}
x[i] = y[i] / U[i][i];
}
}