// ************************************************************************** //
//                                                                            //
//    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 algorithm solves a linear equation system based on LU decomposition.  //
// ************************************************************************** //

import Matrices;
macro N = 3;

module LinearEquLU([N][N]real ?a,[N]real ?b,[N]real x,event !rdy) {
    [N][N]real c,d;
    [N]real y;

    // compute LU decomposition in matrix c of matrix a (d stores a)
    DecompLU2(a,d,c,rdy);
    
    // having the LU decomposition, we first compute y as solution of L*y = b
    y[0] = b[0];
    for(i=1..N-1)
        y[i] = b[i] - sum(j=0..i-1) (c[i][j] * y[j]);

    // having y, we now compute the solution x of U*x = y, i.e., L*U*x = L*y = b
    for(i=0..N-1)
        x[i] = (y[i] - sum(j=i+1..N-1) (c[i][j] * x[j])) / c[i][i];
}
drivenby {
    a[0][0] = 1.0;
    a[0][1] = 2.0;
    a[0][2] = 0.0;
    a[1][0] = 3.0;
    a[1][1] = 4.0;
    a[1][2] = 4.0;
    a[2][0] = 5.0;
    a[2][1] = 6.0;
    a[2][2] = 3.0;
    b[0] = 3.0;
    b[1] = 7.0;
    b[2] = 8.0;
    await(rdy);
    // check the result, i.e. that a*x=b holds up to numeric inaccuracies
    for(i=0..N-1) {
        assert(-0.00001<=b[i] - sum(k=0..N-1) (a[i][k] * x[k]));
        assert(b[i] - sum(k=0..N-1) (a[i][k] * x[k])<=0.00001);
    }
}