Choleski decomposition

A simple description of the algorithm in a sequential code may be,

  on(Adlib.global) {
    Range x = new CyclicRange(size, Adlib.global.dim(0));

    float a[[,#]] = new float [[size, x]];
    // initialize the array here;

    float b[[]] = new float [[size]]; // used as a buffer

    Location j;

    at(j=x|0)
      a[0,j]=sqrt(a[0,j]);

    for(int k=0; k<size-1; k++) {
      for(int s=k+1; s<size; s++)
        at(j=x|k)
          a[s,j]/=a[j,j];

      Adlib.remap(b[[k+1:]],a[[k,x|k+1:]]);

      over (j=x|k+1:)
        for (int i=x; i<size; i++)
          a[i,j]-=b[i]*b[j];

      at(j=x|k+1)
        a[k+1,j]=sqrt(a[k+1,j]);
    }
  }