// Instructions:
// Compile and execute

public class QRExample {

    static double[][] Q, R;

    public static void main (String[] argv)
    {
	double[][] A = {
	    {2,1},
	    {1,2}
	};

	int numIterations = 10;

	for (int k=0; k<numIterations; k++) {

	    // Step 1: QR decomposition
	    computeQR (A);

	    // Step 2: switch order and multiply:
	    A = MatrixTool.matrixMult (R, Q);

	    MatrixTool.print (A);
	}

    }

    static void computeQR (double[][] A)
    {
	// For a 2x2, this is just a single projection.
	// We'll use the same notation as in Module 7, using
	// indexing from 0.

	// Make space:
        Q = new double [2][2];
	R = new double [2][2];

	// First column of Q is the first column of A, normalized.
	double[] w0 = new double [2];
	w0[0] = A[0][0];  w0[1] = A[1][0];
	double[] v0 = w0;
	double[] u0 = normalize (v0);
	Q[0][0] = u0[0];
	Q[1][0] = u0[1];
	R[0][0] = MatrixTool.norm (v0);

	// Now for step 2 of the modified Gram-Schmidt algorithm:
	double[] w1 = new double [2];
	w1[0] = A[0][1];  w1[1] = A[1][1];
	double c01 = MatrixTool.dotProduct(w1,v0) / MatrixTool.dotProduct(v0,v0);
	double[] temp = MatrixTool.scalarMult (c01, v0);
	double[] v1 = MatrixTool.sub (w1, temp);
	double[] u1 = normalize (v1);
	Q[0][1] = u1[0];  Q[1][1] = u1[1];
	R[1][1] = MatrixTool.norm (v1);
	R[0][1] = c01 * MatrixTool.norm (v0);
    }

    static double[] normalize (double[] x)
    {
	double normX = MatrixTool.norm (x);
	double[] y = new double [x.length];
	for (int i=0; i<x.length; i++) {
	    y[i] = (1.0 / normX) * x[i];
	}
	return y;
    }

}