I'm working on a small fun project on the side for a numerical method I've been working a bit with. Roughly, the computational problem I have to solve is the following: Assume you have a collection of $N$ points, call this collection $X$. For each $x_i \in X$, the following steps occur:
1: Find $k$ nearest neighbors to $x_i$, call this set $X_i$.
2: Build a $k \times k$ matrix out of these points with some simple arithmetic. Call this matrix $A^i$. $A^i(j,k)$ depends only on points $y_j,y_k \in X_i$.
3: Build a vector $b^i$ who is 1 only in the row corresponding to the originally chosen point $x^i$.
4: Solve the linear system $A^if^i =b^i$ with whatever method is best.
5: Return the $k \times 1$ vector $f^i$ from Step 4 as well as the indices which map points in $X^i$ back to $X$ (since $X^i$ is a subcollection of $X$).
So we have those 5 steps to run through for each point $x_i \in X$. So, if $N = 30,000$, for example, then we have to run through Steps 1-5 30,000 times. Each iteration is independent: there needs to be no communication between running iteration $i$ and iteration $j$. They all rely on the same common dataset, which is $X$.
I've done this in Matlab before and it was pretty easy to code up, but even with vectorizing as much as I could, it got a bit slow when N = 30,000 and k = 400. I decided as a good lesson for myself, to try to write it in C++.
My thoughts on the C++ code were the following:
1: Step 1 requires computing nearest neighbors. Perhaps it would be helpful to use a tree structure to keep in mind what points are nearby. I used Boost's spatial library to set up an R Tree and I query that for step 1 to build the set $X_i$. Please see http://www.boost.org/doc/libs/1_58_0/libs/geometry/doc/html/geometry/spatial_indexes/introduction.html
2: Since i needed linear algebra, but my matrices were "small", I figured something like Trilinos or Petsc would be overkill, so I went with Armadillo. I use armadillo to build the matrices $A^i$, the vectors, and handle the solve step.
3: I'm not very experience with OOP, so I may have made a major faux pas here. I made a class Assembler that has the dataset $X$ and the Rtree as member variables (it builds them), and then it has a method I can call that will instantiate an object L, which L then runs steps 1-5 and stores the solution $f^i$ and the indices as member variables.
4: I have what looks to me to be an embarassingly parallel problem, so I call my Assembler class's method to build up the individual objects L and stick that in a for loop for i=1 to i= N. I then just stuck an OpenMP pragma around the for loop. This gave me pretty good speedups (better than Matlab by a lot, depending on what machine I ran it on).
My questions are the following: Am I on the right track? Anyone have any suggestions or advice on how to improve things? Any suggestions on speed-ups or improved parallelism? If it would be helpful, I can get a link to the code so people can take a look at it.
EDIT: Another question I forgot to ask: each iteration of building object L_i runts through Steps 1-5, but at the end of the day, all I want is to store (for later use) the coefficient vector $f^i$ and the index vector mapping $X_i$ back to $X$ (both of which are member variables object $L_i$ finds during Steps 1-5). Naively, I'm inclined to write these to a csv file or binary file or something, and have one file per object $L_i$, but I'm not exactly sure if that's smart. In MATLAB, I had a sparse matrix $S$ where column $i$ of $S$ was the vector $f^i$ indexed properly with the map that sends $X_i \to X$. Since $k << N$, $S$ was a sparse matrix. This got super slow in MATLAB, as I guess I wasn't allocating $S$ properly and it was extending itself or something.