🚀 go-pugleaf

I encountered the following annoying problem while using gcc for
numeric computations and am curious if other people have experienced
similar difficulties.


I wrote a conjugate gradient algorithm to compute the extrema of a
certain function defined on the compact Stiefel manifold V(n,k)
(n-by-k matrices with orthonormal columns).  This algorithm is mostly
QR decompositions and matrix-vector multiplications.

The algorithm converges as predicted, but for moderately sized n and
small k (n = 100, k = 3), I noticed that after about 80 steps or so,
the algorithm slowed down considerably -- from about 0.01 seconds per
step to over 1 second per step.  This was very perplexing because the
steps (should) have a constant flop count.

I assumed that I had made some coding error that resulted in memory
allocation at every step, or had done something that increased the
flop count at every step.  My memory requirements (about 3000 floats)
were much too small for swapping to play a role.  My timing
experiments indicated that the slower running speed was not isolated
in one part of the algorithm.  After ruling out all these
possibilities, I asked a friend for help.

After some puzzling, he suggested that I use SUN's cc compiler instead
of GNU's gcc compiler (v. 2.2.2).  The problem vanished immediately
after recompiling.  Ironically, using gcc's optimization flag -O made
the program run even slower.


I am left wondering if I can trust gcc's reliability for numerical
work.  Has anyone else experienced similar difficulties with gcc, or C
compilers in general, when doing intensive numerical work?  I wish
that I could report a specific bug, but I haven't taken the time (or
have the time) to pare down my algorithm to some small kernel of code
that exhibits this apparently subtle problem.

Steven Smith