\begin{page}{manpageXXc02}{NAG On-line Documentation: c02} \beginscroll \begin{verbatim} C02(3NAG) Foundation Library (12/10/92) C02(3NAG) C02 -- Zeros of Polynomials Introduction -- C02 Chapter C02 Zeros of Polynomials 1. Scope of the Chapter This chapter is concerned with computing the zeros of a polynomial with real or complex coefficients. 2. Background to the Problems Let f(z) be a polynomial of degree n with complex coefficients a : i n n-1 n-2 f(z)==a z +a z +a z +...+a z+a , a /=0. 0 1 2 n-1 n 0 A complex number z is called a zero of f(z) (or equivalently a 1 root of the equation f(z)=0), if: f(z )=0. 1 If z is a zero, then f(z) can be divided by a factor (z-z ): 1 1 f(z)=(z-z )f (z) (1) 1 1 where f (z) is a polynomial of degree n-1. By the Fundamental 1 Theorem of Algebra, a polynomial f(z) always has a zero, and so the process of dividing out factors (z-z ) can be continued until i we have a complete factorization of f(z) f(z)==a (z-z )(z-z )...(z-z ). 0 1 2 n Here the complex numbers z ,z ,...,z are the zeros of f(z); they 1 2 n may not all be distinct, so it is sometimes more convenient to write: m m m 1 2 k f(z)==a (z-z ) (z-z ) ...(z-z ) , k<=n, 0 1 2 k with distinct zeros z ,z ,...,z and multiplicities m >=1. If 1 2 k i m =1, z is called a single zero, if m >1, z is called a i i i i multiple or repeated zero; a multiple zero is also a zero of the derivative of f(z). If the coefficients of f(z) are all real, then the zeros of f(z) are either real or else occur as pairs of conjugate complex numbers x+iy and x-iy. A pair of complex conjugate zeros are the 2 zeros of a quadratic factor of f(z), (z +rz+s), with real coefficients r and s. Mathematicians are accustomed to thinking of polynomials as pleasantly simple functions to work with. However the problem of numerically computing the zeros of an arbitrary polynomial is far from simple. A great variety of algorithms have been proposed, of which a number have been widely used in practice; for a fairly comprehensive survey, see Householder [1]. All general algorithms are iterative. Most converge to one zero at a time; the corresponding factor can then be divided out as in equation (1) above - this process is called deflation or, loosely, dividing out the zero - and the algorithm can be applied again to the polynomial f (z). A pair of complex conjugate zeros can be 1 divided out together - this corresponds to dividing f(z) by a quadratic factor. Whatever the theoretical basis of the algorithm, a number of practical problems arise: for a thorough discussion of some of them see Peters and Wilkinson [2] and Wilkinson [3]. The most elementary point is that, even if z is mathematically an exact 1 zero of f(z), because of the fundamental limitations of computer arithmetic the computed value of f(z ) will not necessarily be 1 exactly 0.0. In practice there is usually a small region of values of z about the exact zero at which the computed value of f(z) becomes swamped by rounding errors. Moreover in many algorithms this inaccuracy in the computed value of f(z) results in a similar inaccuracy in the computed step from one iterate to the next. This limits the precision with which any zero can be computed. Deflation is another potential cause of trouble, since, in the notation of equation (1), the computed coefficients of f (z) will not be completely accurate, especially if z is not an 1 1 exact zero of f(z); so the zeros of the computed f (z) will 1 deviate from the zeros of f(z). A zero is called ill-conditioned if it is sensitive to small changes in the coefficients of the polynomial. An ill-conditioned zero is likewise sensitive to the computational inaccuracies just mentioned. Conversely a zero is called well-conditioned if it is comparatively insensitive to such perturbations. Roughly speaking a zero which is well separated from other zeros is well- conditioned, while zeros which are close together are ill- conditioned, but in talking about 'closeness' the decisive factor is not the absolute distance between neighbouring zeros but their ratio: if the ratio is close to 1 the zeros are ill-conditioned. In particular, multiple zeros are ill-conditioned. A multiple zero is usually split into a cluster of zeros by perturbations in the polynomial or computational inaccuracies. 2.1. References [1] Householder A S (1970) The Numerical Treatment of a Single Nonlinear Equation. McGraw-Hill. [2] Peters G and Wilkinson J H (1971) Practical Problems Arising in the Solution of Polynomial Equations. J. Inst. Maths Applics. 8 16--35. [3] Wilkinson J H (1963) Rounding Errors in Algebraic Processes, Chapter 2. HMSO. 3. Recommendations on Choice and Use of Routines 3.1. Discussion Two routines are available: C02AFF for polynomials with complex coefficients and C02AGF for polynomials with real coefficients. C02AFF and C02AGF both use a variant of Laguerre's Method due to BT Smith to calculate each zero until the degree of the deflated polynomial is less than 3, whereupon the remaining zeros are obtained using the 'standard' closed formulae for a quadratic or linear equation. The accuracy of the roots will depend on how ill-conditioned they are. Peters and Wilkinson [2] describe techniques for estimating the errors in the zeros after they have been computed. 3.2. Index Zeros of a complex polynomial C02AFF Zeros of a real polynomial C02AGF C02 -- Zeros of Polynomials Contents -- C02 Chapter C02 Zeros of Polynomials C02AFF All zeros of complex polynomial, modified Laguerre method C02AGF All zeros of real polynomial, modified Laguerre method \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc02aff}{NAG On-line Documentation: c02aff} \beginscroll \begin{verbatim} C02AFF(3NAG) Foundation Library (12/10/92) C02AFF(3NAG) C02 -- Zeros of Polynomials C02AFF C02AFF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C02AFF finds all the roots of a complex polynomial equation, using a variant of Laguerre's Method. 2. Specification SUBROUTINE C02AFF (A, N, SCALE, Z, W, IFAIL) INTEGER N, IFAIL DOUBLE PRECISION A(2,N+1), Z(2,N), W(4*(N+1)) LOGICAL SCALE 3. Description The routine attempts to find all the roots of the nth degree complex polynomial equation n n-1 n-2 P(z)=a z +a z +a z +...+a z+a =0. 0 1 2 n-1 n The roots are located using a modified form of Laguerre's Method, originally proposed by Smith [2]. The method of Laguerre [3] can be described by the iterative scheme -n*P(z ) k L(z )=z -z = ----------------, k k+1 k P'(z )+- /H(z ) k \/ k 2 where H(z )=(n-1)*[(n-1)*(P'(z )) -n*P(z )P''(z )], and z is k k k k 0 specified. The sign in the denominator is chosen so that the modulus of the Laguerre step at z , viz. |L(z )|, is as small as possible. The k k method can be shown to be cubically convergent for isolated roots (real or complex) and linearly convergent for multiple roots. The routine generates a sequence of iterates z , z , z ,..., such 1 2 3 that |P(z )|<|P(z )| and ensures that z +L(z ) 'roughly' k+1 k k+1 k+1 lies inside a circular region of radius |F| about z known to k contain a zero of P(z); that is, |L(z )|<=|F|, where F denotes k+1 the Fejer bound (see Marden [1]) at the point z . Following Smith k [2], F is taken to be min(B,1.445*n*R), where B is an upper bound for the magnitude of the smallest zero given by 1/n B=1.0001*min(\/n*L(z ),|r |,|a /a | ), k 1 n 0 r is the zero X of smaller magnitude of the quadratic equation 1 2 2(P''(z )/(2*n*(n-1)))X +2(P'(z )/n)X+P(z )=0 k k k and the Cauchy lower bound R for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation n n-1 n-2 |a |z +|a |z +|a |z +...+|a |z-|a |=0. 0 1 2 n-1 n Starting from the origin, successive iterates are generated according to the rule z =z +L(z ) for k = 1,2,3,... and L(z ) k+1 k k k is 'adjusted' so that |P(z )|<|P(z )| and |L(z )|<=|F|. The k+1 k k+1 iterative procedure terminates if P(z ) is smaller in absolute k+1 value than the bound on the rounding error in P(z ) and the k+1 current iterate z =z is taken to be a zero of P(z). The p k+1 ~ deflated polynomial P(z)=P(z)/(z-z ) of degree n-1 is then p formed, and the above procedure is repeated on the deflated polynomial until n<3, whereupon the remaining roots are obtained via the 'standard' closed formulae for a linear (n = 1) or quadratic (n = 2) equation. To obtain the roots of a quadratic polynomial, C02AHF(*) can be used. 4. References [1] Marden M (1966) Geometry of Polynomials. Mathematical Surveys. 3 Am. Math. Soc., Providence, RI. [2] Smith B T (1967) ZERPOL: A Zero Finding Algorithm for Polynomials Using Laguerre's Method. Technical Report. Department of Computer Science, University of Toronto, Canada. [3] Wilkinson J H (1965) The Algebraic Eigenvalue Problem. Clarendon Press. 5. Parameters 1: A(2,N+1) -- DOUBLE PRECISION array Input On entry: if A is declared with bounds (2,0:N), then A(1,i) and A(2,i) must contain the real and imaginary parts of a i n-i (i.e., the coefficient of z ), for i=0,1,...,n. Constraint: A(1,0) /= 0.0 or A(2,0) /= 0.0. 2: N -- INTEGER Input On entry: the degree of the polynomial, n. Constraint: N >= 1. 3: SCALE -- LOGICAL Input On entry: indicates whether or not the polynomial is to be scaled. See Section 8 for advice on when it may be preferable to set SCALE = .FALSE. and for a description of the scaling strategy. Suggested value: SCALE = .TRUE.. 4: Z(2,N) -- DOUBLE PRECISION array Output On exit: the real and imaginary parts of the roots are stored in Z(1,i) and Z(2,i) respectively, for i=1,2,...,n. 5: W(4*(N+1)) -- DOUBLE PRECISION array Workspace 6: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry A(1,0) = 0.0 and A(2,0) = 0.0, or N < 1. IFAIL= 2 The iterative procedure has failed to converge. This error is very unlikely to occur. If it does, please contact NAG immediately, as some basic assumption for the arithmetic has been violated. See also Section 8. IFAIL= 3 Either overflow or underflow prevents the evaluation of P(z) near some of its zeros. This error is very unlikely to occur. If it does, please contact NAG immediately. See also Section 8. 7. Accuracy All roots are evaluated as accurately as possible, but because of the inherent nature of the problem complete accuracy cannot be guaranteed. 8. Further Comments If SCALE = .TRUE., then a scaling factor for the coefficients is chosen as a power of the base B of the machine so that the EMAX-P largest coefficient in magnitude approaches THRESH = B . Users should note that no scaling is performed if the largest coefficient in magnitude exceeds THRESH, even if SCALE = .TRUE.. (For definition of B, EMAX and P see the Chapter Introduction X02.) However, with SCALE = .TRUE., overflow may be encountered when the input coefficients a ,a ,a ,...,a vary widely in magnitude, 0 1 2 n (4*P) particularly on those machines for which B overflows. In such cases, SCALE should be set to .FALSE. and the coefficients scaled so that the largest coefficient in magnitude does not (EMAX-2*P) exceed B . Even so, the scaling strategy used in C02AFF is sometimes insufficient to avoid overflow and/or underflow conditions. In such cases, the user is recommended to scale the independent variable (z) so that the disparity between the largest and smallest coefficient in magnitude is reduced. That is, use the routine to locate the zeros of the polynomial d*P(cz) for some suitable values of c and d. For example, if the original -100 100 20 -10 polynomial was P(z)=2 i+2 z , then choosing c=2 and 100 20 d=2 , for instance, would yield the scaled polynomial i+z , which is well-behaved relative to overflow and underflow and has 10 zeros which are 2 times those of P(z). If the routine fails with IFAIL = 2 or 3, then the real and imaginary parts of any roots obtained before the failure occurred are stored in Z in the reverse order in which they were found. Let n denote the number of roots found before the failure R occurred. Then Z(1,n) and Z(2,n) contain the real and imaginary parts of the 1st root found, Z(1,n-1) and Z(2,n-1) contain the real and imaginary parts of the 2nd root found, ..., Z(1,n ) and R Z(2,n ) contain the real and imaginary parts of the n th root R R found. After the failure has occurred, the remaining 2*(n-n ) R elements of Z contain a large negative number (equal to -1/(X02AMF().\/2)). 9. Example 5 4 3 2 To find the roots of the polynomial a z +a z +a z +a z +a z+a =0, 0 1 2 3 4 5 where a =(5.0+6.0i), a =(30.0+20.0i), a =-(0.2+6.0i), 0 1 2 a =(50.0+100000.0i), a =-(2.0-40.0i) and a =(10.0+1.0i). 3 4 5 The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc02agf}{NAG On-line Documentation: c02agf} \beginscroll \begin{verbatim} C02AGF(3NAG) Foundation Library (12/10/92) C02AGF(3NAG) C02 -- Zeros of Polynomials C02AGF C02AGF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C02AGF finds all the roots of a real polynomial equation, using a variant of Laguerre's Method. 2. Specification SUBROUTINE C02AGF (A, N, SCALE, Z, W, IFAIL) INTEGER N, IFAIL DOUBLE PRECISION A(N+1), Z(2,N), W(2*(N+1)) LOGICAL SCALE 3. Description The routine attempts to find all the roots of the nth degree real polynomial equation n n-1 n-2 P(z)=a z +a z +a z +...+a z+a =0. 0 1 2 n-1 n The roots are located using a modified form of Laguerre's Method, originally proposed by Smith [2]. The method of Laguerre [3] can be described by the iterative scheme -n*P(z ) k L(z )=z -z = ----------------, k k+1 k P'(z )+- /H(z ) k \/ k 2 where H(z )=(n-1)*[(n-1)*(P'(z )) -n*P(z )P''(z )], and z is k k k k 0 specified. The sign in the denominator is chosen so that the modulus of the Laguerre step at z , viz. |L(z )|, is as small as possible. The k k method can be shown to be cubically convergent for isolated roots (real or complex) and linearly convergent for multiple roots. The routine generates a sequence of iterates z , z , z ,..., such 1 2 3 that |P(z +1)|<|P(z )| and ensures that z +L(z ) 'roughly' k k k+1 k+1 lies inside a circular region of radius |F| about z known to k contain a zero of P(z); that is, |L(z )|<=|F|, where F denotes k+1 the Fejer bound (see Marden [1]) at the point z . Following Smith k [2], F is taken to be min(B,1.445*n*R), where B is an upper bound for the magnitude of the smallest zero given by 1/n B=1.0001*min(\/n*L(z ),|r |,|a /a | ), k 1 n 0 r is the zero X of smaller magnitude of the quadratic equation 1 2 2(P''(z )/(2*n*(n-1)))X +2(P'(z )/n)X+P(z )=0 k k k and the Cauchy lower bound R for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation n n-1 n-2 |a |z +|a |z +|a |z +...+|a |z-|a |=0. 0 1 2 n-1 n Starting from the origin, successive iterates are generated according to the rule z =z +L(z ) for k=1,2,3,... and L(z ) is k+1 k k k k+1 k k+1 iterative procedure terminates if P(z ) is smaller in absolute k+1 value than the bound on the rounding error in P(z ) and the k+1 current iterate z =z is taken to be a zero of P(z) (as is its p k-1 conjugate z if z is complex). The deflated polynomial p p ~ P(z)=P(z)/(z-z ) of degree n-1 if z is real p p ~ (P(z)=P(z)/((z-z )(z-z )) of degree n-2 if z is complex) is then p p p formed, and the above procedure is repeated on the deflated polynomial until n<3, whereupon the remaining roots are obtained via the 'standard' closed formulae for a linear (n = 1) or quadratic (n = 2) equation. To obtain the roots of a quadratic polynomial, C02AJF(*) can be used. 4. References [1] Marden M (1966) Geometry of Polynomials. Mathematical Surveys. 3 Am. Math. Soc., Providence, RI. [2] Smith B T (1967) ZERPOL: A Zero Finding Algorithm for Polynomials Using Laguerre's Method. Technical Report. Department of Computer Science, University of Toronto, Canada. [3] Wilkinson J H (1965) The Algebraic Eigenvalue Problem. Clarendon Press. 5. Parameters 1: A(N+1) -- DOUBLE PRECISION array Input On entry: if A is declared with bounds (0:N), then A(i) n-i must contain a (i.e., the coefficient of z ), for i i=0,1,...,n. Constraint: A(0) /= 0.0. 2: N -- INTEGER Input On entry: the degree of the polynomial, n. Constraint: N >= 1. 3: SCALE -- LOGICAL Input On entry: indicates whether or not the polynomial is to be scaled. See Section 8 for advice on when it may be preferable to set SCALE = .FALSE. and for a description of the scaling strategy. Suggested value: SCALE = .TRUE.. 4: Z(2,N) -- DOUBLE PRECISION array Output On exit: the real and imaginary parts of the roots are stored in Z(1,i) and Z(2,i) respectively, for i=1,2,...,n. Complex conjugate pairs of roots are stored in consecutive pairs of elements of Z; that is, Z(1,i+1) = Z(1,i) and Z(2,i+1)=-Z(2,i). 5: W(2*(N+1)) -- DOUBLE PRECISION array Workspace 6: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry A(0) = 0.0, or N < 1. IFAIL= 2 The iterative procedure has failed to converge. This error is very unlikely to occur. If it does, please contact NAG immediately, as some basic assumption for the arithmetic has been violated. See also Section 8. IFAIL= 3 Either overflow or underflow prevents the evaluation of P(z) near some of its zeros. This error is very unlikely to occur. If it does, please contact NAG immediately. See also Section 8. 7. Accuracy All roots are evaluated as accurately as possible, but because of the inherent nature of the problem complete accuracy cannot be guaranteed. 8. Further Comments If SCALE = .TRUE., then a scaling factor for the coefficients is chosen as a power of the base B of the machine so that the EMAX-P largest coefficient in magnitude approaches THRESH = B . Users should note that no scaling is performed if the largest coefficient in magnitude exceeds THRESH, even if SCALE = .TRUE.. (For definition of B, EMAX and P see the Chapter Introduction X02.) However, with SCALE = .TRUE., overflow may be encountered when the input coefficients a ,a ,a ,...,a vary widely in magnitude, 0 1 2 n (4*P) particularly on those machines for which B overflows. In such cases, SCALE should be set to .FALSE. and the coefficients scaled so that the largest coefficient in magnitude does not (EMAX-2*P) exceed B . Even so, the scaling strategy used in C02AGF is sometimes insufficient to avoid overflow and/or underflow conditions. In such cases, the user is recommended to scale the independent variable (z) so that the disparity between the largest and smallest coefficient in magnitude is reduced. That is, use the routine to locate the zeros of the polynomial d*P(cz) for some suitable values of c and d. For example, if the original -100 100 20 -10 polynomial was P(z)=2 +2 z , then choosing c=2 and 100 20 d=2 , for instance, would yield the scaled polynomial 1+z , which is well-behaved relative to overflow and underflow and has 10 zeros which are 2 times those of P(z). If the routine fails with IFAIL = 2 or 3, then the real and imaginary parts of any roots obtained before the failure occurred are stored in Z in the reverse order in which they were found. Let n denote the number of roots found before the failure R occurred. Then Z(1,n) and Z(2,n) contain the real and imaginary parts of the 1st root found, Z(1,n-1) and Z(2,n-1) contain the real and imaginary parts of the 2nd root found, ..., Z(1,n ) and R Z(2,n ) contain the real and imaginary parts of the n th root R R found. After the failure has occurred, the remaining 2*(n-n ) R elements of Z contain a large negative number (equal to -1/(X02AMF().\/2)). 9. Example To find the roots of the 5th degree polynomial 5 4 3 2 z +2z +3z +4z +5z+6=0. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc05}{NAG On-line Documentation: c05} \beginscroll \begin{verbatim} C05(3NAG) Foundation Library (12/10/92) C05(3NAG) C05 -- Roots of One or More Transcendental Equations Introduction -- C05 Chapter C05 Roots of One or More Transcendental Equations 1. Scope of the Chapter This chapter is concerned with the calculation of real zeros of continuous real functions of one or more variables. (Complex equations must be expressed in terms of the equivalent larger system of real equations.) 2. Background to the Problems The chapter divides naturally into two parts. 2.1. A Single Equation The first deals with the real zeros of a real function of a single variable f(x). At present, there is only one routine with a simple calling sequence. This routine assumes that the user can determine an initial interval [a,b] within which the desired zero lies, that is f(a)*f(b)<0, and outside which all other zeros lie. The routine then systematically subdivides the interval to produce a final interval containing the zero. This final interval has a length bounded by the user's specified error requirements; the end of the interval where the function has smallest magnitude is returned as the zero. This routine is guaranteed to converge to a simple zero of the function. (Here we define a simple zero as a zero corresponding to a sign-change of the function.) The algorithm used is due to Bus and Dekker. 2.2. Systems of Equations The routines in the second part of this chapter are designed to solve a set of nonlinear equations in n unknowns T f (x)=0, i=1,2,...,n, x=(x ,x ,...,x ) (1) i 1 2 n where T stands for transpose. It is assumed that the functions are continuous and differentiable so that the matrix of first partial derivatives of the functions, the Jacobian matrix J (x)=ddf /ddx evaluated at ij i j the point x, exists, though it may not be possible to calculate it directly. The functions f must be independent, otherwise there will be an i infinity of solutions and the methods will fail. However, even when the functions are independent the solutions may not be unique. Since the methods are iterative, an initial guess at the solution has to be supplied, and the solution located will usually be the one closest to this initial guess. 2.3. References [1] Gill P E and Murray W (1976) Algorithms for the Solution of the Nonlinear Least-squares Problem. NAC 71 National Physical Laboratory. [2] More J J, Garbow B S and Hillstrom K E (1974) User Guide for Minpack-1. ANL-80-74 Argonne National Laboratory. [3] Ortega J M and Rheinboldt W C (1970) Iterative Solution of Nonlinear Equations in Several Variables. Academic Press. [4] Rabinowitz P (1970) Numerical Methods for Nonlinear Algebraic Equations. Gordon and Breach. 3. Recommendations on Choice and Use of Routines 3.1. Zeros of Functions of One Variable There is only one routine (C05ADF) for solving a single nonlinear equation. This routine is designed for solving problems where the function f(x) whose zero is to be calculated, can be coded as a user-supplied routine. C05ADF may only be used when the user can supply an interval [a,b] containing the zero, that is f(a)*f(b)<0. 3.2. Solution of Sets of Nonlinear Equations The solution of a set of nonlinear equations f (x ,x ,...,x )=0, i=1,2,...,n (2) i 1 2 n can be regarded as a special case of the problem of finding a minimum of a sum of squares m / 2 s(x)= | [f (x ,x ,...,x )] (m>=n). (3) / i 1 2 n i=1 So the routines in Chapter E04 of the Library are relevant as well as the special nonlinear equations routines. There are two routines (C05NBF and C05PBF) for solving a set of nonlinear equations. These routines require the f (and possibly i their derivatives) to be calculated in user-supplied routines. These should be set up carefully so the Library routines can work as efficiently as possible. The main decision which has to be made by the user is whether to ddf i supply the derivatives ----. It is advisable to do so if ddx j possible, since the results obtained by algorithms which use derivatives are generally more reliable than those obtained by algorithms which do not use derivatives. C05PBF requires the user to provide the derivatives, whilst C05NBF does not. C05NBF and C05PBF are easy-to-use routines. A routine, C05ZAF, is provided for use in conjunction with C05PBF to check the user-provided derivatives for consistency with the functions themselves. The user is strongly advised to make use of this routine whenever C05PBF is used. Firstly, the calculation of the functions and their derivatives should be ordered so that cancellation errors are avoided. This is particularly important in a routine that uses these quantities to build up estimates of higher derivatives. Secondly, scaling of the variables has a considerable effect on the efficiency of a routine. The problem should be designed so that the elements of x are of similar magnitude. The same comment applies to the functions, all the f should be of comparable i size. The accuracy is usually determined by the accuracy parameters of the routines, but the following points may be useful: (i) Greater accuracy in the solution may be requested by choosing smaller input values for the accuracy parameters. However, if unreasonable accuracy is demanded, rounding errors may become important and cause a failure. (ii) Some idea of the accuracies of the x may be obtained by i monitoring the progress of the routine to see how many figures remain unchanged during the last few iterations. (iii) An approximation to the error in the solution x, given by e where e is the solution to the set of linear equations J(x)e=-f(x) T where f(x)=(f (x),f (x),...,f (x)) (see Chapter F04). 1 2 n (iv) If the functions f (x) are changed by small amounts i (epsilon) , for i=1,2,...,n, then the corresponding change i in the solution x is given approximately by (sigma), where (sigma) is the solution of the set of linear equations J(x)(sigma)=-(epsilon), (see Chapter F04). Thus one can estimate the sensitivity of x to any uncertainties in the specification of f (x), for i i=1,2,...,n. 3.3. Index Zeros of functions of one variable: Bus and Dekker algorithm C05ADF Zeros of functions of several variables: easy-to-use C05NBF easy-to-use, derivatives required C05PBF Checking Routine: Checks user-supplied Jacobian C05ZAF C05 -- Roots of One or More Transcendental Equations Contents -- C05 Chapter C05 Roots of One or More Transcendental Equations C05ADF Zero of continuous function in given interval, Bus and Dekker algorithm C05NBF Solution of system of nonlinear equations using function values only C05PBF Solution of system of nonlinear equations using 1st derivatives C05ZAF Check user's routine for calculating 1st derivatives \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc05adf}{NAG On-line Documentation: c05adf} \beginscroll \begin{verbatim} C05ADF(3NAG) Foundation Library (12/10/92) C05ADF(3NAG) C05 -- Roots of One or More Transcendental Equations C05ADF C05ADF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C05ADF locates a zero of a continuous function in a given interval by a combination of the methods of linear interpolation, extrapolation and bisection. 2. Specification SUBROUTINE C05ADF (A, B, EPS, ETA, F, X, IFAIL) INTEGER IFAIL DOUBLE PRECISION A, B, EPS, ETA, F, X EXTERNAL F 3. Description The routine attempts to obtain an approximation to a simple zero of the function f(x) given an initial interval [a,b] such that f(a)*f(b)<=0. The zero is found by calls to C05AZF(*) whose specification should be consulted for details of the method used. The approximation x to the zero (alpha) is determined so that one or both of the following criteria are satisfied: (i) |x-(alpha)| 0.0. 4: ETA -- DOUBLE PRECISION Input On entry: a value such that if |f(x)|0.0. IFAIL= 2 Too much accuracy has been requested in the computation, that is, EPS is too small for the computer being used. The final value of X is an accurate approximation to the zero. IFAIL= 3 A change in sign of f(x) has been determined as occurring near the point defined by the final value of X. However, there is some evidence that this sign-change corresponds to a pole of f(x). IFAIL= 4 Indicates that a serious error has occurred in C05AZF(*). Check all routine calls. Seek expert help. 7. Accuracy This depends on the value of EPS and ETA. If full machine accuracy is required, they may be set very small, resulting in an error exit with IFAIL = 2, although this may involve more iterations than a lesser accuracy. The user is recommended to set ETA = 0.0 and to use EPS to control the accuracy, unless he has considerable knowledge of the size of f(x) for values of x near the zero. 8. Further Comments The time taken by the routine depends primarily on the time spent evaluating F (see Section 5). If it is important to determine an interval of length less than EPS containing the zero, or if the function F is expensive to evaluate and the number of calls to F is to be restricted, then use of C05AZF(*) is recommended. Use of C05AZF(*) is also recommended when the structure of the problem to be solved does not permit a simple function F to be written: the reverse communication facilities of C05AZF(*) are more flexible than the direct communication of F required by C05ADF. 9. Example -x The example program below calculates the zero of e -x within the interval [0,1] to approximately 5 decimal places. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc05nbf}{NAG On-line Documentation: c05nbf} \beginscroll \begin{verbatim} C05NBF(3NAG) Foundation Library (12/10/92) C05NBF(3NAG) C05 -- Roots of One or More Transcendental Equations C05NBF C05NBF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C05NBF is an easy-to-use routine to find a solution of a system of nonlinear equations by a modification of the Powell hybrid method. 2. Specification SUBROUTINE C05NBF (FCN, N, X, FVEC, XTOL, WA, LWA, IFAIL) INTEGER N, LWA, IFAIL DOUBLE PRECISION X(N), FVEC(N), XTOL, WA(LWA) EXTERNAL FCN 3. Description The system of equations is defined as: f (x ,x ,...,x )=0, for i=1,2,...,n. i 1 2 n C05NBF is based upon the MINPACK routine HYBRD1 (More et al [1]). It chooses the correction at each step as a convex combination of the Newton and scaled gradient directions. Under reasonable conditions this guarantees global convergence for starting points far from the solution and a fast rate of convergence. The Jacobian is updated by the rank-1 method of Broyden. At the starting point the Jacobian is approximated by forward differences, but these are not used again until the rank-1 method fails to produce satisfactory progress. For more details see Powell [2]. 4. References [1] More J J, Garbow B S and Hillstrom K E User Guide for MINPACK-1. Technical Report ANL-80-74. Argonne National Laboratory. [2] Powell M J D (1970) A Hybrid Method for Nonlinear Algebraic Equations. Numerical Methods for Nonlinear Algebraic Equations. (ed P Rabinowitz) Gordon and Breach. 5. Parameters 1: FCN -- SUBROUTINE, supplied by the user. External Procedure FCN must return the values of the functions f at a point x. i Its specification is: SUBROUTINE FCN (N, X, FVEC, IFLAG) INTEGER N, IFLAG DOUBLE PRECISION X(N), FVEC(N) 1: N -- INTEGER Input On entry: the number of equations, n. 2: X(N) -- DOUBLE PRECISION array Input On entry: the components of the point x at which the functions must be evaluated. 3: FVEC(N) -- DOUBLE PRECISION array Output On exit: the function values f (x) (unless IFLAG is i set to a negative value by FCN). 4: IFLAG -- INTEGER Input/Output On entry: IFLAG > 0. On exit: in general, IFLAG should not be reset by FCN. If, however, the user wishes to terminate execution (perhaps because some illegal point X has been reached), then IFLAG should be set to a negative integer. This value will be returned through IFAIL. FCN must be declared as EXTERNAL in the (sub)program from which C05NBF is called. Parameters denoted as Input must not be changed by this procedure. 2: N -- INTEGER Input On entry: the number of equations, n. Constraint: N > 0. 3: X(N) -- DOUBLE PRECISION array Input/Output On entry: an initial guess at the solution vector. On exit: the final estimate of the solution vector. 4: FVEC(N) -- DOUBLE PRECISION array Output On exit: the function values at the final point, X. 5: XTOL -- DOUBLE PRECISION Input On entry: the accuracy in X to which the solution is required. Suggested value: the square root of the machine precision. Constraint: XTOL >= 0.0. 6: WA(LWA) -- DOUBLE PRECISION array Workspace 7: LWA -- INTEGER Input On entry: the dimension of the array WA. Constraint: LWA>=N*(3*N+13)/2. 8: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL< 0 The user has set IFLAG negative in FCN. The value of IFAIL will be the same as the user's setting of IFLAG. IFAIL= 1 On entry N <= 0, or XTOL < 0.0, or LWA 0. 3: X(N) -- DOUBLE PRECISION array Input/Output On entry: an initial guess at the solution vector. On exit: the final estimate of the solution vector. 4: FVEC(N) -- DOUBLE PRECISION array Output On exit: the function values at the final point, X. 5: FJAC(LDFJAC,N) -- DOUBLE PRECISION array Output On exit: the orthogonal matrix Q produced by the QR factorization of the final approximate Jacobian. 6: LDFJAC -- INTEGER Input On entry: the first dimension of the array FJAC as declared in the (sub)program from which C05PBF is called. Constraint: LDFJAC >= N. 7: XTOL -- DOUBLE PRECISION Input On entry: the accuracy in X to which the solution is required. Suggested value: the square root of the machine precision. Constraint: XTOL >= 0.0. 8: WA(LWA) -- DOUBLE PRECISION array Workspace 9: LWA -- INTEGER Input On entry: the dimension of the array WA. Constraint: LWA>=N*(N+13)/2. 10: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL< 0 A negative value of IFAIL indicates an exit from C05PBF because the user has set IFLAG negative in FCN. The value of IFAIL will be the same as the user's setting of IFLAG. IFAIL= 1 On entry N <= 0, or LDFJAC < N, or XTOL < 0.0, or LWA= M. 7: XP(N) -- DOUBLE PRECISION array Output On exit: when MODE = 1, XP is set to a neighbouring point to X. 8: FVECP(M) -- DOUBLE PRECISION array Input On entry: when MODE = 2, FVECP must contain the functions evaluated at XP. 9: MODE -- INTEGER Input On entry: the value 1 on the first call and the value 2 on the second call of C05ZAF. 10: ERR(M) -- DOUBLE PRECISION array Output On exit: when MODE = 2, ERR contains measures of correctness of the respective gradients. If there is no loss of significance (see Section 8), then if ERR(i) is 1.0 the i th user-supplied gradient is correct, whilst if ERR(i) is 0. 0 the ith gradient is incorrect. For values of ERR(i) between 0.0 and 1.0 the categorisation is less certain. In general, a value of ERR(i)>0.5 indicates that the ith gradient is probably correct. 6. Error Indicators and Warnings None. 7. Accuracy See below. 8. Further Comments The time required by C05ZAF increases with M and N. C05ZAF does not perform reliably if cancellation or rounding errors cause a severe loss of significance in the evaluation of a function. Therefore, none of the components of x should be unusually small (in particular, zero) or any other value which may cause loss of significance. The relative differences between corresponding elements of FVECP and FVEC should be at least two orders of magnitude greater than the machine precision. 9. Example This example checks the Jacobian matrix for a problem with 15 functions of 3 variables. The results indicate that the first 7 gradients are probably incorrect (this is caused by a deliberate error in the code to calculate the Jacobian). The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06}{NAG On-line Documentation: c06} \beginscroll \begin{verbatim} C06(3NAG) Foundation Library (12/10/92) C06(3NAG) C06 -- Summation of Series Introduction -- C06 Chapter C06 Summation of Series 1. Scope of the Chapter This chapter is concerned with calculating the discrete Fourier transform of a sequence of real or complex data values, and applying it to calculate convolutions and correlations. 2. Background to the Problems 2.1. Discrete Fourier Transforms 2.1.1. Complex transforms Most of the routines in this chapter calculate the finite discrete Fourier transform (DFT) of a sequence of n complex numbers z , for j=0,1,...,n-1. The transform is defined by: j n-1 ^ 1 -- ( 2(pi)jk) z = --- > z exp(-i -------) (1) k -- j ( n ) \/n j=0 for k=0,1,...,n-1. Note that equation (1) makes sense for all ^ integral k and with this extension z is periodic with period n, k ^ ^ ^ ^ i.e. z =z , and in particular z =z . k k+-n -k n-k ^ ^ If we write z =x +iy and z =a +ib , then the definition of z j j j k k k k may be written in terms of sines and cosines as: n-1 1 -- ( ( 2(pi)jk) ( 2(pi)jk)) a = --- > (x cos( -------)+y sin( -------)) k -- ( j ( n ) j ( n )) \/n j=0 n-1 1 -- ( ( 2(pi)jk) ( 2(pi)jk)) b = --- > (y cos( -------)-x sin( -------)). k -- ( j ( n ) j ( n )) \/n j=0 The original data values z may conversely be recovered from the j ^ transform z by an inverse discrete Fourier transform: k n-1 1 -- ^ ( 2(pi)jk) z = --- > z exp(+i -------) (2) j -- k ( n ) \/n k=0 for j=0,1,...,n-1. If we take the complex conjugate of (2), we ^ find that the sequence z is the DFT of the sequence z . Hence j k ^ the inverse DFT of the sequence z may be obtained by: taking the k ^ complex conjugates of the z ; performing a DFT; and taking the k complex conjugates of the result. Notes: definitions of the discrete Fourier transform vary. Sometimes (2) is used as the definition of the DFT, and (1) as the definition of the inverse. Also the scale-factor of 1/\/n may be omitted in the definition of the DFT, and replaced by 1/n in the definition of the inverse. 2.1.2. Real transforms If the original sequence is purely real valued, i.e. z =x , then j j n-1 ^ 1 -- ( 2(pi)jk) z =a +ib = --- > x exp(-i -------) k k k -- j ( n ) \/n j=0 ^ ^ and z is the complex conjugate of z . Thus the DFT of a real n-k k sequence is a particular type of complex sequence, called a Hermitian sequence, or half-complex or conjugate symmetric with the properties: a =a b =-b b =0 and, if n is even, b =0. n-k k n-k k 0 n/2 Thus a Hermitian sequence of n complex data values can be represented by only n, rather than 2n, independent real values. This can obviously lead to economies in storage, the following scheme being used in this chapter: the real parts a for k 0<=k<=n/2 are stored in normal order in the first n/2+1 locations of an array X of length n; the corresponding non-zero imaginary parts are stored in reverse order in the remaining locations of X. In other words, if X is declared with bounds (0:n-1) in the ^ user's (sub)program, the real and imaginary parts of z are k stored as follows: if n=2s if n=2s-1 X(0) a a 0 0 X(1) a a 1 1 X(2) a a 2 2 . . . . . . . . . X(s-1) a a s-1 s-1 X(s) a b s s-1 X(s+1) b b s-1 s-2 . . . . . . . . . X(n-2) b b 2 2 X(n-1) b b 1 1 ( n/2-1 ) 1 ( -- ( ( 2(pi)jk) ( 2(pi)jk)) ) Hence x = ---(a +2 > (a cos( -------)-b sin( -------))+a ) j ( 0 -- ( k ( n ) k ( n )) n/2) \/n( k=0 ) where a = 0 if n is odd. n/2 2.1.3. Fourier integral transforms The usual application of the discrete Fourier transform is that of obtaining an approximation of the Fourier integral transform +infty / F(s)= | f(t)exp(-i2(pi)st)dt / -infty when f(t) is negligible outside some region (0,c). Dividing the region into n equal intervals we have n-1 c -- F(s)~= - > f exp(-i2(pi)sjc/n) n -- j j=0 and so n-1 c -- F ~= - > f exp(-i2(pi)jk/n) k n -- j j=0 for k=0,1,...,n-1, where f =f(jc/n) and F =F(k/c). j k Hence the discrete Fourier transform gives an approximation to the Fourier integral transform in the region s=0 to s=n/c. If the function f(t) is defined over some more general interval (a,b), then the integral transform can still be approximated by the discrete transform provided a shift is applied to move the point a to the origin. 2.1.4. Convolutions and correlations One of the most important applications of the discrete Fourier transform is to the computation of the discrete convolution or correlation of two vectors x and y defined (as in Brigham [1]) by: n-1 -- convolution: z = > x y k -- j k-j j=0 n-1 -- correlation: w = > x y k -- j k+j j=0 (Here x and y are assumed to be periodic with period n.) Under certain circumstances (see Brigham [1]) these can be used as approximations to the convolution or correlation integrals defined by: +infty / z(s)= | x(t)y(s-t)dt / -infty and +infty / w(s)= | x(t)y(s+t)dt, -infty10 ). Group 2 routines are designed to perform several transforms in a single call, all with the same value of n. They do however require more working storage. Even on scalar processors, they may be somewhat faster than repeated calls to Group 1 routines because of reduced overheads and because they pre-compute and store the required values of trigonometric functions. Group 2 routines impose no practical restrictions on the value of n; however the fast Fourier transform algorithm ceases to be 'fast' if applied to values of n which cannot be expressed as a product of small prime factors. All the above routines are particularly efficient if the only prime factors of n are 2, 3 or 5. If extensive use is to be made of these routines, users who are concerned about efficiency are advised to conduct their own timing tests. To compute inverse discrete Fourier transforms the above routines should be used in conjunction with the utility routines C06GBF, C06GCF and C06GQF which form the complex conjugate of a Hermitian or general sequence of complex data values. 3.2. Multi-dimensional Fourier Transforms C06FUF computes a 2-dimensional discrete Fourier transform of a 2-dimensional sequence of complex data values. This is defined by n -1 n -1 1 2 ( 2(pi)j k ) ( 2(pi)j k ) ^ 1 -- -- ( 1 1) ( 2 2) z = ------- > > z exp(-i ---------)exp(-i ---------). -- -- ( n ) ( n ) k k /n n j =0 j =0 j j ( 1 ) ( 2 ) 1 2 \/ 1 2 1 2 1 2 3.3. Convolution and Correlation C06EKF computes either the discrete convolution or the discrete correlation of two real vectors. 3.4. Index Complex conjugate, complex sequence C06GCF Hermitian sequence C06GBF multiple Hermitian sequences C06GQF Complex sequence from Hermitian sequences C06GSF Convolution or Correlation real vectors C06EKF Discrete Fourier Transform two-dimensional complex sequence C06FUF one-dimensional, multiple transforms complex sequence C06FRF Hermitian sequence C06FQF real sequence C06FPF one-dimensional, single transforms complex sequence C06ECF Hermitian sequence C06EBF real sequence C06EAF C06 -- Summation of Series Contents -- C06 Chapter C06 Summation of Series C06EAF Single 1-D real discrete Fourier transform, no extra workspace C06EBF Single 1-D Hermitian discrete Fourier transform, no extra workspace C06ECF Single 1-D complex discrete Fourier transform, no extra workspace C06EKF Circular convolution or correlation of two real vectors, no extra workspace C06FPF Multiple 1-D real discrete Fourier transforms C06FQF Multiple 1-D Hermitian discrete Fourier transforms C06FRF Multiple 1-D complex discrete Fourier transforms C06FUF 2-D complex discrete Fourier transform C06GBF Complex conjugate of Hermitian sequence C06GCF Complex conjugate of complex sequence C06GQF Complex conjugate of multiple Hermitian sequences C06GSF Convert Hermitian sequences to general complex sequences \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06eaf}{NAG On-line Documentation: c06eaf} \beginscroll \begin{verbatim} C06EAF(3NAG) Foundation Library (12/10/92) C06EAF(3NAG) C06 -- Summation of Series C06EAF C06EAF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06EAF calculates the discrete Fourier transform of a sequence of n real data values. (No extra workspace required.) 2. Specification SUBROUTINE C06EAF (X, N, IFAIL) INTEGER N, IFAIL DOUBLE PRECISION X(N) 3. Description Given a sequence of n real data values x , for j=0,1,...,n-1, j this routine calculates their discrete Fourier transform defined by: n-1 ^ 1 -- ( 2(pi)jk) z = --- > x *exp(-i -------), k=0,1,...,n-1. k -- j ( n ) \/n j=0 1 (Note the scale factor of --- in this definition.) The \/n ^ transformed values z are complex, but they form a Hermitian k ^ ^ sequence (i.e., z is the complex conjugate of z ), so they are n-k k completely determined by n real numbers (see also the Chapter Introduction). To compute the inverse discrete Fourier transform defined by: n-1 ^ 1 -- ( 2(pi)jk) w = --- > x *exp(+i -------), k -- j ( n ) \/n j=0 this routine should be followed by a call of C06GBF to form the ^ complex conjugates of the z . k The routine uses the fast Fourier transform (FFT) algorithm (Brigham [1]). There are some restrictions on the value of n (see Section 5). 4. References [1] Brigham E O (1973) The Fast Fourier Transform. Prentice- Hall. 5. Parameters 1: X(N) -- DOUBLE PRECISION array Input/Output On entry: if X is declared with bounds (0:N-1) in the (sub) program from which C06EAF is called, then X(j) must contain x , for j=0,1,...,n-1. On exit: the discrete Fourier j transform stored in Hermitian form. If the components of the ^ transform z are written as a +ib , and if X is declared k k k with bounds (0:N-1) in the (sub)program from which C06EAF is called, then for 0<=k<=n/2, a is contained in X(k), and for k 1<=k<=(n-1)/2, b is contained in X(n-k). (See also Section k 2.1.2 of the Chapter Introduction, and the Example Program.) 2: N -- INTEGER Input On entry: the number of data values, n. The largest prime factor of N must not exceed 19, and the total number of prime factors of N, counting repetitions, must not exceed 20. Constraint: N > 1. 3: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 At least one of the prime factors of N is greater than 19. IFAIL= 2 N has more than 20 prime factors. IFAIL= 3 N <= 1. 7. Accuracy Some indication of accuracy can be obtained by performing a subsequent inverse transform and comparing the results with the original sequence (in exact arithmetic they would be identical). 8. Further Comments The time taken by the routine is approximately proportional to n*logn, but also depends on the factorization of n. The routine is somewhat faster than average if the only prime factors of n are 2, 3 or 5; and fastest of all if n is a power of 2. On the other hand, the routine is particularly slow if n has several unpaired prime factors, i.e., if the 'square-free' part of n has several factors. For such values of n, routine C06FAF(*) (which requires an additional n elements of workspace) is considerably faster. 9. Example This program reads in a sequence of real data values, and prints their discrete Fourier transform (as computed by C06EAF), after expanding it from Hermitian form into a full complex sequence. It then performs an inverse transform using C06GBF and C06EBF, and prints the sequence so obtained alongside the original data values. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06ebf}{NAG On-line Documentation: c06ebf} \beginscroll \begin{verbatim} C06EBF(3NAG) Foundation Library (12/10/92) C06EBF(3NAG) C06 -- Summation of Series C06EBF C06EBF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06EBF calculates the discrete Fourier transform of a Hermitian sequence of n complex data values. (No extra workspace required.) 2. Specification SUBROUTINE C06EBF (X, N, IFAIL) INTEGER N, IFAIL DOUBLE PRECISION X(N) 3. Description Given a Hermitian sequence of n complex data values z (i.e., a j sequence such that z is real and z is the complex conjugate 0 n-j of z , for j=1,2,...,n-1) this routine calculates their discrete j Fourier transform defined by: n-1 ^ 1 -- ( 2(pi)jk) x = --- > z *exp(-i -------), k=0,1,...,n-1. k -- j ( n ) \/n j=0 1 (Note the scale factor of --- in this definition.) The \/n ^ transformed values x are purely real (see also the the Chapter k Introduction). To compute the inverse discrete Fourier transform defined by: n-1 ^ 1 -- ( 2(pi)jk) y = --- > z *exp(+i -------), k -- j ( n ) \/n j=0 this routine should be preceded by a call of C06GBF to form the complex conjugates of the z . j The routine uses the fast Fourier transform (FFT) algorithm (Brigham [1]). There are some restrictions on the value of n (see Section 5). 4. References [1] Brigham E O (1973) The Fast Fourier Transform. Prentice- Hall. 5. Parameters 1: X(N) -- DOUBLE PRECISION array Input/Output On entry: the sequence to be transformed stored in Hermitian form. If the data values z are written as x +iy , j j j and if X is declared with bounds (0:N-1) in the subroutine from which C06EBF is called, then for 0<=j<=n/2, x is j contained in X(j), and for 1<=j<=(n-1)/2, y is contained in j X(n-j). (See also Section 2.1.2 of the Chapter Introduction and the Example Program.) On exit: the components of the ^ discrete Fourier transform x . If X is declared with bounds k (0:N-1) in the (sub)program from which C06EBF is called, ^ then x is stored in X(k), for k=0,1,...,n-1. k 2: N -- INTEGER Input On entry: the number of data values, n. The largest prime factor of N must not exceed 19, and the total number of prime factors of N, counting repetitions, must not exceed 20. Constraint: N > 1. 3: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 At least one of the prime factors of N is greater than 19. IFAIL= 2 N has more than 20 prime factors. IFAIL= 3 N <= 1. 7. Accuracy Some indication of accuracy can be obtained by performing a subsequent inverse transform and comparing the results with the original sequence (in exact arithmetic they would be identical). 8. Further Comments The time taken by the routine is approximately proportional to n*logn, but also depends on the factorization of n. The routine is somewhat faster than average if the only prime factors of n are 2, 3 or 5; and fastest of all if n is a power of 2. On the other hand, the routine is particularly slow if n has several unpaired prime factors, i.e., if the 'square-free' part of n has several factors. For such values of n, routine C06FBF(*) (which requires an additional n elements of workspace) is considerably faster. 9. Example This program reads in a sequence of real data values which is assumed to be a Hermitian sequence of complex data values stored in Hermitian form. The input sequence is expanded into a full complex sequence and printed alongside the original sequence. The discrete Fourier transform (as computed by C06EBF) is printed out. The program then performs an inverse transform using C06EAF and C06GBF, and prints the sequence so obtained alongside the original data values. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06ecf}{NAG On-line Documentation: c06ecf} \beginscroll \begin{verbatim} C06ECF(3NAG) Foundation Library (12/10/92) C06ECF(3NAG) C06 -- Summation of Series C06ECF C06ECF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06ECF calculates the discrete Fourier transform of a sequence of n complex data values. (No extra workspace required.) 2. Specification SUBROUTINE C06ECF (X, Y, N, IFAIL) INTEGER N, IFAIL DOUBLE PRECISION X(N), Y(N) 3. Description Given a sequence of n complex data values z , for j=0,1,...,n-1, j this routine calculates their discrete Fourier transform defined by: n-1 ^ 1 -- ( 2(pi)jk) z = --- > z *exp(-i -------), k=0,1,...,n-1. k -- j ( n ) \/n j=0 1 (Note the scale factor of --- in this definition.) \/n To compute the inverse discrete Fourier transform defined by: n-1 ^ 1 -- ( 2(pi)jk) w = --- > z *exp(+i -------), k -- j ( n ) \/n j=0 this routine should be preceded and followed by calls of C06GCF ^ to form the complex conjugates of the z and the z . j k The routine uses the fast Fourier transform (FFT) algorithm (Brigham [1]). There are some restrictions on the value of n (see Section 5). 4. References [1] Brigham E O (1973) The Fast Fourier Transform. Prentice- Hall. 5. Parameters 1: X(N) -- DOUBLE PRECISION array Input/Output On entry: if X is declared with bounds (0:N-1) in the (sub) program from which C06ECF is called, then X(j) must contain x , the real part of z , for j=0,1,...,n-1. On exit: the j j real parts a of the components of the discrete Fourier k transform. If X is declared with bounds (0:N-1) in the (sub) program from which C06ECF is called, then a is contained in k X(k), for k=0,1,...,n-1. 2: Y(N) -- DOUBLE PRECISION array Input/Output On entry: if Y is declared with bounds (0:N-1) in the (sub) program from which C06ECF is called, then Y(j) must contain y , the imaginary part of z , for j=0,1,...,n-1. On exit: j j the imaginary parts b of the components of the discrete k Fourier transform. If Y is declared with bounds (0:N-1) in the (sub)program from which C06ECF is called, then b is k contained in Y(k), for k=0,1,...,n-1. 3: N -- INTEGER Input On entry: the number of data values, n. The largest prime factor of N must not exceed 19, and the total number of prime factors of N, counting repetitions, must not exceed 20. Constraint: N > 1. 4: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 At least one of the prime factors of N is greater than 19. IFAIL= 2 N has more than 20 prime factors. IFAIL= 3 N <= 1. 7. Accuracy Some indication of accuracy can be obtained by performing a subsequent inverse transform and comparing the results with the original sequence (in exact arithmetic they would be identical). 8. Further Comments The time taken by the routine is approximately proportional to n*logn, but also depends on the factorization of n. The routine is somewhat faster than average if the only prime factors of n are 2, 3 or 5; and fastest of all if n is a power of 2. On the other hand, the routine is particularly slow if n has several unpaired prime factors, i.e., if the 'square-free' part of n has several factors. For such values of n, routine C06FCF(*) (which requires an additional n real elements of workspace) is considerably faster. 9. Example This program reads in a sequence of complex data values and prints their discrete Fourier transform. It then performs an inverse transform using C06GCF and C06ECF, and prints the sequence so obtained alongside the original data values. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06ekf}{NAG On-line Documentation: c06ekf} \beginscroll \begin{verbatim} C06EKF(3NAG) Foundation Library (12/10/92) C06EKF(3NAG) C06 -- Summation of Series C06EKF C06EKF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06EKF calculates the circular convolution or correlation of two real vectors of period n. No extra workspace is required. 2. Specification SUBROUTINE C06EKF (JOB, X, Y, N, IFAIL) INTEGER JOB, N, IFAIL DOUBLE PRECISION X(N), Y(N) 3. Description This routine computes: if JOB =1, the discrete convolution of x and y, defined by: n-1 n-1 -- -- z = > x y = > x y ; k -- j k-j -- k-j j j=0 j=0 if JOB =2, the discrete correlation of x and y defined by: n-1 -- w = > x y . k -- j k+j j=0 Here x and y are real vectors, assumed to be periodic, with period n, i.e., x =x =x =...; z and w are then also j j+-n j+-2n periodic with period n. Note: this usage of the terms 'convolution' and 'correlation' is taken from Brigham [1]. The term 'convolution' is sometimes used to denote both these computations. ^ ^ ^ ^ If x, y, z and w are the discrete Fourier transforms of these sequences, n-1 ^ 1 -- ( 2(pi)jk) i.e., x = --- > x *exp(-i -------), etc, k -- j ( n ) \/n j=0 ^ ^ ^ then z =\/n.x y k k k ^ ^ ^ and w =\/n.x y k k k (the bar denoting complex conjugate). This routine calls the same auxiliary routines as C06EAF and C06EBF to compute discrete Fourier transforms, and there are some restrictions on the value of n. 4. References [1] Brigham E O (1973) The Fast Fourier Transform. Prentice- Hall. 5. Parameters 1: JOB -- INTEGER Input On entry: the computation to be performed: n-1 -- if JOB = 1, z = > x y (convolution); k -- j k-j j=0 n-1 -- if JOB = 2, w = > x y (correlation). k -- j k+j j=0 Constraint: JOB = 1 or 2. 2: X(N) -- DOUBLE PRECISION array Input/Output On entry: the elements of one period of the vector x. If X is declared with bounds (0:N-1) in the (sub)program from which C06EKF is called, then X(j) must contain x , for j j=0,1,...,n-1. On exit: the corresponding elements of the discrete convolution or correlation. 3: Y(N) -- DOUBLE PRECISION array Input/Output On entry: the elements of one period of the vector y. If Y is declared with bounds (0:N-1) in the (sub)program from which C06EKF is called, then Y(j) must contain y , for j j=0,1,...,n-1. On exit: the discrete Fourier transform of the convolution or correlation returned in the array X; the transform is stored in Hermitian form, exactly as described in the document C06EAF. 4: N -- INTEGER Input On entry: the number of values, n, in one period of the vectors X and Y. The largest prime factor of N must not exceed 19, and the total number of prime factors of N, counting repetitions, must not exceed 20. Constraint: N > 1. 5: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 At least one of the prime factors of N is greater than 19. IFAIL= 2 N has more than 20 prime factors. IFAIL= 3 N <= 1. IFAIL= 4 JOB /= 1 or 2. 7. Accuracy The results should be accurate to within a small multiple of the machine precision. 8. Further Comments The time taken by the routine is approximately proportional to n*logn, but also depends on the factorization of n. The routine is faster than average if the only prime factors are 2, 3 or 5; and fastest of all if n is a power of 2. The routine is particularly slow if n has several unpaired prime factors, i.e., if the 'square free' part of n has several factors. For such values of n, routine C06FKF(*) is considerably faster (but requires an additional workspace of n elements). 9. Example This program reads in the elements of one period of two real vectors x and y and prints their discrete convolution and correlation (as computed by C06EKF). In realistic computations the number of data values would be much larger. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06fpf}{NAG On-line Documentation: c06fpf} \beginscroll \begin{verbatim} C06FPF(3NAG) Foundation Library (12/10/92) C06FPF(3NAG) C06 -- Summation of Series C06FPF C06FPF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06FPF computes the discrete Fourier transforms of m sequences, each containing n real data values. This routine is designed to be particularly efficient on vector processors. 2. Specification SUBROUTINE C06FPF (M, N, X, INIT, TRIG, WORK, IFAIL) INTEGER M, N, IFAIL DOUBLE PRECISION X(M*N), TRIG(2*N), WORK(M*N) CHARACTER*1 INIT 3. Description p Given m sequences of n real data values x , for j=0,1,...,n-1; j p=1,2,...,m, this routine simultaneously calculates the Fourier transforms of all the sequences defined by: n-1 ^p 1 -- p ( 2(pi)jk) z = --- > x *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m. k -- j ( n ) \/n j=0 1 (Note the scale factor --- in this definition.) \/n ^p The transformed values z are complex, but for each value of p k ^p ^p the z form a Hermitian sequence (i.e.,z is the complex k n-k ^p conjugate of z ), so they are completely determined by mn real k numbers (see also the Chapter Introduction). The discrete Fourier transform is sometimes defined using a positive sign in the exponential term: n-1 ^p 1 -- p ( 2(pi)jk) z = --- > x *exp(+i -------). k -- j ( n ) \/n j=0 To compute this form, this routine should be followed by a call ^p to C06GQF to form the complex conjugates of the z . k The routine uses a variant of the fast Fourier transform (FFT) algorithm (Brigham [1]) known as the Stockham self-sorting algorithm, which is described in Temperton [2]. Special coding is provided for the factors 2, 3, 4, 5 and 6. This routine is designed to be particularly efficient on vector processors, and it becomes especially fast as M, the number of transforms to be computed in parallel, increases. 4. References [1] Brigham E O (1973) The Fast Fourier Transform. Prentice- Hall. [2] Temperton C (1983) Fast Mixed-Radix Real Fourier Transforms. J. Comput. Phys. 52 340--350. 5. Parameters 1: M -- INTEGER Input On entry: the number of sequences to be transformed, m. Constraint: M >= 1. 2: N -- INTEGER Input On entry: the number of real values in each sequence, n. Constraint: N >= 1. 3: X(M,N) -- DOUBLE PRECISION array Input/Output On entry: the data must be stored in X as if in a two- dimensional array of dimension (1:M,0:N-1); each of the m sequences is stored in a row of the array. In other words, if the data values of the pth sequence to be transformed are p denoted by x , for j=0,1,...,n-1, then the mn elements of j the array X must contain the values 1 2 m 1 2 m 1 2 m x ,x ,...,x , x ,x ,..., x ,...,x ,x ,...,x . 0 0 0 1 1 1 n-1 n-1 n-1 On exit: the m discrete Fourier transforms stored as if in a two-dimensional array of dimension (1:M,0:N-1). Each of the m transforms is stored in a row of the array in Hermitian form, overwriting the corresponding original sequence. If the n components of the discrete Fourier ^p p p p transform z are written as a +ib , then for 0<=k<=n/2, a k k k k p is contained in X(p,k), and for 1<=k<=(n-1)/2, b is k contained in X(p,n-k). (See also Section 2.1.2 of the Chapter Introduction.) 4: INIT -- CHARACTER*1 Input On entry: if the trigonometric coefficients required to compute the transforms are to be calculated by the routine and stored in the array TRIG, then INIT must be set equal to 'I' (Initial call). If INIT contains 'S' (Subsequent call), then the routine assumes that trigonometric coefficients for the specified value of n are supplied in the array TRIG, having been calculated in a previous call to one of C06FPF, C06FQF or C06FRF. If INIT contains 'R' (Restart then the routine assumes that trigonometric coefficients for the particular value of n are supplied in the array TRIG, but does not check that C06FPF, C06FQF or C06FRF have previously been called. This option allows the TRIG array to be stored in an external file, read in and re-used without the need for a call with INIT equal to 'I'. The routine carries out a simple test to check that the current value of n is consistent with the array TRIG. Constraint: INIT = 'I', 'S' or 'R'. 5: TRIG(2*N) -- DOUBLE PRECISION array Input/Output On entry: if INIT = 'S' or 'R', TRIG must contain the required coefficients calculated in a previous call of the routine. Otherwise TRIG need not be set. On exit: TRIG contains the required coefficients (computed by the routine if INIT = 'I'). 6: WORK(M*N) -- DOUBLE PRECISION array Workspace 7: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry M < 1. IFAIL= 2 N < 1. IFAIL= 3 INIT is not one of 'I', 'S' or 'R'. IFAIL= 4 INIT = 'S', but none of C06FPF, C06FQF or C06FRF has previously been called. IFAIL= 5 INIT = 'S' or 'R', but the array TRIG and the current value of N are inconsistent. 7. Accuracy Some indication of accuracy can be obtained by performing a subsequent inverse transform and comparing the results with the original sequence (in exact arithmetic they would be identical). 8. Further Comments The time taken by the routine is approximately proportional to nm*logn, but also depends on the factors of n. The routine is fastest if the only prime factors of n are 2, 3 and 5, and is particularly slow if n is a large prime, or has large prime factors. 9. Example This program reads in sequences of real data values and prints their discrete Fourier transforms (as computed by C06FPF). The Fourier transforms are expanded into full complex form using C06GSF and printed. Inverse transforms are then calculated by calling C06GQF followed by C06FQF showing that the original sequences are restored. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06fqf}{NAG On-line Documentation: c06fqf} \beginscroll \begin{verbatim} C06FQF(3NAG) Foundation Library (12/10/92) C06FQF(3NAG) C06 -- Summation of Series C06FQF C06FQF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06FQF computes the discrete Fourier transforms of m Hermitian sequences, each containing n complex data values. This routine is designed to be particularly efficient on vector processors. 2. Specification SUBROUTINE C06FQF (M, N, X, INIT, TRIG, WORK, IFAIL) INTEGER M, N, IFAIL DOUBLE PRECISION X(M*N), TRIG(2*N), WORK(M*N) CHARACTER*1 INIT 3. Description p Given m Hermitian sequences of n complex data values z , for j j=0,1,...,n-1; p=1,2,...,m, this routine simultaneously calculates the Fourier transforms of all the sequences defined by: n-1 ^p 1 -- p ( 2(pi)jk) x = --- > z *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m. k -- j ( n ) \/n j=0 1 (Note the scale factor --- in this definition.) \/n The transformed values are purely real (see also the Chapter Introduction). The discrete Fourier transform is sometimes defined using a positive sign in the exponential term n-1 ^p 1 -- p ( 2(pi)jk) x = --- > z *exp(+i -------). k -- j ( n ) \/n j=0 To compute this form, this routine should be preceded by a call ^p to C06GQF to form the complex conjugates of the z . j The routine uses a variant of the fast Fourier transform (FFT) algorithm (Brigham [1]) known as the Stockham self-sorting algorithm, which is described in Temperton [2]. Special code is included for the factors 2, 3, 4, 5 and 6. This routine is designed to be particularly efficient on vector processors, and it becomes especially fast as m, the number of transforms to be computed in parallel, increases. 4. References [1] Brigham E O (1973) The Fast Fourier Transform. Prentice- Hall. [2] Temperton C (1983) Fast Mixed-Radix Real Fourier Transforms. J. Comput. Phys. 52 340--350. 5. Parameters 1: M -- INTEGER Input On entry: the number of sequences to be transformed, m. Constraint: M >= 1. 2: N -- INTEGER Input On entry: the number of data values in each sequence, n. Constraint: N >= 1. 3: X(M,N) -- DOUBLE PRECISION array Input/Output On entry: the data must be stored in X as if in a two- dimensional array of dimension (1:M,0:N-1); each of the m sequences is stored in a row of the array in Hermitian form. p p p If the n data values z are written as x +iy , then for j j j p 0<=j<=n/2, x is contained in X(p,j), and for 1<=j<=(n-1)/2, j p y is contained in X(p,n-j). (See also Section 2.1.2 of the j Chapter Introduction.) On exit: the components of the m discrete Fourier transforms, stored as if in a two- dimensional array of dimension (1:M,0:N-1). Each of the m transforms is stored as a row of the array, overwriting the corresponding original sequence. If the n components of the ^p discrete Fourier transform are denoted by x , for k k=0,1,...,n-1, then the mn elements of the array X contain the values ^1 ^2 ^m ^1 ^2 ^m ^1 ^2 ^m x ,x ,...,x , x ,x ,..., x ,...,x ,x ,...,x . 0 0 0 1 1 1 n-1 n-1 n-1 4: INIT -- CHARACTER*1 Input On entry: if the trigonometric coefficients required to compute the transforms are to be calculated by the routine and stored in the array TRIG, then INIT must be set equal to 'I' (Initial call). If INIT contains 'S' (Subsequent call), then the routine assumes that trigonometric coefficients for the specified value of n are supplied in the array TRIG, having been calculated in a previous call to one of C06FPF, C06FQF or C06FRF. If INIT contains 'R' (Restart), then the routine assumes that trigonometric coefficients for the particular value of N are supplied in the array TRIG, but does not check that C06FPF, C06FQF or C06FRF have previously been called. This option allows the TRIG array to be stored in an external file, read in and re-used without the need for a call with INIT equal to 'I'. The routine carries out a simple test to check that the current value of n is compatible with the array TRIG. Constraint: INIT = 'I', 'S' or 'R'. 5: TRIG(2*N) -- DOUBLE PRECISION array Input/Output On entry: if INIT = 'S' or 'R', TRIG must contain the required coefficients calculated in a previous call of the routine. Otherwise TRIG need not be set. On exit: TRIG contains the required coefficients (computed by the routine if INIT = 'I'). 6: WORK(M*N) -- DOUBLE PRECISION array Workspace 7: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry M < 1. IFAIL= 2 On entry N < 1. IFAIL= 3 On entry INIT is not one of 'I', 'S' or 'R'. IFAIL= 4 On entry INIT = 'S', but none of C06FPF, C06FQF and C06FRF has previously been called. IFAIL= 5 On entry INIT = 'S' or 'R', but the array TRIG and the current value of n are inconsistent. 7. Accuracy Some indication of accuracy can be obtained by performing a subsequent inverse transform and comparing the results with the original sequence (in exact arithmetic they would be identical). 8. Further Comments The time taken by the routine is approximately proportional to nm*logn, but also depends on the factors of n. The routine is fastest if the only prime factors of n are 2, 3 and 5, and is particularly slow if n is a large prime, or has large prime factors. 9. Example This program reads in sequences of real data values which are assumed to be Hermitian sequences of complex data stored in Hermitian form. The sequences are expanded into full complex form using C06GSF and printed. The discrete Fourier transforms are then computed (using C06FQF) and printed out. Inverse transforms are then calculated by calling C06FPF followed by C06GQF showing that the original sequences are restored. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06frf}{NAG On-line Documentation: c06frf} \beginscroll \begin{verbatim} C06FRF(3NAG) Foundation Library (12/10/92) C06FRF(3NAG) C06 -- Summation of Series C06FRF C06FRF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06FRF computes the discrete Fourier transforms of m sequences, each containing n complex data values. This routine is designed to be particularly efficient on vector processors. 2. Specification SUBROUTINE C06FRF (M, N, X, Y, INIT, TRIG, WORK, IFAIL) INTEGER M, N, IFAIL DOUBLE PRECISION X(M*N), Y(M*N), TRIG(2*N), WORK(2*M*N) CHARACTER*1 INIT 3. Description p Given m sequences of n complex data values z , for j=0,1,...,n-1; j p=1,2,...,m, this routine simultaneously calculates the Fourier transforms of all the sequences defined by: n-1 ^p 1 -- p ( 2(pi)jk) z = --- > z *exp(-i -------), k=0,1,...,n-1; p=1,2,...,m. k -- j ( n ) \/n j=0 1 (Note the scale factor --- in this definition.) \/n The discrete Fourier transform is sometimes defined using a positive sign in the exponential term n-1 ^p 1 -- p ( 2(pi)jk) z = --- > z *exp(+i -------). k -- j ( n ) \/n j=0 To compute this form, this routine should be preceded and followed by a call of C06GCF to form the complex conjugates of p ^p the z and the z . j k The routine uses a variant of the fast Fourier transform (FFT) algorithm (Brigham [1]) known as the Stockham self-sorting algorithm, which is described in Temperton [2]. Special code is provided for the factors 2, 3, 4, 5 and 6. This routine is designed to be particularly efficient on vector processors, and it becomes especially fast as m, the number of transforms to be computed in parallel, increases. 4. References [1] Brigham E O (1973) The Fast Fourier Transform. Prentice- Hall. [2] Temperton C (1983) Self-sorting Mixed-radix Fast Fourier Transforms. J. Comput. Phys. 52 1--23. 5. Parameters 1: M -- INTEGER Input On entry: the number of sequences to be transformed, m. Constraint: M >= 1. 2: N -- INTEGER Input On entry: the number of complex values in each sequence, n. Constraint: N >= 1. 3: X(M,N) -- DOUBLE PRECISION array Input/Output 4: Y(M,N) -- DOUBLE PRECISION array Input/Output On entry: the real and imaginary parts of the complex data must be stored in X and Y respectively as if in a two- dimensional array of dimension (1:M,0:N-1); each of the m sequences is stored in a row of each array. In other words, if the real parts of the pth sequence to be transformed are p denoted by x , for j=0,1,...,n-1, then the mn elements of j the array X must contain the values 1 2 m 1 2 m 1 2 m x ,x ,...,x , x ,x ,...,x ,..., x ,x ,...,x . 0 0 0 1 1 1 n-1 n-1 n-1 On exit: X and Y are overwritten by the real and imaginary parts of the complex transforms. 5: INIT -- CHARACTER*1 Input On entry: if the trigonometric coefficients required to compute the transforms are to be calculated by the routine and stored in the array TRIG, then INIT must be set equal to 'I' (Initial call). If INIT contains 'S' (Subsequent call), then the routine assumes that trigonometric coefficients for the specified value of n are supplied in the array TRIG, having been calculated in a previous call to one of C06FPF, C06FQF or C06FRF. If INIT contains 'R' (Restart) then the routine assumes that trigonometric coefficients for the particular value of n are supplied in the array TRIG, but does not check that C06FPF, C06FQF or C06FRF have previously been called. This option allows the TRIG array to be stored in an external file, read in and re-used without the need for a call with INIT equal to 'I'. The routine carries out a simple test to check that the current value of n is compatible with the array TRIG. Constraint: INIT = 'I', 'S' or 'R'. 6: TRIG(2*N) -- DOUBLE PRECISION array Input/Output On entry: if INIT = 'S' or 'R', TRIG must contain the required coefficients calculated in a previous call of the routine. Otherwise TRIG need not be set. On exit: TRIG contains the required coefficients (computed by the routine if INIT = 'I'). 7: WORK(2*M*N) -- DOUBLE PRECISION array Workspace 8: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry M < 1. IFAIL= 2 On entry N < 1. IFAIL= 3 On entry INIT is not one of 'I', 'S' or 'R'. IFAIL= 4 On entry INIT = 'S', but none of C06FPF, C06FQF and C06FRF has previously been called. IFAIL= 5 On entry INIT = 'S' or 'R', but the array TRIG and the current value of n are inconsistent. 7. Accuracy Some indication of accuracy can be obtained by performing a subsequent inverse transform and comparing the results with the original sequence (in exact arithmetic they would be identical). 8. Further Comments The time taken by the routine is approximately proportional to nm*logn, but also depends on the factors of n. The routine is fastest if the only prime factors of n are 2, 3 and 5, and is particularly slow if n is a large prime, or has large prime factors. 9. Example This program reads in sequences of complex data values and prints their discrete Fourier transforms (as computed by C06FRF). Inverse transforms are then calculated using C06GCF and C06FRF and printed out, showing that the original sequences are restored. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06fuf}{NAG On-line Documentation: c06fuf} \beginscroll \begin{verbatim} C06FUF(3NAG) Foundation Library (12/10/92) C06FUF(3NAG) C06 -- Summation of Series C06FUF C06FUF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06FUF computes the two-dimensional discrete Fourier transform of a bivariate sequence of complex data values. This routine is designed to be particularly efficient on vector processors. 2. Specification SUBROUTINE C06FUF (M, N, X, Y, INIT, TRIGM, TRIGN, WORK, 1 IFAIL) INTEGER M, N, IFAIL DOUBLE PRECISION X(M*N), Y(M*N), TRIGM(2*M), TRIGN(2*N), 1 WORK(2*M*N) CHARACTER*1 INIT 3. Description This routine computes the two-dimensional discrete Fourier transform of a bivariate sequence of complex data values z , j j 1 2 where j =0,1,...,m-1, j =0,1,...,n-1. 1 2 The discrete Fourier transform is here defined by: m-1 n-1 ( ( j k j k )) ^ 1 -- -- ( ( 1 1 2 2)) z = ---- > > z *exp(-2(pi)i( ----+ ----)), -- -- j j ( ( m n )) k k \/mn j =0 j =0 1 2 1 2 1 2 where k =0,1,...,m-1, k =0,1,...,n-1. 1 2 1 (Note the scale factor of ---- in this definition.) \/mn To compute the inverse discrete Fourier transform, defined with exp(+2(pi)i(...)) in the above formula instead of exp(- 2(pi)i(...)), this routine should be preceded and followed by calls of C06GCF to form the complex conjugates of the data values and the transform. This routine calls C06FRF to perform multiple one-dimensional discrete Fourier transforms by the fast Fourier transform (FFT) algorithm in Brigham [1]. It is designed to be particularly efficient on vector processors. 4. References [1] Brigham E O (1973) The Fast Fourier Transform. Prentice- Hall. [2] Temperton C (1983) Self-sorting Mixed-radix Fast Fourier Transforms. J. Comput. Phys. 52 1--23. 5. Parameters 1: M -- INTEGER Input On entry: the number of rows, m, of the arrays X and Y. Constraint: M >= 1. 2: N -- INTEGER Input On entry: the number of columns, n, of the arrays X and Y. Constraint: N >= 1. 3: X(M,N) -- DOUBLE PRECISION array Input/Output 4: Y(M,N) -- DOUBLE PRECISION array Input/Output On entry: the real and imaginary parts of the complex data values must be stored in arrays X and Y respectively. If X and Y are regarded as two-dimensional arrays of dimension (0:M-1,0:N-1), then X(j ,j ) and Y(j ,j ) must contain the 1 2 1 2 real and imaginary parts of z . On exit: the real and j j 1 2 imaginary parts respectively of the corresponding elements of the computed transform. 5: INIT -- CHARACTER*1 Input On entry: if the trigonometric coefficients required to compute the transforms are to be calculated by the routine and stored in the arrays TRIGM and TRIGN, then INIT must be set equal to 'I', (Initial call). If INIT contains 'S', (Subsequent call), then the routine assumes that trigonometric coefficients for the specified values of m and n are supplied in the arrays TRIGM and TRIGN, having been calculated in a previous call to the routine. If INIT contains 'R', (Restart), then the routine assumes that trigonometric coefficients for the particular values of m and n are supplied in the arrays TRIGM and TRIGN, but does not check that the routine has previously been called. This option allows the TRIGM and TRIGN arrays to be stored in an external file, read in and re-used without the need for a call with INIT equal to 'I'. The routine carries out a simple test to check that the current values of m and n are compatible with the arrays TRIGM and TRIGN. Constraint: INIT = 'I', 'S' or 'R'. 6: TRIGM(2*M) -- DOUBLE PRECISION array Input/Output 7: TRIGN(2*N) -- DOUBLE PRECISION array Input/Output On entry: if INIT = 'S' or 'R',TRIGM and TRIGN must contain the required coefficients calculated in a previous call of the routine. Otherwise TRIGM and TRIGN need not be set. If m=n the same array may be supplied for TRIGM and TRIGN. On exit: TRIGM and TRIGN contain the required coefficients (computed by the routine if INIT = 'I'). 8: WORK(2*M*N) -- DOUBLE PRECISION array Workspace 9: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry M < 1. IFAIL= 2 On entry N < 1. IFAIL= 3 On entry INIT is not one of 'I', 'S' or 'R'. IFAIL= 4 On entry INIT = 'S', but C06FUF has not previously been called. IFAIL= 5 On entry INIT = 'S' or 'R', but at least one of the arrays TRIGM and TRIGN is inconsistent with the current value of M or N. 7. Accuracy Some indication of accuracy can be obtained by performing a subsequent inverse transform and comparing the results with the original sequence (in exact arithmetic they would be identical). 8. Further Comments The time taken by the routine is approximately proportional to mn*log(mn), but also depends on the factorization of the individual dimensions m and n. The routine is somewhat faster than average if their only prime factors are 2, 3 or 5; and fastest of all if they are powers of 2. 9. Example This program reads in a bivariate sequence of complex data values and prints the two-dimensional Fourier transform. It then performs an inverse transform and prints the sequence so obtained, which may be compared to the original data values. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06gbf}{NAG On-line Documentation: c06gbf} \beginscroll \begin{verbatim} C06GBF(3NAG) Foundation Library (12/10/92) C06GBF(3NAG) C06 -- Summation of Series C06GBF C06GBF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06GBF forms the complex conjugate of a Hermitian sequence of n data values. 2. Specification SUBROUTINE C06GBF (X, N, IFAIL) INTEGER N, IFAIL DOUBLE PRECISION X(N) 3. Description This is a utility routine for use in conjunction with C06EAF, C06EBF, C06FAF(*) or C06FBF(*) to calculate inverse discrete Fourier transforms (see the Chapter Introduction). 4. References None. 5. Parameters 1: X(N) -- DOUBLE PRECISION array Input/Output On entry: if the data values z are written as x +iy and j j j if X is declared with bounds (0:N-1) in the (sub)program from which C06GBF is called, then for 0<=j<=n/2, X(j) must contain x (=x ), while for n/2= 1. 3: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 N < 1. 7. Accuracy Exact. 8. Further Comments The time taken by the routine is negligible. 9. Example This program reads in a sequence of real data values, calls C06EAF followed by C06GBF to compute their inverse discrete Fourier transform, and prints this after expanding it from Hermitian form into a full complex sequence. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06gcf}{NAG On-line Documentation: c06gcf} \beginscroll \begin{verbatim} C06GCF(3NAG) Foundation Library (12/10/92) C06GCF(3NAG) C06 -- Summation of Series C06GCF C06GCF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06GCF forms the complex conjugate of a sequence of n data values. 2. Specification SUBROUTINE C06GCF (Y, N, IFAIL) INTEGER N, IFAIL DOUBLE PRECISION Y(N) 3. Description This is a utility routine for use in conjunction with C06ECF or C06FCF(*) to calculate inverse discrete Fourier transforms (see the Chapter Introduction). 4. References None. 5. Parameters 1: Y(N) -- DOUBLE PRECISION array Input/Output On entry: if Y is declared with bounds (0:N-1) in the (sub) program which C06GCF is called, then Y(j) must contain the imaginary part of the jth data value, for 0<=j<=n-1. On exit: these values are negated. 2: N -- INTEGER Input On entry: the number of data values, n. Constraint: N >= 1. 3: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: IFAIL= 1 N < 1. 7. Accuracy Exact. 8. Further Comments The time taken by the routine is negligible. 9. Example This program reads in a sequence of complex data values and prints their inverse discrete Fourier transform as computed by calling C06GCF, followed by C06ECF and C06GCF again. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06gqf}{NAG On-line Documentation: c06gqf} \beginscroll \begin{verbatim} C06GQF(3NAG) Foundation Library (12/10/92) C06GQF(3NAG) C06 -- Summation of Series C06GQF C06GQF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06GQF forms the complex conjugates of m Hermitian sequences, each containing n data values. 2. Specification SUBROUTINE C06GQF (M, N, X, IFAIL) INTEGER M, N, IFAIL DOUBLE PRECISION X(M*N) 3. Description This is a utility routine for use in conjunction with C06FPF and C06FQF to calculate inverse discrete Fourier transforms (see the Chapter Introduction). 4. References None. 5. Parameters 1: M -- INTEGER Input On entry: the number of Hermitian sequences to be conjugated, m. Constraint: M >= 1. 2: N -- INTEGER Input On entry: the number of data values in each Hermitian sequence, n. Constraint: N >= 1. 3: X(M,N) -- DOUBLE PRECISION array Input/Output On entry: the data must be stored in array X as if in a two-dimensional array of dimension (1:M,0:N-1); each of the m sequences is stored in a row of the array in Hermitian p p p form. If the n data values z are written as x +iy , then j j j p for 0<=j<=n/2, x is contained in X(p,j), and for 1<=j<=(n- j p 1)/2, y is contained in X(p,n-j). (See also Section 2.1.2 j of the Chapter Introduction.) On exit: the imaginary parts p p y are negated. The real parts x are not referenced. j j 4: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry M < 1. IFAIL= 2 On entry N < 1. 7. Accuracy Exact. 8. Further Comments None. 9. Example This program reads in sequences of real data values which are assumed to be Hermitian sequences of complex data stored in Hermitian form. The sequences are expanded into full complex form using C06GSF and printed. The sequences are then conjugated (using C06GQF) and the conjugated sequences are expanded into complex form using C06GSF and printed out. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page} \begin{page}{manpageXXc06gsf}{NAG On-line Documentation: c06gsf} \beginscroll \begin{verbatim} C06GSF(3NAG) Foundation Library (12/10/92) C06GSF(3NAG) C06 -- Summation of Series C06GSF C06GSF -- NAG Foundation Library Routine Document Note: Before using this routine, please read the Users' Note for your implementation to check implementation-dependent details. The symbol (*) after a NAG routine name denotes a routine that is not included in the Foundation Library. 1. Purpose C06GSF takes m Hermitian sequences, each containing n data values, and forms the real and imaginary parts of the m corresponding complex sequences. 2. Specification SUBROUTINE C06GSF (M, N, X, U, V, IFAIL) INTEGER M, N, IFAIL DOUBLE PRECISION X(M*N), U(M*N), V(M*N) 3. Description This is a utility routine for use in conjunction with C06FPF and C06FQF (see the Chapter Introduction). 4. References None. 5. Parameters 1: M -- INTEGER Input On entry: the number of Hermitian sequences, m, to be converted into complex form. Constraint: M >= 1. 2: N -- INTEGER Input On entry: the number of data values, n, in each sequence. Constraint: N >= 1. 3: X(M,N) -- DOUBLE PRECISION array Input On entry: the data must be stored in X as if in a two- dimensional array of dimension (1:M,0:N-1); each of the m sequences is stored in a row of the array in Hermitian form. p p p If the n data values z are written as x +iy , then for j j j p 0<=j<=n/2, x is contained in X(p,j), and for 1<=j<=(n-1)/2, j p y is contained in X(p,n-j). (See also Section 2.1.2 of the j Chapter Introduction.) 4: U(M,N) -- DOUBLE PRECISION array Output 5: V(M,N) -- DOUBLE PRECISION array Output On exit: the real and imaginary parts of the m sequences of length n, are stored in U and V respectively, as if in two- dimensional arrays of dimension (1:M,0:N-1); each of the m sequences is stored as if in a row of each array. In other words, if the real parts of the pth sequence are denoted by p x , for j=0,1,...,n-1 then the mn elements of the array U j contain the values 1 2 m 1 2 m 1 2 m x ,x ,...,x , x ,x ,...,x ,..., x ,x ,...,x 0 0 0 1 1 1 n-1 n-1 n-1 6: IFAIL -- INTEGER Input/Output On entry: IFAIL must be set to 0, -1 or 1. For users not familiar with this parameter (described in the Essential Introduction) the recommended value is 0. On exit: IFAIL = 0 unless the routine detects an error (see Section 6). 6. Error Indicators and Warnings Errors detected by the routine: If on entry IFAIL = 0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF). IFAIL= 1 On entry M < 1. IFAIL= 2 On entry N < 1. 7. Accuracy Exact. 8. Further Comments None. 9. Example This program reads in sequences of real data values which are assumed to be Hermitian sequences of complex data stored in Hermitian form. The sequences are then expanded into full complex form using C06GSF and printed. The example program is not reproduced here. The source code for all example programs is distributed with the NAG Foundation Library software and should be available on-line. \end{verbatim} \endscroll \end{page}