Givens Rotation in old SLATEC Routine

Open discussion regarding features, bugs, issues, vendors, etc.

Givens Rotation in old SLATEC Routine

Postby DavidB » Thu Dec 13, 2018 11:19 pm

I am translating an old Fortran sub-routine (R1UPDT) from the SLATEC library into C++ and came across a section of code that puzzles me.

A Givens Rotation is pretty straightforward; it comes down to computing r = 1.0/sqrt(a^2 + b^2)
Usually, this equation is put in a slightly different form, depending upon whether abs(a) < abs( b) or abs(a) > abs(b). A temporary variable, c, is used such that c is smaller than 1.

So, if abs(a) < abs(b), c = a/b and the equation is written r = 1.0 / sqrt(c^2 + 1)

On the other hand, if abs(a) > abs(b), c = b/a and the equation is written r = 1.0 / sqrt(1 + c^2)

Either way, the potential risk of overflow is avoided because the term squared is less than 1.
So far, so good.

However, the Fortan code in R1UPDT is written such that r = 0.5 / sqrt(0.25 + 0.25 * c^2)

Code: Select all
SIN = 0.5 / SQRT(0.25 + 0.25 * COTAN**2)

I can't help but ask, why leave the code like this? The equation could be factored by hand simply enough: take the 0.5 factor out of the radicand in the denominator, cancel the common factor from numerator and denominator, and avoid a multiplication process:

Code: Select all
SIN = 1.0 / SQRT(1.0 + COTAN**2)

Since this code was originally written in the 80's, I am wondering if keeping the 0.25 factor in the radicand in the denominator is one more technique for avoiding overflow (or something else.)

Any ideas or suggestions? This routine is actually a helper routine for SNSQE, which is a multi-dimensional root-finding routine. Anybody here use SNSQE?
Posts: 38
Joined: Thu Dec 15, 2005 1:49 pm
Location: Vancouver, B.C. Canada

Return to User Discussion

Who is online

Users browsing this forum: No registered users and 6 guests