rotmg(3F)
SROTMG, DROTMG - Constructs a modified Givens plane rotation
As shipped in IRIX 6.5.19. Last changed in IRIX 6.5.19.
NAME SROTMG, DROTMG - Constructs a modified Givens plane rotation SYNOPSIS Real CALL SROTMG (d1, d2, b1, b2, rparam) Double precision CALL DROTMG (d1, d2, b1, b2, rparam) DESCRIPTION These routines compute the elements of a modified Givens plane rotation matrix. These routines have the following arguments: d1 First diagonal element. (input and output) SROTMG: Real. DROTMG: Double precision. On input, this value is the first diagonal element of the scaling matrix D. On the first call to the routine, this value is typically 1.0. Subsequent calls typically use the value from the previous call. On output, this value is the first diagonal element of the updated scaling matrix D'. d2 First diagonal element. (input and output) SROTMG: Real. DROTMG: Double precision. On input, this is the first diagonal element of the scaling matrix D. On the first call to the routine, this value is typically 1.0. Subsequent calls typically use the value from the previous call. On output, this value is the first diagonal element of the updated scaling matrix D'. b1 x-coordinate. (input and output) SROTMG: Real. DROTMG: Double precision. On input, this value is the x-coordinate of the vector used to define the angle of rotation, before scaling (multiplying by the matrix D). On output, this value is the x-coordinate of the rotated vector, before scaling (multiplying by the matrix D'). b2 y-coordinate. (input) SROTMG: Real. DROTMG: Double precision. On input, this value is the y-coordinate of the vector used to define the angle of rotation, before scaling (multiplying by the matrix D). It is unchanged on output. rparam Real array of dimension 5. (output) SROTMG: Real array. DROTMG: Double precision array. This array contains rotation matrix information. The routine sets up the computed elements in rparam from inputs d1, d2, b1, and b2. Standard Givens Rotation A standard Givens rotation (see SROTG(3F)) is based on an orthogonal matrix G that rotates points on a Cartesian xy-coordinate plane. To calculate the rotation matrix, you must provide the angle of rotation desired, or, equivalently, a vector (point) that lies along the angle of rotation. For a given planar point (xr, yr), G is formed so that: _ _ _ _ _ _ _ _ | x' | | c s | | x | G | x | | | | | | r | | r | | 0 | = |-s c | | y | = | y | | | | | | r | | r | - - - - - - - - where x' = sqrt (xr2 + yr2) . With this rotation matrix G, you can then convert any number of existing planar points to the new (rotated) xy-coordinate system. For n points, the rotations would be as follows: _ _ _ _ _ _ | x | | c s | | x | | i | | | | i | | 0 |<- |-s c | | y | | i | | | | i | - - - - - - for i = 1, 2, ..., n Modified Givens Rotation The algorithm for these routines is based on the following observation. The rotation matrix G can be factored into a scaling matrix (diagonal matrix) and modified rotation matrix H, for which either the diagonal or the off-diagonal elements are units (that is, +1 or -1). Thus, to perform m modified (scaled) rotations on n planar points, requires only 2nm, rather than 4nm multiplications for the standard rotation. Because you may want to perform several successive rotations, this routine assumes that you have leftover scaling factors from your previous modified Givens rotation; that is, the routine requires you to input not only a planar rotation vector (b1, b2) but also the squares of the diagonal elements of the scaling matrix, d1 and d2. The actual rotation vector is specified as follows: _ _ _ _ _ _ 1/2_ _ | xr | | sqrt(d1) 0 | | b1 | D | b1 | | yr | = | 0 sqrt(d2) | | b2 | = | b2 | - - - - - - - - where d1 and d2 are the input scaling factors. Given these inputs, the routine generates a new modified rotation matrix H with units for either the diagonal or off-diagonal elements, and new elements d1' and d2' for the new scaling factors, and rotated but unscaled vector (b1', 0), with the following results: _ _ 1/2_ _ 1/2 _ _ 1/2_ _ _ _ G | xr | G D | b1 | D' H | b1 | = D' | h11 h12 | | b1 | | yr | = | b2 | = | b2 | | h21 h22 | | b2 | - - - - - - - - - - 1/2_ _ _ _ D' | b1' | | x' | = | 0 | = | 0 | - - - - where: 1/2 D' ( = diag {sqrt(d1'), sqrt(d2')} ) uses the updated scaling factors d1' and d2', which are d1 and d2 on output. H is stored in the output array argument rparam; b1' is stored as b1 on output. D'1/2H equals G D 1/2, not G, as implied earlier. You must account for the old scaling factors when calculating the new scaling factors. After calculating the matrix H by using SROTMG/DROTMG, you can then use it in SROTM/DROTM to convert points to the new coordinate system. NOTES Meaning of the Output Values The output values are returned through arguments d1, d2, b1, and rparam. The scaling factors d1 and d2 are updated with each call to the routine. Although SROTM/DROTM does not need the updated factors, they are needed in two other important contexts: * As input for subsequent calls to this routine. * As scaling factors for rotated but unscaled points (xi, yi), which are output from SROTM./DROTM. In this second usage, the actual (scaled) points would be given by (sqrt(d1)*x(i), sqrt(d2)*y(i)). Doing this operation frequently on all your points is counterproductive. The main advantage of the modified rotation algorithm is to reduce the number of operations. If you fold in the scaling factors after each rotation, you are performing the same number of operations as in the standard Givens rotation. These two uses for the scaling factors are mutually exclusive; that is, if you fold the scaling factors back into all your points, you no longer need those factors for these routines. After folding in the scaling factors, and before the next call to SROTMG/DROTMG, reset d1 and d2 to 1.0. On output, b1 represents the new x-coordinate after rotating (but before scaling) the rotation vector. Although the y-coordinate of this vector is 0.0 (see the the previous discussion), the corresponding value b2 is unchanged on output. The output array argument rparam specifies the format of matrix H, and it holds the nonunit values of H. This is the only output of these routines that SROTM/DROTM requires. Each element of rparam has a specific meaning, as follows: rparam(1) a flag parameter that specifies how the matrix is stored = 0.0 Off-diagonal elements of H are units. = 1.0 Diagonal elements of H are units. = -1.0 Rescaling case (see the following subsection). = -2.0 H is the identity matrix; no rotation needed. rparam(2) = h1,1 if needed rparam(3) = h2,1 if needed rparam(4) = h1,2 if needed rparam(5) = h2,2 if needed Calculating the Output Values The following presents the algorithm for calculating the output values d1, d2, b1, and rparam, based on the input values d1, d2, b1, and b2. This algorithm is presented without explanation or proof. For a more complete discussion of the modified Givens algorithm, see the papers listed in the SEE ALSO section. Case 1: b2 = 0 (trivial case) In this case, no rotation is needed. The flag value, rparam(1), is set to -2.0. When passed to routine SROTM/DROTM, this flag indicates that it should not do any rotations. On output from SROTMG/DROTMG, d1, d2, and b1 are unchanged. Case 2: | sqrt(d1)*b1 | > | sqrt(d2)*b2 | (|xr| > |yr|) In this case, the diagonal elements of H (h1,1 and h2,2) are set to 1.0 . Thus, the rparam values set on output are as follows: rparam(1) <-- 0.0 rparam(3) <-- h21 = -b2/b1 rparam(4) <-- h1,2 = (d2*b2)/(d1*b1) The output values of the scaling factors are as follows: d1 <-- d1' = d1/u d2 <-- d2' = d2/u where u = det(H) = 1 + (d2*b2*b2)/(d1*b1*b1) . The output value of b1 is as follows: b1 <-- b1' = b1*u If rescaling is needed, SROTMG/DROTMG will further modify these output values before the end of the routine. See case 4 later in this subsection. Case 3: | sqrt(d1)*b1 | <= | sqrt(d2)*b2 | (|xr| <= |yr|) In this case, the off-diagonal elements of H are units (to be specific, h2,1 = -1 and h1,2 = 1). Thus, the rparam values set on output are as follows: rparam(1) <-- 1.0 rparam(2) <-- h1,1 = (d1*b1)/(d2*b2) rparam(5) <-- h2,2 = b1/b2 The output values of the scaling factors are as follows: d1 <-- d1' = d2/u d2 <-- d2' = d1/u where u = det(H) = 1 + (d1*b1*b1)/(d2*b2*b2) . The output value of b1 is as follows: b1 <-- b1' = b2*u If rescaling is needed, SROTMG/DROTMG will further modify these output values before the end of the routine. See case 4. Case 4: Rescaling If the scaling factors become either very large or very small, the scaling and rotation operations may lose a lot of accuracy; therefore, each scaling factor from case 2 or case 3 is kept within the range: gamma**(-2) <= | d(i)' | <= gamma**2 , for i = 1, 2 where d(1)' = d1', d(2)' = d2', and gamma = 2.0**12 = 4096.0. At the end of case 2 or case 3 assignments, if either of the scaling factors falls outside this range, this routine must rescale that factor (and the corresponding elements of H) to bring its size back within the specified range. 48 - log2(gamma) = 48 - 12 = 36 bits (~ 10 decimal digits) Rescaling is performed as follows: If either d1 or d2 is 0, no rescaling is done; otherwise, let q(i) = int(log_base_gamma(sqrt(|d(i)'|))) = int(log2(|d(i)'|)/24) , for i = 1 , 2. Then the following is true: q(i) < 0 , if |d(i)'| < gamma**2 q(i) = 0 , if gamma**(-2) <= |d(i)'| <= gamma**2 q(i) > 0 , if |d(i)'| > gamma**2 Furthermore, |q(i)| represents the number of times d(i)' must be multiplied (or divided) by gamma**2 to return it to the proper range of values. In this case, the rparam values set on output are as follows: rparam(1) <-- -1.0 rparam(2) <-- h1,1' = h1,1*gamma**q(1) rparam(3) <-- h2,1' = h2,1*gamma**q(2) rparam(4) <-- h1,2' = h1,2*gamma**q(1) rparam(5) <-- h2,2' = h2,2*gamma**q(2) The output values of the scaling factors are as follows: d1 <-- d1'' = d1'*gamma**(-2*q(1)) d2 <-- d2'' = d2'*gamma**(-2*q(2)) The output value of b1 is as follows: b1 <-- b1'' = b1'*gamma**q(1) SEE ALSO ROT(3F), ROTG(3F), ROTM(3F) Gentleman, W. M., "Least Squares Computations by Givens Transformations Without Square Roots," Journal of the Institute for Mathematical Applications 12 (1973), Lawson, C., Hanson, R., Kincaid, D., and Krogh, F., "Basic Linear Algebra Subprograms for Fortran Usage," ACM Transactions on Mathematical Software, 5 (1979),