srotmg(3S)
SROTMG, DROTMG - Constructs a modified Givens plane rotation
Showing IRIX 6.5.30 (default release). Added in IRIX 6.5.15.
NAME SROTMG, DROTMG - Constructs a modified Givens plane rotation SYNOPSIS Single precision Fortran: CALL SROTMG (d1, d2, b1, b2, rparam) C/C++: #include <scsl_blas.h> void srotmg (float *d1, float *d2, float *b1, float b2, float *rparam); Double precision Fortran: CALL DROTMG (d1, d2, b1, b2, rparam) C/C++: #include <scsl_blas.h> void drotmg (double *d1, double *d2, double *b1, double b2, double *rparam); IMPLEMENTATION These routines are part of the SCSL Scientific Library and can be loaded using either the -lscs or the -lscs_mp option. The -lscs_mp option directs the linker to use the multi-processor version of the library. When linking to SCSL with -lscs or -lscs_mp, the default integer size is 4 bytes (32 bits). Another version of SCSL is available in which integers are 8 bytes (64 bits). This version allows the user access to larger memory sizes and helps when porting legacy Cray codes. It can be loaded by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use only one of the two versions; 4-byte integer and 8-byte integer library calls cannot be mixed. The C and C++ prototypes shown above are appropriate for the 4-byte integer version of SCSL. When using the 8-byte integer version, the variables of type int become long long and the <scsl_blas_i8.h> header file should be included. DESCRIPTION These routines compute the elements of a modified Givens plane rotation matrix. See the NOTES section of this man page for information about the interpretation of the data types described in the following arguments. These routines have the following arguments: d1 First diagonal element. (input and output) SROTMG: Single precision. 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'. For C/C++, a pointer to this value is passed. d2 Second diagonal element. (input and output) SROTMG: Single precision. DROTMG: Double precision. On input, this is the second diagonal element of the scaling matrix D. On the second 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 second diagonal element of the updated scaling matrix D'. For C/C++, a pointer to this value is passed. b1 x-coordinate. (input and output) SROTMG: Single precision. 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'). For C/C++, a pointer to this value is passed. b2 y-coordinate. (input) SROTMG: Single precision. 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: Single precision 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(3S)) 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' * 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/2 H 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. 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) <-- h2,1 = -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) NOTES The following data types are described in this documentation: Term Used Data type Fortran: Array dimensioned n x(n) Integer INTEGER (INTEGER*8 for -lscs_i8[_mp]) Single precision REAL Double precision DOUBLE PRECISION C/C++: Array dimensioned n x[n] Integer int (long long for -lscs_i8[_mp]) Single precision float Double precision double SEE ALSO INTRO_SCSL(3S), INTRO_BLAS1(3S) INTRO_CBLAS(3S) for information about using the C interface to Fortran 77 Basic Linear Algebra Subprograms (legacy BLAS) set forth by the Basic Linear Algebra Subprograms Technical Forum. SROT(3S), SROTG(3S), SROTM(3S) Gentleman, W. M., "Least Squares Computations by Givens Transformations Without Square Roots," Journal of the Institute for Mathematical Applications 12 (1973), pp. 329 - 336. Lawson, C., Hanson, R., Kincaid, D., and Krogh, F., "Basic Linear Algebra Subprograms for Fortran Usage," ACM Transactions on Mathematical Software, 5 (1979), pp. 308 - 325.