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