// ************************************************************************** //
//                                                                            //
//    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 approximates the eigenvalues of a given NxN matrix by  //
// means of the QR-decomposition: Given a matrix A, we define the iteration   //
// 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. Since Q_k is orthogonal, its inverse matrix is its//
// transposed matrix, and therefore A_{k+1} = (Q_k^t*Q_k)*(R_k*Q_k) holds,    //
// which is equivalent to A_{k+1} = Q_k^t*A_k*Q_k. For this reason, all A_k   //
// have the same eigenvalues, since e*x = A_{k+1}*x = Q_k^t*A_k*Q_k*x, i.e.   //
// equivalent to e*(Q_k*x) = A_k*(Q_k*x). In case, the sequence A_k converges //
// to a matrix A', it then follows that also Q_k and R_k converge to matrices //
// Q' and R' so that A' = Q' * R' = R' * Q' holds, where Q' is orthogonal and //
// R' is an upper triangular matrix. A' is the a Schur matrix, i.e. a matrix  //
// whose entries are all zero with the exception of diagonal elements A[i][i] //
// and some of the elements A[i][i+1] and A[i+1][i] so that the diagonal of A //
// consists of the eigenvalues and 2x2 matrices that represent two complex    //
// eigenvalues. The simple QR-iteration as implemented here does however not  //
// converge if there are complex eigenvalues, and if there are eigenvalues of //
// the same absolute value. To solve the latter, we would have to deflate the //
// matrix, and to deal with complex eigenvalues, we have to use Wilkinson     //
// shifts with complex numbers (or use the double step QR-iteration).         //
// ************************************************************************** //

macro N = 5;

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


module EigenvaluesQR([N][N]real ?m,a,event !rdy) {
    [N][N]real q,r;
    event QRrdy;
    event iterate;

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

    do {
        // decompose A_k = Q_k * R_k
        DecompQR(a,q,r,QRrdy);
        // determine A_{k+1} = R_k * Q_k and check convergence
        for(i=0..N-1) {
            for(j=0..N-1) {
                let(new_aij = sum(k=0..N-1) (r[i][k]*q[k][j])) {
                next(a[i][j]) = new_aij;
                if(i==j & rabs(a[i][j]-new_aij)>eps)
                    next(iterate) = true;
                }
            }
        }
        pause;
        next(iterate) = false;
    } while(iterate);

    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);
}
drivenby {
    // 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);
}