@uhoh I pasted the HTML code into an online HTML viewer and got something that looks a bit neater/nicer:
@uhoh , @Tyberius when I profile MATLAB code with the default profiler that comes with MATLAB, if I were to see something like 35% │ 102 │ScalarFunction.mini_me , I can then click on ScalarFunction.mini_me and it will give me a new HTML file that shows the amount of time the computer spent on each line of ScalarFunction.mini_me, and I can continue going deeper and deeper inside until I find the precisely line that is the bigger culprit.
Maybe there's a better Python profiler out there. @uhoh have you used MATLAB?
@NikeDattani Hi, technically yes I have used Matlab, but these days the only thing I can use is Python, Math, and English. I do feel that at least for my simple calculation a more microscopic analysis is necessary, the heavy lifting is really concentrated in just a few lines of NumPy and is really being done in complied code I'll never see (mostly by choice!)
@uhoh The MATLAB profiler would tell you which lines of compiled code are taking up the most time. Those compiled codes could then be replaced by compiled CUDA codes to run on the GPU. ScalarFunction.mini_me is not compiled code right?
@NikeDattani Yes I believe that statement is correct.
But Python is a scripting language, and when we use numpy arrays we are aware that everything we do is running in compiled code; numpy is just the easy-to-use packaging for it
@uhoh MATLAB is also a scripting language. I don't think every use of numpy arrays will be running in compiled code, and even if so, we would want to know which of the compiled codes called within ScalarFunction.mini_me are taking the longest to run.
@NikeDattani yes, it's mostly line 185 vectors = array[self.atoms_2] - array[self.atoms_1] which is where the indexing happens which is why "indexing" is in the title of my question
With some wasted space the lattice can be mapped to a square with nearest-neighbor only interactions; the problem should be mappable without indexing to a 2D array of processors like a GPU or even an FPGA. But I'm not a programmer so I rely on scripting languages to figure out how to do that
A and B are lists of (i, j) indices in a 2D array of atoms. Array(i, j) is the 3D positions of atoms (i, j). Since it's a hexagonal rather than a square lattice.
Subtracting returns all the r(i, j) vectors representing bonds.
Since it's a hexagonal rather than a square lattice... an atom has SIX nearest neighbors, and the edges of the domain have to be treated differently
So using indexing (a big list of all possible bonds) is slow because it moves all over the place in the array. Basically pointing "here" then "here" then "here"
If the array completely fits into a CPU cache, then it's pretty fast, but as soon as there is cache swapping it slows way down.
My hope with CuPy is that it handles indexing smartly and knows how to distribute the calculation among the GPU cores
So A is a list of all the (i, j) coordinates of the "left hand atoms" of each of the bonds (the color line segments) and B is a list of all the (i, j) coordinates of the atoms on the other sides of the bonds. Subtracting array(A) - array(B) returns the vectors representing the bond lengths/directions
If there are 1000 atoms then there are roughly 3,000 bonds, but less since we have to omit the edges. And in the future I'll add defects (adatoms, vacancies) so the indexing is a very nice way to stick defects into the problem.
These are defects in a honeycomb lattice (I'll do that later) ncbi.nlm.nih.gov/pmc/articles/PMC6190370 and they are easy to generate in an indexing scheme and probably hard to do directly mapping atoms to points in a 2D array.
@ParmeetSinghEP066 Do you need access to more HPC resources? Let me know: nike@hpqc.org
@uhoh so for N atoms, A and B each contain 2N elements which are the (i,j) coordinates in Angstroms, and array[self.A] - array[self.B] contains the Euclidean distances: sqrt( (A(i) - B(i))^2 + (A(j) - B(j))^2)?
I guess you're saying that for N =1000, there's roughly 3000 bonds because 3000 is roughly what you'd get when calculating 3N-6? How does Python know that though? If it calculates all distances between atoms in A and atoms in B, then it would be more like N^2 distances.
@NikeDattani wait, I said something wrong, it's not (i, j) it's just (i).
Let me start again
array is N x 3 where N is the number of atoms and 3 is the spatial dimensions
A and B are long lists of indices (i) A might be (1,2 3, 4...) and B might be (4, 3, 2, 8...)
array(A) - array(B) will be a new array of shape N x 3
it represents the vectors between atoms in the A list and the atoms in the B list.
So r_14, r_23, r_32, r_48, etc.
then the variable lengths is the lengths of those bonds, used to calculate energy based on a simple harmonic oscillator potential (deviation from nominal distance)
The number of bonds is a little less than 3N. Each atom has 6 bonds (except at the edges) and each bond is shared by two atoms.
Python doesn't know that, I have to work out the details in the function define_nearest_neighbor_bonds() starting on line 55. I make those detailed lists of indices A and B by thinking about it and scripting an algorithm.
Once those lists of indices are generated it defines the problem.
Yea that's just all 0s. It might be faster to do A - A than to do A - B.
you're getting the same time with A -B ?
If Python was smarter, it would know that A - A = 0, and not worry about doing the subtraction (which takes time)!
Another thing, are all the A(i) and B(i) values supposed to be integers? It seems to be that case in your code, but I'd imagine that distances would be floating-point numbers (the arithmetic speed would be different).
I tried it and it's not any different and that would be shocking if it was. Numpy does exactly what you tell it to do, there is no hidden cleverness.
The positions of the N atoms are stored in the (N x 3) floating point array.
A and B are indices of that array
this is what "indexing" is, a big blob of indices of another array
So 3000 x 3 = 9000 subtractions in 130 us is 14 nanoseconds per subtraction. Doing strait subtraction of two arrays (i.e. array1 - array2) on my laptop would take only about 1 nanosecond per subtraction.
The bad news I think is that these arrays are too small to gain any benefit from a GPU, but the good news is that I think the code can be sped up a lot.
The problem with small arrays, is that the CPU might take only 0.3 seconds to do the calculations, and even if the GPU takes 0.02 seconds, you lose more time with the GPU because there's the time that it takes to move the data from the CPU to the GPU and then to transfer the final result from the GPU to CPU. This is a "bandwidth limited calculation", as opposed to a "compute limited calculation".
I still worry that 100x larger, means only 100,000 x 3 elements in the array.
If each of the 300,000 numbers is an 8-byte (double-precision floating point) number, then it's only 2.4 MB.
That's a table from this paper of mine it shows that major gains from the GPU don't kick in until the data is a few megabytes in size. If the data is only 1MB in size, the CPU takes 0.64 seconds vs 2.2 seconds on the GPU due to the time that's lost with transferring the data between the CPU and GPU.
Oh that's really helpful to hear and think about. If several iterations of the minimization could be done in some parallel way at once in the GPU, instead of moving it back to the CPU each time, that might take more advantage.
When atoms actually relax on a surface, they do it locally, they don't have a master program calculating the energy and the next step.
I wonder if the problem could just run in GPU with only infrequent transfers somehow?
Can each node just be an atom and move by itself somehow?
btw I just moved a sofa up three flights of stairs with a friend so I will take a look at your paper a little later when I've got my brain oxygenated again.
I think we need to work with bigger arrays in order to get any meaningful timing/benchmarking results. At the microsecond scale, there's a lot of noise.
Even on the laptop, if the array sizes and number of repetitions of the whole process, could be increased, then it would be helpful.
@uhoh You could maybe include a link to the HTML since it should give more navigation options. It looks like scalene doesn't go into libraries with the defaults settings, but will if you specify --profile-all. I don't think it will go all the way down to C code, but it should at least go into the internal numpy/scipy functions.
@Tyberius Thanks I will dig in to scalene's options today and try that. I just got news that a coworker (may have) succeeded in firing up a GPU with CuPy which means I'll hopefully be able to get my hands dirty soon, and that I can apply scalene to a CPU/GPU comparison, yay!
I've added two links to the original results in note 3 of the question but will now go back and chekc the options as soon as this coffee kicks in.
@Tyberius @uhoh I moved the messages to the HPQC room since we're more talking about "high-performance computing" than about "what's going on elsewhere"
