/*  -*- MaTX -*-
 *
 *  【名前】
 *      chol() - コレスキー分解
 *
 *  【形式】
 *      R = chol(X)
 *         Matrix R;
 *         Matrix X;
 *
 *  【機能説明】
 *      R = chol(X)は，対称行列 X のコレスキー解を求める。もし，X が
 *      正定なら，R' * R = X を満たす上三角行列 R を求める。もし，X
 *      が正定でなければ，エラーになる。
 *
 *  【例題】
 *      X = [[2 2 4][2 5 6][4 6 10];
 *      R = chol(X)
 *      === [R] : (  3,  3) ===
 *                 (  1)           (  2)           (  3)     
 *      (  1)  1.41421356E+00  1.41421356E+00  2.82842712E+00
 *      (  2)  0.00000000E+00  1.73205081E+00  1.15470054E+00
 *      (  3)  0.00000000E+00  0.00000000E+00  8.16496581E-01
 *
 *      print R'*R
 *      === [ans] : (  3,  3) ===
 *                 (  1)           (  2)           (  3)     
 *      (  1)  2.00000000E+00  2.00000000E+00  4.00000000E+00
 *      (  2)  2.00000000E+00  5.00000000E+00  6.00000000E+00
 *      (  3)  4.00000000E+00  6.00000000E+00  1.00000000E+01
 *
 *  【関連項目】
 *     
 */

Func Matrix chol(X)
	Matrix X;
{
	Integer n, j;
	Matrix R;
	Real tol;

	tol = frobnorm(X)*EPS;
	n = Cols(X);

	if (max(abs(X - X#)) > tol) {
		printf("max(abs(Array(X - X#))) = %g\n", max(abs(X - X#)));
		warning("chol(): 対称行列かつエルミート行列でない。\n");
	}

	R = X;
	for (j = 1; j <= n; j++) {
		if (j > 1) {
			R(j:n,j) = R(j:n,j) - R(j:n,1:j-1)*R(j,1:j-1)#;
		}
		if (R(j,j) < 0.0) {
			error("chol(): 正定行列でない。\n");
		}
		R(j:n,j) = R(j:n,j)/sqrt(R(j,j));
	}

	for (j = 1; j <= n-1; j++) {
		R(j,j+1:) = Z(1,n-j);
	}
	
	if (max(abs(R*R# - X)) > tol) {
		printf("max(abs(R*R# - X)) = %g\n", max(abs(R*R# - X)));
		warning("chol(): 分解は正確でないかもしれない。\n");
	}

	return R#;
}