#include <numeric.h> int main() { int m = 3; array double complex a[3][3] = { complex(2,0), complex(0,-1), complex(0,0), complex(0,1), complex(2,0), complex(0,0), complex(0,0), complex(0,0), complex(3,0)}; int status; array double complex l[m][m]; status = choldecomp(a, l); if (status == 0) { printf("upper triangle l =\n%5.3f\n", l); printf("transpose(l)*l - a =\n%5.3f\n", conj(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%5.3f\n", l); printf("l*transpose(l) - a =\n%5.3f\n", l*conj(transpose(l))-a); } else printf("error: Matrix must be positive definite.\n" "The %d order of matrix is not positive definite.\n",status); }
upper triangle l = complex(1.414,0.000) complex(0.000,-0.707) complex(0.000,0.000) complex(0.000,0.000) complex(1.225,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(1.732,0.000) transpose(l)*l - a = complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(-0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(-0.000,0.000) lower triangle l = complex(1.414,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.707) complex(1.225,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(1.732,0.000) l*transpose(l) - a = complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(-0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(0.000,0.000) complex(-0.000,0.000)