zvookin
commented
8 years ago I don't have time to go into detail, but my guess is the best approach is to implement the Jacobi iteration as an external array function (the thing you get from define_extern) using the C code you already have. The inputs and outputs are small, but large enough to need to be stored in memory anyway, so the overhead of an external array function won't be too bad and you'll still have options re: scheduling the stages around this. It won't speed up the Jacobi iteration and it won't get you on a GPU however. (I might look at using Halide to implement a routine that does one step of the refinement on many matrices at once using vectorization and then calling it inside a loop from C++...)
Halide doesn't really have a way to iterate until convergence at present. Other folks may have some ideas, but it is a case where I really expect greater language power is needed to do this.
-Z-
dsharletg
SihamTabik
abadams
Dear all,
I would like to implement a 4-stage pipeline using Halide. The first two stages and the last stage can be expressed as pure functions of the input images. However, the middle stage needs to solve an eigenvalues & eigenvectors problem for each output-pixel using the jacobi iterative method (I copied the C-code below) i.e., To update each output pixel I need to calcute the eigenvalues of a symetric 3X3 matrix. The jacobi iterative method can take between 1 to 50 iterations to converge depending on the values of the 3x3 matrices. This means that the number of the needed arithmetic opertations varies drastically between different pixels.
My intention is to explore the full fusion of the 4-stages with Halide. Is it possible to express this stage only using Halide func an Expr?
Best regards,
the C-code of the jacobi iterative method is shown below: /----------------------------------------------------------------------------/ void jacobi3(float a[3][3], float d[3], float v[3][3]) /* This Jacobi routine Accepts and Returns the matrices in C convention. It also includes the ordering of the eigensystem, embedding the eigsrt routine from Numerical Recipes. This ordering is located just before the return. In principle, this routine only works for 3x3 matrices, but it can easily extended to any value of n if:
{
define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s_(h+g_tau);\
int j,iq,ip,i,k; float tresh,theta,tau,t,sm,s,h,g,c,p,b[3],z[3]; int n=3; int nrot;
/* Initialize to the identity matrix. */ for (ip=0;ip<n;ip++) { for (iq=0;iq<n;iq++) v[iq][ip]=0.0; v[ip][ip]=1.0; }
/* Initialize b and d to the diagonal of a. */ for (ip=0;ip<n;ip++) { b[ip]=d[ip]=a[ip][ip]; }
/* This vector accumulates terms tapq as in Eq.(11.1.14) in Numerical Recipes.*/ for (ip=0;ip<n;ip++) z[ip]=0.0;
/* Iterative procedure. */ nrot=0; for (i=0;i<50;i++) {
fprintf(stderr,"\nToo many iterations in routine jacobi: %i\n", nrot);
/* This is not the way by which the routine should exit. / / --- Reordering the eigenvalues/vectors. --- Taken from eigsrt routine from Numerical Recipes. It sorts the eigenvalues in d into descending order, and rearranges the eigenvectors v correspondingly. The method is straight insertion.*/ for (i=0;i<n;i++) { p=d[k=i]; for (j=i+1;j<n;j++) if (d[j] >= p) p=d[k=j]; if (k != i) { d[k]=d[i]; d[i]=p; for (j=0;j<n;j++) { p=v[i][j]; v[i][j]=v[k][j]; v[k][j]=p; } } }
undef ROTATE
} /---------------------------------------------------------------------------/