//
// CholeskyDemo 2014: DJ Greaves. 
// (Based on http://rosettacode.org/wiki/Cholesky_decomposition).
// 
public class CholeskyDemo
{
    

    public static double[][] Cholesky(double[][] Adata)
    {
	int m = Adata.length;
	double[][] R = new double[m][m];//Default 0 for unassigned upper triangle.
	for(int i = 0; i<m ;i++)
	    for(int j = 0; j < i+1; j++)//Fill in diagonal and below only.
		{
		    double sum = 0.0;
		    for(int k = 0; k < j; k++) sum += R[i][k] * R[j][k];
		    R[i][j] = (i == j) ? Math.sqrt(Adata[j][j] - sum):
                                         (Adata[i][j] - sum)/R[j][j];
		}
	return R;
    }
    

    public static void test(double [][] Adata)
    {
	System.out.printf("Input data:\n");
	printa(Adata);
	double[][] c = Cholesky(Adata);
	System.out.printf("Cholesky:\n");
	printa(c);
	System.out.printf("Back subst:\n");
	printa(mpx(c,transpose(c)));

    }

    public static void main(String[] args)
    {
	double[][] test1 = {{25, 15, -5},
			    {15, 18, 0},
			    {-5, 0, 11}};

	double[][] test2 = {{18, 22, 54, 42},
			    {22, 70, 86, 62},
			    {54, 86, 174, 134},
			    {42, 62, 134, 106}};

	test(test1);
	test(test2);
    }





    // Matrix support functions now follow (there is similar code in java.utils).

    public static void printa (double [][] A)
    {

	for (int i=0; i<A.length; i++)
	    {
		for (int j=0; j<A[i].length; j++)
		    {
			if (A[i][j]==0.0) System.out.printf("-.----- ", A[i][j]);
			else System.out.printf("%1.5f ", A[i][j]);
		    }
		System.out.printf("\n");
	    }
	System.out.printf("\n");
    }

    public static double[][] mpx(double[][] AA, double[][] BB)
    {
	double[][] Ans = new double[AA.length][BB[0].length];
	//System.out.printf(" mpx %d,%d with %d,%d\n", AA.length, AA[0].length, BB.length, BB[0].length);	
	assert(AA[0].length == BB.length);
	for (int i=0; i<AA.length; i++)
	    for (int k=0; k<BB[0].length; k++)
		{
		    double sum = 0.0;
		    for (int j=0; j<AA[0].length; j++) sum += AA[i][j] * BB[j][k];
		    Ans[i][k] = sum;
		}
	return Ans;
    }

    public static double[][] transpose(double[][] AA)
    {
	double[][] Ans = new double[AA[0].length][AA.length];
	for (int i=0; i<AA.length; i++)
	    for (int j=0; j<AA[0].length; j++) Ans[j][i] = AA[i][j];
	return Ans;
    }

    public static double [] [] copy2d(double [] [] matrix)
    {
	double [][] copy = new  double[matrix.length][];
	for(int i = 0; i < matrix.length; i++) copy[i] = matrix[i].clone();
	return copy;
    }

}

// eof