// ************************************************************************** //
//                                                                            //
//    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 Cholesky decomposition of a given symmetric //
// and positive definite NxN matrix. Matrix a is symmetric iff a[i][j]=a[j][i]//
// holds, and it is positive definite iff x*a*x>0 for all vectors x.          //
// The algorithm does neither check whether a is symmetric nor does it check  //
// whether it is positive definite. However, if a should be symmetric, but    //
// not positive definite, then the algorithm will fail since it will have     //
// to compute the square root of a negative number. Thus, the algorithm can   //
// also check that a given symmetric matrix is positive definite this way.    //
// Remark: The Cholesky decomposition is numerically more stable than LU      //
// or LUP decompositions. It computes for a given matrix a, the upper triangle//
// matrix l such that a=l^t*l holds. Using forward and backward substitutions,//
// one can thus also use it to solve linear equation systems.                 //
// ************************************************************************** //

macro N = 4;
macro sqrt(x) = exp(x,0.5);
macro equal(a,b) = (-0.0000001 <= a-b) & (a-b<=0.0000001);

module DecompCholesky([N][N]real ?a,l) {
    // compute Cholesky decomposition in one step
    l[0][0] = sqrt(a[0][0]);
    for(j=1..N-1)
        l[0][j] = a[0][j]/l[0][0];
    for(i=1..N-1) {
        l[i][i] = sqrt(a[i][i] - sum(k=0..i-1) (l[k][i] * l[k][i]));
        for(j=i+1..N-1)
            l[i][j] = (a[i][j] - sum(k=0..i-1) (l[k][i] * l[k][j]))/l[i][i]; 
        }
}
drivenby {
    a[0][0] = 9.0;
    a[0][1] = 3.0;
    a[0][2] = -6.0;
    a[0][3] = 12.0;
    a[1][0] = 3.0;
    a[1][1] = 26.0;
    a[1][2] = -7.0;
    a[1][3] = -11.0;
    a[2][0] = -6.0;
    a[2][1] = -7.0;
    a[2][2] = 9.0;
    a[2][3] = 7.0;
    a[3][0] = 12.0;
    a[3][1] = -11.0;
    a[3][2] = 7.0;
    a[3][3] = 65.0;
    // check correctness: a = transpose(l)*l
    for(i=0..N-1)
        for(j=0..N-1)
            assert(equal(a[i][j],sum(k=0 .. N-1) (l[k][i] * l[k][j])));
}