LU-разложение с выбором ведущего элемента ⇐ C++

[Anonymous]

Это мой код для реализации LU-декомпозиции с выбором ведущего элемента. Скажите пожалуйста, почему у меня не работает выделение ведущего элемента?
Ввод:
Матрица А:

















42
68
35
1
70
25
79
59
63
65


6
46
82
28< /td>
62
92
96
43
28
37


92
5
3
54
93
83
22
17
19
96


48
2772
39
70
13
68
100
36
95
< /tr>

4
12
23
34
74
65
42
12
54
69


< td>48
45
63
58
38
60
24
42
30< /td>
79


17
36
91
43
89
7
41
43
65
49


47
6
91
30
71
51
72
94
49


30
24
85
55
57
41
67
77
32
9


45
40
27
24
38
< td>39
19
83
30
42



Выход:
Матрица L:< /p>












< /th>




1
< td>0
0
0
0
0
0
0
0
0< /td>


0.142857
1
0
0
0
0
0
0
0
0


2,19048
-3,96719
1
0
0
0
0
0
0
0


1,14286
-1,39764
0,602303
1
0
0
0
0
0
0


0.0952381
0.1522310.0342737
-1.14891
1
0
0
0
0 0
69


1,14286
-0,901575
0,398699-0,822519
-2,91362
1
00
0
0


0,404762
0,233596
0,25386
0,244059
0.69331
0.0456395
1
0
0
0


1,11905
-1,93176
0,865283
2,75359
1,42826
-0,785902
2,01127
1
0
< td>0


0.714286
-0.677165
0,483772
0,256388
-0,862612
0,361139
-0.213311
-0.333943
1
0


1,07143
-0,905512
0,25549
-0.318508
-1.38632
0.422558
0.189914
1.21794
-3.58096
1



Матрица U:

















42
68
35
1
70
25
79
59< /td>
63
65


0
36,2857
77
27,8571
52
88,4286
84,7143
34,5714
19
27,7143


0
0
231,807
162,324
145,961
379,051185,03
24,9134
-43,6234
63,5669


0
0
0-20,977
-25,2354
-120,284
-15,3301
65,8844
16,8296
21,1623

< tr>
0
0
0
0
25,4216
-102,029
-2,37443
75,9592
65,9384
80,7253


0
0
0
0
0-436,184
-83,2082
271,315
198,485256,966


0
0
0
0
0
0
-48,5515
-76,4068
-12,746
-72,7809


0
0
0
0
0
0
058.1954
139.056
149.559


0
0
0
0
0
0
0
0
45,5715
-43,5858


0
0
0
0
0
0
0
0
0
-333,13



Общее время = 0 сек.
Из сравнения исходной матрицы и матрицы U видно, что первая и третья строки должны были поменяться местами, но этого не произошло.
Мой код:
#include
#include
#include
#include

using namespace std;

void LU_Decomposition(double * A, double * L, double * U, int N)
{
for (int i = 0; i < N; i++)
{
for (int j = 0; j < N; j++)
{
U[N * i + j] = A[N * i + j];
}
}
for (int i = 0; i < N; i++)
{
L[i * N + i] = 1;
for (int k = i + 1; k < N; k++)
{
double mu = U[k * N + i] / U[N * i + i];
for (int j = i; j < N; j++)
{
U[k * N + j] -= mu * U[i * N + j];
}
L[k * N + i] = mu;
L[i * N + k] = 0;
}
}
for (int i = 1; i < N; i++)
{
for (int j = 0; j < i; j++)
{
U[i * N + j] = 0;
}
}
}

void LU_DecompositionP(double* A, double* L, double* U, int* p, int N)
{
for (int i = 0; i < N; i++)
{
p = i;
for (int j = 0; j < N; j++)
{
U[N * i + j] = A[N * i + j];
L[N * i + j] = 0;
}
L[i * N + i] = 1;
}

for (int i = 0; i < N; i++)
{
int pos = i;
for (int k = i + 1; k < N; k++)
{
if (fabs(U[p[k] * N + i]) > fabs(U[p[pos] * N + i]))
{
pos = k;
}
}
int tmp = p;
p = p[pos];
p[pos] = tmp;

for (int k = i + 1; k < N; k++)
{
double mu = U[p[k] * N + i] / U[p * N + i];
for (int j = i; j < N; j++)
{
U[p[k] * N + j] -= mu * U[p * N + j];
}
L[p[k] * N + i] = mu;
}
}

for (int i = 1; i < N; i++)
{
for (int j = 0; j < i; j++)
{
U[p * N + j] = 0;
}
}
}

void PrintResult(double* A, double* L, double* U, int N)
{
cout

Подробнее здесь: https://stackoverflow.com/questions/791 ... -selection