// ************************************************************************** //
//                                                                            //
//    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. In contrast to module DecompQR that works with reals, the//
// version below deals with complex numbers.                                  //
// ************************************************************************** //

macro N = 5;

// -------------------------------------------------------------------------- //
// numeric precision
// -------------------------------------------------------------------------- //
macro eps = 1.0e-8;
// -------------------------------------------------------------------------- //
// macros on real numbers
// -------------------------------------------------------------------------- //
macro rroot(x) = exp(x,0.5);
macro rsign(x) = (x<0 ? -1 : +1);
// -------------------------------------------------------------------------- //
// arithmetic on complex numbers
// -------------------------------------------------------------------------- //
macro cadd(x,y) = (x.0+y.0,x.1+y.1);
macro csub(x,y) = (x.0-y.0,x.1-y.1);
macro cmul(x,y) = (x.0*y.0-x.1*y.1,x.1*y.0+x.0*y.1);
macro cdiv(x,y) = ((x.0*y.0+x.1*y.1)/(y.0*y.0+y.1*y.1),
                   (x.1*y.0-x.0*y.1)/(y.0*y.0+y.1*y.1));
macro sqabs(x) = (x.0*x.0 + x.1*x.1);
macro cabs(x)  = rroot(x.0*x.0 + x.1*x.1);
macro croot(x) = (rsign(x.1) * rroot(0.5 * (cabs(x) + x.0)),
                               rroot(0.5 * (cabs(x) - x.0)));
macro Re(x) = x.0;
macro Im(x) = x.1;
macro unit = (1.0,0.0);
macro zero = (0.0,0.0);
// -------------------------------------------------------------------------- //


module DecompQRC([N][N](real * real) ?m,q,r,event !rdy) {
    (real * real) a,b,Rq,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 ? unit : zero);
            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]!=zero) {
                i = j-1;
                a = r[i][k];
                b = r[j][k];
                Rq = cadd(cmul(a,a),cmul(b,b));
                R = (Rq.1==0 ? (Rq.0>=0.0 ? (rroot(Rq.0),0.0) : (0.0,rroot(Rq.0)))
                             : croot(Rq));
                c = cdiv(a,R);
                s = cdiv(b,R);
                // apply Givens rotation to q and r
                for(l=0..N-1) {
                    next(q[i][l]) = cadd(cmul(c,q[i][l]),cmul(s,q[j][l]));
                    next(q[j][l]) = csub(cmul(c,q[j][l]),cmul(s,q[i][l]));
                    next(r[i][l]) = cadd(cmul(c,r[i][l]),cmul(s,r[j][l]));
                    next(r[j][l]) = (l==k ? zero : csub(cmul(c,r[j][l]),cmul(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,0.0),(2.0,0.0),(3.0,0.0),(4.0,0.0),(5.0,0.0)];
    m[1] = [(2.0,0.0),(3.0,0.0),(4.0,0.0),(5.0,0.0),(1.0,0.0)];
    m[2] = [(3.0,0.0),(4.0,0.0),(5.0,0.0),(1.0,0.0),(2.0,0.0)];
    m[3] = [(4.0,0.0),(5.0,0.0),(1.0,0.0),(2.0,0.0),(3.0,0.0)];
    m[4] = [(5.0,0.0),(1.0,0.0),(2.0,0.0),(3.0,0.0),(4.0,0.0)];
    await(rdy);
}
drivenby {
    m[0] = [ (2.0,0.0), (3.0,0.0), (-6.0,0.0), (-2.0,0.0), (4.0,0.0)];
    m[1] = [ unit,  zero,  zero,  zero,  zero];
    m[2] = [ zero,  unit,  zero,  zero,  zero];
    m[3] = [ zero,  zero,  unit,  zero,  zero];
    m[4] = [ zero,  zero,  zero,  unit,  zero];
    await(rdy);
}