Profile to Find the Bottleneck
Others have already made this point, but you need to find where your program is spending all its time, and focus on that.
Look for Ways to Speed up the I/O
This is going to be implementation-dependent, but in many cases,
fscanf(), which doesn’t need to copy or heap-allocate, can be faster than
getline(). In fact, you currently allocate and free a buffer from the heap for every data point.
Not in the standard library, your OS, or a third-party library like glib, probably lets you memory-map a regular file. That lets you scan for tokens and pass the beginning and end of each token to
strtod. This can be extremely fast because there is no buffering or copying, just mapping the pages of the file on disk to memory similarly to the virtual memory manager. (G. Sliepen had an excellent series of comments, now removed, pointing out the ways that might not be true.)
Store the Data Points in a Dynamic Array
Instead of reading each value in one at a time, you want to store all the scanned values in memory (maximizing the I/O throughput) and then process them at once (allowing parallelization and SIMD).
It’s possible to implement a resizing vector of
double in C, and that’s probably the data structure you want for this, although you might also make two passes through the file in memory, one to count the items, and then another to fill the array.
Parallelizing the Calculation
The version of the calculation that you’re using is incremental, and therefore inherently serial. If you have all the data in memory, however, there is no need to use an incremental formula. You can make two passes, one reduction to calculate the mean, and then another to calculate the variance as a map-reduction (the sum of (xᵢ-μ)² for all xᵢ), and finally calculate σ from the variance. You will then be adding the squares of small differences, which all have the same sign, and not get any subtractive cancellation from this step (although there could be large data points with opposite signs making the calculation of the mean numerically unstable). This version of the calculation sacrifices some stability for being able to use SIMD.
You can additionally split the both passes between threads, with OpenMP, or calculate the variance of each subrange of the input and then the total variance from those. Since the mean and variance can be calculated from simple sums of terms, you could use the reduction operation of OpenMP for this (or in the C++ standard library,
std::transform_reduce with a parallel execution policy), rather than a weighted-mean and weighted-variance formula. Otherwise, you would want to partition the input array into one slice per worker thread.
While the way I outlined above sacrifices some numerical stability for the speed of SIMD, you might figure out an incremental formula that adds four terms at once, or however many your hardware can efficiently vectorize.