Three ways to make your code even faster:
- you can eliminate almost all divisions by working with arrays that are numerically equal to 1/v[n] and v[n-1]/v[n].
- you don’t need to renormalize so often. Using an if statement, you can accumulate the square of the 2-norm, and if gets too close to overflow, only then must you renormalize.
- you can vectorize the code across (neighbouring) eigenvectors. This requires a bit more work to tease Julia’s compiler to do this for you.
To give you some ideas, the following C code (probably you’re not looking for a C solution, but nevertheless…) does all of these things when performing a transposed eigenmatrix-vector product (of symmetric tridiagonal eigenvectors), V^\top c. The default version is arguably more legible https://github.com/MikaelSlevinsky/FastTransforms/blob/5ee9c149067fa1d1484ea4cf291f9e8d8f1e0bb6/src/recurrence/recurrence_default.c#L113-#L144
but the macro gets used to generate vectorized versions
https://github.com/MikaelSlevinsky/FastTransforms/blob/master/src/recurrence/recurrence.h#L103-#L176