// ************************************************************************** //
//                                                                            //
//    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 EigenvaluesQR in that shift operations //
// are used to speed-up the QR-iteration. The idea is as follows: instead of  //
// using A_0 := A and A_{k+1} := R_k*Q_k where A_k = Q_k * R_k is obtained by //
// the QR-decomposition of A_k, we use A_0 := A and A_{k+1} := R_k*Q_k+p_k*E  //
// where A_k-p_k*E=Q_k*R_k is obtained by the QR-decomposition of A_k-p_k*E.  //
// It is then easily seen that A_{k+1} = Q_k^T*A_k*Q_k so that all A_k have   //
// the same eigenvalues, regardless which value is chosen for p_k. It can be  //
// proved that the convergence of A[i][j]_k is linear with |e[i]/e[j]|^k. Now,//
// note that A_k-p_k*E has the eigenvalues e[i]-p_k so that A[i][j]_k will    //
// converge as |(e[i]-p_k)/(e[j]-p_k)|^k. Thus, one should choose p_k close to//
// e[i+1], but since the latter is not yet known, one has to guess.           //
// In particular, one often chooses A[i][i], but it is known that this does   //
// not work with complex roots. Wilkinson therefore proposed to use the       //
// eigenvalue of the 2x2 matrix [[A[i][i],A[i][i+1]],[A[i][i+1],A[i+1][i+1]]] //
// that is closer to A[i+1][i+1] (note that this eigenvalue may be complex).  //
// ************************************************************************** //

macro N = 5;

// numeric precision
macro rabs(x) = (x<0 ? -x : +x);
macro sqrt(x) = exp(x,0.5);
macro eps = 1.0e-10;


module EigenvaluesShiftQR([N][N]real ?m,a,[N]real e,event !rdy) {
    [N][N]real q,r; // matrices q,r of QR-decomposition
    real p;         // shift parameter
    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 shift p: note that this simple approximation to e[N-1]
            // works well in case the latter is a real number; otherwise the 
            // QR-iteration will not converge
            p = a[n-1][n-1];
      
            // shift matrix a
            for(i=0..N-1)
                next(a[i][i]) = a[i][i]-p;
            pause;
      
            // decompose shifted matrix A_k-p*E into factors Q_k and R_k
            DecompQR(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 : 0.0))
                    next(a[i][j]) = p0 + sum(k=0..N-1) (r[i][k]*q[k][j]);
                }
            }
            pause;
        } while(rabs(a[n-1][n-2])>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]) = 0.0;
            next(a[n-1][i]) = 0.0;
        }
        pause;
    } while(n>1);

    e[n-1] = a[n-1][n-1];
    emit(rdy);
}
drivenby m0 {
    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);
}
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, -85.0, 225.0, -274.0, 120.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);
}
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,   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);
}
// m3 does not converge since there are complex eigenvalues
//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, -14.0, -24.0,  -33.0, -18.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);
//}
// m4 does not converge since there are complex eigenvalues
//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,  -2.0,   4.0,    3.0,  -6.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);
//}