int main() {
int n = 3;
array double a[3][3] = {2,1,-2,
4,-1,2,
2,-1,1} ; /* n-by-n matrix */
array double a2[3][3] = {-1,5,6,
3,-6,1,
6,8, 9} ; /* n-by-n matrix */
array double l[n][n], u[n][n];
array int p[n][n];
int status;
status = ludecomp(a, l, u);
printf("l =\n%f\n", l);
printf("u =\n%f\n", u);
printf("l*u =\n%f\n", l*u);
/* pass the address of p[0][0]
or use ludecomp(a, l, u, &p);
or use ludecomp(a, l, u, &p[0][0]); */
status = ludecomp(a, l, u, p);
printf("l =\n%f\n", l);
printf("u =\n%f\n", u);
printf("p =\n%d\n", p);
printf("p*l*u =\n%f\n", p*l*u);
status = ludecomp(a2, l, u);
if(status == 0) {
printf("l =\n%f\n", l);
printf("u =\n%f\n", u);
printf("l*u =\n%f\n", l*u);
}
else
printf("error: numerical error in ludecomp()\n");
status = ludecomp(a2, l, u, p);
printf("l =\n%f\n", l);
printf("u =\n%f\n", u);
printf("p =\n%d\n", p);
printf("p*l*u =\n%f\n", p*l*u);
}
l =
0.500000 1.000000 0.000000
1.000000 0.000000 0.000000
0.500000 -0.333333 1.000000
u =
4.000000 -1.000000 2.000000
0.000000 1.500000 -3.000000
0.000000 0.000000 -1.000000
l*u =
2.000000 1.000000 -2.000000
4.000000 -1.000000 2.000000
2.000000 -1.000000 1.000000
l =
1.000000 0.000000 0.000000
0.500000 1.000000 0.000000
0.500000 -0.333333 1.000000
u =
4.000000 -1.000000 2.000000
0.000000 1.500000 -3.000000
0.000000 0.000000 -1.000000
p =
0 1 0
1 0 0
0 0 1
p*l*u =
2.000000 1.000000 -2.000000
4.000000 -1.000000 2.000000
2.000000 -1.000000 1.000000
l =
-0.166667 -0.633333 1.000000
0.500000 1.000000 0.000000
1.000000 0.000000 0.000000
u =
6.000000 8.000000 9.000000
0.000000 -10.000000 -3.500000
0.000000 0.000000 5.283333
l*u =
-1.000000 5.000000 6.000000
3.000000 -6.000000 1.000000
6.000000 8.000000 9.000000
l =
1.000000 0.000000 0.000000
0.500000 1.000000 0.000000
-0.166667 -0.633333 1.000000
u =
6.000000 8.000000 9.000000
0.000000 -10.000000 -3.500000
0.000000 0.000000 5.283333
p =
0 0 1
0 1 0
1 0 0
p*l*u =
-1.000000 5.000000 6.000000
3.000000 -6.000000 1.000000
6.000000 8.000000 9.000000