// ************************************************************************** //
//                                                                            //
//    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                                             //
//                                                                            //
// ************************************************************************** //
// The algorithm below improves module is a version of EigenvaluesShiftQR that//
// works on complex numbers.                                                  //
// ************************************************************************** //

macro N = 5;

// -------------------------------------------------------------------------- //
// numeric precision
// -------------------------------------------------------------------------- //
macro eps = 1.0e-8;
// -------------------------------------------------------------------------- //
// macros on real numbers
// -------------------------------------------------------------------------- //
macro rsign(x) = (x<0.0 ? -1.0 : +1.0);
macro rabs(x)  = (x<0.0 ? -x : +x);
macro rroot(x) = exp(x,0.5);
// -------------------------------------------------------------------------- //
// 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.1*x.1 + x.0*x.0);
macro cabs(x)  = rroot(x.1*x.1 + x.0*x.0);
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 EigenvaluesShiftQRC([N][N](real * real) ?m,a,[N](real * real) e,event !rdy) {
    [N][N](real * real) q,r;                    // matrices q,r of QR-decomposition
    (real * real) trc2,det,dis,rtdis,e1,e2,p;   // to determine shift parameter p
    real dist1,dist2;   // distances of eigenvectors for Wilkinson shift
    nat{N+1} n;         // current size of the matrix (due to deflation)
    event QRrdy;        // whether QR-decomposition terminated

    // initialize
    n = N;
    for(i=0..N-1)
        for(j=0..N-1)
            a[i][j] = m[i][j];

    do {
        do {

            // determine Wilkinson shift p: the idea is to compute the eigenvalues
            // of the matrix [[a[N-2][N-2],a[N-2][N-1]],[a[N-1][N-2],a[N-1][N-1]]]
            // and to select that eigenvalue that is closer to a[N-1][N-1] as shift p
            // note that this is the only step where complex numbers may be introduced
            // for real matrices, otherwise the QR-iteration would fail for real matrices
            // that have complex eigenvalues

            let(a00 = a[n-2][n-2])
            let(a01 = a[n-2][n-1])
            let(a10 = a[n-1][n-2])
            let(a11 = a[n-1][n-1]) {
            trc2 = cdiv(cadd(a00,a11),(2.0,0.0));
            det = csub(cmul(a00,a11),cmul(a01,a10));
            dis = csub(cmul(trc2,trc2),det);
            // calling croot directly caused inaccuracies, so the following is better
            rtdis = (dis.1==0.0 ? (dis.0<0.0 ? (0.0,rroot(-dis.0)) : (rroot(dis.0),0.0))
                                : croot(dis));
            e1 = cadd(trc2,rtdis);
            e2 = csub(trc2,rtdis);
            dist1 = sqabs(csub(e1,a11));
            dist2 = sqabs(csub(e2,a11));
            p = (dist1<dist2 ? e1 : e2);
            }

            // shift matrix a
            for(i=0..N-1)
                next(a[i][i]) = csub(a[i][i],p);
            pause;

            // decompose shifted matrix A_k-p*E into factors Q_k and R_k
            DecompQRC(a,q,r,QRrdy);

            // compute A_{k+1} := R_k*Q_k+p_k*E and check for convergence
            for(i=0..N-1) {
                for(j=0..N-1) {
                    let(p0 = (i==j ? p : zero))
                    let(new_aij_re = Re(p0) + sum(k=0..N-1) Re(cmul(r[i][k],q[k][j])))
                    let(new_aij_im = Im(p0) + sum(k=0..N-1) Im(cmul(r[i][k],q[k][j])))
                    next(a[i][j]) = (new_aij_re,new_aij_im);
                }
            }
            pause;
        } while(sqabs(a[n-1][n-2])>eps*eps);
        next(n) = n-1;
        // eigenvalue e[n-1] has been found, we can deflate the matrix
        e[n-1] = a[n-1][n-1];
        for(i=0..N-1) {
            next(a[i][n-1]) = zero;
            next(a[n-1][i]) = zero;
        }
        pause;
    } while(n>1);

    e[n-1] = a[n-1][n-1];
    emit(rdy);
}
drivenby m0 {
    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 m1 {
    // note: x^5 - 15x^4 + 85x^3 - 225x^2 + 274x - 120 = (x-1)(x-2)(x-3)(x-4)(x-5)
    m[0] = [ (15.0,0.0), (-85.0,0.0), (225.0,0.0), (-274.0,0.0), (120.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);
}
drivenby m2 {
    // note: x^5-2x^4-3x^3+6x^2+2x-4 = (x-2)(x-1)(x+1)(x-sqrt(2))(x+sqrt(2))
    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);
}
drivenby m3 {
    // note: x^5+6x^4+14x^3+24x^2+33x+18 = (x+3)(x+2)(x+1)(x^2+3)
    m[0] = [ (-6.0,0.0), (-14.0,0.0), (-24.0,0.0), (-33.0,0.0), (-18.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);
}
drivenby m4 {
    // note: x^5-2x^4+2x^3-4x^2-3x+6 = (x-2)(x-1)(x+1)(x^2+3)
    m[0] = [ (2.0,0.0), (-2.0,0.0), (4.0,0.0), (3.0,0.0), (-6.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);
}