#include <numeric.h> int main() { int m = 5; array double a[5][5] = { 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 1, 3, 6,10,15, 1, 4, 10,20,35, 1, 5, 15,35,70}; array double a1[5][5] = { -1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 1, 3, 6,10,15, 1, 4, 10,20,35, 1, 5, 15,35,70}; int status; array double l[m][m]; status = choldecomp(a, l); if (status == 0) { printf("upper triangle l =\n%f\n", l); printf("transpose(l)*l - a =\n%f\n", transpose(l)*l-a); } else printf("error: Matrix must be positive definite.\n" "The %d order of matrix is not positive definite.\n",status); status = choldecomp(a, l, 'l'); if (status == 0) { printf("lower triangle l =\n%f\n", l); printf("l*transpose(l) - a =\n%f\n", l*transpose(l)-a); } else printf("error: Matrix must be positive definite.\n" "The %d order of matrix is not positive definite.\n",status); status = choldecomp(a1, l); if (status == 0) { printf("l =\n%f\n", l); } else { printf("a =\n%f", a1); printf("error: Matrix must be positive definite.\n" "The %d order of matrix is not positive definite.\n",status); } }
upper triangle l = 1.000000 1.000000 1.000000 1.000000 1.000000 0.000000 1.000000 2.000000 3.000000 4.000000 0.000000 0.000000 1.000000 3.000000 6.000000 0.000000 0.000000 0.000000 1.000000 4.000000 0.000000 0.000000 0.000000 0.000000 1.000000 transpose(l)*l - a = 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 lower triangle l = 1.000000 0.000000 0.000000 0.000000 0.000000 1.000000 1.000000 0.000000 0.000000 0.000000 1.000000 2.000000 1.000000 0.000000 0.000000 1.000000 3.000000 3.000000 1.000000 0.000000 1.000000 4.000000 6.000000 4.000000 1.000000 l*transpose(l) - a = 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 a = -1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 2.000000 3.000000 4.000000 5.000000 1.000000 3.000000 6.000000 10.000000 15.000000 1.000000 4.000000 10.000000 20.000000 35.000000 1.000000 5.000000 15.000000 35.000000 70.000000 error: Matrix must be positive definite. The 1 order of matrix is not positive definite.