/**
 * Gram-Schmidt-Verfahren zur Orthogonalisierung und QR-Zerlegung
 */
public class GS {

    /*
     * Testmatrizen
     */
    private static double[][] A = {
        { 0, 1 },
	{ 1, 1 },
	{ 2, 1 },
	{ 3, 1 },
	{ 4, 1 },
	{ 5, 1 }
    };

    private static double[][] B = {
        {  2.0,   4.0, -9.0 },
        { -2.0, -10.0,  3.0 },
        {  1.0,  -1.0,  6.0 }
    };

    /*
     * Hole Spaltenvektor j aus Matrix a
     */
    public static double[] colGet(double[][] a, int j) {
        int n = a.length;
	double[] x = new double[n];

	for (int i=0 ; i<n ; i++) {
	    x[i] = a[i][j]; 
	}

	return x;
    }

    /*
     * Lege Vektor x in Spalte j von Matrix a ab.
     */
    public static void colSet(double[][] a, int j, double[] x) {
        int n = x.length;

	for (int i=0 ; i<n ; i++) {
	    a[i][j] = x[i]; 
	}
    }

    /*
     * Skaliere Vektor x mit Skalar alpha
     */
    public static double[] scale(double[] x, double alpha) {
        int n = x.length;
	double[] y = new double[n];

	for (int i=0 ; i<n ; i++) {
	    y[i] = alpha * x[i];
	}

	return y;
    }

    /*
     * Vektorsubtraktion
     */
    public static double[] minus(double[] x, double[] y) {
        int n = x.length;
	double[] z = new double[n];

	for (int i=0 ; i<n ; i++) {
	    z[i] = x[i] - y[i];
	}

	return z;
    }

    /*
     * Berechne die euklidischen Norm fuer einen Vektor x
     */
    public static double norm2(double[] x) {
        int n = x.length;
	double norm = 0.0;

	for (int i=0 ; i<n ; i++) {
	    norm = norm + x[i]*x[i];
	}

	return Math.sqrt(norm);
    }

    /*
     * Berechne Skalarprodukt von zwei Vektoren x und y
     */
    public static double dot(double[] x, double[] y) {
        int n = x.length;
	double dp = 0.0;

	for (int i=0 ; i<n ; i++) {
	    dp = dp + x[i]*y[i];
	}

	return dp;
    }

    /*
     * Ausgabe einer n x m-Matrix
     */
    public static void print(double[][] a) {
        int n = a.length;
	int m = a[0].length;

	for (int i=0 ; i<n ; i++) {
	    for (int j=0 ; j<m ; j++) {
	        System.out.printf("%10.6f ", a[i][j]);
	    }
	    System.out.println();
	}
	System.out.println();
    }

    /*
     * Berechnung einer reduzierten QR-Zerlegung einer Matrix A mit dem
     * Gram-Schmidt-Verfahren
     *
     * Input: n x m-Matrix A = (a_ij)
     *
     * Output: n x m- Matrix Q = (q_ij) mit orthonormierten Spalten
     *         im x m-obere Dreicksmatrix R = (r_ij)
     *
     * Hinweise: 
     *     - q und r muessen in passenden Groessen vorliegen.
     *     - Die untere Haelfte von r wird nicht veraendert, muss also
     *       mit 0en vorbelegt sein.
     *     - Wenn A eine quadratische Matrix ist, bilden Q und R auch eine
     *       volle QR-Zerlegung.
     */
    public static void gramSchmidt(double[][] a, double [][] q, double[][] r) {
	int m = a[0].length;			// m = #Spalten von A
        double[] a0 = colGet(a,0);		// hole erste Spalte von A
        r[0][0] = norm2(a0);			// r_1,1 := ||a^(1)||
	double[] q0 = scale(a0,1/r[0][0]);	// q^(1) := 1/r_1,1 * a^(1)
	colSet(q,0,q0);				// setze erste Spalte in Q
	for (int j=1 ; j<m ; j++) {		// fuer j:=2 bis m
	    double[] aj = colGet(a,j);		// hole j-te Spalte aus A
	    double[] qj = colGet(a,j);		// q^(j) := a^(j)
	    for (int i=0 ; i<j ; i++) {		// fuer i:=1 bis j-1
		double[] qi = colGet(q,i);	// r_i,j := <q^(i),a^(j)>
	        r[i][j] = dot(qi,aj);		// q^(j) := q^(j) - r_i,j*q^(i)
		qj = minus(qj,scale(qi,r[i][j]));
	    }
	    r[j][j] = norm2(qj);		// r_j,j := ||q^(j)||
	    qj = scale(qj,1/r[j][j]);		// q^(j) := 1/r_j,j * q^(j)
	    colSet(q,j,qj);			// kopiere Vektor in die Matrix q
	}
    }

    public static void main(String[] arg) {
        double[][] a = B;			// hier evtl. Testmatrix austauschen
	int n = a.length;			// Anzahl Zeilen von A
	int m = a[0].length;			// Anzahl Spalten von A
	double[][] q = new double[n][m];	// Matrizen in passender Groesse anlegen
	double[][] r = new double[m][m];

	gramSchmidt(a,q,r);			// QR-Zerlegung berechnen
	print(q);				// und Matrizen ausgeben
	print(r);
    }
}