// ************************************************************************** //
//                                                                            //
//    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 computes a LU decomposition of a given NxN matrix as   //
// described in Cormen, Leiserson, Rivest and Stein's textbook on algorithms. //
// Hence, for the given matrix a, the matrices l and u are computed so that   //
// the equation a = l*u holds, l is a lower left triangular, and u is a upper //
// right triangular matrix.                                                   //
// ************************************************************************** //

macro N = 4;

module DecompLU([N][N]real ?a,b,l,u,event !rdy) {
    
    // store input matrix a in matrix b
    for(i=0..N-1)
        for(j=0..N-1)
            b[i][j] = a[i][j];
    
    // compute LU decomposition of b
    for(k=0..N-1) {
        let(piv=b[k][k]) {
            u[k][k] = piv;
            for(i=k+1..N-1) {
                l[i][k] = b[i][k]/piv;
                u[k][i] = b[k][i];
            }
        }
        for(i=k+1..N-1) {
            for(j=k+1..N-1)
                next(b[i][j]) = b[i][j] - l[i][k] * u[k][j];
        }
        cp: pause;
    }
    for(i=0..N-1)
        l[i][i] = 1.0;
    emit(rdy);
}
drivenby {
    a[0][0] = 2.0;
    a[0][1] = 3.0;
    a[0][2] = 1.0;
    a[0][3] = 5.0;
    a[1][0] = 6.0;
    a[1][1] = 13.0;
    a[1][2] = 5.0;
    a[1][3] = 19.0;
    a[2][0] = 2.0;
    a[2][1] = 19.0;
    a[2][2] = 10.0;
    a[2][3] = 23.0;
    a[3][0] = 4.0;
    a[3][1] = 10.0;
    a[3][2] = 11.0;
    a[3][3] = 31.0;
    await(rdy);
    // check the result, i.e. that a = L * U holds, where L and U are stored in b   
    for(i=0..N-1)
        for(j=0..N-1)
            assert(a[i][j] == sum(k=0..N-1) (l[i][k] * u[k][j]));
}