Tempo2 unnecessarily computes "Uinv" when doing generalised least squares

[Michael Keith]

The “cholesky” code first does a cholesky decomposition and then inverts the lower triangular matrix U (don’t ask) in order to whiten the data, however it is much faster to directly solve U.x=b rather than invert U.

This issue was raised by Bill Coles (and George) via email, so this issue is more just for other authors to see what’s happening.

From Bill:

Finding Uinv is not necessary. If we think about the prewhitening as follows:

Dat = M * Params * Err                          where Dat and Err contain red noise, goes to

Dwhite = Mwhite * Params + Ewhite       where Dwhite = Uinv * Dat and Mwhite = Uinv * M

However  Dat = U * Dwhite can be solved for Dwhite in Nobs^2 time using back substitution and
                M = U * Mwhite can be solved in Nobs^2 * Npar the same way.

So this runs  Nobs/Npar faster than using Uinv. That can be a big help. It makes finding U the limiting
factor. I tested it successfully simulating T2 with Matlab.

‌

[Michael Keith]

In practice because “U” is lower triangular, we actually do forward substitution. This has been implemented “by hand” in TKforwardSub and TKforwardSubVec (for matrix and vector) in TKmatrix.C. Seems to work as of f9fd985.

Note that TKforwardSub actually transposes the matricies (other than U), because of the way that we create/store the design matrix in a transposed state. This is quite normal for tempo2 routines though!

Finally, this also uses F77_dpotrf rather than F77_dpotf2 which I believe will be faster for large problems, but I’ve not tested it.

This means that the “uinv” variable now stores “U” rather than “Uinv” - we should probably rename this. For weighted least-squares this still contains 1/sigma at the moment, but I suggest also should be changed to sigma before a merge to master.

Currently has an compile-time option, but I think should be made a permanent change once tested.