<![CDATA[
#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);
}
}
]]>
<![CDATA[
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.
]]>