// ************************************************************************** //
//                                                                            //
//    eses                   eses                                             //
//   eses                     eses                                            //
//  eses    eseses  esesese    eses   Embedded Systems Group                  //
//  ese    ese  ese ese         ese                                           //
//  ese    eseseses eseseses    ese   Department of Computer Science          //
//  eses   eses          ese   eses                                           //
//   eses   eseses  eseseses  eses    University of Kaiserslautern            //
//    eses                   eses                                             //
//                                                                            //
// ************************************************************************** //
// This module computes a QR-decomposition of a given NxN matrix by means of  //
// Givens rotations. Alternative methods are based on Householder reductions  //
// or the Gram-Schmidt procedure to determine orthogonal matrices. In general,//
// every matrix M can be written as M = Q*R with a orthogonal matrix Q and a  //
// upper right triangular matrix R. Since Q is orthogonal, all its columns are//
// orthogonal to each other and their lengths are 1. Therefore, the transpose //
// of Q is the inverse of Q.                                                  //
//   A Givens rotation on i,j with angle phi is a linear function f(x) = G*x  //
// that maps a N-vector x to a N-vector f(x) such that x is rotated in the    //
// (i,j)-plane with angle phi. The rotation can be described by a NxN matrix  //
// G that is the unit matrix with the following four exceptions:              //
//  * G[i][i] = +cos(phi)                                                     //
//  * G[i][j] = +sin(phi)                                                     //
//  * G[j][i] = -sin(phi)                                                     //
//  * G[j][j] = +cos(phi).                                                    //
//                                                                            //
// These matrices are also called Jacobi matrices. The vector y := f(x) is    //
// therefore only changed at positions i and j as follows:                    //
//  * y[i] = cos(phi)*x[i] + sin(phi)*x[j]                                    //
//  * y[j] = cos(phi)*x[i] - sin(phi)*x[j]                                    //
//  * y[k] = x[k] for k!=i and k!=j                                           //
//                                                                            //
// It is not necessary to compute phi and then cos(phi) and sin(phi), instead //
// it is sufficient to directly compute cos(phi) and sin(phi). If the goal is //
// to zero out an element b of a vector, then we can reduce the consideration //
// to two rows:                                                               //
//                                                                            //
//  | c   s|   |a|   |R|                                                      //
//  |      | * | | = | |                                                      //
//  |-s   c|   |b|   |0|                                                      //
//                                                                            //
// Clearly, R = sqrt(a*a+b*b) is the length of the vector (a,b) and we have   //
// c = cos(phi) = a/R and s = sin(phi) = b/R.                                //
//                                                                            //
// A QR-decomposition of a given matrix m is then obtained by Givens rotations//
// R := G_0 * G_1 * ... * G_L * m such that all sub-diagonal elements of m are//
// zeroed out by the rotations. The matrix Q is the inverse of (G_0 *…*G_L),  //
// which can be computed as Q = transpose(G_0 * ... * G_L).                   //
// ************************************************************************** //

macro N = 5;
macro sqrt(x) = exp(x,0.5);
macro rabs(x) = (x<0 ? -x : +x);
macro eps = 1.0e-8;


module DecompQR([N][N]real ?m,q,r,event !rdy) {
    real a,b,R,s,c;
    nat{N} i;

    // initially, q is the unit matrix and r is m
    for(i=0..N-1)
        for(j=0..N-1) {
            q[i][j] = (i==j ? 1.0 : 0.0);
            r[i][j] = m[i][j];
            }
    pause;

    // next apply Givens rotations to zero out sub-diagonal elements of r
    for(k=0..N-2) // columns
        for(jup=0..N-2-k)
        let(j=N-1-jup) { // rows j=N-1,N-2,…,k+1
            if(r[j][k]!=0.0) {
                i = j-1;
                a = r[i][k];
                b = r[j][k];
                R = sqrt(a*a+b*b);
                c = a/R;
                s = b/R;
                // apply Givens rotation to q and r
                for(l=0..N-1) {
                    next(q[i][l]) = c * q[i][l] + s * q[j][l];
                    next(q[j][l]) = c * q[j][l] - s * q[i][l];
                    next(r[i][l]) = c * r[i][l] + s * r[j][l];
                    next(r[j][l]) = (l==k ? 0.0 : c * r[j][l] - s * r[i][l]);
                }
            pause;
            }
        }

    // finally, compute the transpose of q
    for(i=0..N-1)
        for(j=0..N-1)
            next(q[i][j]) = q[j][i];
    pause;
    emit(rdy);
}
drivenby {
    m[0] = [1.0,2.0,3.0,4.0,5.0];
    m[1] = [2.0,3.0,4.0,5.0,1.0];
    m[2] = [3.0,4.0,5.0,1.0,2.0];
    m[3] = [4.0,5.0,1.0,2.0,3.0];
    m[4] = [5.0,1.0,2.0,3.0,4.0];
    await(rdy);
    // check correctness:
    for(i=0..N-1)
        for(j=0..N-1)
            assert(rabs(m[i][j]-sum(k=0..N-1)(q[i][k]*r[k][j]))<eps);
}
drivenby {
    m[0] = [ 2.0,  3.0, -6.0, -2.0, +4.0];
    m[1] = [ 1.0,  0.0,  0.0,  0.0,  0.0];
    m[2] = [ 0.0,  1.0,  0.0,  0.0,  0.0];
    m[3] = [ 0.0,  0.0,  1.0,  0.0,  0.0];
    m[4] = [ 0.0,  0.0,  0.0,  1.0,  0.0];
    await(rdy);
    // check correctness:
    for(i=0..N-1)
        for(j=0..N-1)
            assert(rabs(m[i][j]-sum(k=0..N-1)(q[i][k]*r[k][j]))<eps);
}